Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
21 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
22 : cp_dbcsr_sm_fm_multiply
23 : USE cp_fm_struct, ONLY: cp_fm_struct_type
24 : USE cp_fm_types, ONLY: cp_fm_create,&
25 : cp_fm_release,&
26 : cp_fm_type
27 : USE dbcsr_api, ONLY: &
28 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_block_p, &
29 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
30 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_type, &
31 : dbcsr_type_no_symmetry
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 1222 : 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 1222 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
78 1222 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
79 1222 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat, rho_ao, rho_ao_im, rho_new, &
80 1222 : S_der, SinvB, SinvH, SinvH_imag
81 : TYPE(dbcsr_type), POINTER :: S_inv, tmp
82 : TYPE(dft_control_type), POINTER :: dft_control
83 1222 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
84 1222 : 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 1222 : CALL timeset(routineN, handle)
90 1222 : 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 1222 : dft_control=dft_control)
99 :
100 1222 : 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 1222 : SinvH=SinvH, SinvB=SinvB)
103 :
104 1222 : 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 1222 : ALLOCATE (tmp)
108 1222 : CALL dbcsr_create(tmp, template=SinvB(1)%matrix)
109 :
110 1222 : IF (rtp%linear_scaling) THEN
111 352 : CALL get_rtp(rtp=rtp, rho_new=rho_new)
112 : ELSE
113 870 : 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 1222 : IF (rtp%propagate_complex_ks) CALL get_rtp(rtp=rtp, SinvH_imag=SinvH_imag)
118 :
119 2686 : DO ispin = 1, SIZE(SinvH)
120 1464 : re = 2*ispin - 1
121 1464 : im = 2*ispin
122 2686 : IF (rtp%linear_scaling) THEN
123 : CALL dbcsr_multiply("N", "N", one, SinvH(ispin)%matrix, rho_new(re)%matrix, zero, tmp, &
124 352 : filter_eps=rtp%filter_eps)
125 352 : IF (rtp%propagate_complex_ks) &
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 : CALL dbcsr_multiply("N", "N", one, SinvB(ispin)%matrix, rho_new(im)%matrix, one, tmp, &
129 352 : filter_eps=rtp%filter_eps)
130 352 : CALL compute_forces(force, tmp, S_der, rho_new(im)%matrix, C_mat, kind_of, atom_of_kind)
131 : ELSE
132 1112 : CALL dbcsr_multiply("N", "N", one, SinvH(ispin)%matrix, rho_ao(ispin)%matrix, zero, tmp)
133 1112 : IF (rtp%propagate_complex_ks) &
134 162 : CALL dbcsr_multiply("N", "N", one, SinvH_imag(ispin)%matrix, rho_ao_im(ispin)%matrix, one, tmp)
135 1112 : CALL dbcsr_multiply("N", "N", one, SinvB(ispin)%matrix, rho_ao_im(ispin)%matrix, one, tmp)
136 1112 : CALL compute_forces(force, tmp, S_der, rho_ao_im(ispin)%matrix, C_mat, kind_of, atom_of_kind)
137 : END IF
138 : END DO
139 :
140 : ! recall QS forces, at this point have the other sign.
141 3396 : DO i = 1, SIZE(force)
142 17780 : force(i)%ehrenfest(:, :) = -force(i)%ehrenfest(:, :)
143 : END DO
144 :
145 1222 : CALL dbcsr_deallocate_matrix(tmp)
146 :
147 1222 : CALL timestop(handle)
148 :
149 2444 : END SUBROUTINE
150 :
151 : ! **************************************************************************************************
152 : !> \brief ...
153 : !> \param force ...
154 : !> \param tmp ...
155 : !> \param S_der ...
156 : !> \param rho_im ...
157 : !> \param C_mat ...
158 : !> \param kind_of ...
159 : !> \param atom_of_kind ...
160 : ! **************************************************************************************************
161 1464 : SUBROUTINE compute_forces(force, tmp, S_der, rho_im, C_mat, kind_of, atom_of_kind)
162 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
163 : TYPE(dbcsr_type), POINTER :: tmp
164 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: S_der
165 : TYPE(dbcsr_type), POINTER :: rho_im
166 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: C_mat
167 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, atom_of_kind
168 :
169 : INTEGER :: col_atom, i, ikind, kind_atom, row_atom
170 : LOGICAL :: found
171 1464 : REAL(dp), DIMENSION(:), POINTER :: block_values, block_values2
172 : TYPE(dbcsr_iterator_type) :: iter
173 :
174 5856 : DO i = 1, 3
175 : !Calculate the sum over the hadmard product
176 : !S_der part
177 :
178 4392 : CALL dbcsr_iterator_start(iter, tmp)
179 25443 : DO WHILE (dbcsr_iterator_blocks_left(iter))
180 21051 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
181 21051 : CALL dbcsr_get_block_p(S_der(i)%matrix, row_atom, col_atom, block_values2, found=found)
182 25443 : IF (found) THEN
183 19897 : ikind = kind_of(col_atom)
184 19897 : kind_atom = atom_of_kind(col_atom)
185 : !The block_values are in a vector format,
186 : ! so the dot_product is the sum over all elements of the hamand product, that I need
187 : force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + &
188 460457 : 2.0_dp*DOT_PRODUCT(block_values, block_values2)
189 : END IF
190 : END DO
191 4392 : CALL dbcsr_iterator_stop(iter)
192 :
193 : !C_mat part
194 :
195 4392 : CALL dbcsr_iterator_start(iter, rho_im)
196 19485 : DO WHILE (dbcsr_iterator_blocks_left(iter))
197 15093 : CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
198 15093 : CALL dbcsr_get_block_p(C_mat(i)%matrix, row_atom, col_atom, block_values2, found=found)
199 19485 : IF (found) THEN
200 9437 : ikind = kind_of(col_atom)
201 9437 : kind_atom = atom_of_kind(col_atom)
202 : !The block_values are in a vector format, so the dot_product is
203 : ! the sum over all elements of the hamand product, that I need
204 : force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + &
205 222202 : 2.0_dp*DOT_PRODUCT(block_values, block_values2)
206 : END IF
207 : END DO
208 5856 : CALL dbcsr_iterator_stop(iter)
209 : END DO
210 :
211 1464 : END SUBROUTINE compute_forces
212 :
213 : ! **************************************************************************************************
214 : !> \brief ...
215 : !> \param qs_env ...
216 : ! **************************************************************************************************
217 20 : SUBROUTINE rt_admm_force(qs_env)
218 : TYPE(qs_environment_type), POINTER :: qs_env
219 :
220 : TYPE(admm_type), POINTER :: admm_env
221 20 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos, mos_admm
222 20 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: KS_aux_im, KS_aux_re, matrix_s_aux_fit, &
223 20 : matrix_s_aux_fit_vs_orb
224 : TYPE(rt_prop_type), POINTER :: rtp
225 :
226 : CALL get_qs_env(qs_env, &
227 : admm_env=admm_env, &
228 20 : rtp=rtp)
229 : CALL get_admm_env(admm_env, matrix_ks_aux_fit=KS_aux_re, &
230 : matrix_ks_aux_fit_im=KS_aux_im, &
231 : matrix_s_aux_fit=matrix_s_aux_fit, &
232 20 : matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
233 :
234 20 : CALL get_rtp(rtp=rtp, mos_new=mos, admm_mos=mos_admm)
235 :
236 : ! currently only none option
237 : CALL rt_admm_forces_none(qs_env, admm_env, KS_aux_re, KS_aux_im, &
238 20 : matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, mos_admm, mos)
239 :
240 20 : END SUBROUTINE rt_admm_force
241 :
242 : ! **************************************************************************************************
243 : !> \brief ...
244 : !> \param qs_env ...
245 : !> \param admm_env ...
246 : !> \param KS_aux_re ...
247 : !> \param KS_aux_im ...
248 : !> \param matrix_s_aux_fit ...
249 : !> \param matrix_s_aux_fit_vs_orb ...
250 : !> \param mos_admm ...
251 : !> \param mos ...
252 : ! **************************************************************************************************
253 20 : 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)
254 : TYPE(qs_environment_type), POINTER :: qs_env
255 : TYPE(admm_type), POINTER :: admm_env
256 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: KS_aux_re, KS_aux_im, matrix_s_aux_fit, &
257 : matrix_s_aux_fit_vs_orb
258 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_admm, mos
259 :
260 : INTEGER :: im, ispin, jspin, nao, natom, naux, nmo, &
261 : re
262 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: admm_force
263 20 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
264 : TYPE(cp_fm_struct_type), POINTER :: mstruct
265 140 : TYPE(cp_fm_type), DIMENSION(2) :: tmp_aux_aux, tmp_aux_mo, tmp_aux_mo1, &
266 60 : tmp_aux_nao
267 : TYPE(dbcsr_type), POINTER :: matrix_w_q, matrix_w_s
268 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
269 20 : POINTER :: sab_aux_fit_asymm, sab_aux_fit_vs_orb
270 20 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
271 : TYPE(qs_ks_env_type), POINTER :: ks_env
272 :
273 20 : NULLIFY (sab_aux_fit_asymm, sab_aux_fit_vs_orb, ks_env)
274 :
275 20 : CALL get_qs_env(qs_env, ks_env=ks_env)
276 : CALL get_admm_env(admm_env, sab_aux_fit_asymm=sab_aux_fit_asymm, &
277 20 : sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
278 :
279 20 : ALLOCATE (matrix_w_s)
280 : CALL dbcsr_create(matrix_w_s, template=matrix_s_aux_fit(1)%matrix, &
281 20 : name='W MATRIX AUX S', matrix_type=dbcsr_type_no_symmetry)
282 20 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s, sab_aux_fit_asymm)
283 :
284 20 : ALLOCATE (matrix_w_q)
285 : CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
286 20 : "W MATRIX AUX Q")
287 :
288 60 : DO jspin = 1, 2
289 40 : CALL cp_fm_create(tmp_aux_aux(jspin), admm_env%work_aux_aux%matrix_struct, name="taa")
290 60 : CALL cp_fm_create(tmp_aux_nao(jspin), admm_env%work_aux_orb%matrix_struct, name="tao")
291 : END DO
292 :
293 40 : DO ispin = 1, SIZE(KS_aux_re)
294 20 : re = 2*ispin - 1; im = 2*ispin
295 20 : naux = admm_env%nao_aux_fit; nmo = admm_env%nmo(ispin); nao = admm_env%nao_orb
296 :
297 20 : mstruct => admm_env%work_aux_nmo(ispin)%matrix_struct
298 60 : DO jspin = 1, 2
299 40 : CALL cp_fm_create(tmp_aux_mo(jspin), mstruct, name="tam")
300 60 : CALL cp_fm_create(tmp_aux_mo1(jspin), mstruct, name="tam")
301 : END DO
302 :
303 : ! First calculate H=KS_aux*C~, real part ends on work_aux_aux2, imaginary part ends at work_aux_aux3
304 20 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(re), tmp_aux_mo(re), nmo, 4.0_dp, 0.0_dp)
305 20 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_re(ispin)%matrix, mos_admm(im), tmp_aux_mo(im), nmo, 4.0_dp, 0.0_dp)
306 20 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(im), tmp_aux_mo(re), nmo, -4.0_dp, 1.0_dp)
307 20 : CALL cp_dbcsr_sm_fm_multiply(KS_aux_im(ispin)%matrix, mos_admm(re), tmp_aux_mo(im), nmo, 4.0_dp, 1.0_dp)
308 :
309 : ! Next step compute S-1*H
310 20 : CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(re), 0.0_dp, tmp_aux_mo1(re))
311 20 : CALL parallel_gemm('N', 'N', naux, nmo, naux, 1.0_dp, admm_env%S_inv, tmp_aux_mo(im), 0.0_dp, tmp_aux_mo1(im))
312 :
313 : ! Here we go on with Ws=S-1*H * C^H (take care of sign of the imaginary part!!!)
314 :
315 : CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(re), mos(re), 0.0_dp, &
316 20 : tmp_aux_nao(re))
317 : CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(im), mos(im), 1.0_dp, &
318 20 : tmp_aux_nao(re))
319 : CALL parallel_gemm("N", "T", naux, nao, nmo, 1.0_dp, tmp_aux_mo1(re), mos(im), 0.0_dp, &
320 20 : tmp_aux_nao(im))
321 : CALL parallel_gemm("N", "T", naux, nao, nmo, -1.0_dp, tmp_aux_mo1(im), mos(re), 1.0_dp, &
322 20 : tmp_aux_nao(im))
323 :
324 : ! Let's do the final bit Wq=S-1*H * C^H * A^T
325 20 : CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(re), admm_env%A, 0.0_dp, tmp_aux_aux(re))
326 20 : CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao(im), admm_env%A, 0.0_dp, tmp_aux_aux(im))
327 :
328 : ! *** copy to sparse matrix
329 20 : CALL copy_fm_to_dbcsr(tmp_aux_nao(re), matrix_w_q, keep_sparsity=.TRUE.)
330 :
331 : ! *** copy to sparse matrix
332 20 : CALL copy_fm_to_dbcsr(tmp_aux_aux(re), matrix_w_s, keep_sparsity=.TRUE.)
333 :
334 60 : DO jspin = 1, 2
335 40 : CALL cp_fm_release(tmp_aux_mo(jspin))
336 60 : CALL cp_fm_release(tmp_aux_mo1(jspin))
337 : END DO
338 :
339 : ! *** This can be done in one call w_total = w_alpha + w_beta
340 : ! allocate force vector
341 20 : CALL get_qs_env(qs_env=qs_env, natom=natom)
342 60 : ALLOCATE (admm_force(3, natom))
343 260 : admm_force = 0.0_dp
344 : CALL build_overlap_force(ks_env, admm_force, &
345 : basis_type_a="AUX_FIT", basis_type_b="AUX_FIT", &
346 20 : sab_nl=sab_aux_fit_asymm, matrix_p=matrix_w_s)
347 : CALL build_overlap_force(ks_env, admm_force, &
348 : basis_type_a="AUX_FIT", basis_type_b="ORB", &
349 20 : sab_nl=sab_aux_fit_vs_orb, matrix_p=matrix_w_q)
350 : ! add forces
351 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
352 20 : force=force)
353 20 : CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
354 20 : DEALLOCATE (admm_force)
355 :
356 : ! *** Deallocated weighted density matrices
357 20 : CALL dbcsr_deallocate_matrix(matrix_w_s)
358 40 : CALL dbcsr_deallocate_matrix(matrix_w_q)
359 : END DO
360 :
361 60 : DO jspin = 1, 2
362 40 : CALL cp_fm_release(tmp_aux_aux(jspin))
363 60 : CALL cp_fm_release(tmp_aux_nao(jspin))
364 : END DO
365 :
366 40 : END SUBROUTINE rt_admm_forces_none
367 :
368 : END MODULE rt_propagation_forces
|
