(git:374b731)
Loading...
Searching...
No Matches
lri_forces.F
Go to the documentation of this file.
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! **************************************************************************************************
17
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,&
28 USE lri_environment_types, ONLY: &
34 dint_type,&
39 USE qs_ks_types, ONLY: get_ks_env,&
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
67
68! **************************************************************************************************
69
70CONTAINS
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 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 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
93 LOGICAL :: use_virial
94 REAL(kind=dp), DIMENSION(:), POINTER :: v_dadr, v_dfdr
95 TYPE(atomic_kind_type), POINTER :: atomic_kind
96 TYPE(kpoint_type), POINTER :: kpoints
97 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
98 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
99 TYPE(qs_ks_env_type), POINTER :: ks_env
100 TYPE(virial_type), POINTER :: virial
101
102 CALL timeset(routinen, handle)
103 NULLIFY (atomic_kind, force, lri_coef, v_dadr, v_dfdr, virial)
104
105 IF (ASSOCIATED(lri_env%soo_list)) THEN
106
107 nkind = lri_env%lri_ints%nkind
108 nspin = SIZE(pmatrix, 1)
109 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
110 CALL get_qs_env(qs_env=qs_env, force=force, virial=virial)
111 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
112 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
113 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
114 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 use_virial, virial)
119 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 use_virial, virial)
122 END IF
123
124 DO ispin = 1, nspin
125
126 lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
127
128 DO ikind = 1, nkind
129 atomic_kind => atomic_kind_set(ikind)
130 CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom)
131 DO iatom = 1, natom
132 v_dadr => lri_coef(ikind)%v_dadr(iatom, :)
133 v_dfdr => lri_coef(ikind)%v_dfdr(iatom, :)
134
135 force(ikind)%rho_lri_elec(:, iatom) = force(ikind)%rho_lri_elec(:, iatom) &
136 + v_dfdr(:) + v_dadr(:)
137
138 END DO
139 END DO
140 END DO
141
142 END IF
143
144 CALL timestop(handle)
145
146 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 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: vint
179 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 REAL(kind=dp), DIMENSION(:), POINTER :: st, v_dadra, v_dadrb
183 REAL(kind=dp), DIMENSION(:, :), POINTER :: dssn, dsst, dtvec, pbij, sdssn, sdsst, &
184 sdt
185 TYPE(dbcsr_type), POINTER :: pmat
186 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 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
191 TYPE(lri_list_type), POINTER :: lri_rho
192 TYPE(lri_rhoab_type), POINTER :: lrho
194 DIMENSION(:), POINTER :: nl_iterator
195 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
196 POINTER :: soo_list
197
198 CALL timeset(routinen, handle)
199
200 NULLIFY (lri_coef, lri_force, lrii, lri_rho, lrho, &
201 nl_iterator, pbij, pmat, soo_list, v_dadra, v_dadrb)
202 NULLIFY (dssn, dsst, dtvec, sdssn, sdsst, sdt, st)
203
204 IF (ASSOCIATED(lri_env%soo_list)) THEN
205 soo_list => lri_env%soo_list
206
207 nkind = lri_env%lri_ints%nkind
208 nspin = SIZE(pmatrix, 1)
209 nthread = 1
210!$ nthread = omp_get_max_threads()
211 use_cell_mapping = (SIZE(pmatrix, 2) > 1)
212
213 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
214
215 DO ispin = 1, nspin
216
217 lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
218 lri_rho => lri_density%lri_rhos(ispin)%lri_list
219
220 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!$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 CALL neighbor_list_iterator_release(nl_iterator)
404
405 END DO
406
407 END IF
408
409 CALL timestop(handle)
410
411 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 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 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: sta, stb, v_dadra, v_dadrb, vinta, vintb
447 REAL(kind=dp), DIMENSION(:, :), POINTER :: dtveca, dtvecb, pbij
448 TYPE(dbcsr_type), POINTER :: pmat
449 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 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
454 TYPE(lri_list_type), POINTER :: lri_rho
455 TYPE(lri_rhoab_type), POINTER :: lrho
457 DIMENSION(:), POINTER :: nl_iterator
458 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
459 POINTER :: soo_list
460
461 CALL timeset(routinen, handle)
462
463 NULLIFY (lri_coef, lrii, lri_rho, lrho, &
464 nl_iterator, pbij, pmat, soo_list, v_dadra, v_dadrb)
465
466 IF (ASSOCIATED(lri_env%soo_list)) THEN
467 soo_list => lri_env%soo_list
468
469 nkind = lri_env%lri_ints%nkind
470 nspin = SIZE(pmatrix, 1)
471 nthread = 1
472!$ nthread = omp_get_max_threads()
473 use_cell_mapping = (SIZE(pmatrix, 2) > 1)
474
475 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
476
477 DO ispin = 1, nspin
478
479 lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
480 lri_rho => lri_density%lri_rhos(ispin)%lri_list
481
482 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!$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 CALL neighbor_list_iterator_release(nl_iterator)
655
656 END DO
657
658 END IF
659
660 CALL timestop(handle)
661
662 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: v_dadr, v_dfdr
686 TYPE(atomic_kind_type), POINTER :: atomic_kind
687 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
688 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
689 TYPE(virial_type), POINTER :: virial
690
691 CALL timeset(routinen, handle)
692 NULLIFY (atomic_kind, force, lri_coef, v_dadr, v_dfdr, virial)
693
694 nkind = SIZE(atomic_kind_set)
695 nspin = SIZE(pmatrix)
696 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom=natom)
697 CALL get_qs_env(qs_env=qs_env, force=force, virial=virial)
698 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 use_virial, virial)
703
704 DO ispin = 1, nspin
705
706 lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
707
708 DO ikind = 1, nkind
709 atomic_kind => atomic_kind_set(ikind)
710 CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom)
711 DO iatom = 1, natom
712 v_dadr => lri_coef(ikind)%v_dadr(iatom, :)
713 v_dfdr => lri_coef(ikind)%v_dfdr(iatom, :)
714
715 force(ikind)%rho_lri_elec(:, iatom) = force(ikind)%rho_lri_elec(:, iatom) &
716 + v_dfdr(:) + v_dadr(:)
717
718 END DO
719 END DO
720 END DO
721
722 CALL timestop(handle)
723
724 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
755 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 REAL(kind=dp), DIMENSION(:), POINTER :: v_dadra, v_dadrb, v_dadrc
759 REAL(kind=dp), DIMENSION(:, :), POINTER :: fi, fj, fk, fo
760 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef
761 TYPE(o3c_iterator_type) :: o3c_iterator
762
763 CALL timeset(routinen, handle)
764
765 nspin = SIZE(pmatrix)
766 nkind = SIZE(atomic_kind_set)
767 DO ispin = 1, nspin
768 lri_coef => lri_density%lri_coefs(ispin)%lri_kinds
769 DO ikind = 1, nkind
770 lri_coef(ikind)%v_dadr(:, :) = 0.0_dp
771 END DO
772 END DO
773
774 bas_ptr => lri_env%ri_fit%bas_ptr
775
776 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 fo => lri_env%ri_fit%fout
780 nthread = 1
781!$ nthread = omp_get_max_threads()
782
783 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!$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 CALL o3c_iterator_release(o3c_iterator)
839
840 CALL timestop(handle)
841
842 END SUBROUTINE calculate_v_dadr_ri
843
844END MODULE lri_forces
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
subroutine, public deallocate_lri_force_components(lri_force)
releases the given lri_force_type
subroutine, public allocate_lri_force_components(lri_force, nfa, nfb)
creates and initializes lri_force
Calculates forces for LRIGPW method lri : local resolution of the identity.
Definition lri_forces.F:16
subroutine, public calculate_lri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the lri forces
Definition lri_forces.F:81
subroutine, public calculate_ri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the ri forces
Definition lri_forces.F:673
Calculates integrals for LRIGPW method lri : local resolution of the identity.
subroutine, public deallocate_int_type(lriint, lridint)
...
subroutine, public allocate_int_type(lriint, lridint, nba, nbb, nfa, nfb, skip_sab, skip_soo, skip_aba, skip_abb, skip_dsab, skip_dsoo, skip_daba, skip_dabb)
...
subroutine, public lri_dint2(lri_env, lrii, lridint, rab, obasa, obasb, fbasa, fbasb, iatom, jatom, ikind, jkind)
...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
3-center overlap type integrals containers
subroutine, public get_o3c_iterator_info(o3c_iterator, mepos, iatom, jatom, katom, ikind, jkind, kkind, rij, rik, cellj, cellk, integral, tvec, force_i, force_j, force_k)
...
subroutine, public o3c_iterator_create(o3c, o3c_iterator, nthread)
...
subroutine, public o3c_iterator_release(o3c_iterator)
...
integer function, public o3c_iterate(o3c_iterator, mepos)
...
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Provides all information about an atomic kind.
Contains information about kpoints.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...