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_fm_to_dbcsr,&
27 : cp_dbcsr_sm_fm_multiply
28 : USE cp_fm_struct, ONLY: cp_fm_struct_type
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 1218 : 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 1218 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
78 1218 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
79 1218 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat, rho_ao, rho_ao_im, rho_new, &
80 1218 : S_der, SinvB, SinvH, SinvH_imag
81 : TYPE(dbcsr_type), POINTER :: S_inv, tmp
82 : TYPE(dft_control_type), POINTER :: dft_control
83 1218 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
84 1218 : 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 1218 : CALL timeset(routineN, handle)
90 1218 : 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 1218 : dft_control=dft_control)
99 :
100 1218 : 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 1218 : SinvH=SinvH, SinvB=SinvB)
103 :
104 1218 : 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 1218 : ALLOCATE (tmp)
108 1218 : CALL dbcsr_create(tmp, template=SinvB(1)%matrix)
109 :
110 1218 : IF (rtp%linear_scaling) THEN
111 346 : 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 1218 : IF (rtp%propagate_complex_ks) CALL get_rtp(rtp=rtp, SinvH_imag=SinvH_imag)
118 :
119 2688 : DO ispin = 1, SIZE(SinvH)
120 1470 : re = 2*ispin - 1
121 1470 : im = 2*ispin
122 2688 : IF (rtp%linear_scaling) THEN
123 : CALL dbcsr_multiply("N", "N", one, SinvH(ispin)%matrix, rho_new(re)%matrix, zero, tmp, &
124 346 : filter_eps=rtp%filter_eps)
125 346 : IF (rtp%propagate_complex_ks) THEN
126 : CALL dbcsr_multiply("N", "N", one, SinvH_imag(ispin)%matrix, rho_new(im)%matrix, one, tmp, &
127 104 : 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 346 : filter_eps=rtp%filter_eps)
131 346 : 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 3384 : DO i = 1, SIZE(force)
144 17720 : force(i)%ehrenfest(:, :) = -force(i)%ehrenfest(:, :)
145 : END DO
146 :
147 1218 : CALL dbcsr_deallocate_matrix(tmp)
148 :
149 1218 : CALL timestop(handle)
150 :
151 2436 : 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 1470 : 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 1470 : REAL(dp), DIMENSION(:, :), POINTER :: block_values, block_values2
174 : TYPE(dbcsr_iterator_type) :: iter
175 :
176 5880 : DO i = 1, 3
177 : !Calculate the sum over the hadmard product
178 : !S_der part
179 :
180 4410 : CALL dbcsr_iterator_start(iter, tmp)
181 25542 : DO WHILE (dbcsr_iterator_blocks_left(iter))
182 21132 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
183 21132 : CALL dbcsr_get_block_p(S_der(i)%matrix, row_atom, col_atom, block_values2, found=found)
184 25542 : IF (found) THEN
185 19996 : ikind = kind_of(col_atom)
186 19996 : 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 559098 : 2.0_dp*SUM(block_values*block_values2)
191 : END IF
192 : END DO
193 4410 : CALL dbcsr_iterator_stop(iter)
194 :
195 : !C_mat part
196 :
197 4410 : CALL dbcsr_iterator_start(iter, rho_im)
198 19536 : DO WHILE (dbcsr_iterator_blocks_left(iter))
199 15126 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
200 15126 : CALL dbcsr_get_block_p(C_mat(i)%matrix, row_atom, col_atom, block_values2, found=found)
201 19536 : IF (found) THEN
202 9491 : ikind = kind_of(col_atom)
203 9491 : 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 272226 : 2.0_dp*SUM(block_values*block_values2)
208 : END IF
209 : END DO
210 14700 : CALL dbcsr_iterator_stop(iter)
211 : END DO
212 :
213 1470 : END SUBROUTINE compute_forces
214 :
215 : ! **************************************************************************************************
216 : !> \brief ...
217 : !> \param qs_env ...
218 : ! **************************************************************************************************
219 30 : SUBROUTINE rt_admm_force(qs_env)
220 : TYPE(qs_environment_type), POINTER :: qs_env
221 :
222 : TYPE(admm_type), POINTER :: admm_env
223 30 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos, mos_admm
224 30 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: KS_aux_im, KS_aux_re, matrix_s_aux_fit, &
225 30 : matrix_s_aux_fit_vs_orb
226 : TYPE(rt_prop_type), POINTER :: rtp
227 :
228 : CALL get_qs_env(qs_env, &
229 : admm_env=admm_env, &
230 30 : rtp=rtp)
231 : CALL get_admm_env(admm_env, matrix_ks_aux_fit=KS_aux_re, &
232 : matrix_ks_aux_fit_im=KS_aux_im, &
233 : matrix_s_aux_fit=matrix_s_aux_fit, &
234 30 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
235 :
236 30 : CALL get_rtp(rtp=rtp, mos_new=mos, admm_mos=mos_admm)
237 :
238 : ! currently only none option
239 : CALL rt_admm_forces_none(qs_env, admm_env, KS_aux_re, KS_aux_im, &
240 30 : matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, mos_admm, mos)
241 :
242 30 : END SUBROUTINE rt_admm_force
243 :
244 : ! **************************************************************************************************
245 : !> \brief ...
246 : !> \param qs_env ...
247 : !> \param admm_env ...
248 : !> \param KS_aux_re ...
249 : !> \param KS_aux_im ...
250 : !> \param matrix_s_aux_fit ...
251 : !> \param matrix_s_aux_fit_vs_orb ...
252 : !> \param mos_admm ...
253 : !> \param mos ...
254 : ! **************************************************************************************************
255 30 : SUBROUTINE rt_admm_forces_none(qs_env, admm_env, KS_aux_re, KS_aux_im, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, mos_admm, mos)
256 : TYPE(qs_environment_type), POINTER :: qs_env
257 : TYPE(admm_type), POINTER :: admm_env
258 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: KS_aux_re, KS_aux_im, matrix_s_aux_fit, &
259 : matrix_s_aux_fit_vs_orb
260 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_admm, mos
261 :
262 : INTEGER :: im, ispin, jspin, nao, natom, naux, nmo, &
263 : re
264 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force
265 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
266 : TYPE(cp_fm_struct_type), POINTER :: mstruct
267 210 : TYPE(cp_fm_type), DIMENSION(2) :: tmp_aux_aux, tmp_aux_mo, tmp_aux_mo1, &
268 90 : tmp_aux_nao
269 : TYPE(dbcsr_type), POINTER :: matrix_w_q, matrix_w_s
270 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
271 30 : POINTER :: sab_aux_fit_asymm, sab_aux_fit_vs_orb
272 30 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
273 : TYPE(qs_ks_env_type), POINTER :: ks_env
274 :
275 30 : NULLIFY (sab_aux_fit_asymm, sab_aux_fit_vs_orb, ks_env)
276 :
277 30 : CALL get_qs_env(qs_env, ks_env=ks_env)
278 : CALL get_admm_env(admm_env, sab_aux_fit_asymm=sab_aux_fit_asymm, &
279 30 : sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
280 :
281 30 : ALLOCATE (matrix_w_s)
282 : CALL dbcsr_create(matrix_w_s, template=matrix_s_aux_fit(1)%matrix, &
283 30 : name='W MATRIX AUX S', matrix_type=dbcsr_type_no_symmetry)
284 30 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s, sab_aux_fit_asymm)
285 :
286 30 : ALLOCATE (matrix_w_q)
287 : CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
288 30 : "W MATRIX AUX Q")
289 :
290 90 : DO jspin = 1, 2
291 60 : CALL cp_fm_create(tmp_aux_aux(jspin), admm_env%work_aux_aux%matrix_struct, name="taa")
292 90 : CALL cp_fm_create(tmp_aux_nao(jspin), admm_env%work_aux_orb%matrix_struct, name="tao")
293 : END DO
294 :
295 70 : DO ispin = 1, SIZE(KS_aux_re)
296 40 : re = 2*ispin - 1; im = 2*ispin
297 40 : naux = admm_env%nao_aux_fit; nmo = admm_env%nmo(ispin); nao = admm_env%nao_orb
298 :
299 40 : mstruct => admm_env%work_aux_nmo(ispin)%matrix_struct
300 120 : DO jspin = 1, 2
301 80 : CALL cp_fm_create(tmp_aux_mo(jspin), mstruct, name="tam")
302 120 : CALL cp_fm_create(tmp_aux_mo1(jspin), mstruct, name="tam")
303 : END DO
304 :
305 : ! First calculate H=KS_aux*C~, real part ends on work_aux_aux2, imaginary part ends at work_aux_aux3
306 40 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(re), tmp_aux_mo(1), nmo, 4.0_dp, 0.0_dp)
307 40 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(im), tmp_aux_mo(2), nmo, 4.0_dp, 0.0_dp)
308 40 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(im), tmp_aux_mo(1), nmo, -4.0_dp, 1.0_dp)
309 40 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(re), tmp_aux_mo(2), nmo, 4.0_dp, 1.0_dp)
310 :
311 : ! Next step compute S-1*H
312 40 : CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(1), 0.0_dp, tmp_aux_mo1(1))
313 40 : CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(2), 0.0_dp, tmp_aux_mo1(2))
314 :
315 : ! Here we go on with Ws=S-1*H * C^H (take care of sign of the imaginary part!!!)
316 :
317 : CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(1), mos(re), 0.0_dp, &
318 40 : tmp_aux_nao(1))
319 : CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(2), mos(im), 1.0_dp, &
320 40 : tmp_aux_nao(1))
321 : CALL parallel_gemm("N", "T", naux, nao, nmo, 1.0_dp, tmp_aux_mo1(1), mos(im), 0.0_dp, &
322 40 : tmp_aux_nao(2))
323 : CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(2), mos(re), 1.0_dp, &
324 40 : tmp_aux_nao(2))
325 :
326 : ! Let's do the final bit Wq=S-1*H * C^H * A^T
327 40 : CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(1), admm_env%A, 0.0_dp, tmp_aux_aux(1))
328 40 : CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(2), admm_env%A, 0.0_dp, tmp_aux_aux(2))
329 :
330 : ! *** copy to sparse matrix
331 40 : CALL copy_fm_to_dbcsr(tmp_aux_nao(1), matrix_w_q, keep_sparsity=.TRUE.)
332 :
333 : ! *** copy to sparse matrix
334 40 : CALL copy_fm_to_dbcsr(tmp_aux_aux(1), matrix_w_s, keep_sparsity=.TRUE.)
335 :
336 120 : DO jspin = 1, 2
337 80 : CALL cp_fm_release(tmp_aux_mo(jspin))
338 120 : CALL cp_fm_release(tmp_aux_mo1(jspin))
339 : END DO
340 :
341 : ! *** This can be done in one call w_total = w_alpha + w_beta
342 : ! allocate force vector
343 40 : CALL get_qs_env(qs_env=qs_env, natom=natom)
344 120 : ALLOCATE (admm_force(3, natom))
345 40 : admm_force = 0.0_dp
346 : CALL build_overlap_force(ks_env, admm_force, &
347 : basis_type_a="AUX_FIT", basis_type_b="AUX_FIT", &
348 40 : sab_nl=sab_aux_fit_asymm, matrix_p=matrix_w_s)
349 : CALL build_overlap_force(ks_env, admm_force, &
350 : basis_type_a="AUX_FIT", basis_type_b="ORB", &
351 40 : sab_nl=sab_aux_fit_vs_orb, matrix_p=matrix_w_q)
352 : ! add forces
353 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
354 40 : force=force)
355 40 : CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
356 70 : DEALLOCATE (admm_force)
357 :
358 : END DO
359 :
360 : ! *** Deallocated weighted density matrices
361 30 : CALL dbcsr_deallocate_matrix(matrix_w_s)
362 30 : CALL dbcsr_deallocate_matrix(matrix_w_q)
363 :
364 90 : DO jspin = 1, 2
365 60 : CALL cp_fm_release(tmp_aux_aux(jspin))
366 90 : CALL cp_fm_release(tmp_aux_nao(jspin))
367 : END DO
368 :
369 60 : END SUBROUTINE rt_admm_forces_none
370 :
371 : END MODULE rt_propagation_forces
|
