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