(git:ccc2433)
ai_kinetic.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 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_os_rr, ONLY: os_rr_ovlp
22  USE kinds, ONLY: dp
23  USE mathconstants, ONLY: pi
24  USE orbital_pointers, ONLY: coset,&
25  ncoset
26 #include "../base/base_uses.f90"
27 
28  IMPLICIT NONE
29 
30  PRIVATE
31 
32  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_kinetic'
33 
34 ! *** Public subroutines ***
35 
36  PUBLIC :: kinetic
37 
38 CONTAINS
39 
40 ! **************************************************************************************************
41 !> \brief Calculation of the two-center kinetic energy integrals [a|T|b] over
42 !> Cartesian Gaussian-type functions.
43 !> \param la_max Maximum L of basis on A
44 !> \param la_min Minimum L of basis on A
45 !> \param npgfa Number of primitive functions in set of basis on A
46 !> \param rpgfa Range of functions on A (used for prescreening)
47 !> \param zeta Exponents of basis on center A
48 !> \param lb_max Maximum L of basis on A
49 !> \param lb_min Minimum L of basis on A
50 !> \param npgfb Number of primitive functions in set of basis on B
51 !> \param rpgfb Range of functions on B (used for prescreening)
52 !> \param zetb Exponents of basis on center B
53 !> \param rab Distance vector between centers A and B
54 !> \param kab Kinetic energy integrals, optional
55 !> \param dab First derivatives of Kinetic energy integrals, optional
56 !> \date 07.07.2014
57 !> \author JGH
58 ! **************************************************************************************************
59  SUBROUTINE kinetic(la_max, la_min, npgfa, rpgfa, zeta, &
60  lb_max, lb_min, npgfb, rpgfb, zetb, &
61  rab, kab, dab)
62  INTEGER, INTENT(IN) :: la_max, la_min, npgfa
63  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
64  INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
65  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
66  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
67  REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT), &
68  OPTIONAL :: kab
69  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
70  OPTIONAL :: dab
71 
72  INTEGER :: ax, ay, az, bx, by, bz, coa, cob, ia, &
73  ib, idx, idy, idz, ipgf, jpgf, la, lb, &
74  ldrr, lma, lmb, ma, mb, na, nb, ofa, &
75  ofb
76  REAL(kind=dp) :: a, b, dsx, dsy, dsz, dtx, dty, dtz, f0, &
77  rab2, tab, xhi, zet
78  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: rr, tt
79  REAL(kind=dp), DIMENSION(3) :: rap, rbp
80 
81 ! Distance of the centers a and b
82 
83  rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
84  tab = sqrt(rab2)
85 
86  ! Maximum l for auxiliary integrals
87  IF (PRESENT(kab)) THEN
88  lma = la_max + 1
89  lmb = lb_max + 1
90  END IF
91  IF (PRESENT(dab)) THEN
92  lma = la_max + 2
93  lmb = lb_max + 1
94  idx = coset(1, 0, 0) - coset(0, 0, 0)
95  idy = coset(0, 1, 0) - coset(0, 0, 0)
96  idz = coset(0, 0, 1) - coset(0, 0, 0)
97  END IF
98  ldrr = max(lma, lmb) + 1
99 
100  ! Allocate space for auxiliary integrals
101  ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), tt(0:ldrr - 1, 0:ldrr - 1, 3))
102 
103  ! Number of integrals, check size of arrays
104  ofa = ncoset(la_min - 1)
105  ofb = ncoset(lb_min - 1)
106  na = ncoset(la_max) - ofa
107  nb = ncoset(lb_max) - ofb
108  IF (PRESENT(kab)) THEN
109  cpassert((SIZE(kab, 1) >= na*npgfa))
110  cpassert((SIZE(kab, 2) >= nb*npgfb))
111  END IF
112  IF (PRESENT(dab)) THEN
113  cpassert((SIZE(dab, 1) >= na*npgfa))
114  cpassert((SIZE(dab, 2) >= nb*npgfb))
115  cpassert((SIZE(dab, 3) >= 3))
116  END IF
117 
118  ! Loops over all pairs of primitive Gaussian-type functions
119  ma = 0
120  DO ipgf = 1, npgfa
121  mb = 0
122  DO jpgf = 1, npgfb
123  ! Distance Screening
124  IF (rpgfa(ipgf) + rpgfb(jpgf) < tab) THEN
125  IF (PRESENT(kab)) kab(ma + 1:ma + na, mb + 1:mb + nb) = 0.0_dp
126  IF (PRESENT(dab)) dab(ma + 1:ma + na, mb + 1:mb + nb, 1:3) = 0.0_dp
127  mb = mb + nb
128  cycle
129  END IF
130 
131  ! Calculate some prefactors
132  a = zeta(ipgf)
133  b = zetb(jpgf)
134  zet = a + b
135  xhi = a*b/zet
136  rap = b*rab/zet
137  rbp = -a*rab/zet
138 
139  ! [s|s] integral
140  f0 = 0.5_dp*(pi/zet)**(1.5_dp)*exp(-xhi*rab2)
141 
142  ! Calculate the recurrence relation, overlap
143  CALL os_rr_ovlp(rap, lma, rbp, lmb, zet, ldrr, rr)
144 
145  ! kinetic energy auxiliary integrals, overlap of [da/dx|db/dx]
146  DO la = 0, lma - 1
147  DO lb = 0, lmb - 1
148  tt(la, lb, 1) = 4.0_dp*a*b*rr(la + 1, lb + 1, 1)
149  tt(la, lb, 2) = 4.0_dp*a*b*rr(la + 1, lb + 1, 2)
150  tt(la, lb, 3) = 4.0_dp*a*b*rr(la + 1, lb + 1, 3)
151  IF (la > 0 .AND. lb > 0) THEN
152  tt(la, lb, 1) = tt(la, lb, 1) + real(la*lb, dp)*rr(la - 1, lb - 1, 1)
153  tt(la, lb, 2) = tt(la, lb, 2) + real(la*lb, dp)*rr(la - 1, lb - 1, 2)
154  tt(la, lb, 3) = tt(la, lb, 3) + real(la*lb, dp)*rr(la - 1, lb - 1, 3)
155  END IF
156  IF (la > 0) THEN
157  tt(la, lb, 1) = tt(la, lb, 1) - 2.0_dp*real(la, dp)*b*rr(la - 1, lb + 1, 1)
158  tt(la, lb, 2) = tt(la, lb, 2) - 2.0_dp*real(la, dp)*b*rr(la - 1, lb + 1, 2)
159  tt(la, lb, 3) = tt(la, lb, 3) - 2.0_dp*real(la, dp)*b*rr(la - 1, lb + 1, 3)
160  END IF
161  IF (lb > 0) THEN
162  tt(la, lb, 1) = tt(la, lb, 1) - 2.0_dp*real(lb, dp)*a*rr(la + 1, lb - 1, 1)
163  tt(la, lb, 2) = tt(la, lb, 2) - 2.0_dp*real(lb, dp)*a*rr(la + 1, lb - 1, 2)
164  tt(la, lb, 3) = tt(la, lb, 3) - 2.0_dp*real(lb, dp)*a*rr(la + 1, lb - 1, 3)
165  END IF
166  END DO
167  END DO
168 
169  DO lb = lb_min, lb_max
170  DO bx = 0, lb
171  DO by = 0, lb - bx
172  bz = lb - bx - by
173  cob = coset(bx, by, bz) - ofb
174  ib = mb + cob
175  DO la = la_min, la_max
176  DO ax = 0, la
177  DO ay = 0, la - ax
178  az = la - ax - ay
179  coa = coset(ax, ay, az) - ofa
180  ia = ma + coa
181  ! integrals
182  IF (PRESENT(kab)) THEN
183  kab(ia, ib) = f0*(tt(ax, bx, 1)*rr(ay, by, 2)*rr(az, bz, 3) + &
184  rr(ax, bx, 1)*tt(ay, by, 2)*rr(az, bz, 3) + &
185  rr(ax, bx, 1)*rr(ay, by, 2)*tt(az, bz, 3))
186  END IF
187  ! first derivatives
188  IF (PRESENT(dab)) THEN
189  ! dx
190  dsx = 2.0_dp*a*rr(ax + 1, bx, 1)
191  IF (ax > 0) dsx = dsx - real(ax, dp)*rr(ax - 1, bx, 1)
192  dtx = 2.0_dp*a*tt(ax + 1, bx, 1)
193  IF (ax > 0) dtx = dtx - real(ax, dp)*tt(ax - 1, bx, 1)
194  dab(ia, ib, idx) = dtx*rr(ay, by, 2)*rr(az, bz, 3) + &
195  dsx*(tt(ay, by, 2)*rr(az, bz, 3) + rr(ay, by, 2)*tt(az, bz, 3))
196  ! dy
197  dsy = 2.0_dp*a*rr(ay + 1, by, 2)
198  IF (ay > 0) dsy = dsy - real(ay, dp)*rr(ay - 1, by, 2)
199  dty = 2.0_dp*a*tt(ay + 1, by, 2)
200  IF (ay > 0) dty = dty - real(ay, dp)*tt(ay - 1, by, 2)
201  dab(ia, ib, idy) = dty*rr(ax, bx, 1)*rr(az, bz, 3) + &
202  dsy*(tt(ax, bx, 1)*rr(az, bz, 3) + rr(ax, bx, 1)*tt(az, bz, 3))
203  ! dz
204  dsz = 2.0_dp*a*rr(az + 1, bz, 3)
205  IF (az > 0) dsz = dsz - real(az, dp)*rr(az - 1, bz, 3)
206  dtz = 2.0_dp*a*tt(az + 1, bz, 3)
207  IF (az > 0) dtz = dtz - real(az, dp)*tt(az - 1, bz, 3)
208  dab(ia, ib, idz) = dtz*rr(ax, bx, 1)*rr(ay, by, 2) + &
209  dsz*(tt(ax, bx, 1)*rr(ay, by, 2) + rr(ax, bx, 1)*tt(ay, by, 2))
210  ! scale
211  dab(ia, ib, 1:3) = f0*dab(ia, ib, 1:3)
212  END IF
213  !
214  END DO
215  END DO
216  END DO !la
217  END DO
218  END DO
219  END DO !lb
220 
221  mb = mb + nb
222  END DO
223  ma = ma + na
224  END DO
225 
226  DEALLOCATE (rr, tt)
227 
228  END SUBROUTINE kinetic
229 
230 END MODULE ai_kinetic
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition: grid_common.h:153
ai_kinetic
Calculation of the kinetic energy integrals over Cartesian Gaussian-type functions.
Definition: ai_kinetic.F:20
ai_kinetic::kinetic
subroutine, public kinetic(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, kab, dab)
Calculation of the two-center kinetic energy integrals [a|T|b] over Cartesian Gaussian-type functions...
Definition: ai_kinetic.F:62
ai_os_rr
Definition: ai_os_rr.F:8
ai_os_rr::os_rr_ovlp
subroutine, public os_rr_ovlp(rap, la_max, rbp, lb_max, zet, ldrr, rr)
Calculation of the basic Obara-Saika recurrence relation.
Definition: ai_os_rr.F:39
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition: orbital_pointers.F:45