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 general three-center integrals over Cartesian
10 : !> Gaussian-type functions and a spherical operator centered at position C
11 : !>
12 : !> <a|V(local)|b> = <a|F(|r-C|)|b>
13 : !> \par Literature
14 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
15 : !> \par History
16 : !> - Based in part on code by MK
17 : !> \par Parameters
18 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
19 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
20 : !> - coset : Cartesian orbital set pointer.
21 : !> - dab : Distance between the atomic centers a and b.
22 : !> - dac : Distance between the atomic centers a and c.
23 : !> - dbc : Distance between the atomic centers b and c.
24 : !> - l{a,b} : Angular momentum quantum number of shell a or b.
25 : !> - l{a,b}_max : Maximum angular momentum quantum number of shell a or b.
26 : !> - ncoset : Number of Cartesian orbitals up to l.
27 : !> - rab : Distance vector between the atomic centers a and b.
28 : !> - rac : Distance vector between the atomic centers a and c.
29 : !> - rbc : Distance vector between the atomic centers b and c.
30 : !> - rpgf{a,b,c}: Radius of the primitive Gaussian-type function a or b.
31 : !> - zet{a,b} : Exponents of the Gaussian-type functions a or b.
32 : !> \author jhu (05.2011)
33 : ! **************************************************************************************************
34 : MODULE ai_oneelectron
35 :
36 : USE kinds, ONLY: dp
37 : USE orbital_pointers, ONLY: coset,&
38 : ncoset
39 : #include "../base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_oneelectron'
46 :
47 : ! *** Public subroutines ***
48 :
49 : PUBLIC :: os_3center, os_2center
50 :
51 : CONTAINS
52 :
53 : ! **************************************************************************************************
54 : !> \brief Calculation of three-center integrals <a|c|b> over
55 : !> Cartesian Gaussian functions and a spherical potential
56 : !>
57 : !> \param la_max_set ...
58 : !> \param la_min_set ...
59 : !> \param npgfa ...
60 : !> \param rpgfa ...
61 : !> \param zeta ...
62 : !> \param lb_max_set ...
63 : !> \param lb_min_set ...
64 : !> \param npgfb ...
65 : !> \param rpgfb ...
66 : !> \param zetb ...
67 : !> \param auxint ...
68 : !> \param rpgfc ...
69 : !> \param rab ...
70 : !> \param dab ...
71 : !> \param rac ...
72 : !> \param dac ...
73 : !> \param rbc ...
74 : !> \param dbc ...
75 : !> \param vab ...
76 : !> \param s ...
77 : !> \param pab ...
78 : !> \param force_a ...
79 : !> \param force_b ...
80 : !> \param fs ...
81 : !> \param vab2 The derivative of the 3-center integrals according to the weighting factors.
82 : !> \param vab2_work ...
83 : !> \param deltaR DIMENSION(3, natoms), weighting factors of the derivatives for each atom and direction
84 : !> \param iatom ...
85 : !> \param jatom ...
86 : !> \param katom ...
87 : !> \date May 2011
88 : !> \author Juerg Hutter
89 : !> \version 1.0
90 : !> \note Extended by the derivatives for DFPT [Sandra Luber, Edward Ditler, 2021]
91 : ! **************************************************************************************************
92 34945148 : SUBROUTINE os_3center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
93 34945148 : lb_max_set, lb_min_set, npgfb, rpgfb, zetb, auxint, rpgfc, &
94 69890296 : rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, &
95 69890296 : vab2, vab2_work, deltaR, iatom, jatom, katom)
96 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
97 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
98 : INTEGER, INTENT(IN) :: lb_max_set, lb_min_set, npgfb
99 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
100 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: auxint
101 : REAL(KIND=dp), INTENT(IN) :: rpgfc
102 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
103 : REAL(KIND=dp), INTENT(IN) :: dab
104 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
105 : REAL(KIND=dp), INTENT(IN) :: dac
106 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc
107 : REAL(KIND=dp), INTENT(IN) :: dbc
108 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: vab
109 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: s
110 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
111 : OPTIONAL :: pab
112 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: force_a, force_b
113 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
114 : OPTIONAL :: fs, vab2, vab2_work
115 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
116 : OPTIONAL :: deltaR
117 : INTEGER, INTENT(IN), OPTIONAL :: iatom, jatom, katom
118 :
119 : INTEGER :: ax, ay, az, bx, by, bz, cda, cdax, cday, cdaz, cdb, cdbx, cdby, cdbz, coa, coamx, &
120 : coamy, coamz, coapx, coapy, coapz, cob, cobmx, cobmy, cobmz, cobpx, cobpy, cobpz, da, &
121 : da_max, dax, day, daz, db, db_max, dbx, dby, dbz, i, ia, iap, iax, iay, iaz, ib, ibm, &
122 : ibx, iby, ibz, idir, ii(3), iim(3), ij, ipgf, ir, ir1, ir2, irm(3), irr(3), irx, iry, &
123 : irz, ix, ixx(1), j, jj(3), jjp(3), jpgf, la, la_max, la_min, lb, lb_max, lb_min, llr, m, &
124 : ma, mb, mmax, na, nb
125 34945148 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: iiap
126 : LOGICAL :: calculate_force_a, calculate_force_b
127 : REAL(KIND=dp) :: aai, abx, fax, fay, faz, fbx, fby, fbz, &
128 : ftz, orho, rho, s1, s2
129 : REAL(KIND=dp), DIMENSION(3) :: pai, pbi, pci
130 :
131 34945148 : IF (PRESENT(pab)) THEN
132 8963078 : CPASSERT(PRESENT(fs))
133 8963078 : IF (PRESENT(force_a)) THEN
134 : calculate_force_a = .TRUE.
135 : ELSE
136 0 : calculate_force_a = .FALSE.
137 : END IF
138 8963078 : IF (PRESENT(force_b)) THEN
139 : calculate_force_b = .TRUE.
140 : ELSE
141 0 : calculate_force_b = .FALSE.
142 : END IF
143 : ELSE
144 : calculate_force_a = .FALSE.
145 : calculate_force_b = .FALSE.
146 : END IF
147 :
148 8963078 : IF (calculate_force_a) THEN
149 8963078 : da_max = 1
150 8963078 : force_a = 0.0_dp
151 : ELSE
152 : da_max = 0
153 : END IF
154 :
155 34945148 : IF (calculate_force_b) THEN
156 8963078 : db_max = 1
157 8963078 : force_b = 0.0_dp
158 : ELSE
159 : db_max = 0
160 : END IF
161 :
162 34945148 : IF (PRESENT(vab2)) THEN
163 870 : da_max = 1
164 870 : db_max = 1
165 : END IF
166 :
167 34945148 : la_max = la_max_set + da_max
168 34945148 : la_min = MAX(0, la_min_set - da_max)
169 :
170 34945148 : lb_max = lb_max_set + db_max
171 34945148 : lb_min = MAX(0, lb_min_set - db_max)
172 :
173 34945148 : mmax = la_max + lb_max
174 :
175 : ! precalculate indices for horizontal recursion
176 104835444 : ALLOCATE (iiap(ncoset(mmax), 3))
177 159075239 : DO ma = 0, mmax
178 478050123 : DO iax = 0, ma
179 1134730830 : DO iay = 0, ma - iax
180 691625855 : iaz = ma - iax - iay
181 691625855 : ia = coset(iax, iay, iaz)
182 691625855 : jj(1) = iax; jj(2) = iay; jj(3) = iaz
183 691625855 : jjp = jj
184 691625855 : jjp(1) = jjp(1) + 1
185 691625855 : iap = coset(jjp(1), jjp(2), jjp(3))
186 691625855 : iiap(ia, 1) = iap
187 691625855 : jjp = jj
188 691625855 : jjp(2) = jjp(2) + 1
189 691625855 : iap = coset(jjp(1), jjp(2), jjp(3))
190 691625855 : iiap(ia, 2) = iap
191 691625855 : jjp = jj
192 691625855 : jjp(3) = jjp(3) + 1
193 691625855 : iap = coset(jjp(1), jjp(2), jjp(3))
194 1010600739 : iiap(ia, 3) = iap
195 : END DO
196 : END DO
197 : END DO
198 :
199 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
200 :
201 34945148 : na = 0
202 :
203 154913027 : DO ipgf = 1, npgfa
204 :
205 : ! *** Screening ***
206 119967879 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
207 75249771 : na = na + ncoset(la_max_set)
208 75249771 : CYCLE
209 : END IF
210 :
211 44718108 : nb = 0
212 :
213 204272435 : DO jpgf = 1, npgfb
214 :
215 : ! *** Screening ***
216 159554327 : IF ((rpgfb(jpgf) + rpgfc < dbc) .OR. &
217 : (rpgfa(ipgf) + rpgfb(jpgf) < dab)) THEN
218 101232176 : nb = nb + ncoset(lb_max_set)
219 101232176 : CYCLE
220 : END IF
221 :
222 : ! *** Calculate some prefactors ***
223 58322151 : rho = zeta(ipgf) + zetb(jpgf)
224 233288604 : pai(:) = zetb(jpgf)/rho*rab(:)
225 233288604 : pbi(:) = -zeta(ipgf)/rho*rab(:)
226 233288604 : pci(:) = -(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:))/rho
227 58322151 : orho = 0.5_dp/rho
228 :
229 58322151 : ij = (ipgf - 1)*npgfb + jpgf
230 270671515 : s(1, 1, 1:mmax + 1) = auxint(0:mmax, ij)
231 :
232 58322151 : IF (la_max > 0) THEN
233 : ! *** Recurrence steps: [s|c|s] -> [a|c|s] ***
234 : ! *** [a|c|s](m) = (Pi - Ai)*[a-1i|c|s](m) - ***
235 : ! *** (Pi - Ci)*[a-1i|c|s](m+1)) + ***
236 : ! *** Ni(a-1i)/2(a+b)*[a-2i|c|s](m) - ***
237 : ! *** Ni(a-1i)/2(a+b)*[a-2i|c|s](m+1) ***
238 202789724 : DO llr = 1, mmax
239 202789724 : IF (llr == 1) THEN
240 202789724 : DO m = 0, mmax - llr
241 151673210 : s1 = s(1, 1, m + 1)
242 151673210 : s2 = s(1, 1, m + 2)
243 151673210 : s(2, 1, m + 1) = pai(1)*s1 - pci(1)*s2 ! [px|o|s]
244 151673210 : s(3, 1, m + 1) = pai(2)*s1 - pci(2)*s2 ! [py|o|s]
245 202789724 : s(4, 1, m + 1) = pai(3)*s1 - pci(3)*s2 ! [pz|o|s]
246 : END DO
247 100556696 : ELSE IF (llr == 2) THEN
248 149791499 : DO m = 0, mmax - llr
249 100556696 : s1 = s(1, 1, m + 1) - s(1, 1, m + 2)
250 100556696 : s(5, 1, m + 1) = pai(1)*s(2, 1, m + 1) - pci(1)*s(2, 1, m + 2) + orho*s1 ! [dx2|o|s]
251 100556696 : s(6, 1, m + 1) = pai(1)*s(3, 1, m + 1) - pci(1)*s(3, 1, m + 2) ! [dxy|o|s]
252 100556696 : s(7, 1, m + 1) = pai(1)*s(4, 1, m + 1) - pci(1)*s(4, 1, m + 2) ! [dxz|o|s]
253 100556696 : s(8, 1, m + 1) = pai(2)*s(3, 1, m + 1) - pci(2)*s(3, 1, m + 2) + orho*s1 ! [dy2|o|s]
254 100556696 : s(9, 1, m + 1) = pai(2)*s(4, 1, m + 1) - pci(2)*s(4, 1, m + 2) ! [dyz|o|s]
255 149791499 : s(10, 1, m + 1) = pai(3)*s(4, 1, m + 1) - pci(3)*s(4, 1, m + 2) + orho*s1 ! [dz2|o|s]
256 : END DO
257 51321893 : ELSE IF (llr == 3) THEN
258 77230814 : DO m = 0, mmax - llr
259 : s(11, 1, m + 1) = pai(1)*s(5, 1, m + 1) - pci(1)*s(5, 1, m + 2) & ! [fx3 |o|s]
260 51321893 : + 2._dp*orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
261 : s(12, 1, m + 1) = pai(1)*s(6, 1, m + 1) - pci(1)*s(6, 1, m + 2) & ! [fx2y|o|s]
262 51321893 : + orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
263 : s(13, 1, m + 1) = pai(1)*s(7, 1, m + 1) - pci(1)*s(7, 1, m + 2) & ! [fx2z|o|s]
264 51321893 : + orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
265 : s(14, 1, m + 1) = pai(2)*s(6, 1, m + 1) - pci(2)*s(6, 1, m + 2) & ! [fxy2|o|s]
266 51321893 : + orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
267 51321893 : s(15, 1, m + 1) = pai(1)*s(9, 1, m + 1) - pci(1)*s(9, 1, m + 2) ! [fxyz|o|s]
268 : s(16, 1, m + 1) = pai(3)*s(7, 1, m + 1) - pci(3)*s(7, 1, m + 2) & ! [fxz2|o|s]
269 51321893 : + orho*(s(2, 1, m + 1) - s(2, 1, m + 2))
270 : s(17, 1, m + 1) = pai(2)*s(8, 1, m + 1) - pci(2)*s(8, 1, m + 2) & ! [fy3 |o|s]
271 51321893 : + 2._dp*orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
272 : s(18, 1, m + 1) = pai(2)*s(9, 1, m + 1) - pci(2)*s(9, 1, m + 2) & ! [fy2z|o|s]
273 51321893 : + orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
274 : s(19, 1, m + 1) = pai(3)*s(9, 1, m + 1) - pci(3)*s(9, 1, m + 2) & ! [fyz2|o|s]
275 51321893 : + orho*(s(3, 1, m + 1) - s(3, 1, m + 2))
276 : s(20, 1, m + 1) = pai(3)*s(10, 1, m + 1) - pci(3)*s(10, 1, m + 2) & ! [fz3 |o|s]
277 77230814 : + 2._dp*orho*(s(4, 1, m + 1) - s(4, 1, m + 2))
278 : END DO
279 25412972 : ELSE IF (llr == 4) THEN
280 43807676 : DO m = 0, mmax - llr
281 : s(21, 1, m + 1) = pai(1)*s(11, 1, m + 1) - pci(1)*s(11, 1, m + 2) & ! [gx4 |s|s]
282 25412972 : + 3._dp*orho*(s(5, 1, m + 1) - s(5, 1, m + 2))
283 : s(22, 1, m + 1) = pai(1)*s(12, 1, m + 1) - pci(1)*s(12, 1, m + 2) & ! [gx3y |s|s]
284 25412972 : + 2._dp*orho*(s(6, 1, m + 1) - s(6, 1, m + 2))
285 : s(23, 1, m + 1) = pai(1)*s(13, 1, m + 1) - pci(1)*s(13, 1, m + 2) & ! [gx3z |s|s]
286 25412972 : + 2._dp*orho*(s(7, 1, m + 1) - s(7, 1, m + 2))
287 : s(24, 1, m + 1) = pai(1)*s(14, 1, m + 1) - pci(1)*s(14, 1, m + 2) & ! [gx2y2|s|s]
288 25412972 : + orho*(s(8, 1, m + 1) - s(8, 1, m + 2))
289 : s(25, 1, m + 1) = pai(1)*s(15, 1, m + 1) - pci(1)*s(15, 1, m + 2) & ! [gx2yz|s|s]
290 25412972 : + orho*(s(9, 1, m + 1) - s(9, 1, m + 2))
291 : s(26, 1, m + 1) = pai(1)*s(16, 1, m + 1) - pci(1)*s(16, 1, m + 2) & ! [gx2z2|s|s]
292 25412972 : + orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
293 25412972 : s(27, 1, m + 1) = pai(1)*s(17, 1, m + 1) - pci(1)*s(17, 1, m + 2) ! [gxy3 |s|s]
294 25412972 : s(28, 1, m + 1) = pai(1)*s(18, 1, m + 1) - pci(1)*s(18, 1, m + 2) ! [gxy2z|s|s]
295 25412972 : s(29, 1, m + 1) = pai(1)*s(19, 1, m + 1) - pci(1)*s(19, 1, m + 2) ! [gxyz2|s|s]
296 25412972 : s(30, 1, m + 1) = pai(1)*s(20, 1, m + 1) - pci(1)*s(20, 1, m + 2) ! [gxz3 |s|s]
297 : s(31, 1, m + 1) = pai(2)*s(17, 1, m + 1) - pci(2)*s(17, 1, m + 2) & ! [gy4 |s|s]
298 25412972 : + 3._dp*orho*(s(8, 1, m + 1) - s(8, 1, m + 2))
299 : s(32, 1, m + 1) = pai(2)*s(18, 1, m + 1) - pci(2)*s(18, 1, m + 2) & ! [gy3z |s|s]
300 25412972 : + 2._dp*orho*(s(9, 1, m + 1) - s(9, 1, m + 2))
301 : s(33, 1, m + 1) = pai(2)*s(19, 1, m + 1) - pci(2)*s(19, 1, m + 2) & ! [gy2z2|s|s]
302 25412972 : + orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
303 25412972 : s(34, 1, m + 1) = pai(2)*s(20, 1, m + 1) - pci(2)*s(20, 1, m + 2) ! [gyz3 |s|s]
304 : s(35, 1, m + 1) = pai(3)*s(20, 1, m + 1) - pci(3)*s(20, 1, m + 2) & ! [gz4 |s|s]
305 43807676 : + 3._dp*orho*(s(10, 1, m + 1) - s(10, 1, m + 2))
306 : END DO
307 : ELSE
308 51810421 : DO irx = 0, llr
309 218556202 : DO iry = 0, llr - irx
310 166745781 : irz = llr - irx - iry
311 166745781 : irr(1) = irx; irr(2) = iry; irr(3) = irz
312 833728905 : ixx = MAXLOC(irr)
313 166745781 : ix = ixx(1)
314 166745781 : ir = coset(irx, iry, irz)
315 166745781 : irm = irr
316 166745781 : irm(ix) = irm(ix) - 1
317 166745781 : aai = REAL(MAX(irm(ix), 0), dp)*orho
318 166745781 : ir1 = coset(irm(1), irm(2), irm(3))
319 166745781 : irm(ix) = irm(ix) - 1
320 166745781 : ir2 = coset(irm(1), irm(2), irm(3))
321 438853748 : DO m = 0, mmax - llr
322 : s(ir, 1, m + 1) = pai(ix)*s(ir1, 1, m + 1) - pci(ix)*s(ir1, 1, m + 2) &
323 394061595 : + aai*(s(ir2, 1, m + 1) - s(ir2, 1, m + 2))
324 : END DO
325 : END DO
326 : END DO
327 : END IF
328 : END DO
329 :
330 : ! *** Horizontal recurrence steps ***
331 : ! *** [a|c|b+1i] = [a+1i|c|b] + (Ai - Bi)*[a|c|b] ***
332 :
333 125818609 : DO mb = 1, lb_max
334 304621111 : DO ibx = 0, mb
335 569668549 : DO iby = 0, mb - ibx
336 316163952 : ibz = mb - ibx - iby
337 316163952 : ib = coset(ibx, iby, ibz)
338 316163952 : ii(1) = ibx; ii(2) = iby; ii(3) = ibz
339 1580819760 : ixx = MAXLOC(ii)
340 316163952 : ix = ixx(1)
341 316163952 : abx = -rab(ix)
342 316163952 : iim = ii
343 316163952 : iim(ix) = iim(ix) - 1
344 316163952 : ibm = coset(iim(1), iim(2), iim(3))
345 4826264601 : DO ia = 1, ncoset(mmax - mb)
346 4331298147 : iap = iiap(ia, ix)
347 4647462099 : s(ia, ib, 1) = s(iap, ibm, 1) + abx*s(ia, ibm, 1)
348 : END DO
349 : END DO
350 : END DO
351 : END DO
352 :
353 7205637 : ELSE IF (lb_max > 0) THEN
354 :
355 : ! *** Recurrence steps: [s|c|s] -> [s|c|b] ***
356 : ! *** [s|c|b](m) = (Pi - Bi)*[s|c|b-1i](m) - ***
357 : ! *** (Pi - Ci)*[s|c|b-1i](m+1)) + ***
358 : ! *** Ni(b-1i)/2(a+b)*[s|c|b-2i](m) - ***
359 : ! *** Ni(b-1i)/2(a+b)*[s|c|b-2i](m+1) ***
360 4479677 : DO llr = 1, lb_max
361 4479677 : IF (llr == 1) THEN
362 4479677 : DO m = 0, lb_max - llr
363 2354003 : s1 = s(1, 1, m + 1)
364 2354003 : s2 = s(1, 1, m + 2)
365 2354003 : s(1, 2, m + 1) = pbi(1)*s1 - pci(1)*s2 ! [px|o|s]
366 2354003 : s(1, 3, m + 1) = pbi(2)*s1 - pci(2)*s2 ! [py|o|s]
367 4479677 : s(1, 4, m + 1) = pbi(3)*s1 - pci(3)*s2 ! [pz|o|s]
368 : END DO
369 228329 : ELSE IF (llr == 2) THEN
370 445016 : DO m = 0, lb_max - llr
371 228329 : s1 = s(1, 1, m + 1) - s(1, 1, m + 2)
372 228329 : s(1, 5, m + 1) = pbi(1)*s(1, 2, m + 1) - pci(1)*s(1, 2, m + 2) + orho*s1 ! [dx2|o|s]
373 228329 : s(1, 6, m + 1) = pbi(1)*s(1, 3, m + 1) - pci(1)*s(1, 3, m + 2) ! [dxy|o|s]
374 228329 : s(1, 7, m + 1) = pbi(1)*s(1, 4, m + 1) - pci(1)*s(1, 4, m + 2) ! [dxz|o|s]
375 228329 : s(1, 8, m + 1) = pbi(2)*s(1, 3, m + 1) - pci(2)*s(1, 3, m + 2) + orho*s1 ! [dy2|o|s]
376 228329 : s(1, 9, m + 1) = pbi(2)*s(1, 4, m + 1) - pci(2)*s(1, 4, m + 2) ! [dyz|o|s]
377 445016 : s(1, 10, m + 1) = pbi(3)*s(1, 4, m + 1) - pci(3)*s(1, 4, m + 2) + orho*s1 ! [dz2|o|s]
378 : END DO
379 11642 : ELSE IF (llr == 3) THEN
380 23202 : DO m = 0, lb_max - llr
381 : s(1, 11, m + 1) = pbi(1)*s(1, 5, m + 1) - pci(1)*s(1, 5, m + 2) & ! [fx3 |o|s]
382 11642 : + 2._dp*orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
383 : s(1, 12, m + 1) = pbi(1)*s(1, 6, m + 1) - pci(1)*s(1, 6, m + 2) & ! [fx2y|o|s]
384 11642 : + orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
385 : s(1, 13, m + 1) = pbi(1)*s(1, 7, m + 1) - pci(1)*s(1, 7, m + 2) & ! [fx2z|o|s]
386 11642 : + orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
387 : s(1, 14, m + 1) = pbi(2)*s(1, 6, m + 1) - pci(2)*s(1, 6, m + 2) & ! [fxy2|o|s]
388 11642 : + orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
389 11642 : s(1, 15, m + 1) = pbi(1)*s(1, 9, m + 1) - pci(1)*s(1, 9, m + 2) ! [fxyz|o|s]
390 : s(1, 16, m + 1) = pbi(3)*s(1, 7, m + 1) - pci(3)*s(1, 7, m + 2) & ! [fxz2|o|s]
391 11642 : + orho*(s(1, 2, m + 1) - s(1, 2, m + 2))
392 : s(1, 17, m + 1) = pbi(2)*s(1, 8, m + 1) - pci(2)*s(1, 8, m + 2) & ! [fy3 |o|s]
393 11642 : + 2._dp*orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
394 : s(1, 18, m + 1) = pbi(2)*s(1, 9, m + 1) - pci(2)*s(1, 9, m + 2) & ! [fy2z|o|s]
395 11642 : + orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
396 : s(1, 19, m + 1) = pbi(3)*s(1, 9, m + 1) - pci(3)*s(1, 9, m + 2) & ! [fyz2|o|s]
397 11642 : + orho*(s(1, 3, m + 1) - s(1, 3, m + 2))
398 : s(1, 20, m + 1) = pbi(3)*s(1, 10, m + 1) - pci(3)*s(1, 10, m + 2) & ! [fz3 |o|s]
399 23202 : + 2._dp*orho*(s(1, 4, m + 1) - s(1, 4, m + 2))
400 : END DO
401 : ELSE
402 492 : DO irx = 0, llr
403 1722 : DO iry = 0, llr - irx
404 1230 : irz = llr - irx - iry
405 1230 : irr(1) = irx; irr(2) = iry; irr(3) = irz
406 6150 : ixx = MAXLOC(irr)
407 1230 : ix = ixx(1)
408 1230 : ir = coset(irx, iry, irz)
409 1230 : irm = irr
410 1230 : irm(ix) = irm(ix) - 1
411 1230 : aai = REAL(MAX(irm(ix), 0), dp)
412 1230 : ir1 = coset(irm(1), irm(2), irm(3))
413 1230 : irm(ix) = irm(ix) - 1
414 1230 : ir2 = coset(irm(1), irm(2), irm(3))
415 2870 : DO m = 0, lb_max - llr
416 : s(1, ir, m + 1) = pbi(ix)*s(1, ir1, m + 1) - pci(ix)*s(1, ir1, m + 2) &
417 2460 : + aai*orho*(s(1, ir2, m + 1) - s(1, ir2, m + 2))
418 : END DO
419 : END DO
420 : END DO
421 : END IF
422 : END DO
423 :
424 : END IF
425 :
426 : ! *** Store the primitive three-center overlap integrals ***
427 311192419 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
428 1722266346 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
429 1663944195 : vab(na + i, nb + j) = vab(na + i, nb + j) + s(i, j, 1)
430 : END DO
431 : END DO
432 :
433 : ! *** Calculate the requested derivatives with respect ***
434 : ! *** to the nuclear coordinates of the atomic center a ***
435 :
436 74539500 : DO da = 0, da_max - 1
437 16217349 : ftz = 2.0_dp*zeta(ipgf)
438 90756849 : DO dax = 0, da
439 48652047 : DO day = 0, da - dax
440 16217349 : daz = da - dax - day
441 16217349 : cda = coset(dax, day, daz)
442 16217349 : cdax = coset(dax + 1, day, daz)
443 16217349 : cday = coset(dax, day + 1, daz)
444 16217349 : cdaz = coset(dax, day, daz + 1)
445 :
446 : ! *** [da/dAi|c|b] = 2*zeta*[a+1i|c|b] - Ni(a)[a-1i|c|b] ***
447 :
448 63205888 : DO la = la_min_set, la_max - da - 1
449 97594279 : DO ax = 0, la
450 50605740 : fax = REAL(ax, dp)
451 155473798 : DO ay = 0, la - ax
452 74096868 : fay = REAL(ay, dp)
453 74096868 : az = la - ax - ay
454 74096868 : faz = REAL(az, dp)
455 74096868 : coa = coset(ax, ay, az)
456 74096868 : coamx = coset(ax - 1, ay, az)
457 74096868 : coamy = coset(ax, ay - 1, az)
458 74096868 : coamz = coset(ax, ay, az - 1)
459 74096868 : coapx = coset(ax + 1, ay, az)
460 74096868 : coapy = coset(ax, ay + 1, az)
461 74096868 : coapz = coset(ax, ay, az + 1)
462 557702069 : DO cob = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
463 432999461 : fs(coa, cob, cdax) = ftz*s(coapx, cob, cda) - fax*s(coamx, cob, cda)
464 432999461 : fs(coa, cob, cday) = ftz*s(coapy, cob, cda) - fay*s(coamy, cob, cda)
465 507096329 : fs(coa, cob, cdaz) = ftz*s(coapz, cob, cda) - faz*s(coamz, cob, cda)
466 : END DO
467 : END DO
468 : END DO
469 : END DO
470 : END DO
471 :
472 : END DO
473 : END DO
474 :
475 : ! DFPT for APTs
476 58322151 : IF (PRESENT(vab2_work)) THEN
477 28512 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
478 72090 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
479 43578 : vab2_work(na + i, nb + j, 1) = vab2_work(na + i, nb + j, 1) + fs(i, j, 2)
480 43578 : vab2_work(na + i, nb + j, 2) = vab2_work(na + i, nb + j, 2) + fs(i, j, 3)
481 62676 : vab2_work(na + i, nb + j, 3) = vab2_work(na + i, nb + j, 3) + fs(i, j, 4)
482 : END DO
483 : END DO
484 : END IF
485 :
486 : ! *** Calculate the force contribution for the atomic center a ***
487 :
488 58322151 : IF (calculate_force_a) THEN
489 90437902 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
490 523393785 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
491 432955883 : force_a(1) = force_a(1) + pab(na + i, nb + j)*fs(i, j, 2)
492 432955883 : force_a(2) = force_a(2) + pab(na + i, nb + j)*fs(i, j, 3)
493 507185850 : force_a(3) = force_a(3) + pab(na + i, nb + j)*fs(i, j, 4)
494 : END DO
495 : END DO
496 : END IF
497 :
498 : ! *** Calculate the requested derivatives with respect ***
499 : ! *** to the nuclear coordinates of the atomic center b ***
500 :
501 74539500 : DO db = 0, db_max - 1
502 16217349 : ftz = 2.0_dp*zetb(jpgf)
503 90756849 : DO dbx = 0, db
504 48652047 : DO dby = 0, db - dbx
505 16217349 : dbz = db - dbx - dby
506 16217349 : cdb = coset(dbx, dby, dbz)
507 16217349 : cdbx = coset(dbx + 1, dby, dbz)
508 16217349 : cdby = coset(dbx, dby + 1, dbz)
509 16217349 : cdbz = coset(dbx, dby, dbz + 1)
510 :
511 : ! *** [a|c|db/dBi] = 2*zetb*[a|c|b+1i] - Ni(b)[a|c|b-1i] ***
512 :
513 63254184 : DO lb = lb_min_set, lb_max - db - 1
514 97741656 : DO bx = 0, lb
515 50704821 : fbx = REAL(bx, dp)
516 155773372 : DO by = 0, lb - bx
517 74249065 : fby = REAL(by, dp)
518 74249065 : bz = lb - bx - by
519 74249065 : fbz = REAL(bz, dp)
520 74249065 : cob = coset(bx, by, bz)
521 74249065 : cobmx = coset(bx - 1, by, bz)
522 74249065 : cobmy = coset(bx, by - 1, bz)
523 74249065 : cobmz = coset(bx, by, bz - 1)
524 74249065 : cobpx = coset(bx + 1, by, bz)
525 74249065 : cobpy = coset(bx, by + 1, bz)
526 74249065 : cobpz = coset(bx, by, bz + 1)
527 557953347 : DO coa = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
528 432999461 : fs(coa, cob, cdbx) = ftz*s(coa, cobpx, cdb) - fbx*s(coa, cobmx, cdb)
529 432999461 : fs(coa, cob, cdby) = ftz*s(coa, cobpy, cdb) - fby*s(coa, cobmy, cdb)
530 507248526 : fs(coa, cob, cdbz) = ftz*s(coa, cobpz, cdb) - fbz*s(coa, cobmz, cdb)
531 : END DO
532 : END DO
533 : END DO
534 : END DO
535 :
536 : END DO
537 : END DO
538 : END DO
539 :
540 : ! DFPT for APTs
541 58322151 : IF (PRESENT(vab2_work)) THEN
542 28512 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
543 72090 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
544 43578 : vab2_work(na + i, nb + j, 4) = vab2_work(na + i, nb + j, 4) + fs(i, j, 2)
545 43578 : vab2_work(na + i, nb + j, 5) = vab2_work(na + i, nb + j, 5) + fs(i, j, 3)
546 62676 : vab2_work(na + i, nb + j, 6) = vab2_work(na + i, nb + j, 6) + fs(i, j, 4)
547 : END DO
548 : END DO
549 :
550 9414 : CPASSERT(PRESENT(iatom) .AND. PRESENT(jatom) .AND. PRESENT(katom))
551 9414 : CPASSERT(PRESENT(deltaR))
552 37656 : DO idir = 1, 3
553 94950 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
554 216270 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
555 : vab2(na + i, nb + j, idir) = vab2(na + i, nb + j, idir) &
556 : + vab2_work(na + i, nb + j, idir)*deltaR(idir, iatom) &
557 : - vab2_work(na + i, nb + j, idir)*deltaR(idir, katom) &
558 : + vab2_work(na + i, nb + j, idir + 3)*deltaR(idir, jatom) &
559 188028 : - vab2_work(na + i, nb + j, idir + 3)*deltaR(idir, katom)
560 : END DO
561 : END DO
562 : END DO
563 : END IF
564 :
565 : ! *** Calculate the force contribution for the atomic center b ***
566 :
567 58322151 : IF (calculate_force_b) THEN
568 90437902 : DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set)
569 523393785 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
570 432955883 : force_b(1) = force_b(1) + pab(na + i, nb + j)*fs(i, j, 2)
571 432955883 : force_b(2) = force_b(2) + pab(na + i, nb + j)*fs(i, j, 3)
572 507185850 : force_b(3) = force_b(3) + pab(na + i, nb + j)*fs(i, j, 4)
573 : END DO
574 : END DO
575 : END IF
576 :
577 103040259 : nb = nb + ncoset(lb_max_set)
578 :
579 : END DO
580 :
581 79663256 : na = na + ncoset(la_max_set)
582 :
583 : END DO
584 :
585 34945148 : DEALLOCATE (iiap)
586 :
587 34945148 : END SUBROUTINE os_3center
588 : ! **************************************************************************************************
589 : !> \brief Calculation of two-center integrals <a|c> over
590 : !> Cartesian Gaussian functions and a spherical potential
591 : !>
592 : !> \param la_max_set ...
593 : !> \param la_min_set ...
594 : !> \param npgfa ...
595 : !> \param rpgfa ...
596 : !> \param zeta ...
597 : !> \param auxint ...
598 : !> \param rpgfc ...
599 : !> \param rac ...
600 : !> \param dac ...
601 : !> \param va ...
602 : !> \param dva ...
603 : !> \date December 2017
604 : !> \author Juerg Hutter
605 : !> \version 1.0
606 : ! **************************************************************************************************
607 328 : SUBROUTINE os_2center(la_max_set, la_min_set, npgfa, rpgfa, zeta, &
608 656 : auxint, rpgfc, rac, dac, va, dva)
609 : INTEGER, INTENT(IN) :: la_max_set, la_min_set, npgfa
610 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
611 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: auxint
612 : REAL(KIND=dp), INTENT(IN) :: rpgfc
613 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac
614 : REAL(KIND=dp), INTENT(IN) :: dac
615 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: va
616 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
617 : OPTIONAL :: dva
618 :
619 : INTEGER :: ax, ay, az, coa, coamx, coamy, coamz, coapx, coapy, coapz, da_max, i, ipgf, ir, &
620 : ir1, ir2, irm(3), irr(3), irx, iry, irz, ix, ixx(1), la, la_max, la_min, llr, m, mmax, na
621 : REAL(KIND=dp) :: aai, fax, fay, faz, ftz, orho, s1
622 328 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: s
623 :
624 328 : IF (PRESENT(dva)) THEN
625 : da_max = 1
626 : ELSE
627 164 : da_max = 0
628 : END IF
629 :
630 328 : la_max = la_max_set + da_max
631 328 : la_min = MAX(0, la_min_set - da_max)
632 :
633 328 : mmax = la_max
634 :
635 1312 : ALLOCATE (s(ncoset(mmax), mmax + 1))
636 328 : na = 0
637 656 : DO ipgf = 1, npgfa
638 328 : IF (rpgfa(ipgf) + rpgfc < dac) THEN
639 0 : na = na + ncoset(la_max_set)
640 0 : CYCLE
641 : END IF
642 1260 : s(1, 1:mmax + 1) = auxint(0:mmax, ipgf)
643 328 : IF (la_max > 0) THEN
644 : ! Recurrence steps: [s|c] -> [a|c]
645 : ! [a|c](m) = (Ci - Ai)*[a-1i|c](m+1) +
646 : ! Ni(a-1i)/2a*[a-2i|c](m) -
647 : ! Ni(a-1i)/2a*[a-2i|c](m+1)
648 : !
649 :
650 268 : orho = 0.5_dp/zeta(ipgf)
651 :
652 872 : DO llr = 1, mmax
653 872 : IF (llr == 1) THEN
654 872 : DO m = 0, mmax - llr
655 604 : s1 = s(1, m + 2)
656 604 : s(2, m + 1) = -rac(1)*s1 ! [px|o]
657 604 : s(3, m + 1) = -rac(2)*s1 ! [py|o]
658 872 : s(4, m + 1) = -rac(3)*s1 ! [pz|o]
659 : END DO
660 336 : ELSE IF (llr == 2) THEN
661 544 : DO m = 0, mmax - llr
662 336 : s1 = s(1, m + 1) - s(1, m + 2)
663 336 : s(5, m + 1) = -rac(1)*s(2, m + 2) + orho*s1 ! [dx2|o]
664 336 : s(6, m + 1) = -rac(1)*s(3, m + 2) ! [dxy|o]
665 336 : s(7, m + 1) = -rac(1)*s(4, m + 2) ! [dxz|o]
666 336 : s(8, m + 1) = -rac(2)*s(3, m + 2) + orho*s1 ! [dy2|o]
667 336 : s(9, m + 1) = -rac(2)*s(4, m + 2) ! [dyz|o]
668 544 : s(10, m + 1) = -rac(3)*s(4, m + 2) + orho*s1 ! [dz2|o]
669 : END DO
670 128 : ELSE IF (llr == 3) THEN
671 244 : DO m = 0, mmax - llr
672 128 : s(11, m + 1) = -rac(1)*s(5, m + 2) + 2._dp*orho*(s(2, m + 1) - s(2, m + 2)) ! [fx3 |o]
673 128 : s(12, m + 1) = -rac(1)*s(6, m + 2) + orho*(s(3, m + 1) - s(3, m + 2)) ! [fx2y|o]
674 128 : s(13, m + 1) = -rac(1)*s(7, m + 2) + orho*(s(4, m + 1) - s(4, m + 2)) ! [fx2z|o]
675 128 : s(14, m + 1) = -rac(2)*s(6, m + 2) + orho*(s(2, m + 1) - s(2, m + 2)) ! [fxy2|o]
676 128 : s(15, m + 1) = -rac(1)*s(9, m + 2) ! [fxyz|o]
677 128 : s(16, m + 1) = -rac(3)*s(7, m + 2) + orho*(s(2, m + 1) - s(2, m + 2)) ! [fxz2|o]
678 128 : s(17, m + 1) = -rac(2)*s(8, m + 2) + 2._dp*orho*(s(3, m + 1) - s(3, m + 2)) ! [fy3 |o]
679 128 : s(18, m + 1) = -rac(2)*s(9, m + 2) + orho*(s(4, m + 1) - s(4, m + 2)) ! [fy2z|o]
680 128 : s(19, m + 1) = -rac(3)*s(9, m + 2) + orho*(s(3, m + 1) - s(3, m + 2)) ! [fyz2|o]
681 244 : s(20, m + 1) = -rac(3)*s(10, m + 2) + 2._dp*orho*(s(4, m + 1) - s(4, m + 2)) ! [fz3 |o]
682 : END DO
683 12 : ELSE IF (llr == 4) THEN
684 24 : DO m = 0, mmax - llr
685 12 : s(21, m + 1) = -rac(1)*s(11, m + 2) + 3._dp*orho*(s(5, m + 1) - s(5, m + 2)) ! [gx4 |s]
686 12 : s(22, m + 1) = -rac(1)*s(12, m + 2) + 2._dp*orho*(s(6, m + 1) - s(6, m + 2)) ! [gx3y |s]
687 12 : s(23, m + 1) = -rac(1)*s(13, m + 2) + 2._dp*orho*(s(7, m + 1) - s(7, m + 2)) ! [gx3z |s]
688 12 : s(24, m + 1) = -rac(1)*s(14, m + 2) + orho*(s(8, m + 1) - s(8, m + 2)) ! [gx2y2|s]
689 12 : s(25, m + 1) = -rac(1)*s(15, m + 2) + orho*(s(9, m + 1) - s(9, m + 2)) ! [gx2yz|s]
690 12 : s(26, m + 1) = -rac(1)*s(16, m + 2) + orho*(s(10, m + 1) - s(10, m + 2)) ! [gx2z2|s]
691 12 : s(27, m + 1) = -rac(1)*s(17, m + 2) ! [gxy3 |s]
692 12 : s(28, m + 1) = -rac(1)*s(18, m + 2) ! [gxy2z|s]
693 12 : s(29, m + 1) = -rac(1)*s(19, m + 2) ! [gxyz2|s]
694 12 : s(30, m + 1) = -rac(1)*s(20, m + 2) ! [gxz3 |s]
695 12 : s(31, m + 1) = -rac(2)*s(17, m + 2) + 3._dp*orho*(s(8, m + 1) - s(8, m + 2)) ! [gy4 |s]
696 12 : s(32, m + 1) = -rac(2)*s(18, m + 2) + 2._dp*orho*(s(9, m + 1) - s(9, m + 2)) ! [gy3z |s]
697 12 : s(33, m + 1) = -rac(2)*s(19, m + 2) + orho*(s(10, m + 1) - s(10, m + 2)) ! [gy2z2|s]
698 12 : s(34, m + 1) = -rac(2)*s(20, m + 2) ! [gyz3 |s]
699 24 : s(35, m + 1) = -rac(3)*s(20, m + 2) + 3._dp*orho*(s(10, m + 1) - s(10, m + 2)) ! [gz4 |s]
700 : END DO
701 : ELSE
702 0 : DO irx = 0, llr
703 0 : DO iry = 0, llr - irx
704 0 : irz = llr - irx - iry
705 0 : irr(1) = irx; irr(2) = iry; irr(3) = irz
706 0 : ixx = MAXLOC(irr)
707 0 : ix = ixx(1)
708 0 : ir = coset(irx, iry, irz)
709 0 : irm = irr
710 0 : irm(ix) = irm(ix) - 1
711 0 : aai = REAL(MAX(irm(ix), 0), dp)*orho
712 0 : ir1 = coset(irm(1), irm(2), irm(3))
713 0 : irm(ix) = irm(ix) - 1
714 0 : ir2 = coset(irm(1), irm(2), irm(3))
715 0 : DO m = 0, mmax - llr
716 0 : s(ir, m + 1) = -rac(ix)*s(ir1, m + 2) + aai*(s(ir2, m + 1) - s(ir2, m + 2))
717 : END DO
718 : END DO
719 : END DO
720 : END IF
721 : END DO
722 :
723 : END IF
724 :
725 : ! Store the primitive three-center overlap integrals
726 2404 : DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set)
727 2404 : va(na + i) = va(na + i) + s(i, 1)
728 : END DO
729 :
730 : ! Calculate the requested derivatives with respect ***
731 : ! to the nuclear coordinates of the atomic center a ***
732 : ! [da/dAi|c] = 2*zeta*[a+1i|c] - Ni(a)[a-1i|c] ***
733 328 : IF (PRESENT(dva)) THEN
734 164 : ftz = 2.0_dp*zeta(ipgf)
735 450 : DO la = la_min_set, la_max_set
736 1048 : DO ax = 0, la
737 598 : fax = REAL(ax, dp)
738 1922 : DO ay = 0, la - ax
739 1038 : fay = REAL(ay, dp)
740 1038 : az = la - ax - ay
741 1038 : faz = REAL(az, dp)
742 1038 : coa = coset(ax, ay, az)
743 1038 : coamx = coset(ax - 1, ay, az)
744 1038 : coamy = coset(ax, ay - 1, az)
745 1038 : coamz = coset(ax, ay, az - 1)
746 1038 : coapx = coset(ax + 1, ay, az)
747 1038 : coapy = coset(ax, ay + 1, az)
748 1038 : coapz = coset(ax, ay, az + 1)
749 1038 : dva(na + coa, 1) = dva(na + coa, 1) + ftz*s(coapx, 1) - fax*s(coamx, 1)
750 1038 : dva(na + coa, 2) = dva(na + coa, 2) + ftz*s(coapy, 1) - fay*s(coamy, 1)
751 1636 : dva(na + coa, 3) = dva(na + coa, 3) + ftz*s(coapz, 1) - faz*s(coamz, 1)
752 : END DO
753 : END DO
754 : END DO
755 : END IF
756 :
757 656 : na = na + ncoset(la_max_set)
758 :
759 : END DO
760 :
761 328 : DEALLOCATE (s)
762 :
763 328 : END SUBROUTINE os_2center
764 : ! **************************************************************************************************
765 :
766 : END MODULE ai_oneelectron
|
