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 Calculates forces for LRIGPW method
10 : !> lri : local resolution of the identity
11 : !> \par History
12 : !> created Dorothea Golze [03.2014]
13 : !> refactoring and distant pair approximation JGH [08.2017]
14 : !> \authors Dorothea Golze
15 : ! **************************************************************************************************
16 : MODULE lri_forces
17 :
18 : USE atomic_kind_types, ONLY: atomic_kind_type,&
19 : get_atomic_kind,&
20 : get_atomic_kind_set
21 : USE basis_set_types, ONLY: gto_basis_set_type
22 : USE dbcsr_api, ONLY: dbcsr_get_block_p,&
23 : dbcsr_p_type,&
24 : dbcsr_type
25 : USE kinds, ONLY: dp
26 : USE kpoint_types, ONLY: get_kpoint_info,&
27 : kpoint_type
28 : USE lri_environment_types, ONLY: &
29 : allocate_lri_force_components, deallocate_lri_force_components, lri_density_type, &
30 : lri_environment_type, lri_force_type, lri_int_type, lri_kind_type, lri_list_type, &
31 : lri_rhoab_type
32 : USE lri_integrals, ONLY: allocate_int_type,&
33 : deallocate_int_type,&
34 : dint_type,&
35 : lri_dint2
36 : USE qs_environment_types, ONLY: get_qs_env,&
37 : qs_environment_type
38 : USE qs_force_types, ONLY: qs_force_type
39 : USE qs_ks_types, ONLY: get_ks_env,&
40 : qs_ks_env_type
41 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
42 : neighbor_list_iterate,&
43 : neighbor_list_iterator_create,&
44 : neighbor_list_iterator_p_type,&
45 : neighbor_list_iterator_release,&
46 : neighbor_list_set_p_type
47 : USE qs_o3c_types, ONLY: get_o3c_iterator_info,&
48 : o3c_iterate,&
49 : o3c_iterator_create,&
50 : o3c_iterator_release,&
51 : o3c_iterator_type
52 : USE virial_methods, ONLY: virial_pair_force
53 : USE virial_types, ONLY: virial_type
54 :
55 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : ! **************************************************************************************************
63 :
64 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_forces'
65 :
66 : PUBLIC :: calculate_lri_forces, calculate_ri_forces
67 :
68 : ! **************************************************************************************************
69 :
70 : CONTAINS
71 :
72 : ! **************************************************************************************************
73 : !> \brief calculates the lri forces
74 : !> \param lri_env ...
75 : !> \param lri_density ...
76 : !> \param qs_env ...
77 : !> \param pmatrix density matrix
78 : !> \param atomic_kind_set ...
79 : ! **************************************************************************************************
80 26 : SUBROUTINE calculate_lri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
81 :
82 : TYPE(lri_environment_type), POINTER :: lri_env
83 : TYPE(lri_density_type), POINTER :: lri_density
84 : TYPE(qs_environment_type), POINTER :: qs_env
85 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
86 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
87 :
88 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_lri_forces'
89 :
90 : INTEGER :: handle, iatom, ikind, ispin, natom, &
91 : nkind, nspin
92 26 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
93 : LOGICAL :: use_virial
94 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: v_dadr, v_dfdr
95 : TYPE(atomic_kind_type), POINTER :: atomic_kind
96 : TYPE(kpoint_type), POINTER :: kpoints
97 26 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
98 26 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
99 : TYPE(qs_ks_env_type), POINTER :: ks_env
100 : TYPE(virial_type), POINTER :: virial
101 :
102 26 : CALL timeset(routineN, handle)
103 26 : NULLIFY (atomic_kind, force, lri_coef, v_dadr, v_dfdr, virial)
104 :
105 26 : IF (ASSOCIATED(lri_env%soo_list)) THEN
106 :
107 26 : nkind = lri_env%lri_ints%nkind
108 26 : nspin = SIZE(pmatrix, 1)
109 26 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
110 26 : CALL get_qs_env(qs_env=qs_env, force=force, virial=virial)
111 26 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
112 26 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
113 26 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
114 26 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
115 :
116 : !calculate SUM_i integral(V*fbas_i)*davec/dR
117 : CALL calculate_v_dadr_sr(lri_env, lri_density, pmatrix, cell_to_index, atomic_kind_set, &
118 26 : use_virial, virial)
119 26 : IF (lri_env%distant_pair_approximation) THEN
120 : CALL calculate_v_dadr_ff(lri_env, lri_density, pmatrix, cell_to_index, atomic_kind_set, &
121 8 : use_virial, virial)
122 : END IF
123 :
124 52 : DO ispin = 1, nspin
125 :
126 26 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
127 :
128 104 : DO ikind = 1, nkind
129 52 : atomic_kind => atomic_kind_set(ikind)
130 52 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom)
131 200 : DO iatom = 1, natom
132 122 : v_dadr => lri_coef(ikind)%v_dadr(iatom, :)
133 122 : v_dfdr => lri_coef(ikind)%v_dfdr(iatom, :)
134 :
135 : force(ikind)%rho_lri_elec(:, iatom) = force(ikind)%rho_lri_elec(:, iatom) &
136 1028 : + v_dfdr(:) + v_dadr(:)
137 :
138 : END DO
139 : END DO
140 : END DO
141 :
142 : END IF
143 :
144 26 : CALL timestop(handle)
145 :
146 26 : END SUBROUTINE calculate_lri_forces
147 :
148 : ! **************************************************************************************************
149 : !> \brief calculates second term of derivative with respect to R, i.e.
150 : !> SUM_i integral(V * fbas_i)*davec/dR
151 : !> \param lri_env ...
152 : !> \param lri_density ...
153 : !> \param pmatrix ...
154 : !> \param cell_to_index ...
155 : !> \param atomic_kind_set ...
156 : !> \param use_virial ...
157 : !> \param virial ...
158 : ! **************************************************************************************************
159 26 : SUBROUTINE calculate_v_dadr_sr(lri_env, lri_density, pmatrix, cell_to_index, atomic_kind_set, &
160 : use_virial, virial)
161 :
162 : TYPE(lri_environment_type), POINTER :: lri_env
163 : TYPE(lri_density_type), POINTER :: lri_density
164 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
165 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
166 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
167 : LOGICAL, INTENT(IN) :: use_virial
168 : TYPE(virial_type), POINTER :: virial
169 :
170 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_v_dadr_sr'
171 :
172 : INTEGER :: atom_a, atom_b, handle, i, iac, iatom, ic, ikind, ilist, ispin, jatom, jkind, &
173 : jneighbor, k, mepos, nba, nbb, nfa, nfb, nkind, nlist, nn, nspin, nthread
174 26 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
175 : INTEGER, DIMENSION(3) :: cell
176 : LOGICAL :: found, use_cell_mapping
177 : REAL(KIND=dp) :: ai, dab, dfw, fw, isn, threshold
178 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vint
179 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab
180 : REAL(KIND=dp), DIMENSION(3) :: dcharge, dlambda, force_a, force_b, &
181 : idav, nsdssn, nsdsst, nsdt, rab
182 26 : REAL(KIND=dp), DIMENSION(:), POINTER :: st, v_dadra, v_dadrb
183 26 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dssn, dsst, dtvec, pbij, sdssn, sdsst, &
184 26 : sdt
185 : TYPE(dbcsr_type), POINTER :: pmat
186 26 : TYPE(dint_type) :: lridint
187 : TYPE(gto_basis_set_type), POINTER :: fbasa, fbasb, obasa, obasb
188 : TYPE(lri_force_type), POINTER :: lri_force
189 : TYPE(lri_int_type), POINTER :: lrii
190 26 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
191 : TYPE(lri_list_type), POINTER :: lri_rho
192 : TYPE(lri_rhoab_type), POINTER :: lrho
193 : TYPE(neighbor_list_iterator_p_type), &
194 26 : DIMENSION(:), POINTER :: nl_iterator
195 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
196 26 : POINTER :: soo_list
197 :
198 26 : CALL timeset(routineN, handle)
199 :
200 26 : NULLIFY (lri_coef, lri_force, lrii, lri_rho, lrho, &
201 26 : nl_iterator, pbij, pmat, soo_list, v_dadra, v_dadrb)
202 26 : NULLIFY (dssn, dsst, dtvec, sdssn, sdsst, sdt, st)
203 :
204 26 : IF (ASSOCIATED(lri_env%soo_list)) THEN
205 26 : soo_list => lri_env%soo_list
206 :
207 26 : nkind = lri_env%lri_ints%nkind
208 26 : nspin = SIZE(pmatrix, 1)
209 : nthread = 1
210 26 : !$ nthread = omp_get_max_threads()
211 26 : use_cell_mapping = (SIZE(pmatrix, 2) > 1)
212 :
213 26 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
214 :
215 52 : DO ispin = 1, nspin
216 :
217 26 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
218 26 : lri_rho => lri_density%lri_rhos(ispin)%lri_list
219 :
220 26 : CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
221 : !$OMP PARALLEL DEFAULT(NONE)&
222 : !$OMP SHARED (nthread,nl_iterator,pmatrix,nkind,lri_env,lri_rho,lri_coef,&
223 : !$OMP atom_of_kind,virial,use_virial,cell_to_index,use_cell_mapping,ispin)&
224 : !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,jneighbor,rab,iac,lrii,lrho,&
225 : !$OMP lri_force,nfa,nfb,nba,nbb,nn,st,nsdsst,nsdssn,dsst,sdsst,dssn,sdssn,sdt,&
226 : !$OMP nsdt,dtvec,pbij,found,obasa,obasb,fbasa,fbasb,i,k,threshold,&
227 : !$OMP dcharge,dlambda,atom_a,atom_b,v_dadra,v_dadrb,force_a,force_b,&
228 26 : !$OMP lridint,pab,fw,dfw,dab,ai,vint,isn,idav,cell,ic,pmat)
229 :
230 : mepos = 0
231 : !$ mepos = omp_get_thread_num()
232 :
233 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
234 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
235 : iatom=iatom, jatom=jatom, nlist=nlist, ilist=ilist, &
236 : inode=jneighbor, r=rab, cell=cell)
237 :
238 : iac = ikind + nkind*(jkind - 1)
239 :
240 : IF (.NOT. ASSOCIATED(lri_env%lri_ints%lri_atom(iac)%lri_node)) CYCLE
241 :
242 : lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
243 : lrho => lri_rho%lri_atom(iac)%lri_node(ilist)%lri_rhoab(jneighbor)
244 :
245 : ! only proceed if short-range exact LRI is needed
246 : IF (.NOT. lrii%lrisr) CYCLE
247 : fw = lrii%wsr
248 : dfw = lrii%dwsr
249 :
250 : ! zero contribution for pairs aa or bb and periodc pairs aa' or bb'
251 : ! calculate forces for periodic pairs aa' and bb' only for virial
252 : IF (.NOT. lrii%calc_force_pair) CYCLE
253 :
254 : nfa = lrii%nfa
255 : nfb = lrii%nfb
256 : nba = lrii%nba
257 : nbb = lrii%nbb
258 : nn = nfa + nfb
259 :
260 : IF (use_cell_mapping) THEN
261 : ic = cell_to_index(cell(1), cell(2), cell(3))
262 : CPASSERT(ic > 0)
263 : ELSE
264 : ic = 1
265 : END IF
266 : pmat => pmatrix(ispin, ic)%matrix
267 :
268 : ! get the density matrix Pab
269 : ALLOCATE (pab(nba, nbb))
270 : NULLIFY (pbij)
271 : IF (iatom <= jatom) THEN
272 : CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pbij, found=found)
273 : CPASSERT(found)
274 : pab(1:nba, 1:nbb) = pbij(1:nba, 1:nbb)
275 : ELSE
276 : CALL dbcsr_get_block_p(matrix=pmat, row=jatom, col=iatom, block=pbij, found=found)
277 : CPASSERT(found)
278 : pab(1:nba, 1:nbb) = TRANSPOSE(pbij(1:nbb, 1:nba))
279 : END IF
280 :
281 : obasa => lri_env%orb_basis(ikind)%gto_basis_set
282 : obasb => lri_env%orb_basis(jkind)%gto_basis_set
283 : fbasa => lri_env%ri_basis(ikind)%gto_basis_set
284 : fbasb => lri_env%ri_basis(jkind)%gto_basis_set
285 :
286 : ! Calculate derivatives of overlap integrals (a,b), (fa,fb), (a,fa,b), (a,b,fb)
287 : CALL allocate_int_type(lridint=lridint, nba=nba, nbb=nbb, nfa=nfa, nfb=nfb)
288 : CALL lri_dint2(lri_env, lrii, lridint, rab, obasa, obasb, fbasa, fbasb, &
289 : iatom, jatom, ikind, jkind)
290 :
291 : NULLIFY (lri_force)
292 : CALL allocate_lri_force_components(lri_force, nfa, nfb)
293 : st => lri_force%st
294 : dsst => lri_force%dsst
295 : dssn => lri_force%dssn
296 : dtvec => lri_force%dtvec
297 :
298 : ! compute dtvec/dRa = SUM_ab Pab *d(a,b,x)/dRa
299 : DO k = 1, 3
300 : threshold = 0.01_dp*lri_env%eps_o3_int/MAX(SUM(ABS(pab(1:nba, 1:nbb))), 1.0e-14_dp)
301 : dcharge(k) = SUM(pab(1:nba, 1:nbb)*lridint%dsooint(1:nba, 1:nbb, k))
302 : DO i = 1, nfa
303 : IF (lrii%abascr(i) > threshold) THEN
304 : dtvec(i, k) = SUM(pab(1:nba, 1:nbb)*lridint%dabaint(1:nba, 1:nbb, i, k))
305 : END IF
306 : END DO
307 : DO i = 1, nfb
308 : IF (lrii%abbscr(i) > threshold) THEN
309 : dtvec(nfa + i, k) = SUM(pab(1:nba, 1:nbb)*lridint%dabbint(1:nba, 1:nbb, i, k))
310 : END IF
311 : END DO
312 : END DO
313 :
314 : DEALLOCATE (pab)
315 :
316 : atom_a = atom_of_kind(iatom)
317 : atom_b = atom_of_kind(jatom)
318 : ALLOCATE (vint(nn))
319 : vint(1:nfa) = lri_coef(ikind)%v_int(atom_a, 1:nfa)
320 : vint(nfa + 1:nn) = lri_coef(jkind)%v_int(atom_b, 1:nfb)
321 :
322 : isn = SUM(vint(1:nn)*lrii%sn(1:nn))
323 : DO k = 1, 3
324 : ai = isn/lrii%nsn*dcharge(k)
325 : force_a(k) = 2.0_dp*fw*ai
326 : force_b(k) = -2.0_dp*fw*ai
327 : END DO
328 :
329 : ! derivative of S (overlap) matrix dS
330 : !dS: dsaa and dsbb are zero, only work with ab blocks in following
331 : st(1:nn) = MATMUL(lrii%sinv(1:nn, 1:nn), lrho%tvec(1:nn))
332 : DO k = 1, 3
333 : dsst(1:nfa, k) = MATMUL(lridint%dsabint(1:nfa, 1:nfb, k), st(nfa + 1:nn))
334 : dsst(nfa + 1:nn, k) = MATMUL(st(1:nfa), lridint%dsabint(1:nfa, 1:nfb, k))
335 : nsdsst(k) = SUM(lrii%sn(1:nn)*dsst(1:nn, k))
336 : dssn(1:nfa, k) = MATMUL(lridint%dsabint(1:nfa, 1:nfb, k), lrii%sn(nfa + 1:nn))
337 : dssn(nfa + 1:nn, k) = MATMUL(lrii%sn(1:nfa), lridint%dsabint(1:nfa, 1:nfb, k))
338 : nsdssn(k) = SUM(lrii%sn(1:nn)*dssn(1:nn, k))
339 : nsdt(k) = SUM(dtvec(1:nn, k)*lrii%sn(1:nn))
340 : END DO
341 : ! dlambda/dRa
342 : DO k = 1, 3
343 : dlambda(k) = (nsdsst(k) - nsdt(k))/lrii%nsn &
344 : + (lrho%charge - lrho%nst)*nsdssn(k)/(lrii%nsn*lrii%nsn)
345 : END DO
346 : DO k = 1, 3
347 : force_a(k) = force_a(k) + 2.0_dp*fw*isn*dlambda(k)
348 : force_b(k) = force_b(k) - 2.0_dp*fw*isn*dlambda(k)
349 : END DO
350 : DO k = 1, 3
351 : st(1:nn) = dtvec(1:nn, k) - dsst(1:nn, k) - lrho%lambda*dssn(1:nn, k)
352 : idav(k) = SUM(vint(1:nn)*MATMUL(lrii%sinv(1:nn, 1:nn), st(1:nn)))
353 : END DO
354 :
355 : ! deallocate derivative integrals
356 : CALL deallocate_int_type(lridint=lridint)
357 :
358 : ! sum over atom pairs
359 : DO k = 1, 3
360 : ai = 2.0_dp*fw*idav(k)
361 : force_a(k) = force_a(k) + ai
362 : force_b(k) = force_b(k) - ai
363 : END DO
364 : IF (ABS(dfw) > 0.0_dp) THEN
365 : dab = SQRT(SUM(rab(1:3)*rab(1:3)))
366 : ai = 2.0_dp*dfw/dab*SUM(lrho%avec(1:nn)*vint(1:nn))
367 : DO k = 1, 3
368 : force_a(k) = force_a(k) - ai*rab(k)
369 : force_b(k) = force_b(k) + ai*rab(k)
370 : END DO
371 : END IF
372 :
373 : DEALLOCATE (vint)
374 :
375 : v_dadra => lri_coef(ikind)%v_dadr(atom_a, :)
376 : v_dadrb => lri_coef(jkind)%v_dadr(atom_b, :)
377 : !$OMP CRITICAL(addforces)
378 : DO k = 1, 3
379 : v_dadra(k) = v_dadra(k) + force_a(k)
380 : v_dadrb(k) = v_dadrb(k) + force_b(k)
381 : END DO
382 : !$OMP END CRITICAL(addforces)
383 :
384 : ! contribution to virial
385 : IF (use_virial) THEN
386 : !periodic self-pairs aa' contribute only with factor 0.5
387 : !$OMP CRITICAL(addvirial)
388 : IF (iatom == jatom) THEN
389 : CALL virial_pair_force(virial%pv_lrigpw, 0.5_dp, force_a, rab)
390 : CALL virial_pair_force(virial%pv_virial, 0.5_dp, force_a, rab)
391 : ELSE
392 : CALL virial_pair_force(virial%pv_lrigpw, 1.0_dp, force_a, rab)
393 : CALL virial_pair_force(virial%pv_virial, 1.0_dp, force_a, rab)
394 : END IF
395 : !$OMP END CRITICAL(addvirial)
396 : END IF
397 :
398 : CALL deallocate_lri_force_components(lri_force)
399 :
400 : END DO
401 : !$OMP END PARALLEL
402 :
403 52 : CALL neighbor_list_iterator_release(nl_iterator)
404 :
405 : END DO
406 :
407 : END IF
408 :
409 26 : CALL timestop(handle)
410 :
411 52 : END SUBROUTINE calculate_v_dadr_sr
412 :
413 : ! **************************************************************************************************
414 : !> \brief calculates second term of derivative with respect to R, i.e.
415 : !> SUM_i integral(V * fbas_i)*davec/dR for distant pair approximation
416 : !> \param lri_env ...
417 : !> \param lri_density ...
418 : !> \param pmatrix ...
419 : !> \param cell_to_index ...
420 : !> \param atomic_kind_set ...
421 : !> \param use_virial ...
422 : !> \param virial ...
423 : ! **************************************************************************************************
424 8 : SUBROUTINE calculate_v_dadr_ff(lri_env, lri_density, pmatrix, cell_to_index, atomic_kind_set, &
425 : use_virial, virial)
426 :
427 : TYPE(lri_environment_type), POINTER :: lri_env
428 : TYPE(lri_density_type), POINTER :: lri_density
429 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmatrix
430 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
431 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
432 : LOGICAL, INTENT(IN) :: use_virial
433 : TYPE(virial_type), POINTER :: virial
434 :
435 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_v_dadr_ff'
436 :
437 : INTEGER :: atom_a, atom_b, handle, i, iac, iatom, ic, ikind, ilist, ispin, jatom, jkind, &
438 : jneighbor, k, mepos, nba, nbb, nfa, nfb, nkind, nlist, nn, nspin, nthread
439 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
440 : INTEGER, DIMENSION(3) :: cell
441 : LOGICAL :: found, use_cell_mapping
442 : REAL(KIND=dp) :: ai, dab, dfw, fw, isna, isnb, threshold
443 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab, wab, wbb
444 : REAL(KIND=dp), DIMENSION(3) :: dchargea, dchargeb, force_a, force_b, &
445 : idava, idavb, rab
446 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: sta, stb, v_dadra, v_dadrb, vinta, vintb
447 8 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dtveca, dtvecb, pbij
448 : TYPE(dbcsr_type), POINTER :: pmat
449 8 : TYPE(dint_type) :: lridint
450 : TYPE(gto_basis_set_type), POINTER :: fbasa, fbasb, obasa, obasb
451 : TYPE(lri_force_type), POINTER :: lri_force_a, lri_force_b
452 : TYPE(lri_int_type), POINTER :: lrii
453 8 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
454 : TYPE(lri_list_type), POINTER :: lri_rho
455 : TYPE(lri_rhoab_type), POINTER :: lrho
456 : TYPE(neighbor_list_iterator_p_type), &
457 8 : DIMENSION(:), POINTER :: nl_iterator
458 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
459 8 : POINTER :: soo_list
460 :
461 8 : CALL timeset(routineN, handle)
462 :
463 8 : NULLIFY (lri_coef, lrii, lri_rho, lrho, &
464 8 : nl_iterator, pbij, pmat, soo_list, v_dadra, v_dadrb)
465 :
466 8 : IF (ASSOCIATED(lri_env%soo_list)) THEN
467 8 : soo_list => lri_env%soo_list
468 :
469 8 : nkind = lri_env%lri_ints%nkind
470 8 : nspin = SIZE(pmatrix, 1)
471 : nthread = 1
472 8 : !$ nthread = omp_get_max_threads()
473 8 : use_cell_mapping = (SIZE(pmatrix, 2) > 1)
474 :
475 8 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
476 :
477 16 : DO ispin = 1, nspin
478 :
479 8 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
480 8 : lri_rho => lri_density%lri_rhos(ispin)%lri_list
481 :
482 8 : CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
483 : !$OMP PARALLEL DEFAULT(NONE)&
484 : !$OMP SHARED (nthread,nl_iterator,pmatrix,nkind,lri_env,lri_rho,lri_coef,&
485 : !$OMP atom_of_kind,virial,use_virial,cell_to_index,use_cell_mapping,ispin)&
486 : !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,jneighbor,rab,iac,lrii,lrho,&
487 : !$OMP lri_force_a,lri_force_b,nfa,nfb,nba,nbb,nn,sta,stb,&
488 : !$OMP dtveca,dtvecb,pbij,found,obasa,obasb,fbasa,fbasb,i,k,threshold,&
489 : !$OMP dchargea,dchargeb,atom_a,atom_b,v_dadra,v_dadrb,force_a,force_b,&
490 8 : !$OMP lridint,pab,wab,wbb,fw,dfw,dab,ai,vinta,vintb,isna,isnb,idava,idavb,cell,ic,pmat)
491 :
492 : mepos = 0
493 : !$ mepos = omp_get_thread_num()
494 :
495 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
496 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
497 : iatom=iatom, jatom=jatom, nlist=nlist, ilist=ilist, &
498 : inode=jneighbor, r=rab, cell=cell)
499 :
500 : iac = ikind + nkind*(jkind - 1)
501 :
502 : IF (.NOT. ASSOCIATED(lri_env%lri_ints%lri_atom(iac)%lri_node)) CYCLE
503 :
504 : lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
505 : lrho => lri_rho%lri_atom(iac)%lri_node(ilist)%lri_rhoab(jneighbor)
506 :
507 : ! only proceed if short-range exact LRI is needed
508 : IF (.NOT. lrii%lriff) CYCLE
509 : fw = lrii%wff
510 : dfw = lrii%dwff
511 :
512 : ! zero contribution for pairs aa or bb and periodc pairs aa' or bb'
513 : ! calculate forces for periodic pairs aa' and bb' only for virial
514 : IF (.NOT. lrii%calc_force_pair) CYCLE
515 :
516 : nfa = lrii%nfa
517 : nfb = lrii%nfb
518 : nba = lrii%nba
519 : nbb = lrii%nbb
520 : nn = nfa + nfb
521 :
522 : IF (use_cell_mapping) THEN
523 : ic = cell_to_index(cell(1), cell(2), cell(3))
524 : CPASSERT(ic > 0)
525 : ELSE
526 : ic = 1
527 : END IF
528 : pmat => pmatrix(ispin, ic)%matrix
529 :
530 : ! get the density matrix Pab
531 : ALLOCATE (pab(nba, nbb))
532 : NULLIFY (pbij)
533 : IF (iatom <= jatom) THEN
534 : CALL dbcsr_get_block_p(matrix=pmat, row=iatom, col=jatom, block=pbij, found=found)
535 : CPASSERT(found)
536 : pab(1:nba, 1:nbb) = pbij(1:nba, 1:nbb)
537 : ELSE
538 : CALL dbcsr_get_block_p(matrix=pmat, row=jatom, col=iatom, block=pbij, found=found)
539 : CPASSERT(found)
540 : pab(1:nba, 1:nbb) = TRANSPOSE(pbij(1:nbb, 1:nba))
541 : END IF
542 :
543 : ALLOCATE (wab(nba, nbb), wbb(nba, nbb))
544 : wab(1:nba, 1:nbb) = lri_env%wmat(ikind, jkind)%mat(1:nba, 1:nbb)
545 : wbb(1:nba, 1:nbb) = 1.0_dp - lri_env%wmat(ikind, jkind)%mat(1:nba, 1:nbb)
546 :
547 : obasa => lri_env%orb_basis(ikind)%gto_basis_set
548 : obasb => lri_env%orb_basis(jkind)%gto_basis_set
549 : fbasa => lri_env%ri_basis(ikind)%gto_basis_set
550 : fbasb => lri_env%ri_basis(jkind)%gto_basis_set
551 :
552 : ! Calculate derivatives of overlap integrals (a,b), (fa,fb), (a,fa,b), (a,b,fb)
553 : CALL allocate_int_type(lridint=lridint, nba=nba, nbb=nbb, nfa=nfa, nfb=nfb, &
554 : skip_dsab=.TRUE.)
555 : CALL lri_dint2(lri_env, lrii, lridint, rab, obasa, obasb, fbasa, fbasb, &
556 : iatom, jatom, ikind, jkind)
557 :
558 : NULLIFY (lri_force_a, lri_force_b)
559 : CALL allocate_lri_force_components(lri_force_a, nfa, 0)
560 : CALL allocate_lri_force_components(lri_force_b, 0, nfb)
561 : dtveca => lri_force_a%dtvec
562 : dtvecb => lri_force_b%dtvec
563 : sta => lri_force_a%st
564 : stb => lri_force_b%st
565 :
566 : ! compute dtvec/dRa = SUM_ab Pab *d(a,b,x)/dRa
567 : threshold = 0.01_dp*lri_env%eps_o3_int/MAX(SUM(ABS(pab(1:nba, 1:nbb))), 1.0e-14_dp)
568 : DO k = 1, 3
569 : dchargea(k) = SUM(wab(1:nba, 1:nbb)*pab(1:nba, 1:nbb)*lridint%dsooint(1:nba, 1:nbb, k))
570 : DO i = 1, nfa
571 : IF (lrii%abascr(i) > threshold) THEN
572 : dtveca(i, k) = SUM(wab(1:nba, 1:nbb)*pab(1:nba, 1:nbb)*lridint%dabaint(1:nba, 1:nbb, i, k))
573 : END IF
574 : END DO
575 : dchargeb(k) = SUM(wbb(1:nba, 1:nbb)*pab(1:nba, 1:nbb)*lridint%dsooint(1:nba, 1:nbb, k))
576 : DO i = 1, nfb
577 : IF (lrii%abbscr(i) > threshold) THEN
578 : dtvecb(i, k) = SUM(wbb(1:nba, 1:nbb)*pab(1:nba, 1:nbb)*lridint%dabbint(1:nba, 1:nbb, i, k))
579 : END IF
580 : END DO
581 : END DO
582 :
583 : DEALLOCATE (pab, wab, wbb)
584 :
585 : atom_a = atom_of_kind(iatom)
586 : atom_b = atom_of_kind(jatom)
587 : vinta => lri_coef(ikind)%v_int(atom_a, :)
588 : vintb => lri_coef(jkind)%v_int(atom_b, :)
589 :
590 : isna = SUM(vinta(1:nfa)*lrii%sna(1:nfa))
591 : isnb = SUM(vintb(1:nfb)*lrii%snb(1:nfb))
592 : DO k = 1, 3
593 : ai = isna/lrii%nsna*dchargea(k) + isnb/lrii%nsnb*dchargeb(k)
594 : force_a(k) = 2.0_dp*fw*ai
595 : force_b(k) = -2.0_dp*fw*ai
596 : END DO
597 :
598 : DO k = 1, 3
599 : sta(1:nfa) = MATMUL(lrii%asinv(1:nfa, 1:nfa), dtveca(1:nfa, k))
600 : idava(k) = SUM(vinta(1:nfa)*sta(1:nfa)) - isna/lrii%nsna*SUM(lrii%na(1:nfa)*sta(1:nfa))
601 : stb(1:nfb) = MATMUL(lrii%bsinv(1:nfb, 1:nfb), dtvecb(1:nfb, k))
602 : idavb(k) = SUM(vintb(1:nfb)*stb(1:nfb)) - isnb/lrii%nsnb*SUM(lrii%nb(1:nfb)*stb(1:nfb))
603 : END DO
604 :
605 : ! deallocate derivative integrals
606 : CALL deallocate_int_type(lridint=lridint)
607 :
608 : ! sum over atom pairs
609 : DO k = 1, 3
610 : ai = 2.0_dp*fw*(idava(k) + idavb(k))
611 : force_a(k) = force_a(k) + ai
612 : force_b(k) = force_b(k) - ai
613 : END DO
614 : IF (ABS(dfw) > 0.0_dp) THEN
615 : dab = SQRT(SUM(rab(1:3)*rab(1:3)))
616 : ai = 2.0_dp*dfw/dab* &
617 : (SUM(lrho%aveca(1:nfa)*vinta(1:nfa)) + &
618 : SUM(lrho%avecb(1:nfb)*vintb(1:nfb)))
619 : DO k = 1, 3
620 : force_a(k) = force_a(k) - ai*rab(k)
621 : force_b(k) = force_b(k) + ai*rab(k)
622 : END DO
623 : END IF
624 : v_dadra => lri_coef(ikind)%v_dadr(atom_a, :)
625 : v_dadrb => lri_coef(jkind)%v_dadr(atom_b, :)
626 :
627 : !$OMP CRITICAL(addforces)
628 : DO k = 1, 3
629 : v_dadra(k) = v_dadra(k) + force_a(k)
630 : v_dadrb(k) = v_dadrb(k) + force_b(k)
631 : END DO
632 : !$OMP END CRITICAL(addforces)
633 :
634 : ! contribution to virial
635 : IF (use_virial) THEN
636 : !periodic self-pairs aa' contribute only with factor 0.5
637 : !$OMP CRITICAL(addvirial)
638 : IF (iatom == jatom) THEN
639 : CALL virial_pair_force(virial%pv_lrigpw, 0.5_dp, force_a, rab)
640 : CALL virial_pair_force(virial%pv_virial, 0.5_dp, force_a, rab)
641 : ELSE
642 : CALL virial_pair_force(virial%pv_lrigpw, 1.0_dp, force_a, rab)
643 : CALL virial_pair_force(virial%pv_virial, 1.0_dp, force_a, rab)
644 : END IF
645 : !$OMP END CRITICAL(addvirial)
646 : END IF
647 :
648 : CALL deallocate_lri_force_components(lri_force_a)
649 : CALL deallocate_lri_force_components(lri_force_b)
650 :
651 : END DO
652 : !$OMP END PARALLEL
653 :
654 16 : CALL neighbor_list_iterator_release(nl_iterator)
655 :
656 : END DO
657 :
658 : END IF
659 :
660 8 : CALL timestop(handle)
661 :
662 16 : END SUBROUTINE calculate_v_dadr_ff
663 :
664 : ! **************************************************************************************************
665 : !> \brief calculates the ri forces
666 : !> \param lri_env ...
667 : !> \param lri_density ...
668 : !> \param qs_env ...
669 : !> \param pmatrix density matrix
670 : !> \param atomic_kind_set ...
671 : ! **************************************************************************************************
672 0 : SUBROUTINE calculate_ri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
673 :
674 : TYPE(lri_environment_type), POINTER :: lri_env
675 : TYPE(lri_density_type), POINTER :: lri_density
676 : TYPE(qs_environment_type), POINTER :: qs_env
677 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix
678 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
679 :
680 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_ri_forces'
681 :
682 : INTEGER :: handle, iatom, ikind, ispin, natom, &
683 : nkind, nspin
684 : LOGICAL :: use_virial
685 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: v_dadr, v_dfdr
686 : TYPE(atomic_kind_type), POINTER :: atomic_kind
687 0 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
688 0 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
689 : TYPE(virial_type), POINTER :: virial
690 :
691 0 : CALL timeset(routineN, handle)
692 0 : NULLIFY (atomic_kind, force, lri_coef, v_dadr, v_dfdr, virial)
693 :
694 0 : nkind = SIZE(atomic_kind_set)
695 0 : nspin = SIZE(pmatrix)
696 0 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
697 0 : CALL get_qs_env(qs_env=qs_env, force=force, virial=virial)
698 0 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
699 :
700 : !calculate SUM_i integral(V*fbas_i)*davec/dR
701 : CALL calculate_v_dadr_ri(lri_env, lri_density, pmatrix, atomic_kind_set, &
702 0 : use_virial, virial)
703 :
704 0 : DO ispin = 1, nspin
705 :
706 0 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
707 :
708 0 : DO ikind = 1, nkind
709 0 : atomic_kind => atomic_kind_set(ikind)
710 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom)
711 0 : DO iatom = 1, natom
712 0 : v_dadr => lri_coef(ikind)%v_dadr(iatom, :)
713 0 : v_dfdr => lri_coef(ikind)%v_dfdr(iatom, :)
714 :
715 : force(ikind)%rho_lri_elec(:, iatom) = force(ikind)%rho_lri_elec(:, iatom) &
716 0 : + v_dfdr(:) + v_dadr(:)
717 :
718 : END DO
719 : END DO
720 : END DO
721 :
722 0 : CALL timestop(handle)
723 :
724 0 : END SUBROUTINE calculate_ri_forces
725 :
726 : ! **************************************************************************************************
727 : !> \brief calculates second term of derivative with respect to R, i.e.
728 : !> SUM_i integral(V * fbas_i)*davec/dR
729 : !> However contributions from charge constraint and derivative of overlap matrix have been
730 : !> calculated in calculate_ri_integrals
731 : !> \param lri_env ...
732 : !> \param lri_density ...
733 : !> \param pmatrix ...
734 : !> \param atomic_kind_set ...
735 : !> \param use_virial ...
736 : !> \param virial ...
737 : ! **************************************************************************************************
738 0 : SUBROUTINE calculate_v_dadr_ri(lri_env, lri_density, pmatrix, atomic_kind_set, &
739 : use_virial, virial)
740 :
741 : TYPE(lri_environment_type), POINTER :: lri_env
742 : TYPE(lri_density_type), POINTER :: lri_density
743 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix
744 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
745 : LOGICAL, INTENT(IN) :: use_virial
746 : TYPE(virial_type), POINTER :: virial
747 :
748 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_v_dadr_ri'
749 :
750 : INTEGER :: atom_a, atom_b, atom_c, handle, i, i1, &
751 : i2, iatom, ikind, ispin, jatom, jkind, &
752 : katom, kkind, m, mepos, nkind, nspin, &
753 : nthread
754 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
755 0 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
756 : REAL(KIND=dp) :: fscal
757 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, force_c, rij, rik, rjk
758 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: v_dadra, v_dadrb, v_dadrc
759 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fi, fj, fk, fo
760 0 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
761 : TYPE(o3c_iterator_type) :: o3c_iterator
762 :
763 0 : CALL timeset(routineN, handle)
764 :
765 0 : nspin = SIZE(pmatrix)
766 0 : nkind = SIZE(atomic_kind_set)
767 0 : DO ispin = 1, nspin
768 0 : lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
769 0 : DO ikind = 1, nkind
770 0 : lri_coef(ikind)%v_dadr(:, :) = 0.0_dp
771 : END DO
772 : END DO
773 :
774 0 : bas_ptr => lri_env%ri_fit%bas_ptr
775 :
776 0 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
777 :
778 : ! contribution from 3-center overlap integral derivatives
779 0 : fo => lri_env%ri_fit%fout
780 : nthread = 1
781 0 : !$ nthread = omp_get_max_threads()
782 :
783 0 : CALL o3c_iterator_create(lri_env%o3c, o3c_iterator, nthread=nthread)
784 :
785 : !$OMP PARALLEL DEFAULT(NONE)&
786 : !$OMP SHARED (nthread,o3c_iterator,bas_ptr,fo,nspin,atom_of_kind,lri_coef,virial,use_virial)&
787 : !$OMP PRIVATE (mepos,ikind,jkind,kkind,iatom,jatom,katom,fi,fj,fk,fscal,&
788 : !$OMP force_a,force_b,force_c,i1,i2,m,i,ispin,atom_a,atom_b,atom_c,&
789 0 : !$OMP v_dadra,v_dadrb,v_dadrc,rij,rik,rjk)
790 :
791 : mepos = 0
792 : !$ mepos = omp_get_thread_num()
793 :
794 : DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
795 : CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, &
796 : ikind=ikind, jkind=jkind, kkind=kkind, &
797 : iatom=iatom, jatom=jatom, katom=katom, &
798 : rij=rij, rik=rik, force_i=fi, force_j=fj, force_k=fk)
799 : i1 = bas_ptr(1, katom)
800 : i2 = bas_ptr(2, katom)
801 : m = i2 - i1 + 1
802 : DO i = 1, 3
803 : force_a(i) = 0.0_dp
804 : force_b(i) = 0.0_dp
805 : force_c(i) = 0.0_dp
806 : DO ispin = 1, nspin
807 : force_a(i) = force_a(i) + SUM(fi(1:m, i)*fo(i1:i2, ispin))
808 : force_b(i) = force_b(i) + SUM(fj(1:m, i)*fo(i1:i2, ispin))
809 : force_c(i) = force_c(i) + SUM(fk(1:m, i)*fo(i1:i2, ispin))
810 : END DO
811 : END DO
812 : atom_a = atom_of_kind(iatom)
813 : atom_b = atom_of_kind(jatom)
814 : atom_c = atom_of_kind(katom)
815 : !
816 : v_dadra => lri_coef(ikind)%v_dadr(atom_a, :)
817 : v_dadrb => lri_coef(jkind)%v_dadr(atom_b, :)
818 : v_dadrc => lri_coef(kkind)%v_dadr(atom_c, :)
819 : !
820 : !$OMP CRITICAL(addforce)
821 : v_dadra(1:3) = v_dadra(1:3) + force_a(1:3)
822 : v_dadrb(1:3) = v_dadrb(1:3) + force_b(1:3)
823 : v_dadrc(1:3) = v_dadrc(1:3) + force_c(1:3)
824 : !
825 : IF (use_virial) THEN
826 : rjk(1:3) = rik(1:3) - rij(1:3)
827 : ! to be debugged
828 : fscal = 1.0_dp
829 : IF (iatom == jatom) fscal = 1.0_dp
830 : CALL virial_pair_force(virial%pv_lrigpw, 1.0_dp, force_a, rik)
831 : CALL virial_pair_force(virial%pv_lrigpw, 1.0_dp, force_b, rjk)
832 : CALL virial_pair_force(virial%pv_virial, 1.0_dp, force_a, rik)
833 : CALL virial_pair_force(virial%pv_virial, 1.0_dp, force_b, rjk)
834 : END IF
835 : !$OMP END CRITICAL(addforce)
836 : END DO
837 : !$OMP END PARALLEL
838 0 : CALL o3c_iterator_release(o3c_iterator)
839 :
840 0 : CALL timestop(handle)
841 :
842 0 : END SUBROUTINE calculate_v_dadr_ri
843 :
844 : END MODULE lri_forces
|
