MPQC 3.0.0-alpha
Loading...
Searching...
No Matches
obosrr.timpl.h
1//
2// obosrr.timpl.h
3//
4// Copyright (C) 2001 Edward Valeev
5//
6// Author: Edward Valeev <evaleev@vt.edu>
7// Maintainer: EV
8//
9// This file is part of the SC Toolkit.
10//
11// The SC Toolkit is free software; you can redistribute it and/or modify
12// it under the terms of the GNU Library General Public License as published by
13// the Free Software Foundation; either version 2, or (at your option)
14// any later version.
15//
16// The SC Toolkit is distributed in the hope that it will be useful,
17// but WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19// GNU Library General Public License for more details.
20//
21// You should have received a copy of the GNU Library General Public License
22// along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24//
25// The U.S. Government is granted a limited license as per AL 91-7.
26//
27
28/* Recurrence relation are from the Obara-Saika paper - pp. 3971-3972 */
29
30template <int Order>
31void
32Int1eLibint2::AI_OSrecurs_(double PA[3], double PB[3],
33 double PC[3], double gamma, int iang, int jang)
34{
35 int a,b,m;
36 int izm = 1;
37 int iym = iang + 1;
38 int ixm = iym * iym;
39 int jzm = 1;
40 int jym = jang + 1;
41 int jxm = jym * jym;
42 int ix,iy,iz,jx,jy,jz;
43 int iind,jind;
44 double pp = 1/(2*gamma);
45 int mmax = iang+jang + Order;
46 double tmp = sqrt(gamma)*M_2_SQRTPI;
47 double u = gamma*(PC[0]*PC[0] + PC[1]*PC[1] + PC[2]*PC[2]);
48 Fm_Eval_->eval(Fm_table_,u,mmax);
49
50 /* Computing starting integrals for recursion */
51
52 for(m=0;m<=mmax;m++)
53 AI0_[0][0][m] = tmp * Fm_table_[m];
54 if (Order >= 1)
55 for(m=0;m<=mmax-1;m++) {
56 AIX_[0][0][m] = 2*gamma*PC[0]*AI0_[0][0][m+1];
57 AIY_[0][0][m] = 2*gamma*PC[1]*AI0_[0][0][m+1];
58 AIZ_[0][0][m] = 2*gamma*PC[2]*AI0_[0][0][m+1];
59 }
60 if (Order >= 2)
61 for (m = 0; m <= mmax - 2; m++) {
62 AIXX_[0][0][m] = 4 * gamma * gamma * PC[0] * PC[0] * AI0_[0][0][m + 2]
63 - 2 * gamma * AI0_[0][0][m + 1];
64 AIYY_[0][0][m] = 4 * gamma * gamma * PC[1] * PC[1] * AI0_[0][0][m + 2]
65 - 2 * gamma * AI0_[0][0][m + 1];
66 AIZZ_[0][0][m] = 4 * gamma * gamma * PC[2] * PC[2] * AI0_[0][0][m + 2]
67 - 2 * gamma * AI0_[0][0][m + 1];
68 AIXY_[0][0][m] = 4 * gamma * gamma * PC[0] * PC[1] * AI0_[0][0][m + 2];
69 AIXZ_[0][0][m] = 4 * gamma * gamma * PC[0] * PC[2] * AI0_[0][0][m + 2];
70 AIYZ_[0][0][m] = 4 * gamma * gamma * PC[1] * PC[2] * AI0_[0][0][m + 2];
71 }
72
73 /* Upward recursion in j with i=0 */
74
75 for(b=1;b<=jang;b++)
76 for(jx=0;jx<=b;jx++)
77 for(jy=0;jy<=b-jx;jy++) {
78 jz = b-jx-jy;
79 jind = jx*jxm+jy*jym+jz*jzm;
80 if (jz > 0) {
81 for(m=0;m<=mmax-b;m++) /* Electrostatic potential integrals */
82 AI0_[0][jind][m] = PB[2]*AI0_[0][jind-jzm][m] -
83 PC[2]*AI0_[0][jind-jzm][m+1];
84 if (Order >= 1)
85 for(m=0;m<=mmax-b-1;m++) { /* Electric field integrals */
86 AIX_[0][jind][m] = PB[2]*AIX_[0][jind-jzm][m] -
87 PC[2]*AIX_[0][jind-jzm][m+1];
88 AIY_[0][jind][m] = PB[2]*AIY_[0][jind-jzm][m] -
89 PC[2]*AIY_[0][jind-jzm][m+1];
90 AIZ_[0][jind][m] = PB[2]*AIZ_[0][jind-jzm][m] -
91 PC[2]*AIZ_[0][jind-jzm][m+1] +
92 AI0_[0][jind-jzm][m+1];
93 }
94 if (Order >= 2)
95 for(m=0;m<=mmax-b-2;m++) { /* Gradients of the electric field */
96 AIXX_[0][jind][m] = PB[2]*AIXX_[0][jind-jzm][m] -
97 PC[2]*AIXX_[0][jind-jzm][m+1];
98 AIYY_[0][jind][m] = PB[2]*AIYY_[0][jind-jzm][m] -
99 PC[2]*AIYY_[0][jind-jzm][m+1];
100 AIZZ_[0][jind][m] = PB[2]*AIZZ_[0][jind-jzm][m] -
101 PC[2]*AIZZ_[0][jind-jzm][m+1] +
102 2*AIZ_[0][jind-jzm][m+1];
103 AIXY_[0][jind][m] = PB[2]*AIXY_[0][jind-jzm][m] -
104 PC[2]*AIXY_[0][jind-jzm][m+1];
105 AIXZ_[0][jind][m] = PB[2]*AIXZ_[0][jind-jzm][m] -
106 PC[2]*AIXZ_[0][jind-jzm][m+1] +
107 AIX_[0][jind-jzm][m+1];
108 AIYZ_[0][jind][m] = PB[2]*AIYZ_[0][jind-jzm][m] -
109 PC[2]*AIYZ_[0][jind-jzm][m+1] +
110 AIY_[0][jind-jzm][m+1];
111 }
112 if (jz > 1) {
113 for(m=0;m<=mmax-b;m++)
114 AI0_[0][jind][m] += pp*(jz-1)*(AI0_[0][jind-2*jzm][m] -
115 AI0_[0][jind-2*jzm][m+1]);
116 if (Order >= 1)
117 for(m=0;m<=mmax-b-1;m++) {
118 AIX_[0][jind][m] += pp*(jz-1)*(AIX_[0][jind-2*jzm][m] -
119 AIX_[0][jind-2*jzm][m+1]);
120 AIY_[0][jind][m] += pp*(jz-1)*(AIY_[0][jind-2*jzm][m] -
121 AIY_[0][jind-2*jzm][m+1]);
122 AIZ_[0][jind][m] += pp*(jz-1)*(AIZ_[0][jind-2*jzm][m] -
123 AIZ_[0][jind-2*jzm][m+1]);
124 }
125 if (Order >= 2)
126 for(m=0;m<=mmax-b-2;m++) {
127 AIXX_[0][jind][m] += pp*(jz-1)*(AIXX_[0][jind-2*jzm][m] -
128 AIXX_[0][jind-2*jzm][m+1]);
129 AIYY_[0][jind][m] += pp*(jz-1)*(AIYY_[0][jind-2*jzm][m] -
130 AIYY_[0][jind-2*jzm][m+1]);
131 AIZZ_[0][jind][m] += pp*(jz-1)*(AIZZ_[0][jind-2*jzm][m] -
132 AIZZ_[0][jind-2*jzm][m+1]);
133 AIXY_[0][jind][m] += pp*(jz-1)*(AIXY_[0][jind-2*jzm][m] -
134 AIXY_[0][jind-2*jzm][m+1]);
135 AIXZ_[0][jind][m] += pp*(jz-1)*(AIXZ_[0][jind-2*jzm][m] -
136 AIXZ_[0][jind-2*jzm][m+1]);
137 AIYZ_[0][jind][m] += pp*(jz-1)*(AIYZ_[0][jind-2*jzm][m] -
138 AIYZ_[0][jind-2*jzm][m+1]);
139 }
140 }
141 }
142 else
143 if (jy > 0) {
144 for(m=0;m<=mmax-b;m++)
145 AI0_[0][jind][m] = PB[1]*AI0_[0][jind-jym][m] -
146 PC[1]*AI0_[0][jind-jym][m+1];
147 if (Order >= 1)
148 for(m=0;m<=mmax-b-1;m++) {
149 AIX_[0][jind][m] = PB[1]*AIX_[0][jind-jym][m] -
150 PC[1]*AIX_[0][jind-jym][m+1];
151 AIY_[0][jind][m] = PB[1]*AIY_[0][jind-jym][m] -
152 PC[1]*AIY_[0][jind-jym][m+1] +
153 AI0_[0][jind-jym][m+1];
154 AIZ_[0][jind][m] = PB[1]*AIZ_[0][jind-jym][m] -
155 PC[1]*AIZ_[0][jind-jym][m+1];
156 }
157 if (Order >= 2)
158 for(m=0;m<=mmax-b-2;m++) {
159 AIXX_[0][jind][m] = PB[1]*AIXX_[0][jind-jym][m] -
160 PC[1]*AIXX_[0][jind-jym][m+1];
161 AIYY_[0][jind][m] = PB[1]*AIYY_[0][jind-jym][m] -
162 PC[1]*AIYY_[0][jind-jym][m+1] +
163 2*AIY_[0][jind-jym][m+1];
164 AIZZ_[0][jind][m] = PB[1]*AIZZ_[0][jind-jym][m] -
165 PC[1]*AIZZ_[0][jind-jym][m+1];
166 AIXY_[0][jind][m] = PB[1]*AIXY_[0][jind-jym][m] -
167 PC[1]*AIXY_[0][jind-jym][m+1] +
168 AIX_[0][jind-jym][m+1];
169 AIXZ_[0][jind][m] = PB[1]*AIXZ_[0][jind-jym][m] -
170 PC[1]*AIXZ_[0][jind-jym][m+1];
171 AIYZ_[0][jind][m] = PB[1]*AIYZ_[0][jind-jym][m] -
172 PC[1]*AIYZ_[0][jind-jym][m+1] +
173 AIZ_[0][jind-jym][m+1];
174 }
175 if (jy > 1) {
176 for(m=0;m<=mmax-b;m++)
177 AI0_[0][jind][m] += pp*(jy-1)*(AI0_[0][jind-2*jym][m] -
178 AI0_[0][jind-2*jym][m+1]);
179 if (Order >= 1)
180 for(m=0;m<=mmax-b-1;m++) {
181 AIX_[0][jind][m] += pp*(jy-1)*(AIX_[0][jind-2*jym][m] -
182 AIX_[0][jind-2*jym][m+1]);
183 AIY_[0][jind][m] += pp*(jy-1)*(AIY_[0][jind-2*jym][m] -
184 AIY_[0][jind-2*jym][m+1]);
185 AIZ_[0][jind][m] += pp*(jy-1)*(AIZ_[0][jind-2*jym][m] -
186 AIZ_[0][jind-2*jym][m+1]);
187 }
188 if (Order >= 2)
189 for(m=0;m<=mmax-b-2;m++) {
190 AIXX_[0][jind][m] += pp*(jy-1)*(AIXX_[0][jind-2*jym][m] -
191 AIXX_[0][jind-2*jym][m+1]);
192 AIYY_[0][jind][m] += pp*(jy-1)*(AIYY_[0][jind-2*jym][m] -
193 AIYY_[0][jind-2*jym][m+1]);
194 AIZZ_[0][jind][m] += pp*(jy-1)*(AIZZ_[0][jind-2*jym][m] -
195 AIZZ_[0][jind-2*jym][m+1]);
196 AIXY_[0][jind][m] += pp*(jy-1)*(AIXY_[0][jind-2*jym][m] -
197 AIXY_[0][jind-2*jym][m+1]);
198 AIXZ_[0][jind][m] += pp*(jy-1)*(AIXZ_[0][jind-2*jym][m] -
199 AIXZ_[0][jind-2*jym][m+1]);
200 AIYZ_[0][jind][m] += pp*(jy-1)*(AIYZ_[0][jind-2*jym][m] -
201 AIYZ_[0][jind-2*jym][m+1]);
202 }
203 }
204 }
205 else
206 if (jx > 0) {
207 for(m=0;m<=mmax-b;m++)
208 AI0_[0][jind][m] = PB[0]*AI0_[0][jind-jxm][m] -
209 PC[0]*AI0_[0][jind-jxm][m+1];
210 if (Order >= 1)
211 for(m=0;m<=mmax-b-1;m++) {
212 AIX_[0][jind][m] = PB[0]*AIX_[0][jind-jxm][m] -
213 PC[0]*AIX_[0][jind-jxm][m+1] +
214 AI0_[0][jind-jxm][m+1];
215 AIY_[0][jind][m] = PB[0]*AIY_[0][jind-jxm][m] -
216 PC[0]*AIY_[0][jind-jxm][m+1];
217 AIZ_[0][jind][m] = PB[0]*AIZ_[0][jind-jxm][m] -
218 PC[0]*AIZ_[0][jind-jxm][m+1];
219 }
220 if (Order >= 2)
221 for(m=0;m<=mmax-b-2;m++) {
222 AIXX_[0][jind][m] = PB[0]*AIXX_[0][jind-jxm][m] -
223 PC[0]*AIXX_[0][jind-jxm][m+1] +
224 2*AIX_[0][jind-jxm][m+1];
225 AIYY_[0][jind][m] = PB[0]*AIYY_[0][jind-jxm][m] -
226 PC[0]*AIYY_[0][jind-jxm][m+1];
227 AIZZ_[0][jind][m] = PB[0]*AIZZ_[0][jind-jxm][m] -
228 PC[0]*AIZZ_[0][jind-jxm][m+1];
229 AIXY_[0][jind][m] = PB[0]*AIXY_[0][jind-jxm][m] -
230 PC[0]*AIXY_[0][jind-jxm][m+1] +
231 AIY_[0][jind-jxm][m+1];
232 AIXZ_[0][jind][m] = PB[0]*AIXZ_[0][jind-jxm][m] -
233 PC[0]*AIXZ_[0][jind-jxm][m+1] +
234 AIZ_[0][jind-jxm][m+1];
235 AIYZ_[0][jind][m] = PB[0]*AIYZ_[0][jind-jxm][m] -
236 PC[0]*AIYZ_[0][jind-jxm][m+1];
237 }
238 if (jx > 1) {
239 for(m=0;m<=mmax-b;m++)
240 AI0_[0][jind][m] += pp*(jx-1)*(AI0_[0][jind-2*jxm][m] -
241 AI0_[0][jind-2*jxm][m+1]);
242 if (Order >= 1)
243 for(m=0;m<=mmax-b-1;m++) {
244 AIX_[0][jind][m] += pp*(jx-1)*(AIX_[0][jind-2*jxm][m] -
245 AIX_[0][jind-2*jxm][m+1]);
246 AIY_[0][jind][m] += pp*(jx-1)*(AIY_[0][jind-2*jxm][m] -
247 AIY_[0][jind-2*jxm][m+1]);
248 AIZ_[0][jind][m] += pp*(jx-1)*(AIZ_[0][jind-2*jxm][m] -
249 AIZ_[0][jind-2*jxm][m+1]);
250 }
251 if (Order >= 2)
252 for(m=0;m<=mmax-b-2;m++) {
253 AIXX_[0][jind][m] += pp*(jx-1)*(AIXX_[0][jind-2*jxm][m] -
254 AIXX_[0][jind-2*jxm][m+1]);
255 AIYY_[0][jind][m] += pp*(jx-1)*(AIYY_[0][jind-2*jxm][m] -
256 AIYY_[0][jind-2*jxm][m+1]);
257 AIZZ_[0][jind][m] += pp*(jx-1)*(AIZZ_[0][jind-2*jxm][m] -
258 AIZZ_[0][jind-2*jxm][m+1]);
259 AIXY_[0][jind][m] += pp*(jx-1)*(AIXY_[0][jind-2*jxm][m] -
260 AIXY_[0][jind-2*jxm][m+1]);
261 AIXZ_[0][jind][m] += pp*(jx-1)*(AIXZ_[0][jind-2*jxm][m] -
262 AIXZ_[0][jind-2*jxm][m+1]);
263 AIYZ_[0][jind][m] += pp*(jx-1)*(AIYZ_[0][jind-2*jxm][m] -
264 AIYZ_[0][jind-2*jxm][m+1]);
265 }
266 }
267 }
268 else
269 throw ProgrammingError("Logic error in Coulomb 1-e OS RR code",__FILE__,__LINE__);
270 }
271
272
273
274 /* The following fragment cannot be vectorized easily, I guess :-) */
275 /* Upward recursion in i with all possible j's */
276
277 for(b=0;b<=jang;b++)
278 for(jx=0;jx<=b;jx++)
279 for(jy=0;jy<=b-jx;jy++) {
280 jz = b-jx-jy;
281 jind = jx*jxm + jy*jym + jz*jzm;
282 for(a=1;a<=iang;a++)
283 for(ix=0;ix<=a;ix++)
284 for(iy=0;iy<=a-ix;iy++) {
285 iz = a-ix-iy;
286 iind = ix*ixm + iy*iym + iz*izm;
287 if (iz > 0) {
288 for(m=0;m<=mmax-a-b;m++)
289 AI0_[iind][jind][m] = PA[2]*AI0_[iind-izm][jind][m] -
290 PC[2]*AI0_[iind-izm][jind][m+1];
291 if (Order >= 1)
292 for(m=0;m<=mmax-a-b-1;m++) { /* Electric field integrals */
293 AIX_[iind][jind][m] = PA[2]*AIX_[iind-izm][jind][m] -
294 PC[2]*AIX_[iind-izm][jind][m+1];
295 AIY_[iind][jind][m] = PA[2]*AIY_[iind-izm][jind][m] -
296 PC[2]*AIY_[iind-izm][jind][m+1];
297 AIZ_[iind][jind][m] = PA[2]*AIZ_[iind-izm][jind][m] -
298 PC[2]*AIZ_[iind-izm][jind][m+1] +
299 AI0_[iind-izm][jind][m+1];
300 }
301 if (Order >= 2)
302 for(m=0;m<=mmax-a-b-2;m++) { /* Gradients of the electric field */
303 AIXX_[iind][jind][m] = PA[2]*AIXX_[iind-izm][jind][m] -
304 PC[2]*AIXX_[iind-izm][jind][m+1];
305 AIYY_[iind][jind][m] = PA[2]*AIYY_[iind-izm][jind][m] -
306 PC[2]*AIYY_[iind-izm][jind][m+1];
307 AIZZ_[iind][jind][m] = PA[2]*AIZZ_[iind-izm][jind][m] -
308 PC[2]*AIZZ_[iind-izm][jind][m+1] +
309 2*AIZ_[iind-izm][jind][m+1];
310 AIXY_[iind][jind][m] = PA[2]*AIXY_[iind-izm][jind][m] -
311 PC[2]*AIXY_[iind-izm][jind][m+1];
312 AIXZ_[iind][jind][m] = PA[2]*AIXZ_[iind-izm][jind][m] -
313 PC[2]*AIXZ_[iind-izm][jind][m+1] +
314 AIX_[iind-izm][jind][m+1];
315 AIYZ_[iind][jind][m] = PA[2]*AIYZ_[iind-izm][jind][m] -
316 PC[2]*AIYZ_[iind-izm][jind][m+1] +
317 AIY_[iind-izm][jind][m+1];
318 }
319 if (iz > 1) {
320 for(m=0;m<=mmax-a-b;m++)
321 AI0_[iind][jind][m] += pp*(iz-1)*
322 (AI0_[iind-2*izm][jind][m] - AI0_[iind-2*izm][jind][m+1]);
323 if (Order >= 1)
324 for(m=0;m<=mmax-a-b-1;m++) {
325 AIX_[iind][jind][m] += pp*(iz-1)*(AIX_[iind-2*izm][jind][m] -
326 AIX_[iind-2*izm][jind][m+1]);
327 AIY_[iind][jind][m] += pp*(iz-1)*(AIY_[iind-2*izm][jind][m] -
328 AIY_[iind-2*izm][jind][m+1]);
329 AIZ_[iind][jind][m] += pp*(iz-1)*(AIZ_[iind-2*izm][jind][m] -
330 AIZ_[iind-2*izm][jind][m+1]);
331 }
332 if (Order >= 2)
333 for(m=0;m<=mmax-a-b-2;m++) {
334 AIXX_[iind][jind][m] += pp*(iz-1)*(AIXX_[iind-2*izm][jind][m] -
335 AIXX_[iind-2*izm][jind][m+1]);
336 AIYY_[iind][jind][m] += pp*(iz-1)*(AIYY_[iind-2*izm][jind][m] -
337 AIYY_[iind-2*izm][jind][m+1]);
338 AIZZ_[iind][jind][m] += pp*(iz-1)*(AIZZ_[iind-2*izm][jind][m] -
339 AIZZ_[iind-2*izm][jind][m+1]);
340 AIXY_[iind][jind][m] += pp*(iz-1)*(AIXY_[iind-2*izm][jind][m] -
341 AIXY_[iind-2*izm][jind][m+1]);
342 AIXZ_[iind][jind][m] += pp*(iz-1)*(AIXZ_[iind-2*izm][jind][m] -
343 AIXZ_[iind-2*izm][jind][m+1]);
344 AIYZ_[iind][jind][m] += pp*(iz-1)*(AIYZ_[iind-2*izm][jind][m] -
345 AIYZ_[iind-2*izm][jind][m+1]);
346 }
347 }
348 if (jz > 0) {
349 for(m=0;m<=mmax-a-b;m++)
350 AI0_[iind][jind][m] += pp*jz*
351 (AI0_[iind-izm][jind-jzm][m] - AI0_[iind-izm][jind-jzm][m+1]);
352 if (Order >= 1)
353 for(m=0;m<=mmax-a-b-1;m++) {
354 AIX_[iind][jind][m] += pp*jz*(AIX_[iind-izm][jind-jzm][m] -
355 AIX_[iind-izm][jind-jzm][m+1]);
356 AIY_[iind][jind][m] += pp*jz*(AIY_[iind-izm][jind-jzm][m] -
357 AIY_[iind-izm][jind-jzm][m+1]);
358 AIZ_[iind][jind][m] += pp*jz*(AIZ_[iind-izm][jind-jzm][m] -
359 AIZ_[iind-izm][jind-jzm][m+1]);
360 }
361 if (Order >= 2)
362 for(m=0;m<=mmax-a-b-2;m++) {
363 AIXX_[iind][jind][m] += pp*jz*(AIXX_[iind-izm][jind-jzm][m] -
364 AIXX_[iind-izm][jind-jzm][m+1]);
365 AIYY_[iind][jind][m] += pp*jz*(AIYY_[iind-izm][jind-jzm][m] -
366 AIYY_[iind-izm][jind-jzm][m+1]);
367 AIZZ_[iind][jind][m] += pp*jz*(AIZZ_[iind-izm][jind-jzm][m] -
368 AIZZ_[iind-izm][jind-jzm][m+1]);
369 AIXY_[iind][jind][m] += pp*jz*(AIXY_[iind-izm][jind-jzm][m] -
370 AIXY_[iind-izm][jind-jzm][m+1]);
371 AIXZ_[iind][jind][m] += pp*jz*(AIXZ_[iind-izm][jind-jzm][m] -
372 AIXZ_[iind-izm][jind-jzm][m+1]);
373 AIYZ_[iind][jind][m] += pp*jz*(AIYZ_[iind-izm][jind-jzm][m] -
374 AIYZ_[iind-izm][jind-jzm][m+1]);
375 }
376 }
377 }
378 else
379 if (iy > 0) {
380 for(m=0;m<=mmax-a-b;m++)
381 AI0_[iind][jind][m] = PA[1]*AI0_[iind-iym][jind][m] -
382 PC[1]*AI0_[iind-iym][jind][m+1];
383 if (Order >= 1)
384 for(m=0;m<=mmax-a-b-1;m++) {
385 AIX_[iind][jind][m] = PA[1]*AIX_[iind-iym][jind][m] -
386 PC[1]*AIX_[iind-iym][jind][m+1];
387 AIY_[iind][jind][m] = PA[1]*AIY_[iind-iym][jind][m] -
388 PC[1]*AIY_[iind-iym][jind][m+1] +
389 AI0_[iind-iym][jind][m+1];
390 AIZ_[iind][jind][m] = PA[1]*AIZ_[iind-iym][jind][m] -
391 PC[1]*AIZ_[iind-iym][jind][m+1];
392 }
393 if (Order >= 2)
394 for(m=0;m<=mmax-a-b-2;m++) {
395 AIXX_[iind][jind][m] = PA[1]*AIXX_[iind-iym][jind][m] -
396 PC[1]*AIXX_[iind-iym][jind][m+1];
397 AIYY_[iind][jind][m] = PA[1]*AIYY_[iind-iym][jind][m] -
398 PC[1]*AIYY_[iind-iym][jind][m+1] +
399 2*AIY_[iind-iym][jind][m+1];
400 AIZZ_[iind][jind][m] = PA[1]*AIZZ_[iind-iym][jind][m] -
401 PC[1]*AIZZ_[iind-iym][jind][m+1];
402 AIXY_[iind][jind][m] = PA[1]*AIXY_[iind-iym][jind][m] -
403 PC[1]*AIXY_[iind-iym][jind][m+1] +
404 AIX_[iind-iym][jind][m+1];
405 AIXZ_[iind][jind][m] = PA[1]*AIXZ_[iind-iym][jind][m] -
406 PC[1]*AIXZ_[iind-iym][jind][m+1];
407 AIYZ_[iind][jind][m] = PA[1]*AIYZ_[iind-iym][jind][m] -
408 PC[1]*AIYZ_[iind-iym][jind][m+1] +
409 AIZ_[iind-iym][jind][m+1];
410 }
411 if (iy > 1) {
412 for(m=0;m<=mmax-a-b;m++)
413 AI0_[iind][jind][m] += pp*(iy-1)*
414 (AI0_[iind-2*iym][jind][m] - AI0_[iind-2*iym][jind][m+1]);
415 if (Order >= 1)
416 for(m=0;m<=mmax-a-b-1;m++) {
417 AIX_[iind][jind][m] += pp*(iy-1)*(AIX_[iind-2*iym][jind][m] -
418 AIX_[iind-2*iym][jind][m+1]);
419 AIY_[iind][jind][m] += pp*(iy-1)*(AIY_[iind-2*iym][jind][m] -
420 AIY_[iind-2*iym][jind][m+1]);
421 AIZ_[iind][jind][m] += pp*(iy-1)*(AIZ_[iind-2*iym][jind][m] -
422 AIZ_[iind-2*iym][jind][m+1]);
423 }
424 if (Order >= 2)
425 for(m=0;m<=mmax-a-b-2;m++) {
426 AIXX_[iind][jind][m] += pp*(iy-1)*(AIXX_[iind-2*iym][jind][m] -
427 AIXX_[iind-2*iym][jind][m+1]);
428 AIYY_[iind][jind][m] += pp*(iy-1)*(AIYY_[iind-2*iym][jind][m] -
429 AIYY_[iind-2*iym][jind][m+1]);
430 AIZZ_[iind][jind][m] += pp*(iy-1)*(AIZZ_[iind-2*iym][jind][m] -
431 AIZZ_[iind-2*iym][jind][m+1]);
432 AIXY_[iind][jind][m] += pp*(iy-1)*(AIXY_[iind-2*iym][jind][m] -
433 AIXY_[iind-2*iym][jind][m+1]);
434 AIXZ_[iind][jind][m] += pp*(iy-1)*(AIXZ_[iind-2*iym][jind][m] -
435 AIXZ_[iind-2*iym][jind][m+1]);
436 AIYZ_[iind][jind][m] += pp*(iy-1)*(AIYZ_[iind-2*iym][jind][m] -
437 AIYZ_[iind-2*iym][jind][m+1]);
438 }
439 }
440 if (jy > 0) {
441 for(m=0;m<=mmax-a-b;m++)
442 AI0_[iind][jind][m] += pp*jy*
443 (AI0_[iind-iym][jind-jym][m] - AI0_[iind-iym][jind-jym][m+1]);
444 if (Order >= 1)
445 for(m=0;m<=mmax-a-b-1;m++) {
446 AIX_[iind][jind][m] += pp*jy*(AIX_[iind-iym][jind-jym][m] -
447 AIX_[iind-iym][jind-jym][m+1]);
448 AIY_[iind][jind][m] += pp*jy*(AIY_[iind-iym][jind-jym][m] -
449 AIY_[iind-iym][jind-jym][m+1]);
450 AIZ_[iind][jind][m] += pp*jy*(AIZ_[iind-iym][jind-jym][m] -
451 AIZ_[iind-iym][jind-jym][m+1]);
452 }
453 if (Order >= 2)
454 for(m=0;m<=mmax-a-b-2;m++) {
455 AIXX_[iind][jind][m] += pp*jy*(AIXX_[iind-iym][jind-jym][m] -
456 AIXX_[iind-iym][jind-jym][m+1]);
457 AIYY_[iind][jind][m] += pp*jy*(AIYY_[iind-iym][jind-jym][m] -
458 AIYY_[iind-iym][jind-jym][m+1]);
459 AIZZ_[iind][jind][m] += pp*jy*(AIZZ_[iind-iym][jind-jym][m] -
460 AIZZ_[iind-iym][jind-jym][m+1]);
461 AIXY_[iind][jind][m] += pp*jy*(AIXY_[iind-iym][jind-jym][m] -
462 AIXY_[iind-iym][jind-jym][m+1]);
463 AIXZ_[iind][jind][m] += pp*jy*(AIXZ_[iind-iym][jind-jym][m] -
464 AIXZ_[iind-iym][jind-jym][m+1]);
465 AIYZ_[iind][jind][m] += pp*jy*(AIYZ_[iind-iym][jind-jym][m] -
466 AIYZ_[iind-iym][jind-jym][m+1]);
467 }
468 }
469 }
470 else
471 if (ix > 0) {
472 for(m=0;m<=mmax-a-b;m++)
473 AI0_[iind][jind][m] = PA[0]*AI0_[iind-ixm][jind][m] -
474 PC[0]*AI0_[iind-ixm][jind][m+1];
475 if (Order >= 1)
476 for(m=0;m<=mmax-a-b-1;m++) { /* Electric field integrals */
477 AIX_[iind][jind][m] = PA[0]*AIX_[iind-ixm][jind][m] -
478 PC[0]*AIX_[iind-ixm][jind][m+1] +
479 AI0_[iind-ixm][jind][m+1];
480 AIY_[iind][jind][m] = PA[0]*AIY_[iind-ixm][jind][m] -
481 PC[0]*AIY_[iind-ixm][jind][m+1];
482 AIZ_[iind][jind][m] = PA[0]*AIZ_[iind-ixm][jind][m] -
483 PC[0]*AIZ_[iind-ixm][jind][m+1];
484 }
485 if (Order >= 2)
486 for(m=0;m<=mmax-a-b-2;m++) { /* Gradients of the electric field */
487 AIXX_[iind][jind][m] = PA[0]*AIXX_[iind-ixm][jind][m] -
488 PC[0]*AIXX_[iind-ixm][jind][m+1] +
489 2*AIX_[iind-ixm][jind][m+1];
490 AIYY_[iind][jind][m] = PA[0]*AIYY_[iind-ixm][jind][m] -
491 PC[0]*AIYY_[iind-ixm][jind][m+1];
492 AIZZ_[iind][jind][m] = PA[0]*AIZZ_[iind-ixm][jind][m] -
493 PC[0]*AIZZ_[iind-ixm][jind][m+1];
494 AIXY_[iind][jind][m] = PA[0]*AIXY_[iind-ixm][jind][m] -
495 PC[0]*AIXY_[iind-ixm][jind][m+1] +
496 AIY_[iind-ixm][jind][m+1];
497 AIXZ_[iind][jind][m] = PA[0]*AIXZ_[iind-ixm][jind][m] -
498 PC[0]*AIXZ_[iind-ixm][jind][m+1] +
499 AIZ_[iind-ixm][jind][m+1];
500 AIYZ_[iind][jind][m] = PA[0]*AIYZ_[iind-ixm][jind][m] -
501 PC[0]*AIYZ_[iind-ixm][jind][m+1];
502 }
503 if (ix > 1) {
504 for(m=0;m<=mmax-a-b;m++)
505 AI0_[iind][jind][m] += pp*(ix-1)*
506 (AI0_[iind-2*ixm][jind][m] - AI0_[iind-2*ixm][jind][m+1]);
507 if (Order >= 1)
508 for(m=0;m<=mmax-a-b-1;m++) {
509 AIX_[iind][jind][m] += pp*(ix-1)*(AIX_[iind-2*ixm][jind][m] -
510 AIX_[iind-2*ixm][jind][m+1]);
511 AIY_[iind][jind][m] += pp*(ix-1)*(AIY_[iind-2*ixm][jind][m] -
512 AIY_[iind-2*ixm][jind][m+1]);
513 AIZ_[iind][jind][m] += pp*(ix-1)*(AIZ_[iind-2*ixm][jind][m] -
514 AIZ_[iind-2*ixm][jind][m+1]);
515 }
516 if (Order >= 2)
517 for(m=0;m<=mmax-a-b-2;m++) {
518 AIXX_[iind][jind][m] += pp*(ix-1)*(AIXX_[iind-2*ixm][jind][m] -
519 AIXX_[iind-2*ixm][jind][m+1]);
520 AIYY_[iind][jind][m] += pp*(ix-1)*(AIYY_[iind-2*ixm][jind][m] -
521 AIYY_[iind-2*ixm][jind][m+1]);
522 AIZZ_[iind][jind][m] += pp*(ix-1)*(AIZZ_[iind-2*ixm][jind][m] -
523 AIZZ_[iind-2*ixm][jind][m+1]);
524 AIXY_[iind][jind][m] += pp*(ix-1)*(AIXY_[iind-2*ixm][jind][m] -
525 AIXY_[iind-2*ixm][jind][m+1]);
526 AIXZ_[iind][jind][m] += pp*(ix-1)*(AIXZ_[iind-2*ixm][jind][m] -
527 AIXZ_[iind-2*ixm][jind][m+1]);
528 AIYZ_[iind][jind][m] += pp*(ix-1)*(AIYZ_[iind-2*ixm][jind][m] -
529 AIYZ_[iind-2*ixm][jind][m+1]);
530 }
531 }
532 if (jx > 0) {
533 for(m=0;m<=mmax-a-b;m++)
534 AI0_[iind][jind][m] += pp*jx*
535 (AI0_[iind-ixm][jind-jxm][m] - AI0_[iind-ixm][jind-jxm][m+1]);
536 if (Order >= 1)
537 for(m=0;m<=mmax-a-b-1;m++) {
538 AIX_[iind][jind][m] += pp*jx*(AIX_[iind-ixm][jind-jxm][m] -
539 AIX_[iind-ixm][jind-jxm][m+1]);
540 AIY_[iind][jind][m] += pp*jx*(AIY_[iind-ixm][jind-jxm][m] -
541 AIY_[iind-ixm][jind-jxm][m+1]);
542 AIZ_[iind][jind][m] += pp*jx*(AIZ_[iind-ixm][jind-jxm][m] -
543 AIZ_[iind-ixm][jind-jxm][m+1]);
544 }
545 if (Order >= 2)
546 for(m=0;m<=mmax-a-b-2;m++) {
547 AIXX_[iind][jind][m] += pp*jx*(AIXX_[iind-ixm][jind-jxm][m] -
548 AIXX_[iind-ixm][jind-jxm][m+1]);
549 AIYY_[iind][jind][m] += pp*jx*(AIYY_[iind-ixm][jind-jxm][m] -
550 AIYY_[iind-ixm][jind-jxm][m+1]);
551 AIZZ_[iind][jind][m] += pp*jx*(AIZZ_[iind-ixm][jind-jxm][m] -
552 AIZZ_[iind-ixm][jind-jxm][m+1]);
553 AIXY_[iind][jind][m] += pp*jx*(AIXY_[iind-ixm][jind-jxm][m] -
554 AIXY_[iind-ixm][jind-jxm][m+1]);
555 AIXZ_[iind][jind][m] += pp*jx*(AIXZ_[iind-ixm][jind-jxm][m] -
556 AIXZ_[iind-ixm][jind-jxm][m+1]);
557 AIYZ_[iind][jind][m] += pp*jx*(AIYZ_[iind-ixm][jind-jxm][m] -
558 AIYZ_[iind-ixm][jind-jxm][m+1]);
559 }
560 }
561 }
562 else
563 throw ProgrammingError("Logic error in Coulomb 1-e OS RR code",__FILE__,__LINE__);
564 }
565 }
566
567 return;
568}
569
Generated at Wed Sep 25 2024 02:45:29 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.12.0.