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 integrals over Cartesian Gaussian-type functions for [a|(r-Ra)^(2m)|b]
10 : !> Ra is the position of center a
11 : !> \par Literature
12 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 : !> \par History
14 : !> none
15 : !> \par Parameters
16 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
17 : !> - cx,cy,cz : Angular momentum index numbers of orbital c.
18 : !> - coset : Cartesian orbital set pointer.
19 : !> - dac : Distance between the atomic centers a and c.
20 : !> - l{a,c} : Angular momentum quantum number of shell a or c.
21 : !> - l{a,c}_max : Maximum angular momentum quantum number of shell a or c.
22 : !> - l{a,c}_min : Minimum angular momentum quantum number of shell a or c.
23 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
24 : !> - npgf{a,c} : Degree of contraction of shell a or c.
25 : !> - rac : Distance vector between the atomic centers a and c.
26 : !> - rac2 : Square of the distance between the atomic centers a and c.
27 : !> - zet{a,c} : Exponents of the Gaussian-type functions a or c.
28 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
29 : !> - zetw : Reciprocal of the sum of the exponents of orbital a and c.
30 : !> \author Dorothea Golze (08.2016)
31 : ! **************************************************************************************************
32 : MODULE ai_operator_ra2m
33 :
34 : USE ai_cartesian_operators, ONLY: apply_cartesian_1d,&
35 : cartesian_polynomial_1d
36 : USE ai_os_rr, ONLY: os_rr_ovlp
37 : USE kinds, ONLY: dp
38 : USE mathconstants, ONLY: fac,&
39 : pi
40 : USE orbital_pointers, ONLY: indco,&
41 : ncoset
42 : #include "../base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_operator_ra2m'
46 :
47 : PRIVATE
48 :
49 : ! *** Public subroutines ***
50 :
51 : PUBLIC :: operator_ra2m
52 :
53 : ! **************************************************************************************************
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief Calculation of the primitive two-center [a|(r-Ra)^(2m)|b] integrals over Cartesian
59 : !> Gaussian-type functions; operator is here (r-Ra)^(2m)
60 : !> \param la_max ...
61 : !> \param la_min ...
62 : !> \param npgfa ...
63 : !> \param zeta ...
64 : !> \param lb_max ...
65 : !> \param lb_min ...
66 : !> \param npgfb ...
67 : !> \param zetb ...
68 : !> \param m exponent in (r-Ra)^(2m) operator
69 : !> \param rab ...
70 : !> \param sab ...
71 : !> \param dsab ...
72 : !> \param calculate_forces ...
73 : ! **************************************************************************************************
74 14400 : SUBROUTINE operator_ra2m(la_max, la_min, npgfa, zeta, &
75 4800 : lb_max, lb_min, npgfb, zetb, m, rab, sab, dsab, calculate_forces)
76 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
77 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta
78 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
79 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb
80 : INTEGER, INTENT(IN) :: m
81 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
82 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: sab
83 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: dsab
84 : LOGICAL, INTENT(IN) :: calculate_forces
85 :
86 : CHARACTER(LEN=*), PARAMETER :: routineN = 'operator_ra2m'
87 :
88 : INTEGER :: ax, ay, az, bx, by, bz, coa, cob, &
89 : handle, i, ia, ib, ipgf, j, jpgf, k, &
90 : ldrr, lma, ma, mb, na, nb, ofa, ofb, &
91 9600 : powers(0:m)
92 : REAL(KIND=dp) :: a, b, f0, prefac, rab2, rap(3), rbp(3), &
93 9600 : weights(0:m, 0:m), zet
94 4800 : REAL(KIND=dp), ALLOCATABLE :: dpolynomial(:, :, :, :), drr(:, :, :), &
95 4800 : polynomial(:, :, :, :), rr(:, :, :)
96 :
97 4800 : CALL timeset(routineN, handle)
98 6052800 : sab = 0.0_dp
99 18163200 : IF (calculate_forces) dsab = 0.0_dp
100 19200 : rab2 = NORM2(rab)**2
101 4800 : lma = la_max + 2*m
102 4800 : IF (calculate_forces) lma = lma + 1
103 4800 : ldrr = MAX(lma, lb_max) + 1
104 48000 : ALLOCATE (rr(0:ldrr - 1, 0:ldrr - 1, 3), polynomial(0:la_max, 0:lb_max, 0:m, 3))
105 4800 : IF (calculate_forces) THEN
106 28800 : ALLOCATE (drr(0:ldrr - 1, 0:ldrr - 1, 3), dpolynomial(0:la_max, 0:lb_max, 0:m, 3))
107 : END IF
108 43200 : powers = [(2*i, i=0, m)]
109 24000 : DO i = 0, m
110 72000 : DO j = 0, m - i
111 48000 : k = m - i - j
112 67200 : weights(i, j) = fac(m)/fac(i)/fac(j)/fac(k)
113 : END DO
114 : END DO
115 4800 : ofa = ncoset(la_min - 1)
116 4800 : ofb = ncoset(lb_min - 1)
117 4800 : na = ncoset(la_max) - ofa
118 4800 : nb = ncoset(lb_max) - ofb
119 4800 : CPASSERT(SIZE(sab, 1) >= na*npgfa .AND. SIZE(sab, 2) >= nb*npgfb)
120 :
121 9600 : DO ipgf = 1, npgfa
122 4800 : ma = (ipgf - 1)*na
123 14400 : DO jpgf = 1, npgfb
124 4800 : mb = (jpgf - 1)*nb
125 4800 : a = zeta(ipgf)
126 4800 : b = zetb(jpgf)
127 4800 : zet = a + b
128 19200 : rap = b*rab/zet
129 19200 : rbp = -a*rab/zet
130 4800 : f0 = (pi/zet)**1.5_dp*EXP(-a*b*rab2/zet)
131 4800 : CALL os_rr_ovlp(rap, lma, rbp, lb_max, zet, ldrr, rr)
132 4800 : CALL cartesian_polynomial_1d(rr, polynomial, la_max, lb_max, powers, [0.0_dp, 0.0_dp, 0.0_dp])
133 4800 : IF (calculate_forces) THEN
134 : ! The origin moves with A: differentiate the combined Gaussian degree first.
135 : ! Raising that derivative includes both basis and operator-origin contributions.
136 : CALL apply_cartesian_1d(rr, drr, la_max + 2*m, lb_max, .FALSE., .TRUE., a, &
137 4800 : [0.0_dp, 0.0_dp, 0.0_dp])
138 4800 : CALL cartesian_polynomial_1d(drr, dpolynomial, la_max, lb_max, powers, [0.0_dp, 0.0_dp, 0.0_dp])
139 : END IF
140 :
141 118400 : DO cob = ofb + 1, ncoset(lb_max)
142 108800 : bx = indco(1, cob)
143 108800 : by = indco(2, cob)
144 108800 : bz = indco(3, cob)
145 108800 : ib = mb + cob - ofb
146 1691200 : DO coa = ofa + 1, ncoset(la_max)
147 1577600 : ax = indco(1, coa)
148 1577600 : ay = indco(2, coa)
149 1577600 : az = indco(3, coa)
150 1577600 : ia = ma + coa - ofa
151 7996800 : DO i = 0, m
152 23664000 : DO j = 0, m - i
153 15776000 : k = m - i - j
154 15776000 : prefac = weights(i, j)*f0
155 : sab(ia, ib) = sab(ia, ib) + prefac*polynomial(ax, bx, i, 1)* &
156 15776000 : polynomial(ay, by, j, 2)*polynomial(az, bz, k, 3)
157 22086400 : IF (calculate_forces) THEN
158 : dsab(ia, ib, 1) = dsab(ia, ib, 1) + prefac*dpolynomial(ax, bx, i, 1)* &
159 15776000 : polynomial(ay, by, j, 2)*polynomial(az, bz, k, 3)
160 : dsab(ia, ib, 2) = dsab(ia, ib, 2) + prefac*polynomial(ax, bx, i, 1)* &
161 15776000 : dpolynomial(ay, by, j, 2)*polynomial(az, bz, k, 3)
162 : dsab(ia, ib, 3) = dsab(ia, ib, 3) + prefac*polynomial(ax, bx, i, 1)* &
163 15776000 : polynomial(ay, by, j, 2)*dpolynomial(az, bz, k, 3)
164 : END IF
165 : END DO
166 : END DO
167 : END DO
168 : END DO
169 : END DO
170 : END DO
171 4800 : CALL timestop(handle)
172 9600 : END SUBROUTINE operator_ra2m
173 :
174 : END MODULE ai_operator_ra2m
|
