(git:0de0cc2)
lri_integrals.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 integrals for LRIGPW method
10 !> lri : local resolution of the identity
11 !> \par History
12 !> created JGH [08.2012]
13 !> Dorothea Golze [02.2014] (1) extended, re-structured, cleaned
14 !> (2) heavily debugged
15 !> split off JGH [11.2017]
16 !> \authors JGH
17 !> Dorothea Golze
18 ! **************************************************************************************************
20  USE basis_set_types, ONLY: gto_basis_set_type
29  USE kinds, ONLY: dp
30  USE lri_environment_types, ONLY: lri_environment_type,&
31  lri_int_type
32 #include "./base/base_uses.f90"
33 
34  IMPLICIT NONE
35 
36  PRIVATE
37 
38 ! **************************************************************************************************
39  TYPE int_type
40  REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: sabint
41  REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: sooint
42  REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: abaint
43  REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: abbint
44  END TYPE int_type
45  !
46  TYPE dint_type
47  REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dsabint
48  REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dsooint
49  REAL(KIND=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: dabaint
50  REAL(KIND=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: dabbint
51  END TYPE dint_type
52 ! **************************************************************************************************
53 
54  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_integrals'
55 
56  PUBLIC :: lri_int, lri_int2, lri_dint, lri_dint2, int_type, dint_type, allocate_int_type, deallocate_int_type
57 
58 ! **************************************************************************************************
59 
60 CONTAINS
61 
62 ! **************************************************************************************************
63 !> \brief calcuates the lri integrals using solid harmonic Gaussians
64 !> \param lri_env ...
65 !> \param lrii ...
66 !> \param rab distance vector
67 !> \param obasa orb basis on A
68 !> \param obasb orb basis on B
69 !> \param fbasa aux basis on A
70 !> \param fbasb aux basis on B
71 !> \param iatom index atom A
72 !> \param jatom index atom B
73 !> \param ikind kind atom A
74 !> \param jkind kind atom B
75 !> \param calculate_forces ...
76 ! **************************************************************************************************
77  SUBROUTINE lri_int(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
78  iatom, jatom, ikind, jkind, calculate_forces)
79 
80  TYPE(lri_environment_type), POINTER :: lri_env
81  TYPE(lri_int_type), POINTER :: lrii
82  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
83  TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
84  INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
85  LOGICAL, INTENT(IN) :: calculate_forces
86 
87  IF (lri_env%use_shg_integrals) THEN
88  CALL lri_int_shg(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
89  iatom, jatom, ikind, jkind, calculate_forces)
90  ELSE
91  CALL lri_int_os(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
92  iatom, jatom, ikind, calculate_forces)
93  END IF
94 
95  END SUBROUTINE lri_int
96 
97 ! **************************************************************************************************
98 !> \brief calcuates the lri integrals using solid harmonic Gaussians
99 !> \param lri_env ...
100 !> \param lrii ...
101 !> \param lriint ...
102 !> \param rab distance vector
103 !> \param obasa orb basis on A
104 !> \param obasb orb basis on B
105 !> \param fbasa aux basis on A
106 !> \param fbasb aux basis on B
107 !> \param iatom index atom A
108 !> \param jatom index atom B
109 !> \param ikind kind atom A
110 !> \param jkind kind atom B
111 ! **************************************************************************************************
112  SUBROUTINE lri_int2(lri_env, lrii, lriint, rab, obasa, obasb, fbasa, fbasb, iatom, jatom, ikind, jkind)
113 
114  TYPE(lri_environment_type), POINTER :: lri_env
115  TYPE(lri_int_type), POINTER :: lrii
116  TYPE(int_type) :: lriint
117  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
118  TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
119  INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
120 
121  INTEGER :: nba, nfa
122  INTEGER, DIMENSION(:, :, :), POINTER :: fba_index, fbb_index, oba_index, &
123  obb_index
124  REAL(kind=dp) :: dab
125  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dummy1, waux_mat
126  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dummy2, dwaux_mat
127  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: scon_fba, scon_fbb, scon_oba, scon_obb
128  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix
129 
130  dab = sqrt(sum(rab*rab))
131  IF (iatom == jatom .AND. dab < lri_env%delta) THEN
132  cpassert(ALLOCATED(lriint%sabint))
133  cpassert(ALLOCATED(lriint%sooint))
134  cpassert(ALLOCATED(lriint%abaint))
135  nba = obasa%nsgf
136  nfa = fbasa%nsgf
137  lriint%sabint(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp(1:nfa, 1:nfa)
138  lriint%sooint(1:nba, 1:nba) = lri_env%bas_prop(ikind)%orb_ovlp(1:nba, 1:nba)
139  lriint%abaint(1:nba, 1:nba, 1:nfa) = lri_env%bas_prop(ikind)%ovlp3
140  ELSE
141  IF (lri_env%use_shg_integrals) THEN
142  scon_oba => lri_env%bas_prop(ikind)%scon_orb
143  scon_obb => lri_env%bas_prop(jkind)%scon_orb
144  scon_fba => lri_env%bas_prop(ikind)%scon_ri
145  scon_fbb => lri_env%bas_prop(jkind)%scon_ri
146  scona_mix => lri_env%bas_prop(ikind)%scon_mix
147  sconb_mix => lri_env%bas_prop(jkind)%scon_mix
148  oba_index => lri_env%bas_prop(ikind)%orb_index
149  fba_index => lri_env%bas_prop(ikind)%ri_index
150  obb_index => lri_env%bas_prop(jkind)%orb_index
151  fbb_index => lri_env%bas_prop(jkind)%ri_index
152  !
153  CALL lri_precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, waux_mat, dwaux_mat, &
154  .false.)
155  !*** (fa,fb)
156  IF (ALLOCATED(lriint%sabint)) THEN
157  CALL int_overlap_ab_shg_low(lriint%sabint, dummy1, rab, fbasa, fbasb, scon_fba, scon_fbb, &
158  waux_mat, dwaux_mat, .true., .false., contraction_high=.false.)
159  END IF
160  !*** (a,b)
161  IF (ALLOCATED(lriint%sooint)) THEN
162  CALL int_overlap_ab_shg_low(lriint%sooint, dummy1, rab, obasa, obasb, scon_oba, scon_obb, &
163  waux_mat, dwaux_mat, .true., .false., contraction_high=.true.)
164  END IF
165  !*** (a,b,fa)
166  IF (ALLOCATED(lriint%abaint)) THEN
167  CALL int_overlap_aba_shg_low(lriint%abaint, dummy2, rab, obasa, obasb, fbasa, &
168  scon_obb, scona_mix, oba_index, fba_index, &
169  lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
170  lri_env%cg_shg%ncg_none0, &
171  waux_mat, dwaux_mat, .true., .false.)
172  END IF
173  !*** (a,b,fb)
174  IF (ALLOCATED(lriint%abbint)) THEN
175  IF (ikind == jkind) THEN
176  cpassert(ALLOCATED(lriint%abaint))
177  CALL get_abb_same_kind(lriint%abbint, dummy2, lriint%abaint, dummy2, &
178  rab, obasa, fbasa, .true., .false.)
179  ELSE
180  CALL int_overlap_abb_shg_low(lriint%abbint, dummy2, rab, obasa, obasb, fbasb, &
181  scon_oba, sconb_mix, obb_index, fbb_index, &
182  lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
183  lri_env%cg_shg%ncg_none0, &
184  waux_mat, dwaux_mat, .true., .false.)
185  END IF
186  END IF
187  DEALLOCATE (waux_mat, dwaux_mat)
188  ELSE
189  !*** (fa,fb)
190  IF (ALLOCATED(lriint%sabint)) THEN
191  CALL int_overlap_ab_os(lriint%sabint, dummy1, rab, fbasa, fbasb, &
192  .false., lri_env%debug, lrii%dmax_ab)
193  END IF
194  !*** (a,b)
195  IF (ALLOCATED(lriint%sooint)) THEN
196  CALL int_overlap_ab_os(lriint%sooint, dummy1, rab, obasa, obasb, &
197  .false., lri_env%debug, lrii%dmax_oo)
198  END IF
199  !*** (a,b,fa)
200  IF (ALLOCATED(lriint%abaint)) THEN
201  CALL int_overlap_aba_os(lriint%abaint, dummy2, rab, obasa, obasb, fbasa, &
202  .false., lri_env%debug, lrii%dmax_aba)
203  END IF
204  !*** (a,b,fb)
205  IF (ALLOCATED(lriint%abbint)) THEN
206  CALL int_overlap_abb_os(lriint%abbint, dummy2, rab, obasa, obasb, fbasb, &
207  .false., lri_env%debug, lrii%dmax_abb)
208  END IF
209  END IF
210  END IF
211 
212  END SUBROUTINE lri_int2
213 
214 ! **************************************************************************************************
215 
216 ! **************************************************************************************************
217 !> \brief ...
218 !> \param lri_env ...
219 !> \param lrii ...
220 !> \param rab ...
221 !> \param obasa ...
222 !> \param obasb ...
223 !> \param fbasa ...
224 !> \param fbasb ...
225 !> \param iatom ...
226 !> \param jatom ...
227 !> \param ikind ...
228 !> \param jkind ...
229 !> \param calculate_forces ...
230 ! **************************************************************************************************
231  SUBROUTINE lri_dint(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
232  iatom, jatom, ikind, jkind, calculate_forces)
233 
234  TYPE(lri_environment_type), POINTER :: lri_env
235  TYPE(lri_int_type), POINTER :: lrii
236  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
237  TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
238  INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
239  LOGICAL, INTENT(IN) :: calculate_forces
240 
241  IF (lri_env%use_shg_integrals) THEN
242  CALL lri_int_shg(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
243  iatom, jatom, ikind, jkind, calculate_forces)
244  ELSE
245  CALL lri_int_os(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
246  iatom, jatom, ikind, calculate_forces)
247  END IF
248 
249  END SUBROUTINE lri_dint
250 
251 ! **************************************************************************************************
252 !> \brief ...
253 !> \param lri_env ...
254 !> \param lrii ...
255 !> \param lridint ...
256 !> \param rab ...
257 !> \param obasa ...
258 !> \param obasb ...
259 !> \param fbasa ...
260 !> \param fbasb ...
261 !> \param iatom ...
262 !> \param jatom ...
263 !> \param ikind ...
264 !> \param jkind ...
265 ! **************************************************************************************************
266  SUBROUTINE lri_dint2(lri_env, lrii, lridint, rab, obasa, obasb, fbasa, fbasb, &
267  iatom, jatom, ikind, jkind)
268 
269  TYPE(lri_environment_type), POINTER :: lri_env
270  TYPE(lri_int_type), POINTER :: lrii
271  TYPE(dint_type) :: lridint
272  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
273  TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
274  INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
275 
276  INTEGER :: nba, nbb, nfa, nfb
277  INTEGER, DIMENSION(:, :, :), POINTER :: fba_index, fbb_index, oba_index, &
278  obb_index
279  REAL(kind=dp) :: dab
280  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dummy1
281  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dummy2, waux_mat
282  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dwaux_mat
283  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: scon_fba, scon_fbb, scon_oba, scon_obb
284  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix
285 
286  dab = sqrt(sum(rab*rab))
287  nba = obasa%nsgf
288  nbb = obasb%nsgf
289  nfa = fbasa%nsgf
290  nfb = fbasb%nsgf
291 
292  IF (iatom == jatom .AND. dab < lri_env%delta) THEN
293  ! nothing to do
294  ELSE
295  IF (lrii%calc_force_pair) THEN
296  IF (lri_env%use_shg_integrals) THEN
297  scon_oba => lri_env%bas_prop(ikind)%scon_orb
298  scon_obb => lri_env%bas_prop(jkind)%scon_orb
299  scon_fba => lri_env%bas_prop(ikind)%scon_ri
300  scon_fbb => lri_env%bas_prop(jkind)%scon_ri
301  scona_mix => lri_env%bas_prop(ikind)%scon_mix
302  sconb_mix => lri_env%bas_prop(jkind)%scon_mix
303  oba_index => lri_env%bas_prop(ikind)%orb_index
304  fba_index => lri_env%bas_prop(ikind)%ri_index
305  obb_index => lri_env%bas_prop(jkind)%orb_index
306  fbb_index => lri_env%bas_prop(jkind)%ri_index
307 
308  CALL lri_precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, &
309  waux_mat, dwaux_mat, .true.)
310  !*** (fa,fb)
311  IF (ALLOCATED(lridint%dsabint)) THEN
312  CALL int_overlap_ab_shg_low(dummy1, lridint%dsabint, rab, fbasa, fbasb, scon_fba, scon_fbb, &
313  waux_mat, dwaux_mat, .false., .true., contraction_high=.false.)
314  END IF
315  !*** (a,b)
316  IF (ALLOCATED(lridint%dsooint)) THEN
317  CALL int_overlap_ab_shg_low(dummy1, lridint%dsooint, rab, obasa, obasb, scon_oba, scon_obb, &
318  waux_mat, dwaux_mat, .false., .true., contraction_high=.true.)
319  END IF
320  !*** (a,b,fa)
321  IF (ALLOCATED(lridint%dabaint)) THEN
322  CALL int_overlap_aba_shg_low(dummy2, lridint%dabaint, rab, obasa, obasb, fbasa, &
323  scon_obb, scona_mix, oba_index, fba_index, &
324  lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
325  lri_env%cg_shg%ncg_none0, &
326  waux_mat, dwaux_mat, .false., .true.)
327  END IF
328  !*** (a,b,fb)
329  IF (ALLOCATED(lridint%dabbint)) THEN
330  IF (ikind == jkind) THEN
331  cpassert(ALLOCATED(lridint%dabaint))
332  CALL get_abb_same_kind(dummy2, lridint%dabbint, dummy2, lridint%dabaint, &
333  rab, obasa, fbasa, .false., .true.)
334  ELSE
335  CALL int_overlap_abb_shg_low(dummy2, lridint%dabbint, rab, obasa, obasb, fbasb, &
336  scon_oba, sconb_mix, obb_index, fbb_index, &
337  lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
338  lri_env%cg_shg%ncg_none0, &
339  waux_mat, dwaux_mat, .false., .true.)
340  END IF
341  END IF
342  DEALLOCATE (waux_mat, dwaux_mat)
343 
344  ELSE
345 
346  !*** (fa,fb)
347  IF (ALLOCATED(lridint%dsabint)) THEN
348  ALLOCATE (dummy1(nfa, nfb))
349  CALL int_overlap_ab_os(dummy1, lridint%dsabint, rab, fbasa, fbasb, &
350  .true., lri_env%debug, lrii%dmax_ab)
351  DEALLOCATE (dummy1)
352  END IF
353  !*** (a,b)
354  IF (ALLOCATED(lridint%dsooint)) THEN
355  ALLOCATE (dummy1(nba, nbb))
356  CALL int_overlap_ab_os(dummy1, lridint%dsooint, rab, obasa, obasb, &
357  .true., lri_env%debug, lrii%dmax_oo)
358  DEALLOCATE (dummy1)
359  END IF
360  !*** (a,b,fa)
361  IF (ALLOCATED(lridint%dabaint)) THEN
362  ALLOCATE (dummy2(nba, nbb, nfa))
363  CALL int_overlap_aba_os(dummy2, lridint%dabaint, rab, obasa, obasb, fbasa, &
364  .true., lri_env%debug, lrii%dmax_aba)
365  DEALLOCATE (dummy2)
366  END IF
367  !*** (a,b,fb)
368  IF (ALLOCATED(lridint%dabbint)) THEN
369  ALLOCATE (dummy2(nba, nbb, nfb))
370  CALL int_overlap_abb_os(dummy2, lridint%dabbint, rab, obasa, obasb, fbasb, &
371  .true., lri_env%debug, lrii%dmax_abb)
372  DEALLOCATE (dummy2)
373  END IF
374  END IF
375  END IF
376  END IF
377 
378  END SUBROUTINE lri_dint2
379 
380 ! **************************************************************************************************
381 !> \brief calculates the lri intergrals according to the Obara-Saika (OS)
382 !> scheme
383 !> \param lri_env ...
384 !> \param lrii ...
385 !> \param rab distance vector
386 !> \param obasa orb basis on center A
387 !> \param obasb orb basis on center B
388 !> \param fbasa aux basis on center A
389 !> \param fbasb aux basis on center B
390 !> \param iatom index atom A
391 !> \param jatom index atom B
392 !> \param ikind kind atom A
393 !> \param calculate_forces ...
394 ! **************************************************************************************************
395  SUBROUTINE lri_int_os(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
396  iatom, jatom, ikind, calculate_forces)
397 
398  TYPE(lri_environment_type), POINTER :: lri_env
399  TYPE(lri_int_type), POINTER :: lrii
400  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
401  TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
402  INTEGER, INTENT(IN) :: iatom, jatom, ikind
403  LOGICAL, INTENT(IN) :: calculate_forces
404 
405  CHARACTER(LEN=*), PARAMETER :: routinen = 'lri_int_os'
406 
407  INTEGER :: handle, nba, nbb, nfa, nfb
408  LOGICAL :: do_force
409  REAL(kind=dp) :: dab
410 
411  CALL timeset(routinen, handle)
412 
413  dab = sqrt(sum(rab*rab))
414  nba = obasa%nsgf
415  nbb = obasb%nsgf
416  nfa = fbasa%nsgf
417  nfb = fbasb%nsgf
418 
419  !*** calculate overlap integrals; for iatom=jatom this is the self-overlap
420  IF (iatom == jatom .AND. dab < lri_env%delta) THEN
421  !*** (fa,fa)
422  lrii%sab(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp(1:nfa, 1:nfa)
423  !*** (a,a)
424  lrii%soo(1:nba, 1:nba) = lri_env%bas_prop(ikind)%orb_ovlp(1:nba, 1:nba)
425  !*** (a,a,fa)
426  CALL int_overlap_aba_os(lrii%abaint, rab=rab, oba=obasa, obb=obasb, &
427  fba=fbasa, calculate_forces=.false., debug=lri_env%debug, &
428  dmax=lrii%dmax_aba)
429  IF (calculate_forces) THEN
430  lrii%dsab = 0._dp
431  lrii%dsoo = 0._dp
432  lrii%dabdaint = 0.0_dp
433  lrii%dabbint = 0.0_dp
434  END IF
435  ELSE
436  IF (calculate_forces) THEN
437  do_force = lrii%calc_force_pair
438  ELSE
439  do_force = .false.
440  END IF
441  !*** (fa,fb)
442  CALL int_overlap_ab_os(lrii%sab, lrii%dsab, rab, fbasa, fbasb, &
443  do_force, lri_env%debug, lrii%dmax_ab)
444  !*** (a,b)
445  CALL int_overlap_ab_os(lrii%soo, lrii%dsoo, rab, obasa, obasb, &
446  do_force, lri_env%debug, lrii%dmax_oo)
447  !*** (a,b,fa)
448  CALL int_overlap_aba_os(lrii%abaint, lrii%dabdaint, rab, obasa, obasb, fbasa, &
449  do_force, lri_env%debug, lrii%dmax_aba)
450  !*** (a,b,fb)
451  CALL int_overlap_abb_os(lrii%abbint, lrii%dabbint, rab, obasa, obasb, fbasb, &
452  do_force, lri_env%debug, lrii%dmax_abb)
453  END IF
454 
455  CALL timestop(handle)
456 
457  END SUBROUTINE lri_int_os
458 
459 ! **************************************************************************************************
460 !> \brief calcuates the lri integrals using solid harmonic Gaussians
461 !> \param lri_env ...
462 !> \param lrii ...
463 !> \param rab distance vector
464 !> \param obasa orb basis on A
465 !> \param obasb orb basis on B
466 !> \param fbasa aux basis on A
467 !> \param fbasb aux basis on B
468 !> \param iatom index atom A
469 !> \param jatom index atom B
470 !> \param ikind kind atom A
471 !> \param jkind kind atom B
472 !> \param calculate_forces ...
473 ! **************************************************************************************************
474  SUBROUTINE lri_int_shg(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
475  iatom, jatom, ikind, jkind, calculate_forces)
476 
477  TYPE(lri_environment_type), POINTER :: lri_env
478  TYPE(lri_int_type), POINTER :: lrii
479  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
480  TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
481  INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
482  LOGICAL, INTENT(IN) :: calculate_forces
483 
484  CHARACTER(LEN=*), PARAMETER :: routinen = 'lri_int_shg'
485 
486  INTEGER :: handle, nba, nbb, nfa, nfb
487  INTEGER, DIMENSION(:, :, :), POINTER :: fba_index, fbb_index, oba_index, &
488  obb_index
489  LOGICAL :: do_force
490  REAL(kind=dp) :: dab
491  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: waux_mat
492  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dwaux_mat
493  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: scon_fba, scon_fbb, scon_oba, scon_obb
494  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix
495 
496  CALL timeset(routinen, handle)
497  NULLIFY (scon_oba, scon_obb, scon_fba, scon_fbb, scona_mix, sconb_mix, &
498  oba_index, obb_index, fba_index, fbb_index)
499  dab = sqrt(sum(rab*rab))
500  nba = obasa%nsgf
501  nbb = obasb%nsgf
502  nfa = fbasa%nsgf
503  nfb = fbasb%nsgf
504 
505  !*** calculate overlap integrals; for iatom=jatom this is the self-overlap
506  IF (iatom == jatom .AND. dab < lri_env%delta) THEN
507  !*** (fa,fa)
508  lrii%sab(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp(1:nfa, 1:nfa)
509  !*** (a,a)
510  lrii%soo(1:nba, 1:nba) = lri_env%bas_prop(ikind)%orb_ovlp(1:nba, 1:nba)
511  !*** (a,a,fa)
512  lrii%abaint(1:nba, 1:nba, 1:nfa) = lri_env%bas_prop(ikind)%ovlp3
513  IF (calculate_forces) THEN
514  lrii%dsab = 0._dp
515  lrii%dsoo = 0._dp
516  lrii%dabdaint = 0.0_dp
517  lrii%dabbint = 0.0_dp
518  END IF
519  ELSE
520  IF (calculate_forces) THEN
521  do_force = lrii%calc_force_pair
522  ELSE
523  do_force = .false.
524  END IF
525  scon_oba => lri_env%bas_prop(ikind)%scon_orb
526  scon_obb => lri_env%bas_prop(jkind)%scon_orb
527  scon_fba => lri_env%bas_prop(ikind)%scon_ri
528  scon_fbb => lri_env%bas_prop(jkind)%scon_ri
529  scona_mix => lri_env%bas_prop(ikind)%scon_mix
530  sconb_mix => lri_env%bas_prop(jkind)%scon_mix
531  oba_index => lri_env%bas_prop(ikind)%orb_index
532  fba_index => lri_env%bas_prop(ikind)%ri_index
533  obb_index => lri_env%bas_prop(jkind)%orb_index
534  fbb_index => lri_env%bas_prop(jkind)%ri_index
535  CALL lri_precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, waux_mat, dwaux_mat, &
536  do_force)
537  !*** (fa,fb)
538  CALL int_overlap_ab_shg_low(lrii%sab, lrii%dsab, rab, fbasa, fbasb, scon_fba, scon_fbb, &
539  waux_mat, dwaux_mat, .true., do_force, contraction_high=.false.)
540  !*** (a,b)
541  CALL int_overlap_ab_shg_low(lrii%soo, lrii%dsoo, rab, obasa, obasb, scon_oba, scon_obb, &
542  waux_mat, dwaux_mat, .true., do_force, contraction_high=.true.)
543  !*** (a,b,fa)
544  CALL int_overlap_aba_shg_low(lrii%abaint, lrii%dabdaint, rab, obasa, obasb, fbasa, &
545  scon_obb, scona_mix, oba_index, fba_index, &
546  lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
547  lri_env%cg_shg%ncg_none0, &
548  waux_mat, dwaux_mat, .true., do_force)
549  !*** (a,b,fb)
550  IF (ikind == jkind) THEN
551  CALL get_abb_same_kind(lrii%abbint, lrii%dabbint, lrii%abaint, lrii%dabdaint, &
552  rab, obasa, fbasa, .true., do_force)
553  ELSE
554  CALL int_overlap_abb_shg_low(lrii%abbint, lrii%dabbint, rab, obasa, obasb, fbasb, &
555  scon_oba, sconb_mix, obb_index, fbb_index, &
556  lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
557  lri_env%cg_shg%ncg_none0, &
558  waux_mat, dwaux_mat, .true., do_force)
559  END IF
560 
561  DEALLOCATE (waux_mat, dwaux_mat)
562  END IF
563 
564  CALL timestop(handle)
565 
566  END SUBROUTINE lri_int_shg
567 
568 ! **************************************************************************************************
569 !> \brief ...
570 !> \param lriint ...
571 !> \param lridint ...
572 !> \param nba ...
573 !> \param nbb ...
574 !> \param nfa ...
575 !> \param nfb ...
576 !> \param skip_sab ...
577 !> \param skip_soo ...
578 !> \param skip_aba ...
579 !> \param skip_abb ...
580 !> \param skip_dsab ...
581 !> \param skip_dsoo ...
582 !> \param skip_daba ...
583 !> \param skip_dabb ...
584 ! **************************************************************************************************
585  SUBROUTINE allocate_int_type(lriint, lridint, nba, nbb, nfa, nfb, &
586  skip_sab, skip_soo, skip_aba, skip_abb, &
587  skip_dsab, skip_dsoo, skip_daba, skip_dabb)
588 
589  TYPE(int_type), INTENT(INOUT), OPTIONAL :: lriint
590  TYPE(dint_type), INTENT(INOUT), OPTIONAL :: lridint
591  INTEGER, INTENT(IN) :: nba, nbb, nfa, nfb
592  LOGICAL, INTENT(IN), OPTIONAL :: skip_sab, skip_soo, skip_aba, skip_abb, &
593  skip_dsab, skip_dsoo, skip_daba, &
594  skip_dabb
595 
596  LOGICAL :: do_aba, do_abb, do_daba, do_dabb, &
597  do_dsab, do_dsoo, do_sab, do_soo
598 
599  IF (PRESENT(lriint)) THEN
600  do_sab = .true.
601  IF (PRESENT(skip_sab)) do_sab = .NOT. skip_sab
602  do_soo = .true.
603  IF (PRESENT(skip_soo)) do_soo = .NOT. skip_soo
604  do_aba = .true.
605  IF (PRESENT(skip_aba)) do_aba = .NOT. skip_aba
606  do_abb = .true.
607  IF (PRESENT(skip_abb)) do_abb = .NOT. skip_abb
608  !
609  IF (do_sab) THEN
610  IF (ALLOCATED(lriint%sabint)) DEALLOCATE (lriint%sabint)
611  ALLOCATE (lriint%sabint(nfa, nfb))
612  lriint%sabint = 0.0_dp
613  END IF
614  IF (do_soo) THEN
615  IF (ALLOCATED(lriint%sooint)) DEALLOCATE (lriint%sooint)
616  ALLOCATE (lriint%sooint(nba, nbb))
617  lriint%sooint = 0.0_dp
618  END IF
619  IF (do_aba) THEN
620  IF (ALLOCATED(lriint%abaint)) DEALLOCATE (lriint%abaint)
621  ALLOCATE (lriint%abaint(nba, nbb, nfa))
622  lriint%abaint = 0.0_dp
623  END IF
624  IF (do_abb) THEN
625  IF (ALLOCATED(lriint%abbint)) DEALLOCATE (lriint%abbint)
626  ALLOCATE (lriint%abbint(nba, nbb, nfb))
627  lriint%abbint = 0.0_dp
628  END IF
629  END IF
630  !
631  IF (PRESENT(lridint)) THEN
632  do_dsab = .true.
633  IF (PRESENT(skip_dsab)) do_dsab = .NOT. skip_dsab
634  do_dsoo = .true.
635  IF (PRESENT(skip_dsoo)) do_dsoo = .NOT. skip_dsoo
636  do_daba = .true.
637  IF (PRESENT(skip_daba)) do_daba = .NOT. skip_daba
638  do_dabb = .true.
639  IF (PRESENT(skip_dabb)) do_dabb = .NOT. skip_dabb
640  !
641  IF (do_dsab) THEN
642  IF (ALLOCATED(lridint%dsabint)) DEALLOCATE (lridint%dsabint)
643  ALLOCATE (lridint%dsabint(nfa, nfb, 3))
644  lridint%dsabint = 0.0_dp
645  END IF
646  IF (do_dsoo) THEN
647  IF (ALLOCATED(lridint%dsooint)) DEALLOCATE (lridint%dsooint)
648  ALLOCATE (lridint%dsooint(nba, nbb, 3))
649  lridint%dsooint = 0.0_dp
650  END IF
651  IF (do_daba) THEN
652  IF (ALLOCATED(lridint%dabaint)) DEALLOCATE (lridint%dabaint)
653  ALLOCATE (lridint%dabaint(nba, nbb, nfa, 3))
654  lridint%dabaint = 0.0_dp
655  END IF
656  IF (do_dabb) THEN
657  IF (ALLOCATED(lridint%dabbint)) DEALLOCATE (lridint%dabbint)
658  ALLOCATE (lridint%dabbint(nba, nbb, nfb, 3))
659  lridint%dabbint = 0.0_dp
660  END IF
661  END IF
662 
663  END SUBROUTINE allocate_int_type
664 
665 ! **************************************************************************************************
666 !> \brief ...
667 !> \param lriint ...
668 !> \param lridint ...
669 ! **************************************************************************************************
670  SUBROUTINE deallocate_int_type(lriint, lridint)
671 
672  TYPE(int_type), INTENT(INOUT), OPTIONAL :: lriint
673  TYPE(dint_type), INTENT(INOUT), OPTIONAL :: lridint
674 
675  IF (PRESENT(lriint)) THEN
676  IF (ALLOCATED(lriint%sabint)) DEALLOCATE (lriint%sabint)
677  IF (ALLOCATED(lriint%sooint)) DEALLOCATE (lriint%sooint)
678  IF (ALLOCATED(lriint%abaint)) DEALLOCATE (lriint%abaint)
679  IF (ALLOCATED(lriint%abbint)) DEALLOCATE (lriint%abbint)
680  END IF
681  !
682  IF (PRESENT(lridint)) THEN
683  IF (ALLOCATED(lridint%dsabint)) DEALLOCATE (lridint%dsabint)
684  IF (ALLOCATED(lridint%dsooint)) DEALLOCATE (lridint%dsooint)
685  IF (ALLOCATED(lridint%dabaint)) DEALLOCATE (lridint%dabaint)
686  IF (ALLOCATED(lridint%dabbint)) DEALLOCATE (lridint%dabbint)
687  END IF
688 
689  END SUBROUTINE deallocate_int_type
690 
691 END MODULE lri_integrals
Calculation of contracted, spherical Gaussian integrals using the (OS) integral scheme....
subroutine, public int_overlap_aba_os(abaint, dabdaint, rab, oba, obb, fba, calculate_forces, debug, dmax)
calculate integrals (a,b,fa)
subroutine, public int_overlap_ab_os(sab, dsab, rab, fba, fbb, calculate_forces, debug, dmax)
calculate overlap integrals (a,b)
subroutine, public int_overlap_abb_os(abbint, dabbint, rab, oba, obb, fbb, calculate_forces, debug, dmax)
calculate integrals (a,b,fb)
Calculation of contracted, spherical Gaussian integrals using the solid harmonic Gaussian (SHG) integ...
subroutine, public lri_precalc_angular_shg_part(oba, obb, fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces)
precalculates the angular part of the SHG integrals for the matrices (fa,fb), (a,b),...
subroutine, public int_overlap_ab_shg_low(sab, dsab, rab, fba, fbb, scona_shg, sconb_shg, Waux_mat, dWaux_mat, calculate_ints, calculate_forces, contraction_high)
calculate overlap integrals (a,b); requires angular-dependent part as input
subroutine, public get_abb_same_kind(abbint, dabbint, abaint, dabdaint, rab, oba, fba, calculate_ints, calculate_forces)
obtain integrals (a,b,fb) by symmetry relations from (a,b,fa) if basis sets at a and b are of the sam...
subroutine, public int_overlap_abb_shg_low(abbint, dabbint, rab, oba, obb, fbb, scon_oba, sconb_mix, obb_index, fbb_index, cg_coeff, cg_none0_list, ncg_none0, Waux_mat, dWaux_mat, calculate_ints, calculate_forces)
calculate integrals (a,b,fb); requires angular-dependent part as input
subroutine, public int_overlap_aba_shg_low(abaint, dabdaint, rab, oba, obb, fba, scon_obb, scona_mix, oba_index, fba_index, cg_coeff, cg_none0_list, ncg_none0, Waux_mat, dWaux_mat, calculate_ints, calculate_forces)
calculate integrals (a,b,fa); requires angular-dependent part as input
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
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_dint(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, iatom, jatom, ikind, jkind, calculate_forces)
...
subroutine, public lri_int2(lri_env, lrii, lriint, rab, obasa, obasb, fbasa, fbasb, iatom, jatom, ikind, jkind)
calcuates the lri integrals using solid harmonic Gaussians
subroutine, public lri_int(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, iatom, jatom, ikind, jkind, calculate_forces)
calcuates the lri integrals using solid harmonic Gaussians
Definition: lri_integrals.F:79
subroutine, public lri_dint2(lri_env, lrii, lridint, rab, obasa, obasb, fbasa, fbasb, iatom, jatom, ikind, jkind)
...