(git:e7e05ae)
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 ! **************************************************************************************************
16 MODULE lri_forces
17 
18  USE atomic_kind_types, ONLY: atomic_kind_type,&
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: &
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,&
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
44  neighbor_list_iterator_p_type,&
46  neighbor_list_set_p_type
48  o3c_iterate,&
51  o3c_iterator_type
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 
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  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
193  TYPE(neighbor_list_iterator_p_type), &
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
456  TYPE(neighbor_list_iterator_p_type), &
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 
844 END 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.
Definition: kpoint_types.F:15
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.
Definition: kpoint_types.F:333
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.
Definition: lri_integrals.F:19
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)
...
Definition: qs_ks_types.F:330
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
Definition: qs_o3c_types.F:12
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)
...
Definition: qs_o3c_types.F:418
subroutine, public o3c_iterator_create(o3c, o3c_iterator, nthread)
...
Definition: qs_o3c_types.F:357
subroutine, public o3c_iterator_release(o3c_iterator)
...
Definition: qs_o3c_types.F:384
integer function, public o3c_iterate(o3c_iterator, mepos)
...
Definition: qs_o3c_types.F:473
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.