1 #include "simint/boys/boys.h"
2 #include "simint/ostei/gen/ostei_generated.h"
3 #include "simint/vectorization/vectorization.h"
4 #include <math.h>
5 #include <string.h>
6
7
ostei_d_p_p_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__d_p_p_f)8 int ostei_d_p_p_f(struct simint_multi_shellpair const P,
9 struct simint_multi_shellpair const Q,
10 double screen_tol,
11 double * const restrict work,
12 double * const restrict INT__d_p_p_f)
13 {
14
15 SIMINT_ASSUME_ALIGN_DBL(work);
16 SIMINT_ASSUME_ALIGN_DBL(INT__d_p_p_f);
17 int ab, cd, abcd;
18 int istart, jstart;
19 int iprimcd, nprim_icd, icd;
20 const int check_screen = (screen_tol > 0.0);
21 int i, j;
22 int n;
23 int not_screened;
24 int real_abcd;
25 int iket;
26 int ibra;
27
28 // partition workspace
29 double * const INT__d_s_s_f = work + (SIMINT_NSHELL_SIMD * 0);
30 double * const INT__d_s_s_g = work + (SIMINT_NSHELL_SIMD * 60);
31 double * const INT__f_s_s_f = work + (SIMINT_NSHELL_SIMD * 150);
32 double * const INT__f_s_s_g = work + (SIMINT_NSHELL_SIMD * 250);
33 SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*400);
34 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
35 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_p = primwork + 8;
36 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_d = primwork + 29;
37 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_f = primwork + 65;
38 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_g = primwork + 115;
39 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_p = primwork + 175;
40 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_d = primwork + 202;
41 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_f = primwork + 256;
42 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_g = primwork + 346;
43 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_d = primwork + 481;
44 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_f = primwork + 553;
45 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_g = primwork + 673;
46 SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_f = primwork + 853;
47 SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_g = primwork + 953;
48 double * const hrrwork = (double *)(primwork + 1103);
49 double * const HRR_INT__d_p_s_f = hrrwork + 0;
50 double * const HRR_INT__d_p_s_g = hrrwork + 180;
51
52
53 // Create constants
54 const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
55 const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
56 const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
57 const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
58 const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
59
60
61 ////////////////////////////////////////
62 // Loop over shells and primitives
63 ////////////////////////////////////////
64
65 real_abcd = 0;
66 istart = 0;
67 for(ab = 0; ab < P.nshell12_clip; ++ab)
68 {
69 const int iend = istart + P.nprim12[ab];
70
71 cd = 0;
72 jstart = 0;
73
74 for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
75 {
76 const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
77 int jend = jstart;
78 for(i = 0; i < nshellbatch; i++)
79 jend += Q.nprim12[cd+i];
80
81 // Clear the beginning of the workspace (where we are accumulating integrals)
82 memset(work, 0, SIMINT_NSHELL_SIMD * 400 * sizeof(double));
83 abcd = 0;
84
85
86 for(i = istart; i < iend; ++i)
87 {
88 SIMINT_DBLTYPE bra_screen_max; // only used if check_screen
89
90 if(check_screen)
91 {
92 // Skip this whole thing if always insignificant
93 if((P.screen[i] * Q.screen_max) < screen_tol)
94 continue;
95 bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
96 }
97
98 icd = 0;
99 iprimcd = 0;
100 nprim_icd = Q.nprim12[cd];
101 double * restrict PRIM_PTR_INT__d_s_s_f = INT__d_s_s_f + abcd * 60;
102 double * restrict PRIM_PTR_INT__d_s_s_g = INT__d_s_s_g + abcd * 90;
103 double * restrict PRIM_PTR_INT__f_s_s_f = INT__f_s_s_f + abcd * 100;
104 double * restrict PRIM_PTR_INT__f_s_s_g = INT__f_s_s_g + abcd * 150;
105
106
107
108 // Load these one per loop over i
109 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
110 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
111 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
112
113 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
114
115 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
116 {
117 // calculate the shell offsets
118 // these are the offset from the shell pointed to by cd
119 // for each element
120 int shelloffsets[SIMINT_SIMD_LEN] = {0};
121 int lastoffset = 0;
122 const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
123
124 if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
125 {
126 // Handle if the first element of the vector is a new shell
127 if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
128 {
129 nprim_icd += Q.nprim12[cd + (++icd)];
130 PRIM_PTR_INT__d_s_s_f += 60;
131 PRIM_PTR_INT__d_s_s_g += 90;
132 PRIM_PTR_INT__f_s_s_f += 100;
133 PRIM_PTR_INT__f_s_s_g += 150;
134 }
135 iprimcd++;
136 for(n = 1; n < SIMINT_SIMD_LEN; ++n)
137 {
138 if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
139 {
140 shelloffsets[n] = shelloffsets[n-1] + 1;
141 lastoffset++;
142 nprim_icd += Q.nprim12[cd + (++icd)];
143 }
144 else
145 shelloffsets[n] = shelloffsets[n-1];
146 iprimcd++;
147 }
148 }
149 else
150 iprimcd += SIMINT_SIMD_LEN;
151
152 // Do we have to compute this vector (or has it been screened out)?
153 // (not_screened != 0 means we have to do this vector)
154 if(check_screen)
155 {
156 const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
157 if(vmax < screen_tol)
158 {
159 PRIM_PTR_INT__d_s_s_f += lastoffset*60;
160 PRIM_PTR_INT__d_s_s_g += lastoffset*90;
161 PRIM_PTR_INT__f_s_s_f += lastoffset*100;
162 PRIM_PTR_INT__f_s_s_g += lastoffset*150;
163 continue;
164 }
165 }
166
167 const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
168 const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
169 const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
170 const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
171
172
173 /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
174 SIMINT_DBLTYPE PQ[3];
175 PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
176 PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
177 PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
178 SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
179 R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
180 R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
181
182 const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
183 const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
184 const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
185 const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
186 const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
187 const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
188 const SIMINT_DBLTYPE Q_PB[3] = { SIMINT_DBLLOAD(Q.PB_x, j), SIMINT_DBLLOAD(Q.PB_y, j), SIMINT_DBLLOAD(Q.PB_z, j) };
189
190 // NOTE: Minus sign!
191 const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
192 SIMINT_DBLTYPE aop_PQ[3];
193 aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
194 aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
195 aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
196
197 SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
198 SIMINT_DBLTYPE aoq_PQ[3];
199 aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
200 aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
201 aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
202 // Put a minus sign here so we don't have to in RR routines
203 a_over_q = SIMINT_NEG(a_over_q);
204
205
206 //////////////////////////////////////////////
207 // Fjt function section
208 // Maximum v value: 7
209 //////////////////////////////////////////////
210 // The parameter to the Fjt function
211 const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
212
213
214 const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
215
216
217 boys_F_split(PRIM_INT__s_s_s_s, F_x, 7);
218 SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
219 prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
220 for(n = 0; n <= 7; n++)
221 PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
222
223 //////////////////////////////////////////////
224 // Primitive integrals: Vertical recurrance
225 //////////////////////////////////////////////
226
227 const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
228 const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
229 const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
230 const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
231 const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
232 const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
233 const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
234 const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
235 const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
236
237
238
239 // Forming PRIM_INT__s_s_s_p[7 * 3];
240 for(n = 0; n < 7; ++n) // loop over orders of auxiliary function
241 {
242
243 PRIM_INT__s_s_s_p[n * 3 + 0] = SIMINT_MUL(Q_PB[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
244 PRIM_INT__s_s_s_p[n * 3 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_p[n * 3 + 0]);
245
246 PRIM_INT__s_s_s_p[n * 3 + 1] = SIMINT_MUL(Q_PB[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
247 PRIM_INT__s_s_s_p[n * 3 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_p[n * 3 + 1]);
248
249 PRIM_INT__s_s_s_p[n * 3 + 2] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
250 PRIM_INT__s_s_s_p[n * 3 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_p[n * 3 + 2]);
251
252 }
253
254
255
256 // Forming PRIM_INT__s_s_s_d[6 * 6];
257 for(n = 0; n < 6; ++n) // loop over orders of auxiliary function
258 {
259
260 PRIM_INT__s_s_s_d[n * 6 + 0] = SIMINT_MUL(Q_PB[0], PRIM_INT__s_s_s_p[n * 3 + 0]);
261 PRIM_INT__s_s_s_d[n * 6 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__s_s_s_d[n * 6 + 0]);
262 PRIM_INT__s_s_s_d[n * 6 + 0] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_s_d[n * 6 + 0]);
263
264 PRIM_INT__s_s_s_d[n * 6 + 1] = SIMINT_MUL(Q_PB[1], PRIM_INT__s_s_s_p[n * 3 + 0]);
265 PRIM_INT__s_s_s_d[n * 6 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__s_s_s_d[n * 6 + 1]);
266
267 PRIM_INT__s_s_s_d[n * 6 + 2] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_p[n * 3 + 0]);
268 PRIM_INT__s_s_s_d[n * 6 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__s_s_s_d[n * 6 + 2]);
269
270 PRIM_INT__s_s_s_d[n * 6 + 3] = SIMINT_MUL(Q_PB[1], PRIM_INT__s_s_s_p[n * 3 + 1]);
271 PRIM_INT__s_s_s_d[n * 6 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__s_s_s_d[n * 6 + 3]);
272 PRIM_INT__s_s_s_d[n * 6 + 3] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_s_d[n * 6 + 3]);
273
274 PRIM_INT__s_s_s_d[n * 6 + 4] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_p[n * 3 + 1]);
275 PRIM_INT__s_s_s_d[n * 6 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__s_s_s_d[n * 6 + 4]);
276
277 PRIM_INT__s_s_s_d[n * 6 + 5] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_p[n * 3 + 2]);
278 PRIM_INT__s_s_s_d[n * 6 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__s_s_s_d[n * 6 + 5]);
279 PRIM_INT__s_s_s_d[n * 6 + 5] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_s_d[n * 6 + 5]);
280
281 }
282
283
284
285 // Forming PRIM_INT__s_s_s_f[5 * 10];
286 for(n = 0; n < 5; ++n) // loop over orders of auxiliary function
287 {
288
289 PRIM_INT__s_s_s_f[n * 10 + 0] = SIMINT_MUL(Q_PB[0], PRIM_INT__s_s_s_d[n * 6 + 0]);
290 PRIM_INT__s_s_s_f[n * 10 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 0], PRIM_INT__s_s_s_f[n * 10 + 0]);
291 PRIM_INT__s_s_s_f[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__s_s_s_p[n * 3 + 0]), PRIM_INT__s_s_s_f[n * 10 + 0]);
292
293 PRIM_INT__s_s_s_f[n * 10 + 1] = SIMINT_MUL(Q_PB[1], PRIM_INT__s_s_s_d[n * 6 + 0]);
294 PRIM_INT__s_s_s_f[n * 10 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 0], PRIM_INT__s_s_s_f[n * 10 + 1]);
295
296 PRIM_INT__s_s_s_f[n * 10 + 2] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_d[n * 6 + 0]);
297 PRIM_INT__s_s_s_f[n * 10 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 0], PRIM_INT__s_s_s_f[n * 10 + 2]);
298
299 PRIM_INT__s_s_s_f[n * 10 + 3] = SIMINT_MUL(Q_PB[0], PRIM_INT__s_s_s_d[n * 6 + 3]);
300 PRIM_INT__s_s_s_f[n * 10 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 3], PRIM_INT__s_s_s_f[n * 10 + 3]);
301
302 PRIM_INT__s_s_s_f[n * 10 + 4] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_d[n * 6 + 1]);
303 PRIM_INT__s_s_s_f[n * 10 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 1], PRIM_INT__s_s_s_f[n * 10 + 4]);
304
305 PRIM_INT__s_s_s_f[n * 10 + 5] = SIMINT_MUL(Q_PB[0], PRIM_INT__s_s_s_d[n * 6 + 5]);
306 PRIM_INT__s_s_s_f[n * 10 + 5] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 5], PRIM_INT__s_s_s_f[n * 10 + 5]);
307
308 PRIM_INT__s_s_s_f[n * 10 + 6] = SIMINT_MUL(Q_PB[1], PRIM_INT__s_s_s_d[n * 6 + 3]);
309 PRIM_INT__s_s_s_f[n * 10 + 6] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 3], PRIM_INT__s_s_s_f[n * 10 + 6]);
310 PRIM_INT__s_s_s_f[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__s_s_s_p[n * 3 + 1]), PRIM_INT__s_s_s_f[n * 10 + 6]);
311
312 PRIM_INT__s_s_s_f[n * 10 + 7] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_d[n * 6 + 3]);
313 PRIM_INT__s_s_s_f[n * 10 + 7] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 3], PRIM_INT__s_s_s_f[n * 10 + 7]);
314
315 PRIM_INT__s_s_s_f[n * 10 + 8] = SIMINT_MUL(Q_PB[1], PRIM_INT__s_s_s_d[n * 6 + 5]);
316 PRIM_INT__s_s_s_f[n * 10 + 8] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 5], PRIM_INT__s_s_s_f[n * 10 + 8]);
317
318 PRIM_INT__s_s_s_f[n * 10 + 9] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_d[n * 6 + 5]);
319 PRIM_INT__s_s_s_f[n * 10 + 9] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 5], PRIM_INT__s_s_s_f[n * 10 + 9]);
320 PRIM_INT__s_s_s_f[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__s_s_s_p[n * 3 + 2]), PRIM_INT__s_s_s_f[n * 10 + 9]);
321
322 }
323
324
325 VRR_I_p_s_s_f(
326 PRIM_INT__p_s_s_f,
327 PRIM_INT__s_s_s_f,
328 PRIM_INT__s_s_s_d,
329 P_PA,
330 aop_PQ,
331 one_over_2pq,
332 3);
333
334
335
336 // Forming PRIM_INT__p_s_s_d[3 * 18];
337 for(n = 0; n < 3; ++n) // loop over orders of auxiliary function
338 {
339
340 PRIM_INT__p_s_s_d[n * 18 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_d[n * 6 + 0]);
341 PRIM_INT__p_s_s_d[n * 18 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 0], PRIM_INT__p_s_s_d[n * 18 + 0]);
342 PRIM_INT__p_s_s_d[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_d[n * 18 + 0]);
343
344 PRIM_INT__p_s_s_d[n * 18 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_d[n * 6 + 1]);
345 PRIM_INT__p_s_s_d[n * 18 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 1], PRIM_INT__p_s_s_d[n * 18 + 1]);
346 PRIM_INT__p_s_s_d[n * 18 + 1] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__p_s_s_d[n * 18 + 1]);
347
348 PRIM_INT__p_s_s_d[n * 18 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_d[n * 6 + 2]);
349 PRIM_INT__p_s_s_d[n * 18 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 2], PRIM_INT__p_s_s_d[n * 18 + 2]);
350 PRIM_INT__p_s_s_d[n * 18 + 2] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__p_s_s_d[n * 18 + 2]);
351
352 PRIM_INT__p_s_s_d[n * 18 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_d[n * 6 + 3]);
353 PRIM_INT__p_s_s_d[n * 18 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 3], PRIM_INT__p_s_s_d[n * 18 + 3]);
354
355 PRIM_INT__p_s_s_d[n * 18 + 4] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_d[n * 6 + 4]);
356 PRIM_INT__p_s_s_d[n * 18 + 4] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 4], PRIM_INT__p_s_s_d[n * 18 + 4]);
357
358 PRIM_INT__p_s_s_d[n * 18 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_d[n * 6 + 5]);
359 PRIM_INT__p_s_s_d[n * 18 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_d[(n+1) * 6 + 5], PRIM_INT__p_s_s_d[n * 18 + 5]);
360
361 PRIM_INT__p_s_s_d[n * 18 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_d[n * 6 + 0]);
362 PRIM_INT__p_s_s_d[n * 18 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 0], PRIM_INT__p_s_s_d[n * 18 + 6]);
363
364 PRIM_INT__p_s_s_d[n * 18 + 7] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_d[n * 6 + 1]);
365 PRIM_INT__p_s_s_d[n * 18 + 7] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 1], PRIM_INT__p_s_s_d[n * 18 + 7]);
366 PRIM_INT__p_s_s_d[n * 18 + 7] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_d[n * 18 + 7]);
367
368 PRIM_INT__p_s_s_d[n * 18 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_d[n * 6 + 2]);
369 PRIM_INT__p_s_s_d[n * 18 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 2], PRIM_INT__p_s_s_d[n * 18 + 8]);
370
371 PRIM_INT__p_s_s_d[n * 18 + 9] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_d[n * 6 + 3]);
372 PRIM_INT__p_s_s_d[n * 18 + 9] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 3], PRIM_INT__p_s_s_d[n * 18 + 9]);
373 PRIM_INT__p_s_s_d[n * 18 + 9] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__p_s_s_d[n * 18 + 9]);
374
375 PRIM_INT__p_s_s_d[n * 18 + 10] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_d[n * 6 + 4]);
376 PRIM_INT__p_s_s_d[n * 18 + 10] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 4], PRIM_INT__p_s_s_d[n * 18 + 10]);
377 PRIM_INT__p_s_s_d[n * 18 + 10] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__p_s_s_d[n * 18 + 10]);
378
379 PRIM_INT__p_s_s_d[n * 18 + 11] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_d[n * 6 + 5]);
380 PRIM_INT__p_s_s_d[n * 18 + 11] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_d[(n+1) * 6 + 5], PRIM_INT__p_s_s_d[n * 18 + 11]);
381
382 PRIM_INT__p_s_s_d[n * 18 + 12] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_d[n * 6 + 0]);
383 PRIM_INT__p_s_s_d[n * 18 + 12] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 0], PRIM_INT__p_s_s_d[n * 18 + 12]);
384
385 PRIM_INT__p_s_s_d[n * 18 + 13] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_d[n * 6 + 1]);
386 PRIM_INT__p_s_s_d[n * 18 + 13] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 1], PRIM_INT__p_s_s_d[n * 18 + 13]);
387
388 PRIM_INT__p_s_s_d[n * 18 + 14] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_d[n * 6 + 2]);
389 PRIM_INT__p_s_s_d[n * 18 + 14] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 2], PRIM_INT__p_s_s_d[n * 18 + 14]);
390 PRIM_INT__p_s_s_d[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_d[n * 18 + 14]);
391
392 PRIM_INT__p_s_s_d[n * 18 + 15] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_d[n * 6 + 3]);
393 PRIM_INT__p_s_s_d[n * 18 + 15] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 3], PRIM_INT__p_s_s_d[n * 18 + 15]);
394
395 PRIM_INT__p_s_s_d[n * 18 + 16] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_d[n * 6 + 4]);
396 PRIM_INT__p_s_s_d[n * 18 + 16] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 4], PRIM_INT__p_s_s_d[n * 18 + 16]);
397 PRIM_INT__p_s_s_d[n * 18 + 16] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__p_s_s_d[n * 18 + 16]);
398
399 PRIM_INT__p_s_s_d[n * 18 + 17] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_d[n * 6 + 5]);
400 PRIM_INT__p_s_s_d[n * 18 + 17] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_d[(n+1) * 6 + 5], PRIM_INT__p_s_s_d[n * 18 + 17]);
401 PRIM_INT__p_s_s_d[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__p_s_s_d[n * 18 + 17]);
402
403 }
404
405
406 VRR_I_d_s_s_f(
407 PRIM_INT__d_s_s_f,
408 PRIM_INT__p_s_s_f,
409 PRIM_INT__s_s_s_f,
410 PRIM_INT__p_s_s_d,
411 P_PA,
412 a_over_p,
413 aop_PQ,
414 one_over_2p,
415 one_over_2pq,
416 2);
417
418
419 VRR_L_s_s_s_g(
420 PRIM_INT__s_s_s_g,
421 PRIM_INT__s_s_s_f,
422 PRIM_INT__s_s_s_d,
423 Q_PB,
424 a_over_q,
425 aoq_PQ,
426 one_over_2q,
427 4);
428
429
430 VRR_I_p_s_s_g(
431 PRIM_INT__p_s_s_g,
432 PRIM_INT__s_s_s_g,
433 PRIM_INT__s_s_s_f,
434 P_PA,
435 aop_PQ,
436 one_over_2pq,
437 3);
438
439
440 ostei_general_vrr_I(2, 0, 0, 4, 2,
441 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
442 PRIM_INT__p_s_s_g, PRIM_INT__s_s_s_g, NULL, NULL, PRIM_INT__p_s_s_f, PRIM_INT__d_s_s_g);
443
444
445
446 // Forming PRIM_INT__p_s_s_p[3 * 9];
447 for(n = 0; n < 3; ++n) // loop over orders of auxiliary function
448 {
449
450 PRIM_INT__p_s_s_p[n * 9 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_p[n * 3 + 0]);
451 PRIM_INT__p_s_s_p[n * 9 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_p[n * 9 + 0]);
452 PRIM_INT__p_s_s_p[n * 9 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_p[n * 9 + 0]);
453
454 PRIM_INT__p_s_s_p[n * 9 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_p[n * 3 + 1]);
455 PRIM_INT__p_s_s_p[n * 9 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__p_s_s_p[n * 9 + 1]);
456
457 PRIM_INT__p_s_s_p[n * 9 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_p[n * 3 + 2]);
458 PRIM_INT__p_s_s_p[n * 9 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__p_s_s_p[n * 9 + 2]);
459
460 PRIM_INT__p_s_s_p[n * 9 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_p[n * 3 + 0]);
461 PRIM_INT__p_s_s_p[n * 9 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_p[n * 9 + 3]);
462
463 PRIM_INT__p_s_s_p[n * 9 + 4] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_p[n * 3 + 1]);
464 PRIM_INT__p_s_s_p[n * 9 + 4] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__p_s_s_p[n * 9 + 4]);
465 PRIM_INT__p_s_s_p[n * 9 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_p[n * 9 + 4]);
466
467 PRIM_INT__p_s_s_p[n * 9 + 5] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_p[n * 3 + 2]);
468 PRIM_INT__p_s_s_p[n * 9 + 5] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__p_s_s_p[n * 9 + 5]);
469
470 PRIM_INT__p_s_s_p[n * 9 + 6] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_p[n * 3 + 0]);
471 PRIM_INT__p_s_s_p[n * 9 + 6] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_p[n * 9 + 6]);
472
473 PRIM_INT__p_s_s_p[n * 9 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_p[n * 3 + 1]);
474 PRIM_INT__p_s_s_p[n * 9 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__p_s_s_p[n * 9 + 7]);
475
476 PRIM_INT__p_s_s_p[n * 9 + 8] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_p[n * 3 + 2]);
477 PRIM_INT__p_s_s_p[n * 9 + 8] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__p_s_s_p[n * 9 + 8]);
478 PRIM_INT__p_s_s_p[n * 9 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_p[n * 9 + 8]);
479
480 }
481
482
483 VRR_I_d_s_s_d(
484 PRIM_INT__d_s_s_d,
485 PRIM_INT__p_s_s_d,
486 PRIM_INT__s_s_s_d,
487 PRIM_INT__p_s_s_p,
488 P_PA,
489 a_over_p,
490 aop_PQ,
491 one_over_2p,
492 one_over_2pq,
493 2);
494
495
496 ostei_general_vrr_I(3, 0, 0, 3, 1,
497 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
498 PRIM_INT__d_s_s_f, PRIM_INT__p_s_s_f, NULL, NULL, PRIM_INT__d_s_s_d, PRIM_INT__f_s_s_f);
499
500
501 ostei_general_vrr_I(3, 0, 0, 4, 1,
502 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
503 PRIM_INT__d_s_s_g, PRIM_INT__p_s_s_g, NULL, NULL, PRIM_INT__d_s_s_f, PRIM_INT__f_s_s_g);
504
505
506
507
508 ////////////////////////////////////
509 // Accumulate contracted integrals
510 ////////////////////////////////////
511 if(lastoffset == 0)
512 {
513 contract_all(60, PRIM_INT__d_s_s_f, PRIM_PTR_INT__d_s_s_f);
514 contract_all(90, PRIM_INT__d_s_s_g, PRIM_PTR_INT__d_s_s_g);
515 contract_all(100, PRIM_INT__f_s_s_f, PRIM_PTR_INT__f_s_s_f);
516 contract_all(150, PRIM_INT__f_s_s_g, PRIM_PTR_INT__f_s_s_g);
517 }
518 else
519 {
520 contract(60, shelloffsets, PRIM_INT__d_s_s_f, PRIM_PTR_INT__d_s_s_f);
521 contract(90, shelloffsets, PRIM_INT__d_s_s_g, PRIM_PTR_INT__d_s_s_g);
522 contract(100, shelloffsets, PRIM_INT__f_s_s_f, PRIM_PTR_INT__f_s_s_f);
523 contract(150, shelloffsets, PRIM_INT__f_s_s_g, PRIM_PTR_INT__f_s_s_g);
524 PRIM_PTR_INT__d_s_s_f += lastoffset*60;
525 PRIM_PTR_INT__d_s_s_g += lastoffset*90;
526 PRIM_PTR_INT__f_s_s_f += lastoffset*100;
527 PRIM_PTR_INT__f_s_s_g += lastoffset*150;
528 }
529
530 } // close loop over j
531 } // close loop over i
532
533 //Advance to the next batch
534 jstart = SIMINT_SIMD_ROUND(jend);
535
536 //////////////////////////////////////////////
537 // Contracted integrals: Horizontal recurrance
538 //////////////////////////////////////////////
539
540
541 const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
542
543
544 for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
545 {
546 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
547
548 // set up HRR pointers
549 double const * restrict HRR_INT__d_s_s_f = INT__d_s_s_f + abcd * 60;
550 double const * restrict HRR_INT__d_s_s_g = INT__d_s_s_g + abcd * 90;
551 double const * restrict HRR_INT__f_s_s_f = INT__f_s_s_f + abcd * 100;
552 double const * restrict HRR_INT__f_s_s_g = INT__f_s_s_g + abcd * 150;
553 double * restrict HRR_INT__d_p_p_f = INT__d_p_p_f + real_abcd * 540;
554
555 // form INT__d_p_s_f
556 for(iket = 0; iket < 10; ++iket)
557 {
558 HRR_INT__d_p_s_f[0 * 10 + iket] = HRR_INT__f_s_s_f[0 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_s_f[0 * 10 + iket] );
559
560 HRR_INT__d_p_s_f[1 * 10 + iket] = HRR_INT__f_s_s_f[1 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_s_f[0 * 10 + iket] );
561
562 HRR_INT__d_p_s_f[2 * 10 + iket] = HRR_INT__f_s_s_f[2 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_s_f[0 * 10 + iket] );
563
564 HRR_INT__d_p_s_f[3 * 10 + iket] = HRR_INT__f_s_s_f[1 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_s_f[1 * 10 + iket] );
565
566 HRR_INT__d_p_s_f[4 * 10 + iket] = HRR_INT__f_s_s_f[3 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_s_f[1 * 10 + iket] );
567
568 HRR_INT__d_p_s_f[5 * 10 + iket] = HRR_INT__f_s_s_f[4 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_s_f[1 * 10 + iket] );
569
570 HRR_INT__d_p_s_f[6 * 10 + iket] = HRR_INT__f_s_s_f[2 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_s_f[2 * 10 + iket] );
571
572 HRR_INT__d_p_s_f[7 * 10 + iket] = HRR_INT__f_s_s_f[4 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_s_f[2 * 10 + iket] );
573
574 HRR_INT__d_p_s_f[8 * 10 + iket] = HRR_INT__f_s_s_f[5 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_s_f[2 * 10 + iket] );
575
576 HRR_INT__d_p_s_f[9 * 10 + iket] = HRR_INT__f_s_s_f[3 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_s_f[3 * 10 + iket] );
577
578 HRR_INT__d_p_s_f[10 * 10 + iket] = HRR_INT__f_s_s_f[6 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_s_f[3 * 10 + iket] );
579
580 HRR_INT__d_p_s_f[11 * 10 + iket] = HRR_INT__f_s_s_f[7 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_s_f[3 * 10 + iket] );
581
582 HRR_INT__d_p_s_f[12 * 10 + iket] = HRR_INT__f_s_s_f[4 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_s_f[4 * 10 + iket] );
583
584 HRR_INT__d_p_s_f[13 * 10 + iket] = HRR_INT__f_s_s_f[7 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_s_f[4 * 10 + iket] );
585
586 HRR_INT__d_p_s_f[14 * 10 + iket] = HRR_INT__f_s_s_f[8 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_s_f[4 * 10 + iket] );
587
588 HRR_INT__d_p_s_f[15 * 10 + iket] = HRR_INT__f_s_s_f[5 * 10 + iket] + ( hAB[0] * HRR_INT__d_s_s_f[5 * 10 + iket] );
589
590 HRR_INT__d_p_s_f[16 * 10 + iket] = HRR_INT__f_s_s_f[8 * 10 + iket] + ( hAB[1] * HRR_INT__d_s_s_f[5 * 10 + iket] );
591
592 HRR_INT__d_p_s_f[17 * 10 + iket] = HRR_INT__f_s_s_f[9 * 10 + iket] + ( hAB[2] * HRR_INT__d_s_s_f[5 * 10 + iket] );
593
594 }
595
596
597 // form INT__d_p_s_g
598 for(iket = 0; iket < 15; ++iket)
599 {
600 HRR_INT__d_p_s_g[0 * 15 + iket] = HRR_INT__f_s_s_g[0 * 15 + iket] + ( hAB[0] * HRR_INT__d_s_s_g[0 * 15 + iket] );
601
602 HRR_INT__d_p_s_g[1 * 15 + iket] = HRR_INT__f_s_s_g[1 * 15 + iket] + ( hAB[1] * HRR_INT__d_s_s_g[0 * 15 + iket] );
603
604 HRR_INT__d_p_s_g[2 * 15 + iket] = HRR_INT__f_s_s_g[2 * 15 + iket] + ( hAB[2] * HRR_INT__d_s_s_g[0 * 15 + iket] );
605
606 HRR_INT__d_p_s_g[3 * 15 + iket] = HRR_INT__f_s_s_g[1 * 15 + iket] + ( hAB[0] * HRR_INT__d_s_s_g[1 * 15 + iket] );
607
608 HRR_INT__d_p_s_g[4 * 15 + iket] = HRR_INT__f_s_s_g[3 * 15 + iket] + ( hAB[1] * HRR_INT__d_s_s_g[1 * 15 + iket] );
609
610 HRR_INT__d_p_s_g[5 * 15 + iket] = HRR_INT__f_s_s_g[4 * 15 + iket] + ( hAB[2] * HRR_INT__d_s_s_g[1 * 15 + iket] );
611
612 HRR_INT__d_p_s_g[6 * 15 + iket] = HRR_INT__f_s_s_g[2 * 15 + iket] + ( hAB[0] * HRR_INT__d_s_s_g[2 * 15 + iket] );
613
614 HRR_INT__d_p_s_g[7 * 15 + iket] = HRR_INT__f_s_s_g[4 * 15 + iket] + ( hAB[1] * HRR_INT__d_s_s_g[2 * 15 + iket] );
615
616 HRR_INT__d_p_s_g[8 * 15 + iket] = HRR_INT__f_s_s_g[5 * 15 + iket] + ( hAB[2] * HRR_INT__d_s_s_g[2 * 15 + iket] );
617
618 HRR_INT__d_p_s_g[9 * 15 + iket] = HRR_INT__f_s_s_g[3 * 15 + iket] + ( hAB[0] * HRR_INT__d_s_s_g[3 * 15 + iket] );
619
620 HRR_INT__d_p_s_g[10 * 15 + iket] = HRR_INT__f_s_s_g[6 * 15 + iket] + ( hAB[1] * HRR_INT__d_s_s_g[3 * 15 + iket] );
621
622 HRR_INT__d_p_s_g[11 * 15 + iket] = HRR_INT__f_s_s_g[7 * 15 + iket] + ( hAB[2] * HRR_INT__d_s_s_g[3 * 15 + iket] );
623
624 HRR_INT__d_p_s_g[12 * 15 + iket] = HRR_INT__f_s_s_g[4 * 15 + iket] + ( hAB[0] * HRR_INT__d_s_s_g[4 * 15 + iket] );
625
626 HRR_INT__d_p_s_g[13 * 15 + iket] = HRR_INT__f_s_s_g[7 * 15 + iket] + ( hAB[1] * HRR_INT__d_s_s_g[4 * 15 + iket] );
627
628 HRR_INT__d_p_s_g[14 * 15 + iket] = HRR_INT__f_s_s_g[8 * 15 + iket] + ( hAB[2] * HRR_INT__d_s_s_g[4 * 15 + iket] );
629
630 HRR_INT__d_p_s_g[15 * 15 + iket] = HRR_INT__f_s_s_g[5 * 15 + iket] + ( hAB[0] * HRR_INT__d_s_s_g[5 * 15 + iket] );
631
632 HRR_INT__d_p_s_g[16 * 15 + iket] = HRR_INT__f_s_s_g[8 * 15 + iket] + ( hAB[1] * HRR_INT__d_s_s_g[5 * 15 + iket] );
633
634 HRR_INT__d_p_s_g[17 * 15 + iket] = HRR_INT__f_s_s_g[9 * 15 + iket] + ( hAB[2] * HRR_INT__d_s_s_g[5 * 15 + iket] );
635
636 }
637
638
639 // form INT__d_p_p_f
640 HRR_K_p_f(
641 HRR_INT__d_p_p_f,
642 HRR_INT__d_p_s_f,
643 HRR_INT__d_p_s_g,
644 hCD, 18);
645
646
647 } // close HRR loop
648
649
650 } // close loop cdbatch
651
652 istart = iend;
653 } // close loop over ab
654
655 return P.nshell12_clip * Q.nshell12_clip;
656 }
657
658