(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
21MODULE generic_shg_integrals
22 USE basis_set_types, ONLY: gto_basis_set_type
23 USE constants_operator, ONLY: operator_coulomb,&
24 operator_gauss,&
25 operator_verf,&
26 operator_verfc,&
27 operator_vgauss
28 USE construct_shg, ONLY: &
29 construct_dev_shg_ab, construct_int_shg_ab, construct_overlap_shg_aba, &
30 construct_overlap_shg_abb, dev_overlap_shg_aba, dev_overlap_shg_abb, get_w_matrix, &
31 get_dw_matrix, get_real_scaled_solid_harmonic
32 USE kinds, ONLY: dp
33 USE orbital_pointers, ONLY: nsoset
34 USE s_contract_shg, ONLY: &
35 contract_s_overlap_aba, contract_s_overlap_abb, contract_s_ra2m_ab, &
36 contract_sint_ab_chigh, contract_sint_ab_clow, s_coulomb_ab, s_gauss_ab, s_overlap_ab, &
37 s_overlap_abb, s_ra2m_ab, s_verf_ab, s_verfc_ab, s_vgauss_ab
38#include "../base/base_uses.f90"
39
40 IMPLICIT NONE
41
42 PRIVATE
43
44 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals'
45
46 PUBLIC :: int_operators_r12_ab_shg, int_overlap_ab_shg, &
47 int_ra2m_ab_shg, int_overlap_aba_shg, int_overlap_abb_shg, &
48 get_abb_same_kind, lri_precalc_angular_shg_part, &
49 int_overlap_ab_shg_low, int_ra2m_ab_shg_low, int_overlap_aba_shg_low, &
50 int_overlap_abb_shg_low
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
83CONTAINS
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
1326END MODULE generic_shg_integrals
modulo
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
basis_set_types
Definition basis_set_types.F:15
constants_operator
constants for the different operators of the 2c-integrals
Definition constants_operator.F:14
constants_operator::operator_gauss
integer, parameter, public operator_gauss
Definition constants_operator.F:22
constants_operator::operator_verf
integer, parameter, public operator_verf
Definition constants_operator.F:22
constants_operator::operator_vgauss
integer, parameter, public operator_vgauss
Definition constants_operator.F:22
constants_operator::operator_verfc
integer, parameter, public operator_verfc
Definition constants_operator.F:22
constants_operator::operator_coulomb
integer, parameter, public operator_coulomb
Definition constants_operator.F:22
construct_shg
Calculation of the integrals over solid harmonic Gaussian(SHG) functions. Routines for (a|O(r12)|b) a...
Definition construct_shg.F:19
construct_shg::dev_overlap_shg_aba
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
Definition construct_shg.F:692
construct_shg::get_dw_matrix
subroutine, public get_dw_matrix(lamax, lbmax, waux_mat, dwaux_mat)
calculates derivatives of transformation matrix W,
Definition construct_shg.F:314
construct_shg::get_w_matrix
subroutine, public get_w_matrix(lamax, lbmax, lmax, rc, rs, waux_mat)
calculates the angular dependent-part of the SHG integrals, transformation matrix W,...
Definition construct_shg.F:211
construct_shg::construct_dev_shg_ab
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
Definition construct_shg.F:542
construct_shg::construct_overlap_shg_abb
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
Definition construct_shg.F:786
construct_shg::construct_int_shg_ab
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
Definition construct_shg.F:485
construct_shg::dev_overlap_shg_abb
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
Definition construct_shg.F:874
construct_shg::construct_overlap_shg_aba
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
Definition construct_shg.F:607
construct_shg::get_real_scaled_solid_harmonic
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
generic_shg_integrals
Calculation of contracted, spherical Gaussian integrals using the solid harmonic Gaussian (SHG) integ...
Definition generic_shg_integrals.F:21
generic_shg_integrals::int_overlap_aba_shg
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)
Definition generic_shg_integrals.F:251
generic_shg_integrals::lri_precalc_angular_shg_part
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),...
Definition generic_shg_integrals.F:1242
generic_shg_integrals::int_overlap_aba_shg_low
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
Definition generic_shg_integrals.F:1078
generic_shg_integrals::int_overlap_abb_shg_low
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
Definition generic_shg_integrals.F:772
generic_shg_integrals::int_overlap_ab_shg
subroutine, public int_overlap_ab_shg(vab, dvab, rab, fba, fbb, scona_shg, sconb_shg, calculate_forces)
calculate overlap integrals (a,b)
Definition generic_shg_integrals.F:165
generic_shg_integrals::get_abb_same_kind
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...
Definition generic_shg_integrals.F:935
generic_shg_integrals::int_ra2m_ab_shg
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.
Definition generic_shg_integrals.F:205
generic_shg_integrals::int_overlap_ab_shg_low
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
Definition generic_shg_integrals.F:526
generic_shg_integrals::int_overlap_abb_shg
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)
Definition generic_shg_integrals.F:307
generic_shg_integrals::int_ra2m_ab_shg_low
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
Definition generic_shg_integrals.F:654
generic_shg_integrals::int_operators_r12_ab_shg
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.
Definition generic_shg_integrals.F:100
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition orbital_pointers.F:42
s_contract_shg
Routines for calculating the s-integrals and their scalar derivatives with respect to rab2 over solid...
Definition s_contract_shg.F:22
s_contract_shg::contract_s_overlap_aba
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.
Definition s_contract_shg.F:1142
s_contract_shg::s_verf_ab
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
Definition s_contract_shg.F:649
s_contract_shg::s_gauss_ab
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
Definition s_contract_shg.F:843
s_contract_shg::s_overlap_abb
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]
Definition s_contract_shg.F:116
s_contract_shg::s_vgauss_ab
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
Definition s_contract_shg.F:774
s_contract_shg::contract_s_ra2m_ab
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...
Definition s_contract_shg.F:1001
s_contract_shg::contract_sint_ab_chigh
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...
Definition s_contract_shg.F:962
s_contract_shg::s_coulomb_ab
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
Definition s_contract_shg.F:589
s_contract_shg::contract_s_overlap_abb
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.
Definition s_contract_shg.F:1053
s_contract_shg::s_verfc_ab
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
Definition s_contract_shg.F:711
s_contract_shg::contract_sint_ab_clow
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...
Definition s_contract_shg.F:907
s_contract_shg::s_overlap_ab
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
Definition s_contract_shg.F:58
s_contract_shg::s_ra2m_ab
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 ...
Definition s_contract_shg.F:347
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71