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