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