MPQC  3.0.0-alpha
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 
30 template <int Order>
31 void
32 Int1eLibint2::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 Sun Jan 26 2020 23:23:58 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.