(git:0de0cc2)
generic_shg_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 Calculation of contracted, spherical Gaussian integrals using the solid harmonic
10 !> Gaussian (SHG) integral scheme. Routines for the following two-center integrals:
11 !> i) (a|O(r12)|b) where O(r12) is the overlap, coulomb operator etc.
12 !> ii) (aba) and (abb) s-overlaps
13 !> \par Literature
14 !> T.J. Giese and D. M. York, J. Chem. Phys, 128, 064104 (2008)
15 !> T. Helgaker, P Joergensen, J. Olsen, Molecular Electronic-Structure
16 !> Theory, Wiley
17 !> \par History
18 !> created [05.2016]
19 !> \author Dorothea Golze
20 ! **************************************************************************************************
22  USE basis_set_types, ONLY: gto_basis_set_type
28  USE construct_shg, ONLY: &
32  USE kinds, ONLY: dp
33  USE orbital_pointers, ONLY: nsoset
34  USE s_contract_shg, ONLY: &
38 #include "../base/base_uses.f90"
39 
40  IMPLICIT NONE
41 
42  PRIVATE
43 
44  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals'
45 
51 
52  abstract INTERFACE
53 ! **************************************************************************************************
54 !> \brief Interface for the calculation of integrals over s-functions and their scalar derivatives
55 !> with respect to rab2
56 !> \param la_max ...
57 !> \param npgfa ...
58 !> \param zeta ...
59 !> \param lb_max ...
60 !> \param npgfb ...
61 !> \param zetb ...
62 !> \param omega ...
63 !> \param rab ...
64 !> \param v matrix storing the integrals and scalar derivatives
65 !> \param calculate_forces ...
66 ! **************************************************************************************************
67  SUBROUTINE ab_sint_shg(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
68  USE kinds, ONLY: dp
69  INTEGER, INTENT(IN) :: la_max, npgfa
70  REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
71  INTEGER, INTENT(IN) :: lb_max, npgfb
72  REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
73  REAL(KIND=dp), INTENT(IN) :: omega
74  REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
75  REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: v
76  LOGICAL, INTENT(IN) :: calculate_forces
77 
78  END SUBROUTINE ab_sint_shg
79  END INTERFACE
80 
81 ! **************************************************************************************************
82 
83 CONTAINS
84 
85 ! **************************************************************************************************
86 !> \brief Calcululates the two-center integrals of the type (a|O(r12)|b) using the SHG scheme
87 !> \param r12_operator the integral operator, which depends on r12=|r1-r2|
88 !> \param vab integral matrix of spherical contracted Gaussian functions
89 !> \param dvab derivative of the integrals
90 !> \param rab distance vector between center A and B
91 !> \param fba basis at center A
92 !> \param fbb basis at center B
93 !> \param scona_shg SHG contraction matrix for A
94 !> \param sconb_shg SHG contraction matrix for B
95 !> \param omega parameter in the operator
96 !> \param calculate_forces ...
97 ! **************************************************************************************************
98  SUBROUTINE int_operators_r12_ab_shg(r12_operator, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
99  omega, calculate_forces)
100 
101  INTEGER, INTENT(IN) :: r12_operator
102  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
103  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
104  OPTIONAL :: dvab
105  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
106  TYPE(gto_basis_set_type), POINTER :: fba, fbb
107  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: scona_shg, sconb_shg
108  REAL(kind=dp), INTENT(IN), OPTIONAL :: omega
109  LOGICAL, INTENT(IN) :: calculate_forces
110 
111  INTEGER :: la_max, lb_max
112  REAL(kind=dp) :: my_omega
113  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: waux_mat
114  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dwaux_mat
115 
116  PROCEDURE(ab_sint_shg), POINTER :: s_operator_ab
117 
118  NULLIFY (s_operator_ab)
119 
120  la_max = maxval(fba%lmax)
121  lb_max = maxval(fbb%lmax)
122 
123  CALL precalc_angular_shg_part(la_max, lb_max, rab, waux_mat, dwaux_mat, calculate_forces)
124  my_omega = 1.0_dp
125 
126  SELECT CASE (r12_operator)
127  CASE (operator_coulomb)
128  s_operator_ab => s_coulomb_ab
129  CASE (operator_verf)
130  s_operator_ab => s_verf_ab
131  IF (PRESENT(omega)) my_omega = omega
132  CASE (operator_verfc)
133  s_operator_ab => s_verfc_ab
134  IF (PRESENT(omega)) my_omega = omega
135  CASE (operator_vgauss)
136  s_operator_ab => s_vgauss_ab
137  IF (PRESENT(omega)) my_omega = omega
138  CASE (operator_gauss)
139  s_operator_ab => s_gauss_ab
140  IF (PRESENT(omega)) my_omega = omega
141  CASE DEFAULT
142  cpabort("Operator not available")
143  END SELECT
144 
145  CALL int_operator_ab_shg_low(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
146  my_omega, waux_mat, dwaux_mat, calculate_forces)
147 
148  DEALLOCATE (waux_mat, dwaux_mat)
149 
150  END SUBROUTINE int_operators_r12_ab_shg
151 
152 ! **************************************************************************************************
153 !> \brief calculate overlap integrals (a,b)
154 !> \param vab integral (a,b)
155 !> \param dvab derivative of sab
156 !> \param rab distance vector
157 !> \param fba basis at center A
158 !> \param fbb basis at center B
159 !> \param scona_shg contraction matrix A
160 !> \param sconb_shg contraxtion matrix B
161 !> \param calculate_forces ...
162 ! **************************************************************************************************
163  SUBROUTINE int_overlap_ab_shg(vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
164  calculate_forces)
165 
166  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
167  INTENT(INOUT) :: vab
168  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
169  INTENT(INOUT) :: dvab
170  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
171  TYPE(gto_basis_set_type), POINTER :: fba, fbb
172  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg
173  LOGICAL, INTENT(IN) :: calculate_forces
174 
175  INTEGER :: la_max, lb_max
176  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: waux_mat
177  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dwaux_mat
178 
179  la_max = maxval(fba%lmax)
180  lb_max = maxval(fbb%lmax)
181 
182  CALL precalc_angular_shg_part(la_max, lb_max, rab, waux_mat, dwaux_mat, calculate_forces)
183 
184  CALL int_overlap_ab_shg_low(vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
185  waux_mat, dwaux_mat, .true., calculate_forces, contraction_high=.true.)
186 
187  DEALLOCATE (waux_mat, dwaux_mat)
188 
189  END SUBROUTINE int_overlap_ab_shg
190 
191 ! **************************************************************************************************
192 !> \brief Calcululates the two-center integrals of the type (a|(r-Ra)^(2m)|b) using the SHG scheme
193 !> \param vab integral matrix of spherical contracted Gaussian functions
194 !> \param dvab derivative of the integrals
195 !> \param rab distance vector between center A and B
196 !> \param fba basis at center A
197 !> \param fbb basis at center B
198 !> \param scon_ra2m contraction matrix for A including the combinatorial factors
199 !> \param sconb_shg SHG contraction matrix for B
200 !> \param m exponent in (r-Ra)^(2m) operator
201 !> \param calculate_forces ...
202 ! **************************************************************************************************
203  SUBROUTINE int_ra2m_ab_shg(vab, dvab, rab, fba, fbb, scon_ra2m, sconb_shg, &
204  m, calculate_forces)
205 
206  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
207  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dvab
208  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
209  TYPE(gto_basis_set_type), POINTER :: fba, fbb
210  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scon_ra2m
211  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: sconb_shg
212  INTEGER, INTENT(IN) :: m
213  LOGICAL, INTENT(IN) :: calculate_forces
214 
215  INTEGER :: la_max, lb_max
216  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: waux_mat
217  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dwaux_mat
218 
219  la_max = maxval(fba%lmax)
220  lb_max = maxval(fbb%lmax)
221 
222  CALL precalc_angular_shg_part(la_max, lb_max, rab, waux_mat, dwaux_mat, calculate_forces)
223  CALL int_ra2m_ab_shg_low(vab, dvab, rab, fba, fbb, sconb_shg, scon_ra2m, &
224  m, waux_mat, dwaux_mat, calculate_forces)
225 
226  DEALLOCATE (waux_mat, dwaux_mat)
227 
228  END SUBROUTINE int_ra2m_ab_shg
229 
230 ! **************************************************************************************************
231 !> \brief calculate integrals (a,b,fa)
232 !> \param saba integral [aba]
233 !> \param dsaba derivative of [aba]
234 !> \param rab distance vector between A and B
235 !> \param oba orbital basis at center A
236 !> \param obb orbital basis at center B
237 !> \param fba auxiliary basis set at center A
238 !> \param scon_obb contraction matrix for orb bas on B
239 !> \param scona_mix mixed contraction matrix orb + ri basis on A
240 !> \param oba_index orbital basis index for scona_mix
241 !> \param fba_index ri basis index for scona_mix
242 !> \param cg_coeff Clebsch-Gordon coefficients
243 !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
244 !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
245 !> \param calculate_forces ...
246 ! **************************************************************************************************
247  SUBROUTINE int_overlap_aba_shg(saba, dsaba, rab, oba, obb, fba, scon_obb, &
248  scona_mix, oba_index, fba_index, &
249  cg_coeff, cg_none0_list, ncg_none0, &
250  calculate_forces)
251 
252  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
253  INTENT(INOUT) :: saba
254  REAL(kind=dp), ALLOCATABLE, &
255  DIMENSION(:, :, :, :), INTENT(INOUT) :: dsaba
256  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
257  TYPE(gto_basis_set_type), POINTER :: oba, obb, fba
258  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_obb
259  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix
260  INTEGER, DIMENSION(:, :, :), INTENT(IN) :: oba_index, fba_index
261  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
262  INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
263  INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
264  LOGICAL, INTENT(IN) :: calculate_forces
265 
266  INTEGER :: laa_max, lb_max
267  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: waux_mat
268  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dwaux_mat
269 
270  laa_max = maxval(oba%lmax) + maxval(fba%lmax)
271  lb_max = maxval(obb%lmax)
272 
273  saba = 0.0_dp
274  IF (calculate_forces) dsaba = 0.0_dp
275  CALL precalc_angular_shg_part(laa_max, lb_max, rab, waux_mat, dwaux_mat, &
276  calculate_forces)
277  CALL int_overlap_aba_shg_low(saba, dsaba, rab, oba, obb, fba, &
278  scon_obb, scona_mix, oba_index, fba_index, &
279  cg_coeff, cg_none0_list, ncg_none0, &
280  waux_mat, dwaux_mat, .true., calculate_forces)
281 
282  DEALLOCATE (waux_mat, dwaux_mat)
283 
284  END SUBROUTINE int_overlap_aba_shg
285 
286 ! **************************************************************************************************
287 !> \brief calculate integrals (a,b,fb)
288 !> \param sabb integral [abb]
289 !> \param dsabb derivative of [abb]
290 !> \param rab distance vector between A and B
291 !> \param oba orbital basis at center A
292 !> \param obb orbital basis at center B
293 !> \param fbb auxiliary basis set at center B
294 !> \param scon_oba contraction matrix for orb bas on A
295 !> \param sconb_mix mixed contraction matrix orb + ri basis on B
296 !> \param obb_index orbital basis index for sconb_mix
297 !> \param fbb_index ri basis index for sconb_mix
298 !> \param cg_coeff Clebsch-Gordon coefficients
299 !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
300 !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
301 !> \param calculate_forces ...
302 ! **************************************************************************************************
303  SUBROUTINE int_overlap_abb_shg(sabb, dsabb, rab, oba, obb, fbb, scon_oba, &
304  sconb_mix, obb_index, fbb_index, &
305  cg_coeff, cg_none0_list, ncg_none0, &
306  calculate_forces)
307 
308  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
309  INTENT(INOUT) :: sabb
310  REAL(kind=dp), ALLOCATABLE, &
311  DIMENSION(:, :, :, :), INTENT(INOUT) :: dsabb
312  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
313  TYPE(gto_basis_set_type), POINTER :: oba, obb, fbb
314  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_oba
315  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix
316  INTEGER, DIMENSION(:, :, :), INTENT(IN) :: obb_index, fbb_index
317  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
318  INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
319  INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
320  LOGICAL, INTENT(IN) :: calculate_forces
321 
322  INTEGER :: la_max, lbb_max
323  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: waux_mat
324  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dwaux_mat
325 
326  la_max = maxval(oba%lmax)
327  lbb_max = maxval(obb%lmax) + maxval(fbb%lmax)
328 
329  sabb = 0.0_dp
330  IF (calculate_forces) dsabb = 0.0_dp
331  CALL precalc_angular_shg_part(lbb_max, la_max, rab, waux_mat, dwaux_mat, &
332  calculate_forces)
333  CALL int_overlap_abb_shg_low(sabb, dsabb, rab, oba, obb, fbb, &
334  scon_oba, sconb_mix, obb_index, fbb_index, &
335  cg_coeff, cg_none0_list, ncg_none0, &
336  waux_mat, dwaux_mat, .true., calculate_forces)
337 
338  DEALLOCATE (waux_mat, dwaux_mat)
339 
340  END SUBROUTINE int_overlap_abb_shg
341 
342 ! **************************************************************************************************
343 !> \brief precalculates the angular part of the SHG integrals
344 !> \param la_max ...
345 !> \param lb_max ...
346 !> \param rab distance vector between a and b
347 !> \param Waux_mat W matrix that contains the angular-dependent part
348 !> \param dWaux_mat derivative of the W matrix
349 !> \param calculate_forces ...
350 ! **************************************************************************************************
351  SUBROUTINE precalc_angular_shg_part(la_max, lb_max, rab, Waux_mat, dWaux_mat, calculate_forces)
352 
353  INTEGER, INTENT(IN) :: la_max, lb_max
354  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
355  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
356  INTENT(OUT) :: waux_mat
357  REAL(kind=dp), ALLOCATABLE, &
358  DIMENSION(:, :, :, :), INTENT(OUT) :: dwaux_mat
359  LOGICAL, INTENT(IN) :: calculate_forces
360 
361  INTEGER :: lmax, mdim(3)
362  INTEGER, DIMENSION(:), POINTER :: la_max_all
363  REAL(kind=dp) :: rab2
364  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rc, rs
365 
366  NULLIFY (la_max_all)
367  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
368  lmax = max(la_max, lb_max)
369 
370  ALLOCATE (la_max_all(0:lb_max))
371  ALLOCATE (rc(0:lmax, -2*lmax:2*lmax), rs(0:lmax, -2*lmax:2*lmax))
372  rc = 0._dp
373  rs = 0._dp
374  mdim(1) = min(la_max, lb_max) + 1
375  mdim(2) = nsoset(la_max) + 1
376  mdim(3) = nsoset(lb_max) + 1
377  ALLOCATE (waux_mat(mdim(1), mdim(2), mdim(3)))
378  ALLOCATE (dwaux_mat(3, mdim(1), mdim(2), mdim(3)))
379 
380  la_max_all(0:lb_max) = la_max
381  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
382  CALL get_real_scaled_solid_harmonic(rc, rs, lmax, -rab, rab2)
383  CALL get_w_matrix(la_max_all, lb_max, lmax, rc, rs, waux_mat)
384  IF (calculate_forces) THEN
385  CALL get_dw_matrix(la_max_all, lb_max, waux_mat, dwaux_mat)
386  END IF
387 
388  DEALLOCATE (rc, rs, la_max_all)
389 
390  END SUBROUTINE precalc_angular_shg_part
391 
392 ! **************************************************************************************************
393 !> \brief calculate integrals (a|O(r12)|b)
394 !> \param s_operator_ab procedure pointer for the respective operator. The integral evaluation
395 !> differs only in the calculation of the [s|O(r12)|s] integrals and their scalar
396 !> derivatives.
397 !> \param vab integral matrix of spherical contracted Gaussian functions
398 !> \param dvab derivative of the integrals
399 !> \param rab distance vector between center A and B
400 !> \param fba basis at center A
401 !> \param fbb basis at center B
402 !> \param scona_shg SHG contraction matrix for A
403 !> \param sconb_shg SHG contraction matrix for B
404 !> \param omega parameter in the operator
405 !> \param Waux_mat W matrix that contains the angular-dependent part
406 !> \param dWaux_mat derivative of the W matrix
407 !> \param calculate_forces ...
408 ! **************************************************************************************************
409  SUBROUTINE int_operator_ab_shg_low(s_operator_ab, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, &
410  omega, Waux_mat, dWaux_mat, calculate_forces)
411 
412  PROCEDURE(ab_sint_shg), POINTER :: s_operator_ab
413  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
414  REAL(kind=dp), DIMENSION(:, :, :), OPTIONAL, INTENT(INOUT) :: dvab
415  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
416  TYPE(gto_basis_set_type), POINTER :: fba, fbb
417  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: scona_shg, sconb_shg
418  REAL(kind=dp), INTENT(IN) :: omega
419  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: waux_mat
420  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dwaux_mat
421  LOGICAL, INTENT(IN) :: calculate_forces
422 
423  INTEGER :: iset, jset, la_max_set, lb_max_set, ndev, nds, nds_max, npgfa_set, &
424  npgfb_set, nseta, nsetb, nsgfa_set, nsgfb_set, nshella_set, nshellb_set
425  INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
426  nsgfb, nshella, nshellb
427  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb
428  REAL(kind=dp) :: dab
429  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
430  REAL(kind=dp), DIMENSION(:, :), POINTER :: zeta, zetb
431  REAL(kind=dp), DIMENSION(:, :, :), ALLOCATABLE :: swork, swork_cont
432 
433  NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
434  set_radius_b, zeta, zetb)
435 
436  ! basis ikind
437  first_sgfa => fba%first_sgf
438  la_max => fba%lmax
439  la => fba%l
440  npgfa => fba%npgf
441  nsgfa => fba%nsgf_set
442  nseta = fba%nset
443  set_radius_a => fba%set_radius
444  zeta => fba%zet
445  nshella => fba%nshell
446  ! basis jkind
447  first_sgfb => fbb%first_sgf
448  lb_max => fbb%lmax
449  lb => fbb%l
450  npgfb => fbb%npgf
451  nsgfb => fbb%nsgf_set
452  nsetb = fbb%nset
453  set_radius_b => fbb%set_radius
454  zetb => fbb%zet
455  nshellb => fbb%nshell
456 
457  dab = sqrt(sum(rab**2))
458 
459  la_max_set = maxval(la_max)
460  lb_max_set = maxval(lb_max)
461 
462  ! allocate some work matrices
463  npgfa_set = maxval(npgfa)
464  npgfb_set = maxval(npgfb)
465  nshella_set = maxval(nshella)
466  nshellb_set = maxval(nshellb)
467  nsgfa_set = maxval(nsgfa)
468  nsgfb_set = maxval(nsgfb)
469  ndev = 0
470  IF (calculate_forces) ndev = 1
471  nds_max = la_max_set + lb_max_set + ndev + 1
472  ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
473  ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
474 
475  vab = 0.0_dp
476  IF (calculate_forces) dvab = 0.0_dp
477 
478  DO iset = 1, nseta
479 
480  DO jset = 1, nsetb
481 
482  nds = la_max(iset) + lb_max(jset) + ndev + 1
483  swork(1:npgfa(iset), 1:npgfb(jset), 1:nds) = 0.0_dp
484  CALL s_operator_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
485  lb_max(jset), npgfb(jset), zetb(:, jset), &
486  omega, rab, swork, calculate_forces)
487  CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), &
488  scona_shg(1:npgfa(iset), 1:nshella(iset), iset), &
489  npgfb(jset), nshellb(jset), &
490  sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
491  nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), &
492  swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)))
493  CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
494  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
495  swork_cont, waux_mat, vab)
496  IF (calculate_forces) THEN
497  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
498  CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
499  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
500  -rab, swork_cont, waux_mat, dwaux_mat, dvab)
501  END IF
502  END DO
503  END DO
504 
505  DEALLOCATE (swork, swork_cont)
506 
507  END SUBROUTINE int_operator_ab_shg_low
508 
509 ! **************************************************************************************************
510 !> \brief calculate overlap integrals (a,b); requires angular-dependent part as input
511 !> \param sab integral (a,b)
512 !> \param dsab derivative of sab
513 !> \param rab distance vector
514 !> \param fba basis at center A
515 !> \param fbb basis at center B
516 !> \param scona_shg contraction matrix A
517 !> \param sconb_shg contraxtion matrix B
518 !> \param Waux_mat W matrix that contains the angular-dependent part
519 !> \param dWaux_mat derivative of the W matrix
520 !> \param calculate_ints ...
521 !> \param calculate_forces ...
522 !> \param contraction_high ...
523 ! **************************************************************************************************
524  SUBROUTINE int_overlap_ab_shg_low(sab, dsab, rab, fba, fbb, scona_shg, sconb_shg, Waux_mat, dWaux_mat, &
525  calculate_ints, calculate_forces, contraction_high)
526 
527  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
528  INTENT(INOUT) :: sab
529  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
530  INTENT(INOUT) :: dsab
531  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
532  TYPE(gto_basis_set_type), POINTER :: fba, fbb
533  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: scona_shg, sconb_shg, waux_mat
534  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dwaux_mat
535  LOGICAL, INTENT(IN) :: calculate_ints, calculate_forces
536  LOGICAL, INTENT(IN), OPTIONAL :: contraction_high
537 
538  INTEGER :: iset, jset, la_max_set, lb_max_set, &
539  ndev, nds, nds_max, npgfa_set, &
540  npgfb_set, nseta, nsetb, nshella_set, &
541  nshellb_set
542  INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nshella, &
543  nshellb
544  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb
545  LOGICAL :: my_contraction_high
546  REAL(kind=dp) :: dab
547  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: swork, swork_cont
548  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
549  REAL(kind=dp), DIMENSION(:, :), POINTER :: zeta, zetb
550 
551  NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, set_radius_a, &
552  set_radius_b, zeta, zetb)
553 
554  ! basis ikind
555  first_sgfa => fba%first_sgf
556  la_max => fba%lmax
557  la => fba%l
558  npgfa => fba%npgf
559  nseta = fba%nset
560  set_radius_a => fba%set_radius
561  zeta => fba%zet
562  nshella => fba%nshell
563  ! basis jkind
564  first_sgfb => fbb%first_sgf
565  lb_max => fbb%lmax
566  lb => fbb%l
567  npgfb => fbb%npgf
568  nsetb = fbb%nset
569  set_radius_b => fbb%set_radius
570  zetb => fbb%zet
571  nshellb => fbb%nshell
572 
573  dab = sqrt(sum(rab**2))
574 
575  la_max_set = maxval(la_max)
576  lb_max_set = maxval(lb_max)
577 
578  ! allocate some work matrices
579  npgfa_set = maxval(npgfa)
580  npgfb_set = maxval(npgfb)
581  nshella_set = maxval(nshella)
582  nshellb_set = maxval(nshellb)
583  ndev = 0
584  IF (calculate_forces) ndev = 1
585  nds_max = la_max_set + lb_max_set + ndev + 1
586  ALLOCATE (swork(npgfa_set, npgfb_set, nds_max))
587  ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
588 
589  IF (calculate_ints) sab = 0.0_dp
590  IF (calculate_forces) dsab = 0.0_dp
591 
592  my_contraction_high = .true.
593  IF (PRESENT(contraction_high)) my_contraction_high = contraction_high
594 
595  DO iset = 1, nseta
596 
597  DO jset = 1, nsetb
598 
599  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
600 
601  nds = la_max(iset) + lb_max(jset) + ndev + 1
602  swork(1:npgfa(iset), 1:npgfb(jset), 1:nds) = 0.0_dp
603  CALL s_overlap_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
604  lb_max(jset), npgfb(jset), zetb(:, jset), &
605  rab, swork, calculate_forces)
606  IF (my_contraction_high) THEN
607  CALL contract_sint_ab_chigh(npgfa(iset), nshella(iset), &
608  scona_shg(1:npgfa(iset), 1:nshella(iset), iset), &
609  npgfb(jset), nshellb(jset), &
610  sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
611  nds, swork(1:npgfa(iset), 1:npgfb(jset), 1:nds), &
612  swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)))
613  ELSE
614  CALL contract_sint_ab_clow(la(:, iset), npgfa(iset), nshella(iset), &
615  scona_shg(:, :, iset), &
616  lb(:, jset), npgfb(jset), nshellb(jset), &
617  sconb_shg(:, :, jset), &
618  swork, swork_cont, calculate_forces)
619  END IF
620  IF (calculate_ints) THEN
621  CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
622  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
623  swork_cont, waux_mat, sab)
624  END IF
625  IF (calculate_forces) THEN
626  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
627  CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
628  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
629  -rab, swork_cont, waux_mat, dwaux_mat, dsab)
630  END IF
631  END DO
632  END DO
633 
634  DEALLOCATE (swork, swork_cont)
635 
636  END SUBROUTINE int_overlap_ab_shg_low
637 
638 ! **************************************************************************************************
639 !> \brief calculate integrals (a|ra^2m)|b); requires angular-dependent part as input
640 !> \param vab integral matrix of spherical contracted Gaussian functions
641 !> \param dvab derivative of the integrals
642 !> \param rab distance vector between center A and B
643 !> \param fba basis at center A
644 !> \param fbb basis at center B
645 !> \param sconb_shg SHG contraction matrix for B
646 !> \param scon_ra2m contraction matrix for A including the combinatorial factors
647 !> \param m exponent in (r-Ra)^(2m) operator
648 !> \param Waux_mat W matrix that contains the angular-dependent part
649 !> \param dWaux_mat derivative of the W matrix
650 !> \param calculate_forces ...
651 ! **************************************************************************************************
652  SUBROUTINE int_ra2m_ab_shg_low(vab, dvab, rab, fba, fbb, sconb_shg, scon_ra2m, m, Waux_mat, dWaux_mat, &
653  calculate_forces)
654 
655  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
656  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dvab
657  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
658  TYPE(gto_basis_set_type), POINTER :: fba, fbb
659  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: sconb_shg
660  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scon_ra2m
661  INTEGER, INTENT(IN) :: m
662  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: waux_mat
663  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dwaux_mat
664  LOGICAL, INTENT(IN) :: calculate_forces
665 
666  INTEGER :: iset, jset, la_max_set, lb_max_set, &
667  ndev, nds, nds_max, npgfa_set, &
668  npgfb_set, nseta, nsetb, nshella_set, &
669  nshellb_set
670  INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb, nsgfa, &
671  nsgfb, nshella, nshellb
672  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, la, lb
673  REAL(kind=dp) :: dab
674  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: swork_cont
675  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: swork
676  REAL(kind=dp), DIMENSION(:, :), POINTER :: zeta, zetb
677 
678  NULLIFY (la_max, lb_max, npgfa, npgfb, first_sgfa, first_sgfb, zeta, zetb)
679 
680  ! basis ikind
681  first_sgfa => fba%first_sgf
682  la_max => fba%lmax
683  la => fba%l
684  npgfa => fba%npgf
685  nsgfa => fba%nsgf_set
686  nseta = fba%nset
687  zeta => fba%zet
688  nshella => fba%nshell
689  ! basis jkind
690  first_sgfb => fbb%first_sgf
691  lb_max => fbb%lmax
692  lb => fbb%l
693  npgfb => fbb%npgf
694  nsgfb => fbb%nsgf_set
695  nsetb = fbb%nset
696  zetb => fbb%zet
697  nshellb => fbb%nshell
698 
699  dab = sqrt(sum(rab**2))
700 
701  la_max_set = maxval(la_max)
702  lb_max_set = maxval(lb_max)
703 
704  ! allocate some work matrices
705  npgfa_set = maxval(npgfa)
706  npgfb_set = maxval(npgfb)
707  nshella_set = maxval(nshella)
708  nshellb_set = maxval(nshellb)
709  ndev = 0
710  IF (calculate_forces) ndev = 1
711  nds_max = la_max_set + lb_max_set + ndev + 1
712  ALLOCATE (swork(npgfa_set, npgfb_set, 1:m + 1, nds_max))
713  ALLOCATE (swork_cont(nds_max, nshella_set, nshellb_set))
714 
715  vab = 0.0_dp
716  IF (calculate_forces) dvab = 0.0_dp
717 
718  DO iset = 1, nseta
719 
720  DO jset = 1, nsetb
721 
722  nds = la_max(iset) + lb_max(jset) + ndev + 1
723  swork(1:npgfa(iset), 1:npgfb(jset), 1:m + 1, 1:nds) = 0.0_dp
724  CALL s_ra2m_ab(la_max(iset), npgfa(iset), zeta(:, iset), &
725  lb_max(jset), npgfb(jset), zetb(:, jset), &
726  m, rab, swork, calculate_forces)
727  CALL contract_s_ra2m_ab(npgfa(iset), nshella(iset), &
728  scon_ra2m(1:npgfa(iset), 1:m + 1, 1:nshella(iset), iset), &
729  npgfb(jset), nshellb(jset), &
730  sconb_shg(1:npgfb(jset), 1:nshellb(jset), jset), &
731  swork(1:npgfa(iset), 1:npgfb(jset), 1:m + 1, 1:nds), &
732  swork_cont(1:nds, 1:nshella(iset), 1:nshellb(jset)), &
733  m, nds)
734  CALL construct_int_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
735  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
736  swork_cont, waux_mat, vab)
737  IF (calculate_forces) THEN
738  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
739  CALL construct_dev_shg_ab(la(:, iset), first_sgfa(:, iset), nshella(iset), &
740  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
741  -rab, swork_cont, waux_mat, dwaux_mat, dvab)
742  END IF
743  END DO
744  END DO
745 
746  DEALLOCATE (swork, swork_cont)
747 
748  END SUBROUTINE int_ra2m_ab_shg_low
749 ! **************************************************************************************************
750 !> \brief calculate integrals (a,b,fb); requires angular-dependent part as input
751 !> \param abbint integral (a,b,fb)
752 !> \param dabbint derivative of abbint
753 !> \param rab distance vector between A and B
754 !> \param oba orbital basis at center A
755 !> \param obb orbital basis at center B
756 !> \param fbb auxiliary basis set at center B
757 !> \param scon_oba contraction matrix for orb bas on A
758 !> \param sconb_mix mixed contraction matrix orb + ri basis on B
759 !> \param obb_index orbital basis index for sconb_mix
760 !> \param fbb_index ri basis index for sconb_mix
761 !> \param cg_coeff Clebsch-Gordon coefficients
762 !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
763 !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
764 !> \param Waux_mat W matrix that contains the angular-dependent part
765 !> \param dWaux_mat derivative of the W matrix
766 !> \param calculate_ints ...
767 !> \param calculate_forces ...
768 ! **************************************************************************************************
769  SUBROUTINE int_overlap_abb_shg_low(abbint, dabbint, rab, oba, obb, fbb, scon_oba, sconb_mix, &
770  obb_index, fbb_index, cg_coeff, cg_none0_list, &
771  ncg_none0, Waux_mat, dWaux_mat, calculate_ints, calculate_forces)
772 
773  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
774  INTENT(INOUT) :: abbint
775  REAL(kind=dp), ALLOCATABLE, &
776  DIMENSION(:, :, :, :), INTENT(INOUT) :: dabbint
777  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
778  TYPE(gto_basis_set_type), POINTER :: oba, obb, fbb
779  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_oba
780  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: sconb_mix
781  INTEGER, DIMENSION(:, :, :), INTENT(IN) :: obb_index, fbb_index
782  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
783  INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
784  INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
785  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: waux_mat
786  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dwaux_mat
787  LOGICAL, INTENT(IN) :: calculate_ints, calculate_forces
788 
789  INTEGER :: iset, jset, kset, la_max_set, lb_max_set, lbb_max, lbb_max_set, lcb_max_set, na, &
790  nb, ncb, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfcb_set, nseta, nsetb, &
791  nsetcb, nshella_set, nshellb_set, nshellcb_set, sgfa, sgfb, sgfcb
792  INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lcb_max, npgfa, npgfb, &
793  npgfcb, nsgfa, nsgfb, nsgfcb, nshella, &
794  nshellb, nshellcb
795  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfcb, la, &
796  lb, lcb
797  REAL(kind=dp) :: dab, rab2
798  REAL(kind=dp), ALLOCATABLE, &
799  DIMENSION(:, :, :, :, :) :: swork, swork_cont
800  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_cb
801  REAL(kind=dp), DIMENSION(:, :), POINTER :: zeta, zetb, zetcb
802 
803  NULLIFY (la_max, lb_max, lcb_max, npgfa, npgfb, npgfcb)
804  NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, &
805  set_radius_cb, zeta, zetb, zetcb)
806 
807  ! basis ikind
808  first_sgfa => oba%first_sgf
809  la_max => oba%lmax
810  la => oba%l
811  nsgfa => oba%nsgf_set
812  npgfa => oba%npgf
813  nshella => oba%nshell
814  nseta = oba%nset
815  set_radius_a => oba%set_radius
816  zeta => oba%zet
817  ! basis jkind
818  first_sgfb => obb%first_sgf
819  lb_max => obb%lmax
820  lb => obb%l
821  nsgfb => obb%nsgf_set
822  npgfb => obb%npgf
823  nshellb => obb%nshell
824  nsetb = obb%nset
825  set_radius_b => obb%set_radius
826  zetb => obb%zet
827 
828  ! basis RI on B
829  first_sgfcb => fbb%first_sgf
830  lcb_max => fbb%lmax
831  lcb => fbb%l
832  nsgfcb => fbb%nsgf_set
833  npgfcb => fbb%npgf
834  nshellcb => fbb%nshell
835  nsetcb = fbb%nset
836  set_radius_cb => fbb%set_radius
837  zetcb => fbb%zet
838 
839  dab = sqrt(sum(rab**2))
840  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
841 
842  la_max_set = maxval(la_max)
843  lb_max_set = maxval(lb_max)
844  lcb_max_set = maxval(lcb_max)
845  npgfa_set = maxval(npgfa)
846  npgfb_set = maxval(npgfb)
847  npgfcb_set = maxval(npgfcb)
848  nshella_set = maxval(nshella)
849  nshellb_set = maxval(nshellb)
850  nshellcb_set = maxval(nshellcb)
851  !*** for forces: derivative+1 in auxiliary vector required
852  ndev = 0
853  IF (calculate_forces) ndev = 1
854 
855  lbb_max_set = lb_max_set + lcb_max_set
856 
857  ! allocate some work storage....
858  nds_max = la_max_set + lbb_max_set + ndev + 1
859  nl_set = int((lbb_max_set)/2)
860  ALLOCATE (swork(npgfa_set, npgfb_set, npgfcb_set, nl_set + 1, nds_max))
861  ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellcb_set))
862 
863  DO iset = 1, nseta
864 
865  DO jset = 1, nsetb
866 
867  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
868 
869  DO kset = 1, nsetcb
870 
871  IF (set_radius_a(iset) + set_radius_cb(kset) < dab) cycle
872 
873  lbb_max = lb_max(jset) + lcb_max(kset)
874  nds = la_max(iset) + lbb_max + ndev + 1
875  nl = int((lbb_max)/2) + 1
876  swork(1:npgfa(iset), 1:npgfb(jset), 1:npgfcb(kset), 1:nl, 1:nds) = 0.0_dp
877  CALL s_overlap_abb(la_max(iset), npgfa(iset), zeta(:, iset), &
878  lb_max(jset), npgfb(jset), zetb(:, jset), &
879  lcb_max(kset), npgfcb(kset), zetcb(:, kset), &
880  rab, swork, calculate_forces)
881 
882  CALL contract_s_overlap_abb(la(:, iset), npgfa(iset), nshella(iset), &
883  scon_oba(1:npgfa(iset), 1:nshella(iset), iset), &
884  lb(:, jset), npgfb(jset), nshellb(jset), &
885  lcb(:, kset), npgfcb(kset), nshellcb(kset), &
886  obb_index(:, :, jset), fbb_index(:, :, kset), sconb_mix, nl, nds, &
887  swork(1:npgfa(iset), 1:npgfb(jset), 1:npgfcb(kset), 1:nl, 1:nds), &
888  swork_cont, calculate_forces)
889 
890  IF (calculate_ints) THEN
891  CALL construct_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
892  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
893  lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
894  cg_coeff, cg_none0_list, &
895  ncg_none0, swork_cont, waux_mat, abbint)
896  END IF
897  IF (calculate_forces) THEN
898  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
899  CALL dev_overlap_shg_abb(la(:, iset), first_sgfa(:, iset), nshella(iset), &
900  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
901  lcb(:, kset), first_sgfcb(:, kset), nshellcb(kset), &
902  cg_coeff, cg_none0_list, ncg_none0, -rab, swork_cont, &
903  waux_mat, dwaux_mat, dabbint)
904  END IF
905  ! max value of integrals in this set triple
906  sgfa = first_sgfa(1, iset)
907  na = sgfa + nsgfa(iset) - 1
908  sgfb = first_sgfb(1, jset)
909  nb = sgfb + nsgfb(jset) - 1
910  sgfcb = first_sgfcb(1, kset)
911  ncb = sgfcb + nsgfcb(kset) - 1
912  END DO
913  END DO
914  END DO
915 
916  DEALLOCATE (swork_cont)
917  DEALLOCATE (swork)
918 
919  END SUBROUTINE int_overlap_abb_shg_low
920 ! **************************************************************************************************
921 !> \brief obtain integrals (a,b,fb) by symmetry relations from (a,b,fa) if basis sets at a and
922 !> b are of the same kind, i.e. a and b are same atom type
923 !> \param abbint integral (a,b,fb)
924 !> \param dabbint derivative of abbint
925 !> \param abaint integral (a,b,fa)
926 !> \param dabdaint derivative of abaint
927 !> \param rab distance vector between A and B
928 !> \param oba orbital basis at center A
929 !> \param fba auxiliary basis set at center A
930 !> \param calculate_ints ...
931 !> \param calculate_forces ...
932 ! **************************************************************************************************
933  SUBROUTINE get_abb_same_kind(abbint, dabbint, abaint, dabdaint, rab, oba, fba, &
934  calculate_ints, calculate_forces)
935 
936  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
937  INTENT(INOUT) :: abbint
938  REAL(kind=dp), ALLOCATABLE, &
939  DIMENSION(:, :, :, :), INTENT(INOUT) :: dabbint
940  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
941  INTENT(INOUT) :: abaint
942  REAL(kind=dp), ALLOCATABLE, &
943  DIMENSION(:, :, :, :), INTENT(INOUT) :: dabdaint
944  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
945  TYPE(gto_basis_set_type), POINTER :: oba, fba
946  LOGICAL, INTENT(IN) :: calculate_ints, calculate_forces
947 
948  INTEGER :: i, iend, iset, isgfa, ishella, istart, jend, jset, jsgfa, jshella, jstart, kend, &
949  kset, ksgfa, kshella, kstart, lai, laj, lak, nseta, nsetca, nsgfa, nsgfca, sgfa_end, &
950  sgfa_start
951  INTEGER, ALLOCATABLE, DIMENSION(:) :: lai_set, laj_set, lak_set
952  INTEGER, DIMENSION(:), POINTER :: nsgfa_set, nsgfca_set, nshella, nshellca
953  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfca, la, lca
954  REAL(kind=dp) :: dab
955  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_ca
956 
957  NULLIFY (nshellca, first_sgfa, first_sgfca, lca, set_radius_a, &
958  set_radius_ca)
959 
960  ! basis ikind
961  first_sgfa => oba%first_sgf
962  set_radius_a => oba%set_radius
963  nseta = oba%nset
964  nsgfa = oba%nsgf
965  nsgfa_set => oba%nsgf_set
966  nshella => oba%nshell
967  la => oba%l
968 
969  ! basis RI
970  first_sgfca => fba%first_sgf
971  set_radius_ca => fba%set_radius
972  nsetca = fba%nset
973  nshellca => fba%nshell
974  nsgfca = fba%nsgf
975  nsgfca_set => fba%nsgf_set
976  lca => fba%l
977 
978  ALLOCATE (lai_set(nsgfa))
979  ALLOCATE (laj_set(nsgfa))
980  ALLOCATE (lak_set(nsgfca))
981 
982  dab = sqrt(sum(rab**2))
983  DO iset = 1, nseta
984 
985  DO ishella = 1, nshella(iset)
986  sgfa_start = first_sgfa(ishella, iset)
987  sgfa_end = sgfa_start + 2*la(ishella, iset)
988  lai_set(sgfa_start:sgfa_end) = la(ishella, iset)
989  END DO
990  istart = first_sgfa(1, iset)
991  iend = istart + nsgfa_set(iset) - 1
992 
993  DO jset = 1, nseta
994 
995  IF (set_radius_a(iset) + set_radius_a(jset) < dab) cycle
996  DO jshella = 1, nshella(jset)
997  sgfa_start = first_sgfa(jshella, jset)
998  sgfa_end = sgfa_start + 2*la(jshella, jset)
999  laj_set(sgfa_start:sgfa_end) = la(jshella, jset)
1000  END DO
1001  jstart = first_sgfa(1, jset)
1002  jend = jstart + nsgfa_set(jset) - 1
1003 
1004  DO kset = 1, nsetca
1005 
1006  IF (set_radius_a(iset) + set_radius_ca(kset) < dab) cycle
1007 
1008  DO kshella = 1, nshellca(kset)
1009  sgfa_start = first_sgfca(kshella, kset)
1010  sgfa_end = sgfa_start + 2*lca(kshella, kset)
1011  lak_set(sgfa_start:sgfa_end) = lca(kshella, kset)
1012  END DO
1013  kstart = first_sgfca(1, kset)
1014  kend = kstart + nsgfca_set(kset) - 1
1015  DO ksgfa = kstart, kend
1016  lak = lak_set(ksgfa)
1017  DO jsgfa = jstart, jend
1018  laj = laj_set(jsgfa)
1019  DO isgfa = istart, iend
1020  lai = lai_set(isgfa)
1021  IF (modulo((lai + laj + lak), 2) /= 0) THEN
1022  IF (calculate_ints) THEN
1023  abbint(isgfa, jsgfa, ksgfa) = &
1024  -abaint(jsgfa, isgfa, ksgfa)
1025  END IF
1026  IF (calculate_forces) THEN
1027  DO i = 1, 3
1028  dabbint(isgfa, jsgfa, ksgfa, i) = &
1029  -dabdaint(jsgfa, isgfa, ksgfa, i)
1030  END DO
1031  END IF
1032  ELSE
1033  IF (calculate_ints) THEN
1034  abbint(isgfa, jsgfa, ksgfa) = &
1035  abaint(jsgfa, isgfa, ksgfa)
1036  END IF
1037  IF (calculate_forces) THEN
1038  DO i = 1, 3
1039  dabbint(isgfa, jsgfa, ksgfa, i) = &
1040  dabdaint(jsgfa, isgfa, ksgfa, i)
1041  END DO
1042  END IF
1043  END IF
1044  END DO
1045  END DO
1046  END DO
1047  END DO
1048  END DO
1049  END DO
1050 
1051  DEALLOCATE (lai_set, laj_set, lak_set)
1052 
1053  END SUBROUTINE get_abb_same_kind
1054 
1055 ! **************************************************************************************************
1056 !> \brief calculate integrals (a,b,fa); requires angular-dependent part as input
1057 !> \param abaint integral (a,b,fa)
1058 !> \param dabdaint ...
1059 !> \param rab distance vector between A and B
1060 !> \param oba orbital basis at center A
1061 !> \param obb orbital basis at center B
1062 !> \param fba auxiliary basis set at center A
1063 !> \param scon_obb contraction matrix for orb bas on B
1064 !> \param scona_mix mixed contraction matrix orb + ri basis on A
1065 !> \param oba_index orbital basis index for scona_mix
1066 !> \param fba_index ri basis index for scona_mix
1067 !> \param cg_coeff Clebsch-Gordon coefficients
1068 !> \param cg_none0_list list of none-zero Clebsch-Gordon coefficients
1069 !> \param ncg_none0 number of non-zero Clebsch-Gordon coefficients
1070 !> \param Waux_mat W matrix that contains the angular-dependent part
1071 !> \param dWaux_mat derivative of the W matrix
1072 !> \param calculate_ints ...
1073 !> \param calculate_forces ...
1074 ! **************************************************************************************************
1075  SUBROUTINE int_overlap_aba_shg_low(abaint, dabdaint, rab, oba, obb, fba, scon_obb, scona_mix, &
1076  oba_index, fba_index, cg_coeff, cg_none0_list, &
1077  ncg_none0, Waux_mat, dWaux_mat, calculate_ints, calculate_forces)
1078 
1079  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1080  INTENT(INOUT) :: abaint
1081  REAL(kind=dp), ALLOCATABLE, &
1082  DIMENSION(:, :, :, :), INTENT(INOUT) :: dabdaint
1083  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1084  TYPE(gto_basis_set_type), POINTER :: oba, obb, fba
1085  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_obb
1086  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: scona_mix
1087  INTEGER, DIMENSION(:, :, :), INTENT(IN) :: oba_index, fba_index
1088  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: cg_coeff
1089  INTEGER, DIMENSION(:, :, :), INTENT(IN) :: cg_none0_list
1090  INTEGER, DIMENSION(:, :), INTENT(IN) :: ncg_none0
1091  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: waux_mat
1092  REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN) :: dwaux_mat
1093  LOGICAL, INTENT(IN) :: calculate_ints, calculate_forces
1094 
1095  INTEGER :: iset, jset, kset, la_max_set, laa_max, laa_max_set, lb_max_set, lca_max_set, na, &
1096  nb, nca, ndev, nds, nds_max, nl, nl_set, npgfa_set, npgfb_set, npgfca_set, nseta, nsetb, &
1097  nsetca, nshella_set, nshellb_set, nshellca_set, sgfa, sgfb, sgfca
1098  INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, lca_max, npgfa, npgfb, &
1099  npgfca, nsgfa, nsgfb, nsgfca, nshella, &
1100  nshellb, nshellca
1101  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfca, la, &
1102  lb, lca
1103  REAL(kind=dp) :: dab, rab2
1104  REAL(kind=dp), ALLOCATABLE, &
1105  DIMENSION(:, :, :, :, :) :: swork, swork_cont
1106  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_ca
1107  REAL(kind=dp), DIMENSION(:, :), POINTER :: zeta, zetb, zetca
1108 
1109  NULLIFY (la_max, lb_max, lca_max, npgfa, npgfb, npgfca)
1110  NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, &
1111  set_radius_ca, zeta, zetb, zetca)
1112 
1113  ! basis ikind
1114  first_sgfa => oba%first_sgf
1115  la_max => oba%lmax
1116  la => oba%l
1117  nsgfa => oba%nsgf_set
1118  npgfa => oba%npgf
1119  nshella => oba%nshell
1120  nseta = oba%nset
1121  set_radius_a => oba%set_radius
1122  zeta => oba%zet
1123  ! basis jkind
1124  first_sgfb => obb%first_sgf
1125  lb_max => obb%lmax
1126  lb => obb%l
1127  nsgfb => obb%nsgf_set
1128  npgfb => obb%npgf
1129  nshellb => obb%nshell
1130  nsetb = obb%nset
1131  set_radius_b => obb%set_radius
1132  zetb => obb%zet
1133 
1134  ! basis RI A
1135  first_sgfca => fba%first_sgf
1136  lca_max => fba%lmax
1137  lca => fba%l
1138  nsgfca => fba%nsgf_set
1139  npgfca => fba%npgf
1140  nshellca => fba%nshell
1141  nsetca = fba%nset
1142  set_radius_ca => fba%set_radius
1143  zetca => fba%zet
1144 
1145  dab = sqrt(sum(rab**2))
1146  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1147 
1148  la_max_set = maxval(la_max)
1149  lb_max_set = maxval(lb_max)
1150  lca_max_set = maxval(lca_max)
1151  npgfa_set = maxval(npgfa)
1152  npgfb_set = maxval(npgfb)
1153  npgfca_set = maxval(npgfca)
1154  nshella_set = maxval(nshella)
1155  nshellb_set = maxval(nshellb)
1156  nshellca_set = maxval(nshellca)
1157  !*** for forces: derivative+1 in auxiliary vector required
1158  ndev = 0
1159  IF (calculate_forces) ndev = 1
1160 
1161  laa_max_set = la_max_set + lca_max_set
1162 
1163  ! allocate some work storage....
1164  nds_max = laa_max_set + lb_max_set + ndev + 1
1165  nl_set = int((laa_max_set)/2)
1166  ALLOCATE (swork(npgfb_set, npgfa_set, npgfca_set, nl_set + 1, nds_max))
1167  ALLOCATE (swork_cont(nds_max, 0:nl_set, nshella_set, nshellb_set, nshellca_set))
1168 
1169  DO iset = 1, nseta
1170 
1171  DO jset = 1, nsetb
1172 
1173  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1174 
1175  DO kset = 1, nsetca
1176 
1177  IF (set_radius_b(jset) + set_radius_ca(kset) < dab) cycle
1178 
1179  !*** calculate s_baa here
1180  laa_max = la_max(iset) + lca_max(kset)
1181  nds = laa_max + lb_max(jset) + ndev + 1
1182  nl = int(laa_max/2) + 1
1183  swork(1:npgfb(jset), 1:npgfa(iset), 1:npgfca(kset), 1:nl, 1:nds) = 0.0_dp
1184  CALL s_overlap_abb(lb_max(jset), npgfb(jset), zetb(:, jset), &
1185  la_max(iset), npgfa(iset), zeta(:, iset), &
1186  lca_max(kset), npgfca(kset), zetca(:, kset), &
1187  rab, swork, calculate_forces)
1188 
1189  CALL contract_s_overlap_aba(la(:, iset), npgfa(iset), nshella(iset), &
1190  lb(:, jset), npgfb(jset), nshellb(jset), &
1191  scon_obb(1:npgfb(jset), 1:nshellb(jset), jset), &
1192  lca(:, kset), npgfca(kset), nshellca(kset), &
1193  oba_index(:, :, iset), fba_index(:, :, kset), scona_mix, nl, nds, &
1194  swork(1:npgfb(jset), 1:npgfa(iset), 1:npgfca(kset), 1:nl, 1:nds), &
1195  swork_cont, calculate_forces)
1196  IF (calculate_ints) THEN
1197  CALL construct_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
1198  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
1199  lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
1200  cg_coeff, cg_none0_list, ncg_none0, &
1201  swork_cont, waux_mat, abaint)
1202  END IF
1203  IF (calculate_forces) THEN
1204  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
1205  CALL dev_overlap_shg_aba(la(:, iset), first_sgfa(:, iset), nshella(iset), &
1206  lb(:, jset), first_sgfb(:, jset), nshellb(jset), &
1207  lca(:, kset), first_sgfca(:, kset), nshellca(kset), &
1208  cg_coeff, cg_none0_list, ncg_none0, &
1209  -rab, swork_cont, waux_mat, dwaux_mat, dabdaint)
1210  END IF
1211  ! max value of integrals in this set triple
1212  sgfa = first_sgfa(1, iset)
1213  na = sgfa + nsgfa(iset) - 1
1214  sgfb = first_sgfb(1, jset)
1215  nb = sgfb + nsgfb(jset) - 1
1216  sgfca = first_sgfca(1, kset)
1217  nca = sgfca + nsgfca(kset) - 1
1218 
1219  END DO
1220  END DO
1221  END DO
1222 
1223  DEALLOCATE (swork_cont)
1224  DEALLOCATE (swork)
1225 
1226  END SUBROUTINE int_overlap_aba_shg_low
1227 
1228 ! **************************************************************************************************
1229 !> \brief precalculates the angular part of the SHG integrals for the matrices
1230 !> (fa,fb), (a,b), (a,b,fa) and (b,fb,a); the same Waux_mat can then be used for all
1231 !> for integrals; specific for LRIGPW
1232 !> \param oba orbital basis on a
1233 !> \param obb orbital basis on b
1234 !> \param fba aux basis on a
1235 !> \param fbb aux basis on b
1236 !> \param rab distance vector between a and b
1237 !> \param Waux_mat W matrix that contains the angular-dependent part
1238 !> \param dWaux_mat derivative of the W matrix
1239 !> \param calculate_forces ...
1240 ! **************************************************************************************************
1241  SUBROUTINE lri_precalc_angular_shg_part(oba, obb, fba, fbb, rab, Waux_mat, dWaux_mat, calculate_forces)
1242 
1243  TYPE(gto_basis_set_type), POINTER :: oba, obb, fba, fbb
1244  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
1245  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1246  INTENT(OUT) :: waux_mat
1247  REAL(kind=dp), ALLOCATABLE, &
1248  DIMENSION(:, :, :, :), INTENT(OUT) :: dwaux_mat
1249  LOGICAL, INTENT(IN) :: calculate_forces
1250 
1251  INTEGER :: i, isize, j, k, la_max, laa_max, lb_max, &
1252  lbb_max, lca_max, lcb_max, li_max, &
1253  lj_max, lmax, mdim(3), size_int(4, 2), &
1254  temp
1255  INTEGER, DIMENSION(:), POINTER :: li_max_all
1256  REAL(kind=dp) :: rab2
1257  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rc, rs
1258 
1259  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1260 
1261  !*** 1 Waux_mat of size (li_max,lj_max) for elements
1262  ! i j
1263  ! [aab] --> (laa_max, lb_max)
1264  ! [bba] --> (lbb_max, la_max) --> use for [abb]
1265  ! [ab] ri --> (lca_max, lcb_max)
1266  ! [ab] orb --> (la_max , lb_max)
1267 
1268  la_max = maxval(oba%lmax)
1269  lb_max = maxval(obb%lmax)
1270  lca_max = maxval(fba%lmax)
1271  lcb_max = maxval(fbb%lmax)
1272 
1273  laa_max = la_max + lca_max
1274  lbb_max = lb_max + lcb_max
1275  li_max = max(laa_max, lbb_max)
1276  lj_max = max(la_max, lb_max, lcb_max)
1277  lmax = li_max
1278 
1279  ALLOCATE (li_max_all(0:lj_max))
1280  ALLOCATE (rc(0:lmax, -2*lmax:2*lmax), rs(0:lmax, -2*lmax:2*lmax))
1281  rc = 0._dp
1282  rs = 0._dp
1283  mdim(1) = li_max + lj_max + 1
1284  mdim(2) = nsoset(li_max) + 1
1285  mdim(3) = nsoset(lj_max) + 1
1286  ALLOCATE (waux_mat(mdim(1), mdim(2), mdim(3)))
1287  ALLOCATE (dwaux_mat(3, mdim(1), mdim(2), mdim(3)))
1288  !Waux_mat = 0._dp !.. takes time
1289  !dWaux_mat =0._dp !.. takes time
1290 
1291  !*** Waux_mat (li_max,lj_max) contains elements not needed,
1292  !*** make indixing so that only required ones are computed
1293  !*** li_max_all(j) --> li_max dependent on j
1294  size_int(1, :) = (/laa_max, lb_max/)
1295  size_int(2, :) = (/lbb_max, la_max/)
1296  size_int(3, :) = (/lca_max, lcb_max/)
1297  size_int(4, :) = (/la_max, lb_max/)
1298 
1299  li_max_all(:) = 0
1300  DO isize = 1, 4
1301  i = size_int(isize, 1)
1302  j = size_int(isize, 2)
1303  k = li_max_all(j)
1304  IF (k < i) li_max_all(j) = i
1305  END DO
1306  temp = li_max_all(lj_max)
1307  DO j = lj_max, 0, -1
1308  IF (li_max_all(j) < temp) THEN
1309  li_max_all(j) = temp
1310  ELSE
1311  temp = li_max_all(j)
1312  END IF
1313  END DO
1314 
1315  !*** -rab, since Eq. in Ref. use Ra-Rb, not Rb-Ra
1316  CALL get_real_scaled_solid_harmonic(rc, rs, lmax, -rab, rab2)
1317  CALL get_w_matrix(li_max_all, lj_max, lmax, rc, rs, waux_mat)
1318  IF (calculate_forces) THEN
1319  CALL get_dw_matrix(li_max_all, lj_max, waux_mat, dwaux_mat)
1320  END IF
1321 
1322  DEALLOCATE (rc, rs, li_max_all)
1323 
1324  END SUBROUTINE lri_precalc_angular_shg_part
1325 
1326 END MODULE generic_shg_integrals
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
constants for the different operators of the 2c-integrals
integer, parameter, public operator_gauss
integer, parameter, public operator_verf
integer, parameter, public operator_vgauss
integer, parameter, public operator_verfc
integer, parameter, public operator_coulomb
Calculation of the integrals over solid harmonic Gaussian(SHG) functions. Routines for (a|O(r12)|b) a...
Definition: construct_shg.F:19
subroutine, public dev_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, ncg_none0, rab, swork_cont, Waux_mat, dWaux_mat, dsabb)
calculates derivatives of [abb] SHG overlap integrals using precomputed angular-dependent part
subroutine, public get_w_matrix(lamax, lbmax, lmax, Rc, Rs, Waux_mat)
calculates the angular dependent-part of the SHG integrals, transformation matrix W,...
subroutine, public construct_int_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, swork_cont, Waux_mat, sab)
calculates [ab] SHG overlap integrals using precomputed angular- dependent part
subroutine, public construct_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, ncg_none0, swork_cont, Waux_mat, saba)
calculates [aba] SHG overlap integrals using precomputed angular- dependent part
subroutine, public construct_dev_shg_ab(la, first_sgfa, nshella, lb, first_sgfb, nshellb, rab, swork_cont, Waux_mat, dWaux_mat, dsab)
calculates derivatives of [ab] SHG overlap integrals using precomputed angular-dependent part
subroutine, public dev_overlap_shg_aba(la, first_sgfa, nshella, lb, first_sgfb, nshellb, lca, first_sgfca, nshellca, cg_coeff, cg_none0_list, ncg_none0, rab, swork_cont, Waux_mat, dWaux_mat, dsaba)
calculates derivatives of [aba] SHG overlap integrals using precomputed angular-dependent part
subroutine, public get_dw_matrix(lamax, lbmax, Waux_mat, dWaux_mat)
calculates derivatives of transformation matrix W,
subroutine, public get_real_scaled_solid_harmonic(Rlm_c, Rlm_s, l, r, r2)
computes the real scaled solid harmonics Rlm up to a given l
Definition: construct_shg.F:50
subroutine, public construct_overlap_shg_abb(la, first_sgfa, nshella, lb, first_sgfb, nshellb, lcb, first_sgfcb, nshellcb, cg_coeff, cg_none0_list, ncg_none0, swork_cont, Waux_mat, sabb)
calculates [abb] SHG overlap integrals using precomputed angular- dependent part
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_aba_shg(saba, dsaba, rab, oba, obb, fba, scon_obb, scona_mix, oba_index, fba_index, cg_coeff, cg_none0_list, ncg_none0, calculate_forces)
calculate integrals (a,b,fa)
subroutine, public int_ra2m_ab_shg_low(vab, dvab, rab, fba, fbb, sconb_shg, scon_ra2m, m, Waux_mat, dWaux_mat, calculate_forces)
calculate integrals (a|ra^2m)|b); requires angular-dependent part as input
subroutine, public int_overlap_ab_shg(vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, calculate_forces)
calculate overlap integrals (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_ra2m_ab_shg(vab, dvab, rab, fba, fbb, scon_ra2m, sconb_shg, m, calculate_forces)
Calcululates the two-center integrals of the type (a|(r-Ra)^(2m)|b) using the SHG scheme.
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
subroutine, public int_overlap_abb_shg(sabb, dsabb, rab, oba, obb, fbb, scon_oba, sconb_mix, obb_index, fbb_index, cg_coeff, cg_none0_list, ncg_none0, calculate_forces)
calculate integrals (a,b,fb)
subroutine, public int_operators_r12_ab_shg(r12_operator, vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, omega, calculate_forces)
Calcululates the two-center integrals of the type (a|O(r12)|b) using the SHG scheme.
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public nsoset
Routines for calculating the s-integrals and their scalar derivatives with respect to rab2 over solid...
subroutine, public contract_s_overlap_aba(la, npgfa, nshella, lb, npgfb, nshellb, sconb_shg, lca, npgfca, nshellca, orba_index, ria_index, scona_mix, nl_max, nds_max, swork, swork_cont, calculate_forces)
Contraction and normalization of the (aba) overlap.
subroutine, public s_verf_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/erf(omega*r12)/r12|s] two-center integral
subroutine, public s_gauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|exp(-omega*r12**2)|s] two-center integral
subroutine, public s_overlap_abb(la_max, npgfa, zeta, lb_max, npgfb, zetb, lcb_max, npgfcb, zetcb, rab, s, calculate_forces)
calculates [s|ra^n|s] integrals for [aba] and the [s|rb^n|s] integrals for [abb]
subroutine, public s_vgauss_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|exp(-omega*r12**2)/r12|s] two-center integral
subroutine, public contract_s_ra2m_ab(npgfa, nshella, scon_ra2m, npgfb, nshellb, sconb, swork, swork_cont, m, nds_max)
Contraction, normalization and combinatiorial combination of the [0a|(r-Ra)^(2m)|0b] integrals and th...
subroutine, public contract_sint_ab_chigh(npgfa, nshella, scona, npgfb, nshellb, sconb, nds, swork, swork_cont)
Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives; this routin...
subroutine, public s_coulomb_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/r12|s] two-center coulomb integral
subroutine, public contract_s_overlap_abb(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, lcb, npgfcb, nshellcb, orbb_index, rib_index, sconb_mix, nl_max, nds_max, swork, swork_cont, calculate_forces)
Contraction and normalization of the (abb) overlap.
subroutine, public s_verfc_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, omega, rab, v, calculate_forces)
calculates the uncontracted, not normalized [s|1/erfc(omega*r12)/r12|s] two-center integral
subroutine, public contract_sint_ab_clow(la, npgfa, nshella, scona_shg, lb, npgfb, nshellb, sconb_shg, swork, swork_cont, calculate_forces)
Contraction and normalization of the [s|O(r12)|s] integrals and their scalar derivatives; this routin...
subroutine, public s_overlap_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, rab, s, calculate_forces)
calculates the uncontracted, not normalized [s|s] overlap
subroutine, public s_ra2m_ab(la_max, npgfa, zeta, lb_max, npgfb, zetb, m, rab, s, calculate_forces)
calculates the uncontracted, not normalized [0a|ra^(2m)|0b] two-center integral, where ra = r-Ra and ...