Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of spin orbit integrals over Cartesian Gaussian-type functions
10 : !> \par Literature
11 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
12 : !> \par History
13 : !> none
14 : !> \par Parameters
15 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
16 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
17 : !> - coset : Cartesian orbital set pointer.
18 : !> - dab : Distance between the atomic centers a and b.
19 : !> - dac : Distance between the atomic centers a and c.
20 : !> - dbc : Distance between the atomic centers b and c.
21 : !> - l{a,b,c} : Angular momentum quantum number of shell a, b or c.
22 : !> - l{a,b,c}_max: Maximum angular momentum quantum number of shell a, b or c.
23 : !> - l{a,b,c}_min: Minimum angular momentum quantum number of shell a, b or c.
24 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
25 : !> - npgf{a,b} : Degree of contraction of shell a or b.
26 : !> - rab : Distance vector between the atomic centers a and b.
27 : !> - rab2 : Square of the distance between the atomic centers a and b.
28 : !> - rac : Distance vector between the atomic centers a and c.
29 : !> - rac2 : Square of the distance between the atomic centers a and c.
30 : !> - rbc : Distance vector between the atomic centers b and c.
31 : !> - rbc2 : Square of the distance between the atomic centers b and c.
32 : !> - rpgf{a,b,c} : Radius of the primitive Gaussian-type function a, b or c.
33 : !> - zet{a,b,c} : Exponents of the Gaussian-type functions a, b or c.
34 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
35 : !> \author VW
36 : ! **************************************************************************************************
37 : MODULE ai_spin_orbit
38 : USE ai_os_rr, ONLY: os_rr_coul
39 : USE kinds, ONLY: dp
40 : USE mathconstants, ONLY: pi
41 : USE orbital_pointers, ONLY: coset,&
42 : ncoset
43 : #include "../base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_spin_orbit'
47 : PRIVATE
48 :
49 : ! *** Public subroutines ***
50 :
51 : PUBLIC :: pso
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Calculation of the primitive paramagnetic spin orbit integrals over
57 : !> Cartesian Gaussian-type functions.
58 : !> \param la_max ...
59 : !> \param la_min ...
60 : !> \param npgfa ...
61 : !> \param rpgfa ...
62 : !> \param zeta ...
63 : !> \param lb_max ...
64 : !> \param lb_min ...
65 : !> \param npgfb ...
66 : !> \param rpgfb ...
67 : !> \param zetb ...
68 : !> \param rac ...
69 : !> \param rbc ...
70 : !> \param rab ...
71 : !> \param vab ...
72 : !> \param ldrr1 ...
73 : !> \param ldrr2 ...
74 : !> \param rr ...
75 : !> \date 02.03.2009
76 : !> \author VW
77 : !> \version 1.0
78 : ! **************************************************************************************************
79 81476 : SUBROUTINE pso(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, &
80 40738 : rac, rbc, rab, vab, ldrr1, ldrr2, rr)
81 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
82 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
83 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
84 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
85 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc, rab
86 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: vab
87 : INTEGER, INTENT(IN) :: ldrr1, ldrr2
88 : REAL(dp), DIMENSION(0:ldrr1-1, ldrr2, *), &
89 : INTENT(INOUT) :: rr
90 :
91 : INTEGER :: ax, ay, az, bx, by, bz, coa, coam1x, coam1y, coam1z, coap1x, coap1y, coap1z, cob, &
92 : cobm1x, cobm1y, cobm1z, cobp1x, cobp1y, cobp1z, i, ipgf, j, jpgf, la, lb, ma, mb, na, nb
93 : REAL(dp) :: dab, dum1, dum2, f0, rab2, xhi, zet, &
94 : zetab
95 : REAL(dp), DIMENSION(3) :: rap, rbp, rcp
96 :
97 : ! *** Calculate the distance of the centers a and c ***
98 :
99 40738 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
100 40738 : dab = SQRT(rab2)
101 :
102 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
103 :
104 40738 : na = 0
105 :
106 92858 : DO ipgf = 1, npgfa
107 :
108 52120 : nb = 0
109 :
110 118782 : DO jpgf = 1, npgfb
111 :
112 : ! *** Screening ***
113 :
114 66662 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
115 29344 : DO j = nb + 1, nb + ncoset(lb_max)
116 52413 : DO i = na + 1, na + ncoset(la_max)
117 23069 : vab(i, j, 1) = 0.0_dp
118 23069 : vab(i, j, 2) = 0.0_dp
119 40267 : vab(i, j, 3) = 0.0_dp
120 : END DO
121 : END DO
122 12146 : nb = nb + ncoset(lb_max)
123 12146 : CYCLE
124 : END IF
125 :
126 : ! *** Calculate some prefactors ***
127 :
128 54516 : zetab = zeta(ipgf)*zetb(jpgf)
129 54516 : zet = zeta(ipgf) + zetb(jpgf)
130 54516 : xhi = zetab/zet
131 218064 : rap = zetb(jpgf)*rab/zet
132 218064 : rbp = -zeta(ipgf)*rab/zet
133 218064 : rcp = -(zeta(ipgf)*rac + zetb(jpgf)*rbc)/zet
134 :
135 54516 : f0 = 2.0_dp*SQRT(zet/pi)*(pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
136 :
137 : ! *** Calculate the recurrence relation ***
138 :
139 54516 : CALL os_rr_coul(rap, la_max + 1, rbp, lb_max + 1, rcp, zet, ldrr1, ldrr2, rr)
140 :
141 : ! *** Calculate the primitive Fermi contact integrals ***
142 :
143 121388 : DO lb = lb_min, lb_max
144 200627 : DO bx = 0, lb
145 237717 : DO by = 0, lb - bx
146 91606 : bz = lb - bx - by
147 91606 : cob = coset(bx, by, bz)
148 91606 : cobm1x = coset(MAX(bx - 1, 0), by, bz)
149 91606 : cobm1y = coset(bx, MAX(by - 1, 0), bz)
150 91606 : cobm1z = coset(bx, by, MAX(bz - 1, 0))
151 91606 : cobp1x = coset(bx + 1, by, bz)
152 91606 : cobp1y = coset(bx, by + 1, bz)
153 91606 : cobp1z = coset(bx, by, bz + 1)
154 91606 : mb = nb + cob
155 287413 : DO la = la_min, la_max
156 349717 : DO ax = 0, la
157 424629 : DO ay = 0, la - ax
158 166518 : az = la - ax - ay
159 166518 : coa = coset(ax, ay, az)
160 166518 : coam1x = coset(MAX(ax - 1, 0), ay, az)
161 166518 : coam1y = coset(ax, MAX(ay - 1, 0), az)
162 166518 : coam1z = coset(ax, ay, MAX(az - 1, 0))
163 166518 : coap1x = coset(ax + 1, ay, az)
164 166518 : coap1y = coset(ax, ay + 1, az)
165 166518 : coap1z = coset(ax, ay, az + 1)
166 166518 : ma = na + coa
167 : !
168 : !
169 : ! (a|pso_x|b) = (4*zeta*zetb*(a+y||b+z)
170 : ! -2*zeta*Nz(b)*(a+y||b-z)-2*zetb*Ny(a)*(a-y||b+z)
171 : ! +Ny(a)*Nz(b)*(a-y||b-z))
172 : ! -(4*zeta*zetb*(a+z||b+y)
173 : ! -2*zeta*Ny(b)*(a+z||b-y)-2*zetb*Nz(a)*(a-z||b+y)
174 : ! +Nz(a)*Ny(b)*(a-z||b-y))
175 166518 : dum1 = 4.0_dp*zeta(ipgf)*zetb(jpgf)*rr(0, coap1y, cobp1z)
176 166518 : IF (bz > 0) dum1 = dum1 - 2.0_dp*zeta(ipgf)*REAL(bz, dp)*rr(0, coap1y, cobm1z)
177 166518 : IF (ay > 0) dum1 = dum1 - 2.0_dp*zetb(jpgf)*REAL(ay, dp)*rr(0, coam1y, cobp1z)
178 166518 : IF (ay > 0 .AND. bz > 0) dum1 = dum1 + REAL(ay, dp)*REAL(bz, dp)*rr(0, coam1y, cobm1z)
179 : !
180 166518 : dum2 = 4.0_dp*zeta(ipgf)*zetb(jpgf)*rr(0, coap1z, cobp1y)
181 166518 : IF (by > 0) dum2 = dum2 - 2.0_dp*zeta(ipgf)*REAL(by, dp)*rr(0, coap1z, cobm1y)
182 166518 : IF (az > 0) dum2 = dum2 - 2.0_dp*zetb(jpgf)*REAL(az, dp)*rr(0, coam1z, cobp1y)
183 166518 : IF (az > 0 .AND. by > 0) dum2 = dum2 + REAL(az, dp)*REAL(by, dp)*rr(0, coam1z, cobm1y)
184 166518 : vab(ma, mb, 1) = f0*(dum1 - dum2)
185 : !
186 : !
187 : ! (a|pso_y|b) = (4*zeta*zetb*(a+z||b+x)
188 : ! -2*zeta*Nx(b)*(a+z||b-x)-2*zetb*Nz(a)*(a-z||b+x)
189 : ! +Nz(a)*Nx(b)*(a-z||b-x))
190 : ! -(4*zeta*zetb*(a+x||b+z)
191 : ! -2*zeta*Nz(b)*(a+x||b-z)-2*zetb*Nx(a)*(a-x||b+z)
192 : ! +Nx(a)*Nz(b)*(a-x||b-z))
193 166518 : dum1 = 4.0_dp*zeta(ipgf)*zetb(jpgf)*rr(0, coap1z, cobp1x)
194 166518 : IF (bx > 0) dum1 = dum1 - 2.0_dp*zeta(ipgf)*REAL(bx, dp)*rr(0, coap1z, cobm1x)
195 166518 : IF (az > 0) dum1 = dum1 - 2.0_dp*zetb(jpgf)*REAL(az, dp)*rr(0, coam1z, cobp1x)
196 166518 : IF (az > 0 .AND. bx > 0) dum1 = dum1 + REAL(az, dp)*REAL(bx, dp)*rr(0, coam1z, cobm1x)
197 : !
198 166518 : dum2 = 4.0_dp*zeta(ipgf)*zetb(jpgf)*rr(0, coap1x, cobp1z)
199 166518 : IF (bz > 0) dum2 = dum2 - 2.0_dp*zeta(ipgf)*REAL(bz, dp)*rr(0, coap1x, cobm1z)
200 166518 : IF (ax > 0) dum2 = dum2 - 2.0_dp*zetb(jpgf)*REAL(ax, dp)*rr(0, coam1x, cobp1z)
201 166518 : IF (ax > 0 .AND. bz > 0) dum2 = dum2 + REAL(ax, dp)*REAL(bz, dp)*rr(0, coam1x, cobm1z)
202 166518 : vab(ma, mb, 2) = f0*(dum1 - dum2)
203 : !
204 : !
205 : ! (a|pso_z|b) = (4*zeta*zetb*(a+x||b+y)
206 : ! -2*zeta*Ny(b)*(a+x||b-y)-2*zetb*Nx(a)*(a-x||b+y)
207 : ! +Nx(a)*Ny(b)*(a-x||b-y))
208 : ! -(4*zeta*zetb*(a+y||b+x)
209 : ! -2*zeta*Nx(b)*(a+y||b-x)-2*zetb*Ny(a)*(a-y||b+x)
210 : ! +Ny(a)*Nx(b)*(a-y||b-x))
211 166518 : dum1 = 4.0_dp*zeta(ipgf)*zetb(jpgf)*rr(0, coap1x, cobp1y)
212 166518 : IF (by > 0) dum1 = dum1 - 2.0_dp*zeta(ipgf)*REAL(by, dp)*rr(0, coap1x, cobm1y)
213 166518 : IF (ax > 0) dum1 = dum1 - 2.0_dp*zetb(jpgf)*REAL(ax, dp)*rr(0, coam1x, cobp1y)
214 166518 : IF (ax > 0 .AND. by > 0) dum1 = dum1 + REAL(ax, dp)*REAL(by, dp)*rr(0, coam1x, cobm1y)
215 : !
216 166518 : dum2 = 4.0_dp*zeta(ipgf)*zetb(jpgf)*rr(0, coap1y, cobp1x)
217 166518 : IF (bx > 0) dum2 = dum2 - 2.0_dp*zeta(ipgf)*REAL(bx, dp)*rr(0, coap1y, cobm1x)
218 166518 : IF (ay > 0) dum2 = dum2 - 2.0_dp*zetb(jpgf)*REAL(ay, dp)*rr(0, coam1y, cobp1x)
219 166518 : IF (ay > 0 .AND. bx > 0) dum2 = dum2 + REAL(ay, dp)*REAL(bx, dp)*rr(0, coam1y, cobm1x)
220 308061 : vab(ma, mb, 3) = f0*(dum1 - dum2)
221 : !
222 : END DO
223 : END DO
224 : END DO !la
225 :
226 : END DO
227 : END DO
228 : END DO !lb
229 :
230 106636 : nb = nb + ncoset(lb_max)
231 :
232 : END DO
233 :
234 92858 : na = na + ncoset(la_max)
235 :
236 : END DO
237 :
238 40738 : END SUBROUTINE pso
239 :
240 : END MODULE ai_spin_orbit
|
