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 : apply_cartesian_operator_pair,&
23 : cartesian_ladder_map
24 : USE ai_os_rr, ONLY: os_rr_ovlp
25 : USE kinds, ONLY: dp
26 : USE mathconstants, ONLY: pi
27 : USE orbital_pointers, ONLY: indco,&
28 : ncoset
29 : #include "../base/base_uses.f90"
30 :
31 : IMPLICIT NONE
32 :
33 : PRIVATE
34 :
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_kinetic'
36 :
37 : ! *** Public subroutines ***
38 :
39 : PUBLIC :: kinetic
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief Calculation of the two-center kinetic energy integrals [a|T|b] over
45 : !> Cartesian Gaussian-type functions.
46 : !> \param la_max Maximum L of basis on A
47 : !> \param la_min Minimum L of basis on A
48 : !> \param npgfa Number of primitive functions in set of basis on A
49 : !> \param rpgfa Range of functions on A (used for prescreening)
50 : !> \param zeta Exponents of basis on center A
51 : !> \param lb_max Maximum L of basis on A
52 : !> \param lb_min Minimum L of basis on A
53 : !> \param npgfb Number of primitive functions in set of basis on B
54 : !> \param rpgfb Range of functions on B (used for prescreening)
55 : !> \param zetb Exponents of basis on center B
56 : !> \param rab Distance vector between centers A and B
57 : !> \param kab Kinetic energy integrals, optional
58 : !> \param dab First derivatives of Kinetic energy integrals, optional
59 : !> \date 07.07.2014
60 : !> \author JGH
61 : ! **************************************************************************************************
62 4775975 : SUBROUTINE kinetic(la_max, la_min, npgfa, rpgfa, zeta, &
63 4775975 : lb_max, lb_min, npgfb, rpgfb, zetb, &
64 4775975 : rab, kab, dab)
65 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
66 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
67 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
68 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
69 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
70 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
71 : OPTIONAL :: kab
72 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
73 : OPTIONAL :: dab
74 :
75 : INTEGER, PARAMETER :: component(2, 1) = RESHAPE([1, 1], [2, 1])
76 :
77 : INTEGER :: avec(3), bvec(3), coa, cob, ia, ib, &
78 : idir, ipgf, jdir, jpgf, kdir, ldrr, &
79 : lma, lmb, ma, mb, na, nb, ofa, ofb
80 4775975 : INTEGER, ALLOCATABLE :: map_a(:, :, :), map_b(:, :, :)
81 : REAL(KIND=dp) :: a, b, f0, rab2, tab, xhi, zet
82 4775975 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: drr, dtt, rr, tt
83 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp
84 :
85 4775975 : CPASSERT(PRESENT(kab) .OR. PRESENT(dab))
86 :
87 : ! Distance of the centers a and b
88 :
89 4775975 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
90 4775975 : tab = SQRT(rab2)
91 :
92 : ! Maximum l for auxiliary integrals
93 4775975 : IF (PRESENT(kab)) THEN
94 4775975 : lma = la_max + 1
95 4775975 : lmb = lb_max + 1
96 : END IF
97 4775975 : IF (PRESENT(dab)) THEN
98 922371 : lma = la_max + 2
99 922371 : lmb = lb_max + 1
100 : END IF
101 4775975 : ldrr = MAX(lma, lmb) + 1
102 :
103 : ! Allocate space for auxiliary integrals
104 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), tt(0:ldrr - 1, 0:ldrr - 1, 3), &
105 42983775 : drr(0:ldrr - 1, 0:ldrr - 1, 3))
106 6620717 : IF (PRESENT(dab)) ALLOCATE (dtt(0:ldrr - 1, 0:ldrr - 1, 3))
107 23879875 : ALLOCATE (map_a(4, lma, 1), map_b(4, lmb, 1))
108 4775975 : CALL cartesian_ladder_map(lma - 1, 0, 0, map_a(:, :, 1))
109 4775975 : CALL cartesian_ladder_map(lmb - 1, 0, 0, map_b(:, :, 1))
110 :
111 : ! Number of integrals, check size of arrays
112 4775975 : ofa = ncoset(la_min - 1)
113 4775975 : ofb = ncoset(lb_min - 1)
114 4775975 : na = ncoset(la_max) - ofa
115 4775975 : nb = ncoset(lb_max) - ofb
116 4775975 : IF (PRESENT(kab)) THEN
117 4775975 : CPASSERT((SIZE(kab, 1) >= na*npgfa))
118 4775975 : CPASSERT((SIZE(kab, 2) >= nb*npgfb))
119 : END IF
120 4775975 : IF (PRESENT(dab)) THEN
121 922371 : CPASSERT((SIZE(dab, 1) >= na*npgfa))
122 922371 : CPASSERT((SIZE(dab, 2) >= nb*npgfb))
123 922371 : CPASSERT((SIZE(dab, 3) >= 3))
124 : END IF
125 :
126 : ! Loops over all pairs of primitive Gaussian-type functions
127 4775975 : ma = 0
128 18903038 : DO ipgf = 1, npgfa
129 14127063 : mb = 0
130 64467326 : DO jpgf = 1, npgfb
131 : ! Distance Screening
132 50340263 : IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
133 951846492 : IF (PRESENT(kab)) kab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
134 673523634 : IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
135 37034670 : mb = mb + nb
136 37034670 : CYCLE
137 : END IF
138 :
139 : ! Calculate some prefactors
140 13305593 : a = zeta(ipgf)
141 13305593 : b = zetb(jpgf)
142 13305593 : zet = a + b
143 13305593 : xhi = a*b/zet
144 53222372 : rap = b*rab/zet
145 53222372 : rbp = -a*rab/zet
146 :
147 : ! [s|s] integral
148 13305593 : f0 = 0.5_dp*(pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
149 :
150 : ! Calculate the recurrence relation, overlap
151 13305593 : CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
152 :
153 : ! Mixed center derivatives give the overlap of the two electronic gradients.
154 53222372 : DO idir = 1, 3
155 : CALL apply_cartesian_operator_pair(rr(:, :, idir), tt(0:lma - 1, 0:lmb - 1, idir:idir), &
156 : map_a, map_b, component, [.TRUE., .TRUE.], [a, b], &
157 133055930 : [.FALSE., .TRUE.], drr(0:lma - 1, 0:lmb, idir))
158 : END DO
159 13305593 : IF (PRESENT(dab)) THEN
160 2748264 : CALL apply_cartesian_1d(tt, dtt, la_max, lb_max, .FALSE., .TRUE., a, [0.0_dp, 0.0_dp, 0.0_dp])
161 : END IF
162 :
163 61620798 : DO cob = ofb + 1, ncoset(lb_max)
164 193260820 : bvec = indco(:, cob)
165 48315205 : ib = mb + cob - ofb
166 293689017 : DO coa = ofa + 1, ncoset(la_max)
167 928272876 : avec = indco(:, coa)
168 232068219 : ia = ma + coa - ofa
169 232068219 : IF (PRESENT(kab)) kab(ia, ib) = 0.0_dp
170 928272876 : DO idir = 1, 3
171 696204657 : jdir = MOD(idir, 3) + 1
172 696204657 : kdir = MOD(idir + 1, 3) + 1
173 696204657 : IF (PRESENT(kab)) THEN
174 : kab(ia, ib) = kab(ia, ib) + tt(avec(idir), bvec(idir), idir)* &
175 696204657 : rr(avec(jdir), bvec(jdir), jdir)*rr(avec(kdir), bvec(kdir), kdir)
176 : END IF
177 928272876 : IF (PRESENT(dab)) THEN
178 : dab(ia, ib, idir) = dtt(avec(idir), bvec(idir), idir)* &
179 : rr(avec(jdir), bvec(jdir), jdir)*rr(avec(kdir), bvec(kdir), kdir) + &
180 : drr(avec(idir), bvec(idir), idir)* &
181 : (tt(avec(jdir), bvec(jdir), jdir)*rr(avec(kdir), bvec(kdir), kdir) + &
182 151670610 : rr(avec(jdir), bvec(jdir), jdir)*tt(avec(kdir), bvec(kdir), kdir))
183 : END IF
184 : END DO
185 232068219 : IF (PRESENT(kab)) kab(ia, ib) = f0*kab(ia, ib)
186 432054034 : IF (PRESENT(dab)) dab(ia, ib, 1:3) = f0*dab(ia, ib, 1:3)
187 : END DO
188 : END DO
189 :
190 27432656 : mb = mb + nb
191 : END DO
192 18903038 : ma = ma + na
193 : END DO
194 :
195 4775975 : DEALLOCATE (rr, tt)
196 :
197 4775975 : END SUBROUTINE kinetic
198 :
199 : END MODULE ai_kinetic
|
