Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to apply a delta pulse for RTP and EMD
10 : ! **************************************************************************************************
11 :
12 : MODULE rt_delta_pulse
13 : USE bibliography, ONLY: Mattiat2019,&
14 : Mattiat2022,&
15 : cite_reference
16 : USE cell_types, ONLY: cell_type
17 : USE commutator_rpnl, ONLY: build_com_mom_nl,&
18 : build_com_nl_mag
19 : USE cp_blacs_env, ONLY: cp_blacs_env_type
20 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale
21 : USE cp_cfm_diag, ONLY: cp_cfm_heevd
22 : USE cp_cfm_types, ONLY: cp_cfm_create,&
23 : cp_cfm_release,&
24 : cp_cfm_to_cfm,&
25 : cp_cfm_type
26 : USE cp_control_types, ONLY: dft_control_type,&
27 : rtp_control_type
28 : USE cp_dbcsr_api, ONLY: &
29 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_info, &
30 : dbcsr_init_p, dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
31 : dbcsr_type_symmetric
32 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
33 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
34 : cp_dbcsr_sm_fm_multiply,&
35 : dbcsr_allocate_matrix_set,&
36 : dbcsr_deallocate_matrix_set
37 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
38 : cp_fm_triangular_multiply,&
39 : cp_fm_uplo_to_full
40 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
41 : cp_fm_cholesky_invert,&
42 : cp_fm_cholesky_reduce,&
43 : cp_fm_cholesky_restore
44 : USE cp_fm_diag, ONLY: cp_fm_syevd
45 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
46 : cp_fm_struct_release,&
47 : cp_fm_struct_type
48 : USE cp_fm_types, ONLY: cp_fm_create,&
49 : cp_fm_get_info,&
50 : cp_fm_release,&
51 : cp_fm_set_all,&
52 : cp_fm_to_fm,&
53 : cp_fm_type
54 : USE cp_log_handling, ONLY: cp_get_default_logger,&
55 : cp_logger_type
56 : USE cp_output_handling, ONLY: cp_print_key_unit_nr
57 : USE input_section_types, ONLY: section_get_ival,&
58 : section_get_lval,&
59 : section_vals_get_subs_vals,&
60 : section_vals_type,&
61 : section_vals_val_get
62 : USE kinds, ONLY: dp
63 : USE mathconstants, ONLY: one,&
64 : twopi,&
65 : zero
66 : USE message_passing, ONLY: mp_para_env_type
67 : USE moments_utils, ONLY: get_reference_point
68 : USE parallel_gemm_api, ONLY: parallel_gemm
69 : USE particle_types, ONLY: particle_type
70 : USE qs_dftb_matrices, ONLY: build_dftb_overlap
71 : USE qs_environment_types, ONLY: get_qs_env,&
72 : qs_environment_type
73 : USE qs_kind_types, ONLY: qs_kind_type
74 : USE qs_mo_types, ONLY: get_mo_set,&
75 : mo_set_type
76 : USE qs_moments, ONLY: build_berry_moment_matrix,&
77 : build_local_magmom_matrix,&
78 : build_local_moment_matrix
79 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
80 : USE rt_propagation_types, ONLY: get_rtp,&
81 : rt_prop_create_mos,&
82 : rt_prop_type
83 : #include "../base/base_uses.f90"
84 :
85 : IMPLICIT NONE
86 :
87 : PRIVATE
88 :
89 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_delta_pulse'
90 :
91 : PUBLIC :: apply_delta_pulse
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief Interface to call the delta pulse depending on the type of calculation.
97 : !> \param qs_env ...
98 : !> \param rtp ...
99 : !> \param rtp_control ...
100 : !> \author Update: Guillaume Le Breton (2023.01)
101 : ! **************************************************************************************************
102 :
103 58 : SUBROUTINE apply_delta_pulse(qs_env, rtp, rtp_control)
104 : TYPE(qs_environment_type), POINTER :: qs_env
105 : TYPE(rt_prop_type), POINTER :: rtp
106 : TYPE(rtp_control_type), POINTER :: rtp_control
107 :
108 : CHARACTER(LEN=3), DIMENSION(3) :: rlab
109 : INTEGER :: i, output_unit
110 : LOGICAL :: my_apply_pulse, periodic
111 : REAL(KIND=dp), DIMENSION(3) :: kvec
112 : TYPE(cell_type), POINTER :: cell
113 58 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new, mos_old
114 : TYPE(cp_logger_type), POINTER :: logger
115 58 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
116 : TYPE(dft_control_type), POINTER :: dft_control
117 58 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
118 : TYPE(section_vals_type), POINTER :: input, rtp_section
119 :
120 58 : NULLIFY (logger, input, rtp_section)
121 :
122 116 : logger => cp_get_default_logger()
123 : CALL get_qs_env(qs_env, &
124 : cell=cell, &
125 : input=input, &
126 : dft_control=dft_control, &
127 58 : matrix_s=matrix_s)
128 58 : rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
129 : output_unit = cp_print_key_unit_nr(logger, rtp_section, "PRINT%PROGRAM_RUN_INFO", &
130 58 : extension=".scfLog")
131 232 : rlab = [CHARACTER(LEN=3) :: "X", "Y", "Z"]
132 178 : periodic = ANY(cell%perd > 0) ! periodic cell
133 58 : my_apply_pulse = .TRUE.
134 58 : CALL get_qs_env(qs_env, mos=mos)
135 :
136 58 : IF (rtp%linear_scaling) THEN
137 40 : IF (.NOT. ASSOCIATED(mos)) THEN
138 : CALL cp_warn(__LOCATION__, "Delta Pulse not implemented for Linear-Scaling based ground "// &
139 : "state calculation. If you want to perform a Linear-Scaling RTP from a "// &
140 : "Linear-Scaling GS calculation you can do the following: (i) LSCF froms "// &
141 : "scratch, (ii) MO-based SCF (for 1 SCF loop for instance) with the LSCF "// &
142 : "result as a restart and (iii) linear scaling RTP + delta kick (for 1 "// &
143 0 : "SCF loop for instance).")
144 : my_apply_pulse = .FALSE.
145 : ELSE
146 : ! create temporary mos_old and mos_new to use delta kick routine designed for MOs-based RTP
147 : CALL rt_prop_create_mos(rtp, mos, qs_env%mpools, dft_control, &
148 : init_mos_old=.TRUE., init_mos_new=.TRUE., &
149 40 : init_mos_next=.FALSE., init_mos_admn=.FALSE.)
150 : END IF
151 : END IF
152 :
153 : IF (my_apply_pulse) THEN
154 : ! The amplitude of the perturbation for all the method, modulo some prefactor:
155 : kvec(:) = cell%h_inv(1, :)*rtp_control%delta_pulse_direction(1) + &
156 : cell%h_inv(2, :)*rtp_control%delta_pulse_direction(2) + &
157 232 : cell%h_inv(3, :)*rtp_control%delta_pulse_direction(3)
158 232 : kvec = kvec*twopi*rtp_control%delta_pulse_scale
159 :
160 58 : CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
161 58 : IF (rtp_control%apply_delta_pulse) THEN
162 50 : IF (dft_control%qs_control%dftb) THEN
163 0 : CALL build_dftb_overlap(qs_env, 1, matrix_s)
164 : END IF
165 50 : IF (rtp_control%periodic) THEN
166 36 : IF (output_unit > 0) THEN
167 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T40))") &
168 : "An Electric Delta Kick within periodic condition is applied before running RTP. "// &
169 18 : "Its amplitude in atomic unit is:"
170 : WRITE (output_unit, "(T3,3(A,A,E16.8,1X))") &
171 72 : (TRIM(rlab(i)), "=", -kvec(i), i=1, 3)
172 : END IF
173 144 : CALL apply_delta_pulse_electric_periodic(qs_env, mos_old, mos_new, -kvec)
174 : ELSE
175 14 : CPWARN_IF(periodic, "This application of the delta pulse is not compatible with PBC!")
176 14 : IF (output_unit > 0) THEN
177 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T40))") &
178 : "An Electric Delta Kick within the length gauge is applied before running RTP. "// &
179 7 : "Its amplitude in atomic unit is:"
180 : WRITE (output_unit, "(T3,3(A,A,E16.8,1X))") &
181 28 : (TRIM(rlab(i)), "=", -kvec(i), i=1, 3)
182 : END IF
183 56 : CALL apply_delta_pulse_electric(qs_env, mos_old, mos_new, -kvec)
184 : END IF
185 8 : ELSE IF (rtp_control%apply_delta_pulse_mag) THEN
186 8 : CPWARN_IF(periodic, "This application of the delta pulse is not compatible with PBC!")
187 : ! The prefactor (strength of the magnetic field, should be divided by 2c)
188 8 : IF (output_unit > 0) THEN
189 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T40))") &
190 : "A Magnetic Delta Kick is applied before running RTP. "// &
191 4 : "Its amplitude in atomic unit is:"
192 : WRITE (output_unit, "(T3,3(A,A,E16.8,1X))") &
193 16 : (TRIM(rlab(i)), "=", -kvec(i)/2, i=1, 3)
194 : END IF
195 32 : CALL apply_delta_pulse_mag(qs_env, mos_old, mos_new, -kvec(:)/2)
196 : ELSE
197 0 : CPABORT("Code error: this case should not happen!")
198 : END IF
199 : END IF
200 :
201 58 : END SUBROUTINE apply_delta_pulse
202 :
203 : ! **************************************************************************************************
204 : !> \brief uses perturbation theory to get the proper initial conditions
205 : !> The len_rep option is NOT compatible with periodic boundary conditions!
206 : !> \param qs_env ...
207 : !> \param mos_old ...
208 : !> \param mos_new ...
209 : !> \param kvec ...
210 : !> \author Joost & Martin (2011)
211 : ! **************************************************************************************************
212 :
213 180 : SUBROUTINE apply_delta_pulse_electric_periodic(qs_env, mos_old, mos_new, kvec)
214 : TYPE(qs_environment_type), POINTER :: qs_env
215 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old, mos_new
216 : REAL(KIND=dp), DIMENSION(3) :: kvec
217 :
218 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_electric_periodic'
219 :
220 : INTEGER :: handle, icol, idir, irow, ispin, nao, &
221 : ncol_local, nmo, nrow_local, nvirt, &
222 : reference
223 36 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
224 : LOGICAL :: com_nl, len_rep, periodic
225 : REAL(KIND=dp) :: eps_ppnl, factor
226 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
227 36 : POINTER :: local_data
228 : REAL(KIND=dp), DIMENSION(3) :: rcc
229 36 : REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, ref_point
230 : TYPE(cell_type), POINTER :: cell
231 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_tmp
232 : TYPE(cp_fm_type) :: eigenvectors, mat_ks, mat_tmp, momentum, &
233 : S_chol, virtuals
234 36 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_r, matrix_rv, matrix_s
235 : TYPE(dft_control_type), POINTER :: dft_control
236 36 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
237 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
238 36 : POINTER :: sab_orb, sap_ppnl
239 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
240 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
241 : TYPE(rt_prop_type), POINTER :: rtp
242 : TYPE(rtp_control_type), POINTER :: rtp_control
243 : TYPE(section_vals_type), POINTER :: input
244 :
245 36 : CALL timeset(routineN, handle)
246 :
247 36 : NULLIFY (cell, mos, rtp, matrix_s, matrix_ks, input, dft_control, particle_set, fm_struct)
248 : ! we need the overlap and ks matrix for a full diagonalization
249 : CALL get_qs_env(qs_env, &
250 : cell=cell, &
251 : mos=mos, &
252 : rtp=rtp, &
253 : matrix_s=matrix_s, &
254 : matrix_ks=matrix_ks, &
255 : dft_control=dft_control, &
256 : input=input, &
257 36 : particle_set=particle_set)
258 :
259 36 : rtp_control => dft_control%rtp_control
260 102 : periodic = ANY(cell%perd > 0) ! periodic cell
261 :
262 : ! relevant input parameters
263 36 : com_nl = section_get_lval(section_vals=input, keyword_name="DFT%REAL_TIME_PROPAGATION%COM_NL")
264 36 : len_rep = section_get_lval(section_vals=input, keyword_name="DFT%REAL_TIME_PROPAGATION%LEN_REP")
265 :
266 : ! calculate non-local commutator if necessary
267 36 : IF (com_nl) THEN
268 20 : CALL cite_reference(Mattiat2019)
269 20 : NULLIFY (qs_kind_set, sab_orb, sap_ppnl)
270 : CALL get_qs_env(qs_env, &
271 : sap_ppnl=sap_ppnl, &
272 : sab_orb=sab_orb, &
273 20 : qs_kind_set=qs_kind_set)
274 20 : eps_ppnl = dft_control%qs_control%eps_ppnl
275 :
276 20 : NULLIFY (matrix_rv)
277 20 : CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
278 80 : DO idir = 1, 3
279 60 : CALL dbcsr_init_p(matrix_rv(idir)%matrix)
280 : CALL dbcsr_create(matrix_rv(idir)%matrix, template=matrix_s(1)%matrix, &
281 60 : matrix_type=dbcsr_type_antisymmetric)
282 60 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(idir)%matrix, sab_orb)
283 80 : CALL dbcsr_set(matrix_rv(idir)%matrix, 0._dp)
284 : END DO
285 20 : CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv)
286 : END IF
287 :
288 : ! calculate dipole moment matrix if required, NOT for periodic boundary conditions!
289 36 : IF (len_rep) THEN
290 10 : CALL cite_reference(Mattiat2022)
291 10 : CPWARN_IF(periodic, "This application of the delta pulse is not compatible with PBC!")
292 : ! get reference point
293 : reference = section_get_ival(section_vals=input, &
294 10 : keyword_name="DFT%PRINT%MOMENTS%REFERENCE")
295 10 : NULLIFY (ref_point)
296 10 : CALL section_vals_val_get(input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
297 10 : CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
298 :
299 10 : NULLIFY (sab_orb)
300 10 : CALL get_qs_env(qs_env, sab_orb=sab_orb)
301 : ! calculate dipole moment operator
302 10 : NULLIFY (matrix_r)
303 10 : CALL dbcsr_allocate_matrix_set(matrix_r, 3)
304 40 : DO idir = 1, 3
305 30 : CALL dbcsr_init_p(matrix_r(idir)%matrix)
306 30 : CALL dbcsr_create(matrix_r(idir)%matrix, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
307 30 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_r(idir)%matrix, sab_orb)
308 40 : CALL dbcsr_set(matrix_r(idir)%matrix, 0._dp)
309 : END DO
310 10 : CALL build_local_moment_matrix(qs_env, matrix_r, 1, rcc)
311 : END IF
312 :
313 36 : IF (rtp_control%velocity_gauge) THEN
314 0 : rtp_control%vec_pot = rtp_control%vec_pot + kvec
315 : END IF
316 :
317 : ! struct for fm matrices
318 36 : fm_struct => rtp%ao_ao_fmstruct
319 :
320 : ! create matrices and get Cholesky decomposition of S
321 36 : CALL cp_fm_create(mat_ks, matrix_struct=fm_struct, name="mat_ks")
322 36 : CALL cp_fm_create(eigenvectors, matrix_struct=fm_struct, name="eigenvectors")
323 36 : CALL cp_fm_create(S_chol, matrix_struct=fm_struct, name="S_chol")
324 36 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_chol)
325 36 : CALL cp_fm_cholesky_decompose(S_chol)
326 :
327 : ! get number of atomic orbitals
328 36 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
329 :
330 88 : DO ispin = 1, SIZE(matrix_ks)
331 : ! diagonalize KS matrix to get occ and virt mos
332 156 : ALLOCATE (eigenvalues(nao))
333 52 : CALL cp_fm_create(mat_tmp, matrix_struct=fm_struct, name="mat_tmp")
334 52 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, mat_ks)
335 52 : CALL cp_fm_cholesky_reduce(mat_ks, S_chol)
336 52 : CALL cp_fm_syevd(mat_ks, mat_tmp, eigenvalues)
337 52 : CALL cp_fm_cholesky_restore(mat_tmp, nao, S_chol, eigenvectors, "SOLVE")
338 :
339 : ! virtuals
340 52 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
341 52 : nvirt = nao - nmo
342 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=fm_struct%para_env, context=fm_struct%context, &
343 52 : nrow_global=nao, ncol_global=nvirt)
344 52 : CALL cp_fm_create(virtuals, matrix_struct=fm_struct_tmp, name="virtuals")
345 52 : CALL cp_fm_struct_release(fm_struct_tmp)
346 52 : CALL cp_fm_to_fm(eigenvectors, virtuals, nvirt, nmo + 1, 1)
347 :
348 : ! occupied
349 52 : CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1), nmo, 1, 1)
350 :
351 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=fm_struct%para_env, context=fm_struct%context, &
352 52 : nrow_global=nvirt, ncol_global=nmo)
353 52 : CALL cp_fm_create(momentum, matrix_struct=fm_struct_tmp, name="momentum")
354 52 : CALL cp_fm_struct_release(fm_struct_tmp)
355 :
356 : ! the momentum operator (in a given direction)
357 52 : CALL cp_fm_set_all(mos_new(2*ispin - 1), 0.0_dp)
358 :
359 208 : DO idir = 1, 3
360 156 : factor = kvec(idir)
361 208 : IF (factor /= 0.0_dp) THEN
362 56 : IF (.NOT. len_rep) THEN
363 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(idir + 1)%matrix, mos_old(2*ispin - 1), &
364 42 : mos_old(2*ispin), ncol=nmo)
365 : ELSE
366 : CALL cp_dbcsr_sm_fm_multiply(matrix_r(idir)%matrix, mos_old(2*ispin - 1), &
367 14 : mos_old(2*ispin), ncol=nmo)
368 : END IF
369 :
370 56 : CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1), factor, mos_old(2*ispin))
371 56 : IF (com_nl) THEN
372 28 : CALL cp_fm_set_all(mos_old(2*ispin), 0.0_dp)
373 : CALL cp_dbcsr_sm_fm_multiply(matrix_rv(idir)%matrix, mos_old(2*ispin - 1), &
374 28 : mos_old(2*ispin), ncol=nmo)
375 28 : CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1), factor, mos_old(2*ispin))
376 : END IF
377 : END IF
378 : END DO
379 :
380 52 : CALL parallel_gemm('T', 'N', nvirt, nmo, nao, 1.0_dp, virtuals, mos_new(2*ispin - 1), 0.0_dp, momentum)
381 :
382 : ! the tricky bit ... rescale by the eigenvalue difference
383 52 : IF (.NOT. len_rep) THEN
384 : CALL cp_fm_get_info(momentum, nrow_local=nrow_local, ncol_local=ncol_local, &
385 38 : row_indices=row_indices, col_indices=col_indices, local_data=local_data)
386 372 : DO icol = 1, ncol_local
387 6398 : DO irow = 1, nrow_local
388 6026 : factor = 1/(eigenvalues(col_indices(icol)) - eigenvalues(nmo + row_indices(irow)))
389 6360 : local_data(irow, icol) = factor*local_data(irow, icol)
390 : END DO
391 : END DO
392 : END IF
393 52 : CALL cp_fm_release(mat_tmp)
394 52 : DEALLOCATE (eigenvalues)
395 :
396 : ! now obtain the initial condition in mos_old
397 52 : CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1), nmo, 1, 1)
398 52 : CALL parallel_gemm("N", "N", nao, nmo, nvirt, 1.0_dp, virtuals, momentum, 0.0_dp, mos_old(2*ispin))
399 :
400 52 : CALL cp_fm_release(virtuals)
401 296 : CALL cp_fm_release(momentum)
402 : END DO
403 :
404 : ! release matrices
405 36 : CALL cp_fm_release(S_chol)
406 36 : CALL cp_fm_release(mat_ks)
407 36 : CALL cp_fm_release(eigenvectors)
408 36 : IF (com_nl) CALL dbcsr_deallocate_matrix_set(matrix_rv)
409 36 : IF (len_rep) CALL dbcsr_deallocate_matrix_set(matrix_r)
410 :
411 : ! orthonormalize afterwards
412 36 : CALL orthonormalize_complex_mos(qs_env, mos_old)
413 :
414 36 : CALL timestop(handle)
415 :
416 36 : END SUBROUTINE apply_delta_pulse_electric_periodic
417 :
418 : ! **************************************************************************************************
419 : !> \brief applies exp(ikr) to the wavefunction.... stored in mos_old...
420 : !> \param qs_env ...
421 : !> \param mos_old ...
422 : !> \param mos_new ...
423 : !> \param kvec ...
424 : !> \author Joost & Martin (2011)
425 : ! **************************************************************************************************
426 :
427 42 : SUBROUTINE apply_delta_pulse_electric(qs_env, mos_old, mos_new, kvec)
428 : TYPE(qs_environment_type), POINTER :: qs_env
429 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old, mos_new
430 : REAL(KIND=dp), DIMENSION(3) :: kvec
431 :
432 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_electric'
433 :
434 : INTEGER :: handle, i, nao, nmo
435 : TYPE(cell_type), POINTER :: cell
436 : TYPE(cp_fm_type) :: S_inv_fm, tmp
437 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
438 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
439 : TYPE(dft_control_type), POINTER :: dft_control
440 14 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
441 : TYPE(rt_prop_type), POINTER :: rtp
442 : TYPE(rtp_control_type), POINTER :: rtp_control
443 :
444 14 : CALL timeset(routineN, handle)
445 14 : NULLIFY (cell, dft_control, matrix_s, mos, rtp, rtp_control)
446 : CALL get_qs_env(qs_env, &
447 : cell=cell, &
448 : dft_control=dft_control, &
449 : matrix_s=matrix_s, &
450 : mos=mos, &
451 14 : rtp=rtp)
452 14 : rtp_control => dft_control%rtp_control
453 :
454 14 : IF (rtp_control%velocity_gauge) THEN
455 0 : rtp_control%vec_pot = rtp_control%vec_pot + kvec
456 : END IF
457 :
458 : ! calculate exponentials (= Berry moments)
459 14 : NULLIFY (cosmat, sinmat)
460 14 : ALLOCATE (cosmat, sinmat)
461 14 : CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
462 14 : CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
463 14 : CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
464 :
465 : ! need inverse of overlap matrix
466 14 : CALL cp_fm_create(S_inv_fm, matrix_struct=rtp%ao_ao_fmstruct, name="S_inv_fm")
467 14 : CALL cp_fm_create(tmp, matrix_struct=rtp%ao_ao_fmstruct, name="tmp_mat")
468 14 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_inv_fm)
469 14 : CALL cp_fm_cholesky_decompose(S_inv_fm)
470 14 : CALL cp_fm_cholesky_invert(S_inv_fm)
471 14 : CALL cp_fm_uplo_to_full(S_inv_fm, tmp)
472 14 : CALL cp_fm_release(tmp)
473 :
474 38 : DO i = 1, SIZE(mos)
475 : ! apply exponentials to mo coefficients
476 24 : CALL get_mo_set(mos(i), nao=nao, nmo=nmo)
477 24 : CALL cp_dbcsr_sm_fm_multiply(cosmat, mos(i)%mo_coeff, mos_new(2*i - 1), ncol=nmo)
478 24 : CALL cp_dbcsr_sm_fm_multiply(sinmat, mos(i)%mo_coeff, mos_new(2*i), ncol=nmo)
479 :
480 24 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_inv_fm, mos_new(2*i - 1), 0.0_dp, mos_old(2*i - 1))
481 62 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_inv_fm, mos_new(2*i), 0.0_dp, mos_old(2*i))
482 : END DO
483 :
484 14 : CALL cp_fm_release(S_inv_fm)
485 14 : CALL dbcsr_deallocate_matrix(cosmat)
486 14 : CALL dbcsr_deallocate_matrix(sinmat)
487 :
488 : ! orthonormalize afterwards
489 14 : CALL orthonormalize_complex_mos(qs_env, mos_old)
490 :
491 14 : CALL timestop(handle)
492 :
493 14 : END SUBROUTINE apply_delta_pulse_electric
494 :
495 : ! **************************************************************************************************
496 : !> \brief apply magnetic delta pulse to linear order
497 : !> \param qs_env ...
498 : !> \param mos_old ...
499 : !> \param mos_new ...
500 : !> \param kvec ...
501 : ! **************************************************************************************************
502 40 : SUBROUTINE apply_delta_pulse_mag(qs_env, mos_old, mos_new, kvec)
503 : TYPE(qs_environment_type), POINTER :: qs_env
504 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_old, mos_new
505 : REAL(KIND=dp), DIMENSION(3) :: kvec
506 :
507 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_mag'
508 :
509 : INTEGER :: gauge_orig, handle, idir, ispin, nao, &
510 : nmo, nrow_global, nvirt
511 : REAL(KIND=dp) :: eps_ppnl, factor
512 : REAL(KIND=dp), DIMENSION(3) :: rcc
513 8 : REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, ref_point
514 : TYPE(cell_type), POINTER :: cell
515 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
516 : TYPE(cp_fm_type) :: eigenvectors, mat_ks, perturbation, &
517 : S_chol, virtuals
518 8 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_mag, matrix_nl, &
519 8 : matrix_s
520 : TYPE(dft_control_type), POINTER :: dft_control
521 8 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
522 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
523 8 : POINTER :: sab_all, sab_orb, sap_ppnl
524 8 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
525 8 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
526 : TYPE(rt_prop_type), POINTER :: rtp
527 : TYPE(section_vals_type), POINTER :: input
528 :
529 8 : CALL timeset(routineN, handle)
530 :
531 8 : CALL cite_reference(Mattiat2022)
532 :
533 8 : NULLIFY (rtp, dft_control, matrix_ks, matrix_s, input, mos, cell, sab_orb, sab_all, sap_ppnl, &
534 8 : qs_kind_set, particle_set)
535 :
536 : CALL get_qs_env(qs_env, &
537 : rtp=rtp, &
538 : dft_control=dft_control, &
539 : mos=mos, &
540 : matrix_ks=matrix_ks, &
541 : matrix_s=matrix_s, &
542 : input=input, &
543 : cell=cell, &
544 : sab_orb=sab_orb, &
545 : sab_all=sab_all, &
546 8 : sap_ppnl=sap_ppnl)
547 :
548 : gauge_orig = section_get_ival(section_vals=input, &
549 8 : keyword_name="DFT%REAL_TIME_PROPAGATION%GAUGE_ORIG")
550 8 : NULLIFY (ref_point)
551 8 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%GAUGE_ORIG_MANUAL", r_vals=ref_point)
552 8 : CALL get_reference_point(rcc, qs_env=qs_env, reference=gauge_orig, ref_point=ref_point)
553 :
554 : ! Create fm matrices
555 8 : CALL cp_fm_create(S_chol, matrix_struct=rtp%ao_ao_fmstruct, name='Cholesky S')
556 8 : CALL cp_fm_create(eigenvectors, matrix_struct=rtp%ao_ao_fmstruct, name="gs evecs fm")
557 8 : CALL cp_fm_create(mat_ks, matrix_struct=rtp%ao_ao_fmstruct, name='KS matrix')
558 :
559 : ! get nrows_global
560 8 : CALL cp_fm_get_info(mat_ks, nrow_global=nrow_global)
561 :
562 : ! cholesky decomposition of overlap matrix
563 8 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_chol)
564 8 : CALL cp_fm_cholesky_decompose(S_chol)
565 :
566 : ! initiate perturbation matrix
567 8 : NULLIFY (matrix_mag)
568 8 : CALL dbcsr_allocate_matrix_set(matrix_mag, 3)
569 32 : DO idir = 1, 3
570 24 : CALL dbcsr_init_p(matrix_mag(idir)%matrix)
571 : CALL dbcsr_create(matrix_mag(idir)%matrix, template=matrix_s(1)%matrix, &
572 24 : matrix_type=dbcsr_type_antisymmetric)
573 24 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_mag(idir)%matrix, sab_orb)
574 32 : CALL dbcsr_set(matrix_mag(idir)%matrix, 0._dp)
575 : END DO
576 : ! construct magnetic dipole moment matrix
577 8 : CALL build_local_magmom_matrix(qs_env, matrix_mag, 1, ref_point=rcc)
578 :
579 : ! work matrix for non-local potential part if necessary
580 8 : NULLIFY (matrix_nl)
581 8 : IF (ASSOCIATED(sap_ppnl)) THEN
582 8 : CALL dbcsr_allocate_matrix_set(matrix_nl, 3)
583 32 : DO idir = 1, 3
584 24 : CALL dbcsr_init_p(matrix_nl(idir)%matrix)
585 : CALL dbcsr_create(matrix_nl(idir)%matrix, template=matrix_s(1)%matrix, &
586 24 : matrix_type=dbcsr_type_antisymmetric)
587 24 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_nl(idir)%matrix, sab_orb)
588 32 : CALL dbcsr_set(matrix_nl(idir)%matrix, 0._dp)
589 : END DO
590 : ! construct non-local contribution
591 : CALL get_qs_env(qs_env, &
592 : qs_kind_set=qs_kind_set, &
593 8 : particle_set=particle_set)
594 8 : eps_ppnl = dft_control%qs_control%eps_ppnl
595 :
596 8 : CALL build_com_nl_mag(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_nl, rcc, cell)
597 :
598 32 : DO idir = 1, 3
599 32 : CALL dbcsr_add(matrix_mag(idir)%matrix, matrix_nl(idir)%matrix, -one, one)
600 : END DO
601 :
602 8 : CALL dbcsr_deallocate_matrix_set(matrix_nl)
603 : END IF
604 :
605 20 : DO ispin = 1, dft_control%nspins
606 : ! allocate eigenvalues
607 : NULLIFY (eigenvalues)
608 36 : ALLOCATE (eigenvalues(nrow_global))
609 : ! diagonalize KS matrix in AO basis using Cholesky decomp. of S
610 12 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, mat_ks)
611 12 : CALL cp_fm_cholesky_reduce(mat_ks, S_chol)
612 12 : CALL cp_fm_syevd(mat_ks, eigenvectors, eigenvalues)
613 12 : CALL cp_fm_triangular_multiply(S_chol, eigenvectors, invert_tr=.TRUE.)
614 :
615 : ! virtuals
616 12 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
617 12 : nvirt = nao - nmo
618 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=rtp%ao_ao_fmstruct%para_env, context=rtp%ao_ao_fmstruct%context, &
619 12 : nrow_global=nrow_global, ncol_global=nvirt)
620 12 : CALL cp_fm_create(virtuals, matrix_struct=fm_struct_tmp, name="virtuals")
621 12 : CALL cp_fm_struct_release(fm_struct_tmp)
622 12 : CALL cp_fm_to_fm(eigenvectors, virtuals, nvirt, nmo + 1, 1)
623 :
624 : ! occupied
625 12 : CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1), nmo, 1, 1)
626 :
627 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=rtp%ao_ao_fmstruct%para_env, context=rtp%ao_ao_fmstruct%context, &
628 12 : nrow_global=nvirt, ncol_global=nmo)
629 12 : CALL cp_fm_create(perturbation, matrix_struct=fm_struct_tmp, name="perturbation")
630 12 : CALL cp_fm_struct_release(fm_struct_tmp)
631 :
632 : ! apply perturbation
633 12 : CALL cp_fm_set_all(mos_new(2*ispin - 1), 0.0_dp)
634 :
635 48 : DO idir = 1, 3
636 36 : factor = kvec(idir)
637 48 : IF (factor /= 0.0_dp) THEN
638 : CALL cp_dbcsr_sm_fm_multiply(matrix_mag(idir)%matrix, mos_old(2*ispin - 1), &
639 12 : mos_old(2*ispin), ncol=nmo)
640 12 : CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1), factor, mos_old(2*ispin))
641 : END IF
642 : END DO
643 :
644 12 : CALL parallel_gemm('T', 'N', nvirt, nmo, nao, 1.0_dp, virtuals, mos_new(2*ispin - 1), 0.0_dp, perturbation)
645 :
646 12 : DEALLOCATE (eigenvalues)
647 :
648 : ! now obtain the initial condition in mos_old
649 12 : CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1), nmo, 1, 1)
650 12 : CALL parallel_gemm("N", "N", nao, nmo, nvirt, 1.0_dp, virtuals, perturbation, 0.0_dp, mos_old(2*ispin))
651 :
652 12 : CALL cp_fm_release(virtuals)
653 56 : CALL cp_fm_release(perturbation)
654 : END DO
655 :
656 : ! deallocations
657 8 : CALL cp_fm_release(S_chol)
658 8 : CALL cp_fm_release(mat_ks)
659 8 : CALL cp_fm_release(eigenvectors)
660 8 : CALL dbcsr_deallocate_matrix_set(matrix_mag)
661 :
662 : ! orthonormalize afterwards
663 8 : CALL orthonormalize_complex_mos(qs_env, mos_old)
664 :
665 8 : CALL timestop(handle)
666 :
667 8 : END SUBROUTINE apply_delta_pulse_mag
668 :
669 : ! **************************************************************************************************
670 : !> \brief orthonormalize complex mos, e. g. after non-unitary transformations using Löwdin's algorithm
671 : !> \param qs_env ...
672 : !> \param coeffs ...
673 : ! **************************************************************************************************
674 58 : SUBROUTINE orthonormalize_complex_mos(qs_env, coeffs)
675 : TYPE(qs_environment_type), POINTER :: qs_env
676 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT), &
677 : POINTER :: coeffs
678 :
679 58 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_sqrt
680 : INTEGER :: im, ispin, j, nao, nmo, nspins, re
681 58 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
682 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
683 : TYPE(cp_cfm_type) :: oo_c, oo_v, oo_vt
684 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
685 : TYPE(cp_fm_type) :: oo_1, oo_2, S_fm, tmp
686 174 : TYPE(cp_fm_type), DIMENSION(2) :: coeffs_tmp
687 58 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
688 : TYPE(dft_control_type), POINTER :: dft_control
689 58 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
690 : TYPE(mp_para_env_type), POINTER :: para_env
691 :
692 58 : NULLIFY (para_env, blacs_env, dft_control, matrix_s, mos)
693 : CALL get_qs_env(qs_env, &
694 : blacs_env=blacs_env, &
695 : dft_control=dft_control, &
696 : matrix_s=matrix_s, &
697 : mos=mos, &
698 58 : para_env=para_env)
699 58 : nspins = dft_control%nspins
700 58 : CALL cp_fm_get_info(coeffs(1), nrow_global=nao)
701 :
702 : ! get overlap matrix
703 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nao, &
704 58 : context=blacs_env, para_env=para_env)
705 58 : CALL cp_fm_create(S_fm, matrix_struct=fm_struct_tmp, name="overlap fm")
706 58 : CALL cp_fm_struct_release(fm_struct_tmp)
707 : ! copy overlap matrix
708 58 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_fm)
709 :
710 146 : DO ispin = 1, nspins
711 88 : CALL get_mo_set(mos(ispin), nmo=nmo)
712 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
713 88 : nrow_global=nmo, ncol_global=nmo)
714 88 : CALL cp_fm_create(oo_1, matrix_struct=fm_struct_tmp, name="oo_1")
715 88 : CALL cp_fm_create(oo_2, matrix_struct=fm_struct_tmp, name="oo_2")
716 88 : CALL cp_fm_struct_release(fm_struct_tmp)
717 :
718 88 : CALL cp_fm_create(tmp, matrix_struct=coeffs(2*ispin - 1)%matrix_struct, name="tmp_mat")
719 : ! get the complex overlap matrix in MO basis
720 : ! x^T S x + y^T S y + i (-y^TS x+x^T S y)
721 88 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_fm, coeffs(2*ispin - 1), 0.0_dp, tmp)
722 88 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin - 1), tmp, 0.0_dp, oo_1)
723 88 : CALL parallel_gemm("T", "N", nmo, nmo, nao, -1.0_dp, coeffs(2*ispin), tmp, 0.0_dp, oo_2)
724 :
725 88 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_fm, coeffs(2*ispin), 0.0_dp, tmp)
726 88 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin), tmp, 1.0_dp, oo_1)
727 88 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin - 1), tmp, 1.0_dp, oo_2)
728 88 : CALL cp_fm_release(tmp)
729 :
730 : ! complex Löwdin
731 88 : CALL cp_cfm_create(oo_c, oo_1%matrix_struct)
732 88 : CALL cp_cfm_create(oo_v, oo_1%matrix_struct)
733 88 : CALL cp_cfm_create(oo_vt, oo_1%matrix_struct)
734 4455 : oo_c%local_data = CMPLX(oo_1%local_data, oo_2%local_data, KIND=dp)
735 :
736 264 : ALLOCATE (eigenvalues(nmo))
737 264 : ALLOCATE (eigenvalues_sqrt(nmo))
738 88 : CALL cp_cfm_heevd(oo_c, oo_v, eigenvalues)
739 622 : eigenvalues_sqrt(:) = CMPLX(one/SQRT(eigenvalues(:)), zero, dp)
740 88 : CALL cp_cfm_to_cfm(oo_v, oo_vt)
741 88 : CALL cp_cfm_column_scale(oo_v, eigenvalues_sqrt)
742 88 : DEALLOCATE (eigenvalues)
743 88 : DEALLOCATE (eigenvalues_sqrt)
744 : CALL parallel_gemm('N', 'C', nmo, nmo, nmo, (1.0_dp, 0.0_dp), &
745 88 : oo_v, oo_vt, (0.0_dp, 0.0_dp), oo_c)
746 4455 : oo_1%local_data = REAL(oo_c%local_data, KIND=dp)
747 4455 : oo_2%local_data = AIMAG(oo_c%local_data)
748 88 : CALL cp_cfm_release(oo_c)
749 88 : CALL cp_cfm_release(oo_v)
750 88 : CALL cp_cfm_release(oo_vt)
751 :
752 : ! transform coefficients accordingly
753 264 : DO j = 1, 2
754 264 : CALL cp_fm_create(coeffs_tmp(j), matrix_struct=coeffs(2*(ispin - 1) + j)%matrix_struct)
755 : END DO
756 :
757 : ! indices for coeffs_tmp
758 88 : re = 1
759 88 : im = 2
760 88 : CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin - 1), oo_1, zero, coeffs_tmp(re))
761 88 : CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin - 1), oo_2, zero, coeffs_tmp(im))
762 :
763 88 : CALL parallel_gemm("N", "N", nao, nmo, nmo, -one, coeffs(2*ispin), oo_2, zero, coeffs(2*ispin - 1))
764 88 : CALL cp_fm_scale_and_add(one, coeffs(2*ispin - 1), one, coeffs_tmp(re))
765 :
766 88 : CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin), oo_1, one, coeffs_tmp(im))
767 88 : CALL cp_fm_to_fm(coeffs_tmp(im), coeffs(2*ispin))
768 :
769 264 : DO j = 1, 2
770 264 : CALL cp_fm_release(coeffs_tmp(j))
771 : END DO
772 88 : CALL cp_fm_release(oo_1)
773 498 : CALL cp_fm_release(oo_2)
774 : END DO
775 58 : CALL cp_fm_release(S_fm)
776 :
777 116 : END SUBROUTINE orthonormalize_complex_mos
778 :
779 : END MODULE rt_delta_pulse
|
