Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines needed for EMD
10 : !> \author Florian Schiffmann (02.09)
11 : ! **************************************************************************************************
12 :
13 : MODULE rt_propagation_forces
14 : USE admm_types, ONLY: admm_type,&
15 : get_admm_env
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind_set
18 : USE cp_control_types, ONLY: dft_control_type,&
19 : rtp_control_type
20 : USE cp_dbcsr_api, ONLY: &
21 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_block_p, &
22 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
23 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_type, &
24 : dbcsr_type_no_symmetry
25 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
26 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
27 : copy_fm_to_dbcsr,&
28 : cp_dbcsr_sm_fm_multiply
29 : USE cp_fm_types, ONLY: cp_fm_create,&
30 : cp_fm_release,&
31 : cp_fm_type
32 : USE kinds, ONLY: dp
33 : USE mathconstants, ONLY: one,&
34 : zero
35 : USE parallel_gemm_api, ONLY: parallel_gemm
36 : USE particle_types, ONLY: particle_type
37 : USE qs_environment_types, ONLY: get_qs_env,&
38 : qs_environment_type
39 : USE qs_force_types, ONLY: add_qs_force,&
40 : qs_force_type
41 : USE qs_ks_types, ONLY: qs_ks_env_type
42 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
43 : USE qs_overlap, ONLY: build_overlap_force
44 : USE qs_rho_types, ONLY: qs_rho_get,&
45 : qs_rho_type
46 : USE rt_propagation_types, ONLY: get_rtp,&
47 : rt_prop_type
48 : #include "./base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 : PRIVATE
52 :
53 : PUBLIC :: calc_c_mat_force, &
54 : rt_admm_force
55 :
56 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_forces'
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief calculates the three additional force contributions needed in EMD
62 : !> P_imag*C , P_imag*B*S^-1*S_der , P*S^-1*H*S_der
63 : !> \param qs_env ...
64 : !> \par History
65 : !> 02.2014 switched to dbcsr matrices [Samuel Andermatt]
66 : !> 10.2023 merge MO-based and all-atom into one routine [Guillaume Le Breton]
67 : !> \author Florian Schiffmann (02.09)
68 : ! **************************************************************************************************
69 :
70 1228 : SUBROUTINE calc_c_mat_force(qs_env)
71 : TYPE(qs_environment_type), POINTER :: qs_env
72 :
73 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_c_mat_force'
74 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
75 :
76 : INTEGER :: handle, i, im, ispin, re
77 1228 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
78 1228 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
79 1228 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat, rho_ao, rho_ao_im, rho_new, &
80 1228 : S_der, SinvB, SinvH, SinvH_imag
81 : TYPE(dbcsr_type), POINTER :: S_inv, tmp
82 : TYPE(dft_control_type), POINTER :: dft_control
83 1228 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
84 1228 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
85 : TYPE(qs_rho_type), POINTER :: rho
86 : TYPE(rt_prop_type), POINTER :: rtp
87 : TYPE(rtp_control_type), POINTER :: rtp_control
88 :
89 1228 : CALL timeset(routineN, handle)
90 1228 : NULLIFY (rtp, particle_set, atomic_kind_set, dft_control)
91 :
92 : CALL get_qs_env(qs_env, &
93 : rtp=rtp, &
94 : rho=rho, &
95 : particle_set=particle_set, &
96 : atomic_kind_set=atomic_kind_set, &
97 : force=force, &
98 1228 : dft_control=dft_control)
99 :
100 1228 : rtp_control => dft_control%rtp_control
101 : CALL get_rtp(rtp=rtp, C_mat=C_mat, S_der=S_der, S_inv=S_inv, &
102 1228 : SinvH=SinvH, SinvB=SinvB)
103 :
104 1228 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
105 :
106 : NULLIFY (tmp)
107 1228 : ALLOCATE (tmp)
108 1228 : CALL dbcsr_create(tmp, template=SinvB(1)%matrix)
109 :
110 1228 : IF (rtp%linear_scaling) THEN
111 356 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
112 : ELSE
113 872 : CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao, rho_ao_im=rho_ao_im)
114 : END IF
115 :
116 : ! If SinvH has an imaginary part (the minus sign is already in SinvH_imag)
117 1228 : IF (rtp%propagate_complex_ks) CALL get_rtp(rtp=rtp, SinvH_imag=SinvH_imag)
118 :
119 2708 : DO ispin = 1, SIZE(SinvH)
120 1480 : re = 2*ispin - 1
121 1480 : im = 2*ispin
122 2708 : IF (rtp%linear_scaling) THEN
123 : CALL dbcsr_multiply("N", "N", one, SinvH(ispin)%matrix, rho_new(re)%matrix, zero, tmp, &
124 356 : filter_eps=rtp%filter_eps)
125 356 : IF (rtp%propagate_complex_ks) THEN
126 : CALL dbcsr_multiply("N", "N", one, SinvH_imag(ispin)%matrix, rho_new(im)%matrix, one, tmp, &
127 114 : filter_eps=rtp%filter_eps)
128 : END IF
129 : CALL dbcsr_multiply("N", "N", one, SinvB(ispin)%matrix, rho_new(im)%matrix, one, tmp, &
130 356 : filter_eps=rtp%filter_eps)
131 356 : CALL compute_forces(force, tmp, S_der, rho_new(im)%matrix, C_mat, kind_of, atom_of_kind)
132 : ELSE
133 1124 : CALL dbcsr_multiply("N", "N", one, SinvH(ispin)%matrix, rho_ao(ispin)%matrix, zero, tmp)
134 1124 : IF (rtp%propagate_complex_ks) THEN
135 182 : CALL dbcsr_multiply("N", "N", one, SinvH_imag(ispin)%matrix, rho_ao_im(ispin)%matrix, one, tmp)
136 : END IF
137 1124 : CALL dbcsr_multiply("N", "N", one, SinvB(ispin)%matrix, rho_ao_im(ispin)%matrix, one, tmp)
138 1124 : CALL compute_forces(force, tmp, S_der, rho_ao_im(ispin)%matrix, C_mat, kind_of, atom_of_kind)
139 : END IF
140 : END DO
141 :
142 : ! recall QS forces, at this point have the other sign.
143 3414 : DO i = 1, SIZE(force)
144 17990 : force(i)%ehrenfest(:, :) = -force(i)%ehrenfest(:, :)
145 : END DO
146 :
147 1228 : CALL dbcsr_deallocate_matrix(tmp)
148 :
149 1228 : CALL timestop(handle)
150 :
151 2456 : END SUBROUTINE calc_c_mat_force
152 :
153 : ! **************************************************************************************************
154 : !> \brief ...
155 : !> \param force ...
156 : !> \param tmp ...
157 : !> \param S_der ...
158 : !> \param rho_im ...
159 : !> \param C_mat ...
160 : !> \param kind_of ...
161 : !> \param atom_of_kind ...
162 : ! **************************************************************************************************
163 1480 : SUBROUTINE compute_forces(force, tmp, S_der, rho_im, C_mat, kind_of, atom_of_kind)
164 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
165 : TYPE(dbcsr_type), POINTER :: tmp
166 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: S_der
167 : TYPE(dbcsr_type), POINTER :: rho_im
168 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat
169 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, atom_of_kind
170 :
171 : INTEGER :: col_atom, i, ikind, kind_atom, row_atom
172 : LOGICAL :: found
173 1480 : REAL(dp), DIMENSION(:, :), POINTER :: block_values, block_values2
174 : TYPE(dbcsr_iterator_type) :: iter
175 :
176 5920 : DO i = 1, 3
177 : !Calculate the sum over the hadmard product
178 : !S_der part
179 :
180 4440 : CALL dbcsr_iterator_start(iter, tmp)
181 26112 : DO WHILE (dbcsr_iterator_blocks_left(iter))
182 21672 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
183 21672 : CALL dbcsr_get_block_p(S_der(i)%matrix, row_atom, col_atom, block_values2, found=found)
184 26112 : IF (found) THEN
185 20536 : ikind = kind_of(col_atom)
186 20536 : kind_atom = atom_of_kind(col_atom)
187 : !The block_values are in a vector format,
188 : ! so the dot_product is the sum over all elements of the hamand product, that I need
189 : force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + &
190 595518 : 2.0_dp*SUM(block_values*block_values2)
191 : END IF
192 : END DO
193 4440 : CALL dbcsr_iterator_stop(iter)
194 :
195 : !C_mat part
196 :
197 4440 : CALL dbcsr_iterator_start(iter, rho_im)
198 19998 : DO WHILE (dbcsr_iterator_blocks_left(iter))
199 15558 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
200 15558 : CALL dbcsr_get_block_p(C_mat(i)%matrix, row_atom, col_atom, block_values2, found=found)
201 19998 : IF (found) THEN
202 10889 : ikind = kind_of(col_atom)
203 10889 : kind_atom = atom_of_kind(col_atom)
204 : !The block_values are in a vector format, so the dot_product is
205 : ! the sum over all elements of the hamand product, that I need
206 : force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + &
207 302376 : 2.0_dp*SUM(block_values*block_values2)
208 : END IF
209 : END DO
210 14800 : CALL dbcsr_iterator_stop(iter)
211 : END DO
212 :
213 1480 : END SUBROUTINE compute_forces
214 :
215 : ! **************************************************************************************************
216 : !> \brief Force contribution from the dependence of the ADMM projection on the nuclear positions
217 : !> \param qs_env ...
218 : !> \par History
219 : !> 10.2026 use the density matrix instead of the MO coefficients
220 : ! **************************************************************************************************
221 30 : SUBROUTINE rt_admm_force(qs_env)
222 : TYPE(qs_environment_type), POINTER :: qs_env
223 :
224 : TYPE(admm_type), POINTER :: admm_env
225 30 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: KS_aux_im, KS_aux_re, matrix_s_aux_fit, &
226 30 : matrix_s_aux_fit_vs_orb
227 : TYPE(rt_prop_type), POINTER :: rtp
228 :
229 : CALL get_qs_env(qs_env, &
230 : admm_env=admm_env, &
231 30 : rtp=rtp)
232 : CALL get_admm_env(admm_env, matrix_ks_aux_fit=KS_aux_re, &
233 : matrix_ks_aux_fit_im=KS_aux_im, &
234 : matrix_s_aux_fit=matrix_s_aux_fit, &
235 30 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
236 :
237 : ! currently only none option
238 : CALL rt_admm_forces_none(qs_env, admm_env, rtp, KS_aux_re, KS_aux_im, &
239 30 : matrix_s_aux_fit, matrix_s_aux_fit_vs_orb)
240 :
241 30 : END SUBROUTINE rt_admm_force
242 :
243 : ! **************************************************************************************************
244 : !> \brief ADMM projection force without purification.
245 : !> With P~ = A*P*A^T and A = S~^-1*Q, the weighted matrices are
246 : !> W_Q = -2*S~^-1*Re(K~*A*P) and W_S = -W_Q*A^T.
247 : !> The complex density matrix P is used instead of the MO coefficients, so that
248 : !> P is consistent with the auxiliary KS matrix K~ and includes the occupations.
249 : !> \param qs_env ...
250 : !> \param admm_env ...
251 : !> \param rtp ...
252 : !> \param KS_aux_re ...
253 : !> \param KS_aux_im ...
254 : !> \param matrix_s_aux_fit ...
255 : !> \param matrix_s_aux_fit_vs_orb ...
256 : ! **************************************************************************************************
257 30 : SUBROUTINE rt_admm_forces_none(qs_env, admm_env, rtp, KS_aux_re, KS_aux_im, matrix_s_aux_fit, &
258 : matrix_s_aux_fit_vs_orb)
259 : TYPE(qs_environment_type), POINTER :: qs_env
260 : TYPE(admm_type), POINTER :: admm_env
261 : TYPE(rt_prop_type), POINTER :: rtp
262 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: KS_aux_re, KS_aux_im, matrix_s_aux_fit, &
263 : matrix_s_aux_fit_vs_orb
264 :
265 : INTEGER :: ispin, nao, natom, naux
266 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force
267 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
268 : TYPE(cp_fm_type) :: tmp_aux_aux, tmp_aux_kap, tmp_aux_nao, &
269 : tmp_nao_nao
270 90 : TYPE(cp_fm_type), DIMENSION(2) :: tmp_aux_ap
271 30 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_im, rho_new
272 : TYPE(dbcsr_type), POINTER :: matrix_w_q, matrix_w_s, rho_im, rho_re
273 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
274 30 : POINTER :: sab_aux_fit_asymm, sab_aux_fit_vs_orb
275 30 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
276 : TYPE(qs_ks_env_type), POINTER :: ks_env
277 : TYPE(qs_rho_type), POINTER :: rho
278 :
279 30 : NULLIFY (sab_aux_fit_asymm, sab_aux_fit_vs_orb, ks_env, rho, rho_ao, rho_ao_im, rho_new)
280 :
281 : CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho, natom=natom, &
282 30 : atomic_kind_set=atomic_kind_set, force=force)
283 : CALL get_admm_env(admm_env, sab_aux_fit_asymm=sab_aux_fit_asymm, &
284 30 : sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
285 30 : IF (rtp%linear_scaling) THEN
286 0 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
287 : ELSE
288 30 : CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao, rho_ao_im=rho_ao_im)
289 : END IF
290 :
291 30 : naux = admm_env%nao_aux_fit; nao = admm_env%nao_orb
292 :
293 30 : ALLOCATE (matrix_w_s)
294 : CALL dbcsr_create(matrix_w_s, template=matrix_s_aux_fit(1)%matrix, &
295 30 : name='W MATRIX AUX S', matrix_type=dbcsr_type_no_symmetry)
296 30 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s, sab_aux_fit_asymm)
297 :
298 30 : ALLOCATE (matrix_w_q)
299 : CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
300 30 : "W MATRIX AUX Q")
301 :
302 30 : CALL cp_fm_create(tmp_nao_nao, admm_env%work_orb_orb%matrix_struct, name="too")
303 30 : CALL cp_fm_create(tmp_aux_ap(1), admm_env%work_aux_orb%matrix_struct, name="tao")
304 30 : CALL cp_fm_create(tmp_aux_ap(2), admm_env%work_aux_orb%matrix_struct, name="tao")
305 30 : CALL cp_fm_create(tmp_aux_kap, admm_env%work_aux_orb%matrix_struct, name="tao")
306 30 : CALL cp_fm_create(tmp_aux_nao, admm_env%work_aux_orb%matrix_struct, name="tao")
307 30 : CALL cp_fm_create(tmp_aux_aux, admm_env%work_aux_aux%matrix_struct, name="taa")
308 :
309 90 : ALLOCATE (admm_force(3, natom))
310 :
311 70 : DO ispin = 1, SIZE(KS_aux_re)
312 40 : IF (rtp%linear_scaling) THEN
313 0 : rho_re => rho_new(2*ispin - 1)%matrix
314 0 : rho_im => rho_new(2*ispin)%matrix
315 : ELSE
316 40 : rho_re => rho_ao(ispin)%matrix
317 40 : rho_im => rho_ao_im(ispin)%matrix
318 : END IF
319 :
320 : ! First calculate A*P, real part ends on tmp_aux_ap(1), imaginary part on tmp_aux_ap(2)
321 40 : CALL copy_dbcsr_to_fm(rho_re, tmp_nao_nao)
322 40 : CALL parallel_gemm('N', 'N', naux, nao, nao, 1.0_dp, admm_env%A, tmp_nao_nao, 0.0_dp, tmp_aux_ap(1))
323 40 : CALL copy_dbcsr_to_fm(rho_im, tmp_nao_nao)
324 40 : CALL parallel_gemm('N', 'N', naux, nao, nao, 1.0_dp, admm_env%A, tmp_nao_nao, 0.0_dp, tmp_aux_ap(2))
325 :
326 : ! Next step Re(K~*A*P) = K~_re*A*P_re - K~_im*A*P_im
327 40 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, tmp_aux_ap(1), tmp_aux_kap, nao, 1.0_dp, 0.0_dp)
328 40 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, tmp_aux_ap(2), tmp_aux_kap, nao, -1.0_dp, 1.0_dp)
329 :
330 : ! Wq=-2*S~^-1*Re(K~*A*P), minus sign due to force = -dE/dR
331 40 : CALL parallel_gemm('N', 'N', naux, nao, naux, -2.0_dp, admm_env%S_inv, tmp_aux_kap, 0.0_dp, tmp_aux_nao)
332 :
333 : ! Ws=-Wq*A^T
334 40 : CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao, admm_env%A, 0.0_dp, tmp_aux_aux)
335 :
336 : ! *** copy to sparse matrix
337 40 : CALL copy_fm_to_dbcsr(tmp_aux_nao, matrix_w_q, keep_sparsity=.TRUE.)
338 :
339 : ! *** copy to sparse matrix
340 40 : CALL copy_fm_to_dbcsr(tmp_aux_aux, matrix_w_s, keep_sparsity=.TRUE.)
341 :
342 40 : admm_force = 0.0_dp
343 : CALL build_overlap_force(ks_env, admm_force, &
344 : basis_type_a="AUX_FIT", basis_type_b="AUX_FIT", &
345 40 : sab_nl=sab_aux_fit_asymm, matrix_p=matrix_w_s)
346 : CALL build_overlap_force(ks_env, admm_force, &
347 : basis_type_a="AUX_FIT", basis_type_b="ORB", &
348 40 : sab_nl=sab_aux_fit_vs_orb, matrix_p=matrix_w_q)
349 : ! add forces
350 70 : CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
351 :
352 : END DO
353 :
354 30 : DEALLOCATE (admm_force)
355 :
356 : ! *** Deallocated weighted density matrices
357 30 : CALL dbcsr_deallocate_matrix(matrix_w_s)
358 30 : CALL dbcsr_deallocate_matrix(matrix_w_q)
359 :
360 30 : CALL cp_fm_release(tmp_nao_nao)
361 30 : CALL cp_fm_release(tmp_aux_ap(1))
362 30 : CALL cp_fm_release(tmp_aux_ap(2))
363 30 : CALL cp_fm_release(tmp_aux_kap)
364 30 : CALL cp_fm_release(tmp_aux_nao)
365 30 : CALL cp_fm_release(tmp_aux_aux)
366 :
367 120 : END SUBROUTINE rt_admm_forces_none
368 :
369 : END MODULE rt_propagation_forces
|
