Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of the kinetic energy integrals over Cartesian
10 : !> Gaussian-type functions.
11 : !>
12 : !> [a|T|b] = [a|-nabla**2/2|b]
13 : !> \par Literature
14 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
15 : !> \par History
16 : !> - Derivatives added (10.05.2002,MK)
17 : !> - Fully refactored (07.07.2014,JGH)
18 : !> \author Matthias Krack (31.07.2000)
19 : ! **************************************************************************************************
20 : MODULE ai_kinetic
21 : USE ai_cartesian_operators, ONLY: apply_cartesian_1d
22 : USE ai_os_rr, ONLY: os_rr_ovlp
23 : USE kinds, ONLY: dp
24 : USE mathconstants, ONLY: pi
25 : USE orbital_pointers, ONLY: indco,&
26 : ncoset
27 : #include "../base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 :
31 : PRIVATE
32 :
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_kinetic'
34 :
35 : ! *** Public subroutines ***
36 :
37 : PUBLIC :: kinetic
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief Calculation of the two-center kinetic energy integrals [a|T|b] over
43 : !> Cartesian Gaussian-type functions.
44 : !> \param la_max Maximum L of basis on A
45 : !> \param la_min Minimum L of basis on A
46 : !> \param npgfa Number of primitive functions in set of basis on A
47 : !> \param rpgfa Range of functions on A (used for prescreening)
48 : !> \param zeta Exponents of basis on center A
49 : !> \param lb_max Maximum L of basis on A
50 : !> \param lb_min Minimum L of basis on A
51 : !> \param npgfb Number of primitive functions in set of basis on B
52 : !> \param rpgfb Range of functions on B (used for prescreening)
53 : !> \param zetb Exponents of basis on center B
54 : !> \param rab Distance vector between centers A and B
55 : !> \param kab Kinetic energy integrals, optional
56 : !> \param dab First derivatives of Kinetic energy integrals, optional
57 : !> \date 07.07.2014
58 : !> \author JGH
59 : ! **************************************************************************************************
60 4763531 : SUBROUTINE kinetic(la_max, la_min, npgfa, rpgfa, zeta, &
61 4763531 : lb_max, lb_min, npgfb, rpgfb, zetb, &
62 4763531 : rab, kab, dab)
63 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
64 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
65 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
66 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
67 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
68 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
69 : OPTIONAL :: kab
70 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
71 : OPTIONAL :: dab
72 :
73 : INTEGER :: avec(3), bvec(3), coa, cob, ia, ib, &
74 : idir, ipgf, jdir, jpgf, kdir, ldrr, &
75 : lma, lmb, ma, mb, na, nb, ofa, ofb
76 : REAL(KIND=dp) :: a, b, f0, rab2, tab, xhi, zet
77 4763531 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drr, dtt, rr, tt
78 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
79 :
80 4763531 : CPASSERT(PRESENT(kab) .OR. PRESENT(dab))
81 :
82 : ! Distance of the centers a and b
83 :
84 4763531 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
85 4763531 : tab = SQRT(rab2)
86 :
87 : ! Maximum l for auxiliary integrals
88 4763531 : IF (PRESENT(kab)) THEN
89 4763531 : lma = la_max + 1
90 4763531 : lmb = lb_max + 1
91 : END IF
92 4763531 : IF (PRESENT(dab)) THEN
93 922122 : lma = la_max + 2
94 922122 : lmb = lb_max + 1
95 : END IF
96 4763531 : ldrr = MAX(lma, lmb) + 1
97 :
98 : ! Allocate space for auxiliary integrals
99 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), tt(0:ldrr - 1, 0:ldrr - 1, 3), &
100 42871779 : drr(0:ldrr - 1, 0:ldrr - 1, 3))
101 6607775 : IF (PRESENT(dab)) ALLOCATE (dtt(0:ldrr - 1, 0:ldrr - 1, 3))
102 :
103 : ! Number of integrals, check size of arrays
104 4763531 : ofa = ncoset(la_min - 1)
105 4763531 : ofb = ncoset(lb_min - 1)
106 4763531 : na = ncoset(la_max) - ofa
107 4763531 : nb = ncoset(lb_max) - ofb
108 4763531 : IF (PRESENT(kab)) THEN
109 4763531 : CPASSERT((SIZE(kab, 1) >= na*npgfa))
110 4763531 : CPASSERT((SIZE(kab, 2) >= nb*npgfb))
111 : END IF
112 4763531 : IF (PRESENT(dab)) THEN
113 922122 : CPASSERT((SIZE(dab, 1) >= na*npgfa))
114 922122 : CPASSERT((SIZE(dab, 2) >= nb*npgfb))
115 922122 : CPASSERT((SIZE(dab, 3) >= 3))
116 : END IF
117 :
118 : ! Loops over all pairs of primitive Gaussian-type functions
119 4763531 : ma = 0
120 18842644 : DO ipgf = 1, npgfa
121 14079113 : mb = 0
122 64226016 : DO jpgf = 1, npgfb
123 : ! Distance Screening
124 50146903 : IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
125 947984918 : IF (PRESENT(kab)) kab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
126 673299603 : IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
127 36892038 : mb = mb + nb
128 36892038 : CYCLE
129 : END IF
130 :
131 : ! Calculate some prefactors
132 13254865 : a = zeta(ipgf)
133 13254865 : b = zetb(jpgf)
134 13254865 : zet = a + b
135 13254865 : xhi = a*b/zet
136 53019460 : rap = b*rab/zet
137 53019460 : rbp = -a*rab/zet
138 :
139 : ! [s|s] integral
140 13254865 : f0 = 0.5_dp*(pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
141 :
142 : ! Calculate the recurrence relation, overlap
143 13254865 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
144 :
145 : ! Mixed center derivatives give the overlap of the two electronic gradients.
146 13254865 : CALL apply_cartesian_1d(rr, drr, lma - 1, lmb, .FALSE., .TRUE., a, [0.0_dp, 0.0_dp, 0.0_dp])
147 13254865 : CALL apply_cartesian_1d(drr, tt, lma - 1, lmb - 1, .TRUE., .TRUE., b, [0.0_dp, 0.0_dp, 0.0_dp])
148 13254865 : IF (PRESENT(dab)) THEN
149 2747549 : CALL apply_cartesian_1d(tt, dtt, la_max, lb_max, .FALSE., .TRUE., a, [0.0_dp, 0.0_dp, 0.0_dp])
150 : END IF
151 :
152 61348497 : DO cob = ofb + 1, ncoset(lb_max)
153 192374528 : bvec = indco(:, cob)
154 48093632 : ib = mb + cob - ofb
155 292287617 : DO coa = ofa + 1, ncoset(la_max)
156 923756480 : avec = indco(:, coa)
157 230939120 : ia = ma + coa - ofa
158 230939120 : IF (PRESENT(kab)) kab(ia, ib) = 0.0_dp
159 923756480 : DO idir = 1, 3
160 692817360 : jdir = MOD(idir, 3) + 1
161 692817360 : kdir = MOD(idir + 1, 3) + 1
162 692817360 : IF (PRESENT(kab)) THEN
163 : kab(ia, ib) = kab(ia, ib) + tt(avec(idir), bvec(idir), idir)* &
164 692817360 : rr(avec(jdir), bvec(jdir), jdir)*rr(avec(kdir), bvec(kdir), kdir)
165 : END IF
166 923756480 : IF (PRESENT(dab)) THEN
167 : dab(ia, ib, idir) = dtt(avec(idir), bvec(idir), idir)* &
168 : rr(avec(jdir), bvec(jdir), jdir)*rr(avec(kdir), bvec(kdir), kdir) + &
169 : drr(avec(idir), bvec(idir), idir)* &
170 : (tt(avec(jdir), bvec(jdir), jdir)*rr(avec(kdir), bvec(kdir), kdir) + &
171 151655508 : rr(avec(jdir), bvec(jdir), jdir)*tt(avec(kdir), bvec(kdir), kdir))
172 : END IF
173 : END DO
174 230939120 : IF (PRESENT(kab)) kab(ia, ib) = f0*kab(ia, ib)
175 430688260 : IF (PRESENT(dab)) dab(ia, ib, 1:3) = f0*dab(ia, ib, 1:3)
176 : END DO
177 : END DO
178 :
179 27333978 : mb = mb + nb
180 : END DO
181 18842644 : ma = ma + na
182 : END DO
183 :
184 4763531 : DEALLOCATE (rr, tt)
185 :
186 4763531 : END SUBROUTINE kinetic
187 :
188 : END MODULE ai_kinetic
|
