4
#include "prototypes.h"
6
void calc_f(double *, int, double);
9
/* Recursion relations are taken from Obara and Saika paper
10
JCP 84, 3963, 1986. */
12
void MI_OSrecurs(double pax, double pay, double paz,
13
double pbx, double pby, double pbz, double gamma,
14
int lmaxi, int lmaxj, int maxm)
17
double pp = 1/(2*gamma);
19
/* Computing starting integrals for recursive procedure */
22
MIX[0][0][2] = MIY[0][0][2] = MIZ[0][0][2] = pp;
24
/* Upward recursion in j for i=0 */
27
for(k=0;k<=maxm;k++) {
28
MIX[0][j+1][k] = pbx*MIX[0][j][k];
29
MIY[0][j+1][k] = pby*MIY[0][j][k];
30
MIZ[0][j+1][k] = pbz*MIZ[0][j][k];
32
MIX[0][j+1][k] += j*pp*MIX[0][j-1][k];
33
MIY[0][j+1][k] += j*pp*MIY[0][j-1][k];
34
MIZ[0][j+1][k] += j*pp*MIZ[0][j-1][k];
37
MIX[0][j+1][k] += k*pp*MIX[0][j][k-1];
38
MIY[0][j+1][k] += k*pp*MIY[0][j][k-1];
39
MIZ[0][j+1][k] += k*pp*MIZ[0][j][k-1];
43
/* Upward recursion in i for all j's */
47
for(k=0;k<=maxm;k++) {
48
MIX[i+1][j][k] = pax*MIX[i][j][k];
49
MIY[i+1][j][k] = pay*MIY[i][j][k];
50
MIZ[i+1][j][k] = paz*MIZ[i][j][k];
52
MIX[i+1][j][k] += i*pp*MIX[i-1][j][k];
53
MIY[i+1][j][k] += i*pp*MIY[i-1][j][k];
54
MIZ[i+1][j][k] += i*pp*MIZ[i-1][j][k];
57
MIX[i+1][j][k] += j*pp*MIX[i][j-1][k];
58
MIY[i+1][j][k] += j*pp*MIY[i][j-1][k];
59
MIZ[i+1][j][k] += j*pp*MIZ[i][j-1][k];
62
MIX[i+1][j][k] += k*pp*MIX[i][j][k-1];
63
MIY[i+1][j][k] += k*pp*MIY[i][j][k-1];
64
MIZ[i+1][j][k] += k*pp*MIZ[i][j][k-1];
70
/* Recurrence relation are from the same paper - pp. 3971-3972 */
72
void AI_OSrecurs(double pax, double pay, double paz,
73
double pbx, double pby, double pbz,
74
double pcx, double pcy, double pcz,
75
double gamma, int iang, int jang)
84
int ix,iy,iz,jx,jy,jz;
86
double pp = 1/(2*gamma);
87
double mmax = 2*lmax+2;
88
double tmp = sqrt(gamma)*M_2_SQRTPI;
89
double u = gamma*(pcx*pcx + pcy*pcy + pcz*pcz);
92
F = init_array(mmax+1);
96
/* Computing starting integrals for recursion */
99
AI0[0][0][m] = tmp*F[m];
100
for(m=0;m<=mmax-1;m++) {
101
AIX[0][0][m] = 2*gamma*pcx*AI0[0][0][m+1];
102
AIY[0][0][m] = 2*gamma*pcy*AI0[0][0][m+1];
103
AIZ[0][0][m] = 2*gamma*pcz*AI0[0][0][m+1];
105
for(m=0;m<=mmax-2;m++) {
106
AIXX[0][0][m] = 4*gamma*gamma*pcx*pcx*AI0[0][0][m+2] -
107
2*gamma*AI0[0][0][m+1];
108
AIYY[0][0][m] = 4*gamma*gamma*pcy*pcy*AI0[0][0][m+2] -
109
2*gamma*AI0[0][0][m+1];
110
AIZZ[0][0][m] = 4*gamma*gamma*pcz*pcz*AI0[0][0][m+2] -
111
2*gamma*AI0[0][0][m+1];
112
AIXY[0][0][m] = 4*gamma*gamma*pcx*pcy*AI0[0][0][m+2];
113
AIXZ[0][0][m] = 4*gamma*gamma*pcx*pcz*AI0[0][0][m+2];
114
AIYZ[0][0][m] = 4*gamma*gamma*pcy*pcz*AI0[0][0][m+2];
118
/* Upward recursion in j with i=0 */
122
for(jy=0;jy<=b-jx;jy++) {
124
jind = jx*jxm+jy*jym+jz*jzm;
126
for(m=0;m<=mmax-b;m++) /* Electrostatic potential integrals */
127
AI0[0][jind][m] = pbz*AI0[0][jind-jzm][m] -
128
pcz*AI0[0][jind-jzm][m+1];
129
for(m=0;m<=mmax-b-1;m++) { /* Electric field integrals */
130
AIX[0][jind][m] = pbz*AIX[0][jind-jzm][m] -
131
pcz*AIX[0][jind-jzm][m+1];
132
AIY[0][jind][m] = pbz*AIY[0][jind-jzm][m] -
133
pcz*AIY[0][jind-jzm][m+1];
134
AIZ[0][jind][m] = pbz*AIZ[0][jind-jzm][m] -
135
pcz*AIZ[0][jind-jzm][m+1] +
136
AI0[0][jind-jzm][m+1];
138
for(m=0;m<=mmax-b-2;m++) { /* Gradients of the electric field */
139
AIXX[0][jind][m] = pbz*AIXX[0][jind-jzm][m] -
140
pcz*AIXX[0][jind-jzm][m+1];
141
AIYY[0][jind][m] = pbz*AIYY[0][jind-jzm][m] -
142
pcz*AIYY[0][jind-jzm][m+1];
143
AIZZ[0][jind][m] = pbz*AIZZ[0][jind-jzm][m] -
144
pcz*AIZZ[0][jind-jzm][m+1] +
145
2*AIZ[0][jind-jzm][m+1];
146
AIXY[0][jind][m] = pbz*AIXY[0][jind-jzm][m] -
147
pcz*AIXY[0][jind-jzm][m+1];
148
AIXZ[0][jind][m] = pbz*AIXZ[0][jind-jzm][m] -
149
pcz*AIXZ[0][jind-jzm][m+1] +
150
AIX[0][jind-jzm][m+1];
151
AIYZ[0][jind][m] = pbz*AIYZ[0][jind-jzm][m] -
152
pcz*AIYZ[0][jind-jzm][m+1] +
153
AIY[0][jind-jzm][m+1];
156
for(m=0;m<=mmax-b;m++)
157
AI0[0][jind][m] += pp*(jz-1)*(AI0[0][jind-2*jzm][m] -
158
AI0[0][jind-2*jzm][m+1]);
159
for(m=0;m<=mmax-b-1;m++) {
160
AIX[0][jind][m] += pp*(jz-1)*(AIX[0][jind-2*jzm][m] -
161
AIX[0][jind-2*jzm][m+1]);
162
AIY[0][jind][m] += pp*(jz-1)*(AIY[0][jind-2*jzm][m] -
163
AIY[0][jind-2*jzm][m+1]);
164
AIZ[0][jind][m] += pp*(jz-1)*(AIZ[0][jind-2*jzm][m] -
165
AIZ[0][jind-2*jzm][m+1]);
167
for(m=0;m<=mmax-b-2;m++) {
168
AIXX[0][jind][m] += pp*(jz-1)*(AIXX[0][jind-2*jzm][m] -
169
AIXX[0][jind-2*jzm][m+1]);
170
AIYY[0][jind][m] += pp*(jz-1)*(AIYY[0][jind-2*jzm][m] -
171
AIYY[0][jind-2*jzm][m+1]);
172
AIZZ[0][jind][m] += pp*(jz-1)*(AIZZ[0][jind-2*jzm][m] -
173
AIZZ[0][jind-2*jzm][m+1]);
174
AIXY[0][jind][m] += pp*(jz-1)*(AIXY[0][jind-2*jzm][m] -
175
AIXY[0][jind-2*jzm][m+1]);
176
AIXZ[0][jind][m] += pp*(jz-1)*(AIXZ[0][jind-2*jzm][m] -
177
AIXZ[0][jind-2*jzm][m+1]);
178
AIYZ[0][jind][m] += pp*(jz-1)*(AIYZ[0][jind-2*jzm][m] -
179
AIYZ[0][jind-2*jzm][m+1]);
185
for(m=0;m<=mmax-b;m++)
186
AI0[0][jind][m] = pby*AI0[0][jind-jym][m] -
187
pcy*AI0[0][jind-jym][m+1];
188
for(m=0;m<=mmax-b-1;m++) {
189
AIX[0][jind][m] = pby*AIX[0][jind-jym][m] -
190
pcy*AIX[0][jind-jym][m+1];
191
AIY[0][jind][m] = pby*AIY[0][jind-jym][m] -
192
pcy*AIY[0][jind-jym][m+1] +
193
AI0[0][jind-jym][m+1];
194
AIZ[0][jind][m] = pby*AIZ[0][jind-jym][m] -
195
pcy*AIZ[0][jind-jym][m+1];
197
for(m=0;m<=mmax-b-2;m++) {
198
AIXX[0][jind][m] = pby*AIXX[0][jind-jym][m] -
199
pcy*AIXX[0][jind-jym][m+1];
200
AIYY[0][jind][m] = pby*AIYY[0][jind-jym][m] -
201
pcy*AIYY[0][jind-jym][m+1] +
202
2*AIY[0][jind-jym][m+1];
203
AIZZ[0][jind][m] = pby*AIZZ[0][jind-jym][m] -
204
pcy*AIZZ[0][jind-jym][m+1];
205
AIXY[0][jind][m] = pby*AIXY[0][jind-jym][m] -
206
pcy*AIXY[0][jind-jym][m+1] +
207
AIX[0][jind-jym][m+1];
208
AIXZ[0][jind][m] = pby*AIXZ[0][jind-jym][m] -
209
pcy*AIXZ[0][jind-jym][m+1];
210
AIYZ[0][jind][m] = pby*AIYZ[0][jind-jym][m] -
211
pcy*AIYZ[0][jind-jym][m+1] +
212
AIZ[0][jind-jym][m+1];
215
for(m=0;m<=mmax-b;m++)
216
AI0[0][jind][m] += pp*(jy-1)*(AI0[0][jind-2*jym][m] -
217
AI0[0][jind-2*jym][m+1]);
218
for(m=0;m<=mmax-b-1;m++) {
219
AIX[0][jind][m] += pp*(jy-1)*(AIX[0][jind-2*jym][m] -
220
AIX[0][jind-2*jym][m+1]);
221
AIY[0][jind][m] += pp*(jy-1)*(AIY[0][jind-2*jym][m] -
222
AIY[0][jind-2*jym][m+1]);
223
AIZ[0][jind][m] += pp*(jy-1)*(AIZ[0][jind-2*jym][m] -
224
AIZ[0][jind-2*jym][m+1]);
226
for(m=0;m<=mmax-b-2;m++) {
227
AIXX[0][jind][m] += pp*(jy-1)*(AIXX[0][jind-2*jym][m] -
228
AIXX[0][jind-2*jym][m+1]);
229
AIYY[0][jind][m] += pp*(jy-1)*(AIYY[0][jind-2*jym][m] -
230
AIYY[0][jind-2*jym][m+1]);
231
AIZZ[0][jind][m] += pp*(jy-1)*(AIZZ[0][jind-2*jym][m] -
232
AIZZ[0][jind-2*jym][m+1]);
233
AIXY[0][jind][m] += pp*(jy-1)*(AIXY[0][jind-2*jym][m] -
234
AIXY[0][jind-2*jym][m+1]);
235
AIXZ[0][jind][m] += pp*(jy-1)*(AIXZ[0][jind-2*jym][m] -
236
AIXZ[0][jind-2*jym][m+1]);
237
AIYZ[0][jind][m] += pp*(jy-1)*(AIYZ[0][jind-2*jym][m] -
238
AIYZ[0][jind-2*jym][m+1]);
244
for(m=0;m<=mmax-b;m++)
245
AI0[0][jind][m] = pbx*AI0[0][jind-jxm][m] -
246
pcx*AI0[0][jind-jxm][m+1];
247
for(m=0;m<=mmax-b-1;m++) {
248
AIX[0][jind][m] = pbx*AIX[0][jind-jxm][m] -
249
pcx*AIX[0][jind-jxm][m+1] +
250
AI0[0][jind-jxm][m+1];
251
AIY[0][jind][m] = pbx*AIY[0][jind-jxm][m] -
252
pcx*AIY[0][jind-jxm][m+1];
253
AIZ[0][jind][m] = pbx*AIZ[0][jind-jxm][m] -
254
pcx*AIZ[0][jind-jxm][m+1];
256
for(m=0;m<=mmax-b-2;m++) {
257
AIXX[0][jind][m] = pbx*AIXX[0][jind-jxm][m] -
258
pcx*AIXX[0][jind-jxm][m+1] +
259
2*AIX[0][jind-jxm][m+1];
260
AIYY[0][jind][m] = pbx*AIYY[0][jind-jxm][m] -
261
pcx*AIYY[0][jind-jxm][m+1];
262
AIZZ[0][jind][m] = pbx*AIZZ[0][jind-jxm][m] -
263
pcx*AIZZ[0][jind-jxm][m+1];
264
AIXY[0][jind][m] = pbx*AIXY[0][jind-jxm][m] -
265
pcx*AIXY[0][jind-jxm][m+1] +
266
AIY[0][jind-jxm][m+1];
267
AIXZ[0][jind][m] = pbx*AIXZ[0][jind-jxm][m] -
268
pcx*AIXZ[0][jind-jxm][m+1] +
269
AIZ[0][jind-jxm][m+1];
270
AIYZ[0][jind][m] = pbx*AIYZ[0][jind-jxm][m] -
271
pcx*AIYZ[0][jind-jxm][m+1];
274
for(m=0;m<=mmax-b;m++)
275
AI0[0][jind][m] += pp*(jx-1)*(AI0[0][jind-2*jxm][m] -
276
AI0[0][jind-2*jxm][m+1]);
277
for(m=0;m<=mmax-b-1;m++) {
278
AIX[0][jind][m] += pp*(jx-1)*(AIX[0][jind-2*jxm][m] -
279
AIX[0][jind-2*jxm][m+1]);
280
AIY[0][jind][m] += pp*(jx-1)*(AIY[0][jind-2*jxm][m] -
281
AIY[0][jind-2*jxm][m+1]);
282
AIZ[0][jind][m] += pp*(jx-1)*(AIZ[0][jind-2*jxm][m] -
283
AIZ[0][jind-2*jxm][m+1]);
285
for(m=0;m<=mmax-b-2;m++) {
286
AIXX[0][jind][m] += pp*(jx-1)*(AIXX[0][jind-2*jxm][m] -
287
AIXX[0][jind-2*jxm][m+1]);
288
AIYY[0][jind][m] += pp*(jx-1)*(AIYY[0][jind-2*jxm][m] -
289
AIYY[0][jind-2*jxm][m+1]);
290
AIZZ[0][jind][m] += pp*(jx-1)*(AIZZ[0][jind-2*jxm][m] -
291
AIZZ[0][jind-2*jxm][m+1]);
292
AIXY[0][jind][m] += pp*(jx-1)*(AIXY[0][jind-2*jxm][m] -
293
AIXY[0][jind-2*jxm][m+1]);
294
AIXZ[0][jind][m] += pp*(jx-1)*(AIXZ[0][jind-2*jxm][m] -
295
AIXZ[0][jind-2*jxm][m+1]);
296
AIYZ[0][jind][m] += pp*(jx-1)*(AIYZ[0][jind-2*jxm][m] -
297
AIYZ[0][jind-2*jxm][m+1]);
301
else /* This should and will never happen */
302
punt("There's some error in the AI_OSrecurs algorithm");
307
/* The following fragment cannot be vectorized easily, I guess :-) */
308
/* Upward recursion in i with all possible j's */
312
for(jy=0;jy<=b-jx;jy++) {
314
jind = jx*jxm + jy*jym + jz*jzm;
317
for(iy=0;iy<=a-ix;iy++) {
319
iind = ix*ixm + iy*iym + iz*izm;
321
for(m=0;m<=mmax-a-b;m++)
322
AI0[iind][jind][m] = paz*AI0[iind-izm][jind][m] -
323
pcz*AI0[iind-izm][jind][m+1];
324
for(m=0;m<=mmax-a-b-1;m++) { /* Electric field integrals */
325
AIX[iind][jind][m] = paz*AIX[iind-izm][jind][m] -
326
pcz*AIX[iind-izm][jind][m+1];
327
AIY[iind][jind][m] = paz*AIY[iind-izm][jind][m] -
328
pcz*AIY[iind-izm][jind][m+1];
329
AIZ[iind][jind][m] = paz*AIZ[iind-izm][jind][m] -
330
pcz*AIZ[iind-izm][jind][m+1] +
331
AI0[iind-izm][jind][m+1];
333
for(m=0;m<=mmax-a-b-2;m++) { /* Gradients of the electric field */
334
AIXX[iind][jind][m] = paz*AIXX[iind-izm][jind][m] -
335
pcz*AIXX[iind-izm][jind][m+1];
336
AIYY[iind][jind][m] = paz*AIYY[iind-izm][jind][m] -
337
pcz*AIYY[iind-izm][jind][m+1];
338
AIZZ[iind][jind][m] = paz*AIZZ[iind-izm][jind][m] -
339
pcz*AIZZ[iind-izm][jind][m+1] +
340
2*AIZ[iind-izm][jind][m+1];
341
AIXY[iind][jind][m] = paz*AIXY[iind-izm][jind][m] -
342
pcz*AIXY[iind-izm][jind][m+1];
343
AIXZ[iind][jind][m] = paz*AIXZ[iind-izm][jind][m] -
344
pcz*AIXZ[iind-izm][jind][m+1] +
345
AIX[iind-izm][jind][m+1];
346
AIYZ[iind][jind][m] = paz*AIYZ[iind-izm][jind][m] -
347
pcz*AIYZ[iind-izm][jind][m+1] +
348
AIY[iind-izm][jind][m+1];
351
for(m=0;m<=mmax-a-b;m++)
352
AI0[iind][jind][m] += pp*(iz-1)*
353
(AI0[iind-2*izm][jind][m] - AI0[iind-2*izm][jind][m+1]);
354
for(m=0;m<=mmax-a-b-1;m++) {
355
AIX[iind][jind][m] += pp*(iz-1)*(AIX[iind-2*izm][jind][m] -
356
AIX[iind-2*izm][jind][m+1]);
357
AIY[iind][jind][m] += pp*(iz-1)*(AIY[iind-2*izm][jind][m] -
358
AIY[iind-2*izm][jind][m+1]);
359
AIZ[iind][jind][m] += pp*(iz-1)*(AIZ[iind-2*izm][jind][m] -
360
AIZ[iind-2*izm][jind][m+1]);
362
for(m=0;m<=mmax-a-b-2;m++) {
363
AIXX[iind][jind][m] += pp*(iz-1)*(AIXX[iind-2*izm][jind][m] -
364
AIXX[iind-2*izm][jind][m+1]);
365
AIYY[iind][jind][m] += pp*(iz-1)*(AIYY[iind-2*izm][jind][m] -
366
AIYY[iind-2*izm][jind][m+1]);
367
AIZZ[iind][jind][m] += pp*(iz-1)*(AIZZ[iind-2*izm][jind][m] -
368
AIZZ[iind-2*izm][jind][m+1]);
369
AIXY[iind][jind][m] += pp*(iz-1)*(AIXY[iind-2*izm][jind][m] -
370
AIXY[iind-2*izm][jind][m+1]);
371
AIXZ[iind][jind][m] += pp*(iz-1)*(AIXZ[iind-2*izm][jind][m] -
372
AIXZ[iind-2*izm][jind][m+1]);
373
AIYZ[iind][jind][m] += pp*(iz-1)*(AIYZ[iind-2*izm][jind][m] -
374
AIYZ[iind-2*izm][jind][m+1]);
378
for(m=0;m<=mmax-a-b;m++)
379
AI0[iind][jind][m] += pp*jz*
380
(AI0[iind-izm][jind-jzm][m] - AI0[iind-izm][jind-jzm][m+1]);
381
for(m=0;m<=mmax-a-b-1;m++) {
382
AIX[iind][jind][m] += pp*jz*(AIX[iind-izm][jind-jzm][m] -
383
AIX[iind-izm][jind-jzm][m+1]);
384
AIY[iind][jind][m] += pp*jz*(AIY[iind-izm][jind-jzm][m] -
385
AIY[iind-izm][jind-jzm][m+1]);
386
AIZ[iind][jind][m] += pp*jz*(AIZ[iind-izm][jind-jzm][m] -
387
AIZ[iind-izm][jind-jzm][m+1]);
389
for(m=0;m<=mmax-a-b-2;m++) {
390
AIXX[iind][jind][m] += pp*jz*(AIXX[iind-izm][jind-jzm][m] -
391
AIXX[iind-izm][jind-jzm][m+1]);
392
AIYY[iind][jind][m] += pp*jz*(AIYY[iind-izm][jind-jzm][m] -
393
AIYY[iind-izm][jind-jzm][m+1]);
394
AIZZ[iind][jind][m] += pp*jz*(AIZZ[iind-izm][jind-jzm][m] -
395
AIZZ[iind-izm][jind-jzm][m+1]);
396
AIXY[iind][jind][m] += pp*jz*(AIXY[iind-izm][jind-jzm][m] -
397
AIXY[iind-izm][jind-jzm][m+1]);
398
AIXZ[iind][jind][m] += pp*jz*(AIXZ[iind-izm][jind-jzm][m] -
399
AIXZ[iind-izm][jind-jzm][m+1]);
400
AIYZ[iind][jind][m] += pp*jz*(AIYZ[iind-izm][jind-jzm][m] -
401
AIYZ[iind-izm][jind-jzm][m+1]);
407
for(m=0;m<=mmax-a-b;m++)
408
AI0[iind][jind][m] = pay*AI0[iind-iym][jind][m] -
409
pcy*AI0[iind-iym][jind][m+1];
410
for(m=0;m<=mmax-a-b-1;m++) {
411
AIX[iind][jind][m] = pay*AIX[iind-iym][jind][m] -
412
pcy*AIX[iind-iym][jind][m+1];
413
AIY[iind][jind][m] = pay*AIY[iind-iym][jind][m] -
414
pcy*AIY[iind-iym][jind][m+1] +
415
AI0[iind-iym][jind][m+1];
416
AIZ[iind][jind][m] = pay*AIZ[iind-iym][jind][m] -
417
pcy*AIZ[iind-iym][jind][m+1];
419
for(m=0;m<=mmax-a-b-2;m++) {
420
AIXX[iind][jind][m] = pay*AIXX[iind-iym][jind][m] -
421
pcy*AIXX[iind-iym][jind][m+1];
422
AIYY[iind][jind][m] = pay*AIYY[iind-iym][jind][m] -
423
pcy*AIYY[iind-iym][jind][m+1] +
424
2*AIY[iind-iym][jind][m+1];
425
AIZZ[iind][jind][m] = pay*AIZZ[iind-iym][jind][m] -
426
pcy*AIZZ[iind-iym][jind][m+1];
427
AIXY[iind][jind][m] = pay*AIXY[iind-iym][jind][m] -
428
pcy*AIXY[iind-iym][jind][m+1] +
429
AIX[iind-iym][jind][m+1];
430
AIXZ[iind][jind][m] = pay*AIXZ[iind-iym][jind][m] -
431
pcy*AIXZ[iind-iym][jind][m+1];
432
AIYZ[iind][jind][m] = pay*AIYZ[iind-iym][jind][m] -
433
pcy*AIYZ[iind-iym][jind][m+1] +
434
AIZ[iind-iym][jind][m+1];
437
for(m=0;m<=mmax-a-b;m++)
438
AI0[iind][jind][m] += pp*(iy-1)*
439
(AI0[iind-2*iym][jind][m] - AI0[iind-2*iym][jind][m+1]);
440
for(m=0;m<=mmax-a-b-1;m++) {
441
AIX[iind][jind][m] += pp*(iy-1)*(AIX[iind-2*iym][jind][m] -
442
AIX[iind-2*iym][jind][m+1]);
443
AIY[iind][jind][m] += pp*(iy-1)*(AIY[iind-2*iym][jind][m] -
444
AIY[iind-2*iym][jind][m+1]);
445
AIZ[iind][jind][m] += pp*(iy-1)*(AIZ[iind-2*iym][jind][m] -
446
AIZ[iind-2*iym][jind][m+1]);
448
for(m=0;m<=mmax-a-b-2;m++) {
449
AIXX[iind][jind][m] += pp*(iy-1)*(AIXX[iind-2*iym][jind][m] -
450
AIXX[iind-2*iym][jind][m+1]);
451
AIYY[iind][jind][m] += pp*(iy-1)*(AIYY[iind-2*iym][jind][m] -
452
AIYY[iind-2*iym][jind][m+1]);
453
AIZZ[iind][jind][m] += pp*(iy-1)*(AIZZ[iind-2*iym][jind][m] -
454
AIZZ[iind-2*iym][jind][m+1]);
455
AIXY[iind][jind][m] += pp*(iy-1)*(AIXY[iind-2*iym][jind][m] -
456
AIXY[iind-2*iym][jind][m+1]);
457
AIXZ[iind][jind][m] += pp*(iy-1)*(AIXZ[iind-2*iym][jind][m] -
458
AIXZ[iind-2*iym][jind][m+1]);
459
AIYZ[iind][jind][m] += pp*(iy-1)*(AIYZ[iind-2*iym][jind][m] -
460
AIYZ[iind-2*iym][jind][m+1]);
464
for(m=0;m<=mmax-a-b;m++)
465
AI0[iind][jind][m] += pp*jy*
466
(AI0[iind-iym][jind-jym][m] - AI0[iind-iym][jind-jym][m+1]);
467
for(m=0;m<=mmax-a-b-1;m++) {
468
AIX[iind][jind][m] += pp*jy*(AIX[iind-iym][jind-jym][m] -
469
AIX[iind-iym][jind-jym][m+1]);
470
AIY[iind][jind][m] += pp*jy*(AIY[iind-iym][jind-jym][m] -
471
AIY[iind-iym][jind-jym][m+1]);
472
AIZ[iind][jind][m] += pp*jy*(AIZ[iind-iym][jind-jym][m] -
473
AIZ[iind-iym][jind-jym][m+1]);
475
for(m=0;m<=mmax-a-b-2;m++) {
476
AIXX[iind][jind][m] += pp*jy*(AIXX[iind-iym][jind-jym][m] -
477
AIXX[iind-iym][jind-jym][m+1]);
478
AIYY[iind][jind][m] += pp*jy*(AIYY[iind-iym][jind-jym][m] -
479
AIYY[iind-iym][jind-jym][m+1]);
480
AIZZ[iind][jind][m] += pp*jy*(AIZZ[iind-iym][jind-jym][m] -
481
AIZZ[iind-iym][jind-jym][m+1]);
482
AIXY[iind][jind][m] += pp*jy*(AIXY[iind-iym][jind-jym][m] -
483
AIXY[iind-iym][jind-jym][m+1]);
484
AIXZ[iind][jind][m] += pp*jy*(AIXZ[iind-iym][jind-jym][m] -
485
AIXZ[iind-iym][jind-jym][m+1]);
486
AIYZ[iind][jind][m] += pp*jy*(AIYZ[iind-iym][jind-jym][m] -
487
AIYZ[iind-iym][jind-jym][m+1]);
493
for(m=0;m<=mmax-a-b;m++)
494
AI0[iind][jind][m] = pax*AI0[iind-ixm][jind][m] -
495
pcx*AI0[iind-ixm][jind][m+1];
496
for(m=0;m<=mmax-a-b-1;m++) { /* Electric field integrals */
497
AIX[iind][jind][m] = pax*AIX[iind-ixm][jind][m] -
498
pcx*AIX[iind-ixm][jind][m+1] +
499
AI0[iind-ixm][jind][m+1];
500
AIY[iind][jind][m] = pax*AIY[iind-ixm][jind][m] -
501
pcx*AIY[iind-ixm][jind][m+1];
502
AIZ[iind][jind][m] = pax*AIZ[iind-ixm][jind][m] -
503
pcx*AIZ[iind-ixm][jind][m+1];
505
for(m=0;m<=mmax-a-b-2;m++) { /* Gradients of the electric field */
506
AIXX[iind][jind][m] = pax*AIXX[iind-ixm][jind][m] -
507
pcx*AIXX[iind-ixm][jind][m+1] +
508
2*AIX[iind-ixm][jind][m+1];
509
AIYY[iind][jind][m] = pax*AIYY[iind-ixm][jind][m] -
510
pcx*AIYY[iind-ixm][jind][m+1];
511
AIZZ[iind][jind][m] = pax*AIZZ[iind-ixm][jind][m] -
512
pcx*AIZZ[iind-ixm][jind][m+1];
513
AIXY[iind][jind][m] = pax*AIXY[iind-ixm][jind][m] -
514
pcx*AIXY[iind-ixm][jind][m+1] +
515
AIY[iind-ixm][jind][m+1];
516
AIXZ[iind][jind][m] = pax*AIXZ[iind-ixm][jind][m] -
517
pcx*AIXZ[iind-ixm][jind][m+1] +
518
AIZ[iind-ixm][jind][m+1];
519
AIYZ[iind][jind][m] = pax*AIYZ[iind-ixm][jind][m] -
520
pcx*AIYZ[iind-ixm][jind][m+1];
523
for(m=0;m<=mmax-a-b;m++)
524
AI0[iind][jind][m] += pp*(ix-1)*
525
(AI0[iind-2*ixm][jind][m] - AI0[iind-2*ixm][jind][m+1]);
526
for(m=0;m<=mmax-a-b-1;m++) {
527
AIX[iind][jind][m] += pp*(ix-1)*(AIX[iind-2*ixm][jind][m] -
528
AIX[iind-2*ixm][jind][m+1]);
529
AIY[iind][jind][m] += pp*(ix-1)*(AIY[iind-2*ixm][jind][m] -
530
AIY[iind-2*ixm][jind][m+1]);
531
AIZ[iind][jind][m] += pp*(ix-1)*(AIZ[iind-2*ixm][jind][m] -
532
AIZ[iind-2*ixm][jind][m+1]);
534
for(m=0;m<=mmax-a-b-2;m++) {
535
AIXX[iind][jind][m] += pp*(ix-1)*(AIXX[iind-2*ixm][jind][m] -
536
AIXX[iind-2*ixm][jind][m+1]);
537
AIYY[iind][jind][m] += pp*(ix-1)*(AIYY[iind-2*ixm][jind][m] -
538
AIYY[iind-2*ixm][jind][m+1]);
539
AIZZ[iind][jind][m] += pp*(ix-1)*(AIZZ[iind-2*ixm][jind][m] -
540
AIZZ[iind-2*ixm][jind][m+1]);
541
AIXY[iind][jind][m] += pp*(ix-1)*(AIXY[iind-2*ixm][jind][m] -
542
AIXY[iind-2*ixm][jind][m+1]);
543
AIXZ[iind][jind][m] += pp*(ix-1)*(AIXZ[iind-2*ixm][jind][m] -
544
AIXZ[iind-2*ixm][jind][m+1]);
545
AIYZ[iind][jind][m] += pp*(ix-1)*(AIYZ[iind-2*ixm][jind][m] -
546
AIYZ[iind-2*ixm][jind][m+1]);
550
for(m=0;m<=mmax-a-b;m++)
551
AI0[iind][jind][m] += pp*jx*
552
(AI0[iind-ixm][jind-jxm][m] - AI0[iind-ixm][jind-jxm][m+1]);
553
for(m=0;m<=mmax-a-b-1;m++) {
554
AIX[iind][jind][m] += pp*jx*(AIX[iind-ixm][jind-jxm][m] -
555
AIX[iind-ixm][jind-jxm][m+1]);
556
AIY[iind][jind][m] += pp*jx*(AIY[iind-ixm][jind-jxm][m] -
557
AIY[iind-ixm][jind-jxm][m+1]);
558
AIZ[iind][jind][m] += pp*jx*(AIZ[iind-ixm][jind-jxm][m] -
559
AIZ[iind-ixm][jind-jxm][m+1]);
561
for(m=0;m<=mmax-a-b-2;m++) {
562
AIXX[iind][jind][m] += pp*jx*(AIXX[iind-ixm][jind-jxm][m] -
563
AIXX[iind-ixm][jind-jxm][m+1]);
564
AIYY[iind][jind][m] += pp*jx*(AIYY[iind-ixm][jind-jxm][m] -
565
AIYY[iind-ixm][jind-jxm][m+1]);
566
AIZZ[iind][jind][m] += pp*jx*(AIZZ[iind-ixm][jind-jxm][m] -
567
AIZZ[iind-ixm][jind-jxm][m+1]);
568
AIXY[iind][jind][m] += pp*jx*(AIXY[iind-ixm][jind-jxm][m] -
569
AIXY[iind-ixm][jind-jxm][m+1]);
570
AIXZ[iind][jind][m] += pp*jx*(AIXZ[iind-ixm][jind-jxm][m] -
571
AIXZ[iind-ixm][jind-jxm][m+1]);
572
AIYZ[iind][jind][m] += pp*jx*(AIYZ[iind-ixm][jind-jxm][m] -
573
AIYZ[iind-ixm][jind-jxm][m+1]);
577
else /* This should and will never happen */
578
punt("There's some error in the AI_OSrecurs algorithm");
587
/* This function computes infamous integral Fn(t). For its definition
588
see Obara and Saika paper, or Shavitt's chapter in the
589
"Methods in Computational Physics" book (see reference below).
590
This piece of code is from Dr. Justin Fermann's program CINTS */
592
void calc_f(double *F, int n, double t)
600
static double K = 1.0/M_2_SQRTPI;
604
if (t>20.0){ /* For big t's do upward recursion */
609
for(m=0; m<=n-1; m++){
610
F[m+1] = ((2*m + 1)*F[m] - et)/(t2);
613
else { /* For smaller t's compute F with highest n using
614
asymptotic series (see I. Shavitt in
615
"Methods in Computational Physics", ed. B. Alder etal,
616
vol 2, 1963, page 8 */
626
term1 = num/df[m2+2*i+2];
628
} while (fabs(term1) > EPS && i < MAXFACT);
630
for(m=n-1;m>=0;m--){ /* And then do downward recursion */
631
F[m] = (t2*F[m+1] + et)/(2*m+1);
added
removed
src/bin/oeprop/recursion.c
