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_p_p_i_i(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__p_p_i_i)8 int ostei_p_p_i_i(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__p_p_i_i)
13 {
14
15 SIMINT_ASSUME_ALIGN_DBL(work);
16 SIMINT_ASSUME_ALIGN_DBL(INT__p_p_i_i);
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__p_s_i_s = work + (SIMINT_NSHELL_SIMD * 0);
30 double * const INT__p_s_k_s = work + (SIMINT_NSHELL_SIMD * 84);
31 double * const INT__p_s_l_s = work + (SIMINT_NSHELL_SIMD * 192);
32 double * const INT__p_s_m_s = work + (SIMINT_NSHELL_SIMD * 327);
33 double * const INT__p_s_n_s = work + (SIMINT_NSHELL_SIMD * 492);
34 double * const INT__p_s_o_s = work + (SIMINT_NSHELL_SIMD * 690);
35 double * const INT__p_s_q_s = work + (SIMINT_NSHELL_SIMD * 924);
36 double * const INT__d_s_i_s = work + (SIMINT_NSHELL_SIMD * 1197);
37 double * const INT__d_s_k_s = work + (SIMINT_NSHELL_SIMD * 1365);
38 double * const INT__d_s_l_s = work + (SIMINT_NSHELL_SIMD * 1581);
39 double * const INT__d_s_m_s = work + (SIMINT_NSHELL_SIMD * 1851);
40 double * const INT__d_s_n_s = work + (SIMINT_NSHELL_SIMD * 2181);
41 double * const INT__d_s_o_s = work + (SIMINT_NSHELL_SIMD * 2577);
42 double * const INT__d_s_q_s = work + (SIMINT_NSHELL_SIMD * 3045);
43 SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*3591);
44 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
45 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 15;
46 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 57;
47 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_f_s = primwork + 135;
48 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_g_s = primwork + 255;
49 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_h_s = primwork + 420;
50 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_i_s = primwork + 630;
51 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_k_s = primwork + 882;
52 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_l_s = primwork + 1170;
53 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_m_s = primwork + 1485;
54 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_n_s = primwork + 1815;
55 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_o_s = primwork + 2145;
56 SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_q_s = primwork + 2457;
57 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_h_s = primwork + 2730;
58 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_i_s = primwork + 2856;
59 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_k_s = primwork + 3024;
60 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_l_s = primwork + 3240;
61 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_m_s = primwork + 3510;
62 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_n_s = primwork + 3840;
63 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_o_s = primwork + 4236;
64 SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_q_s = primwork + 4704;
65 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_i_s = primwork + 5250;
66 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_k_s = primwork + 5418;
67 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_l_s = primwork + 5634;
68 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_m_s = primwork + 5904;
69 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_n_s = primwork + 6234;
70 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_o_s = primwork + 6630;
71 SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_q_s = primwork + 7098;
72 double * const hrrwork = (double *)(primwork + 7644);
73 double * const HRR_INT__p_p_i_s = hrrwork + 0;
74 double * const HRR_INT__p_p_i_p = hrrwork + 252;
75 double * const HRR_INT__p_p_i_d = hrrwork + 1008;
76 double * const HRR_INT__p_p_i_f = hrrwork + 2520;
77 double * const HRR_INT__p_p_i_g = hrrwork + 5040;
78 double * const HRR_INT__p_p_i_h = hrrwork + 8820;
79 double * const HRR_INT__p_p_k_s = hrrwork + 14112;
80 double * const HRR_INT__p_p_k_p = hrrwork + 14436;
81 double * const HRR_INT__p_p_k_d = hrrwork + 15408;
82 double * const HRR_INT__p_p_k_f = hrrwork + 17352;
83 double * const HRR_INT__p_p_k_g = hrrwork + 20592;
84 double * const HRR_INT__p_p_k_h = hrrwork + 25452;
85 double * const HRR_INT__p_p_l_s = hrrwork + 32256;
86 double * const HRR_INT__p_p_l_p = hrrwork + 32661;
87 double * const HRR_INT__p_p_l_d = hrrwork + 33876;
88 double * const HRR_INT__p_p_l_f = hrrwork + 36306;
89 double * const HRR_INT__p_p_l_g = hrrwork + 40356;
90 double * const HRR_INT__p_p_m_s = hrrwork + 46431;
91 double * const HRR_INT__p_p_m_p = hrrwork + 46926;
92 double * const HRR_INT__p_p_m_d = hrrwork + 48411;
93 double * const HRR_INT__p_p_m_f = hrrwork + 51381;
94 double * const HRR_INT__p_p_n_s = hrrwork + 56331;
95 double * const HRR_INT__p_p_n_p = hrrwork + 56925;
96 double * const HRR_INT__p_p_n_d = hrrwork + 58707;
97 double * const HRR_INT__p_p_o_s = hrrwork + 62271;
98 double * const HRR_INT__p_p_o_p = hrrwork + 62973;
99 double * const HRR_INT__p_p_q_s = hrrwork + 65079;
100
101
102 // Create constants
103 const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
104 const SIMINT_DBLTYPE const_10 = SIMINT_DBLSET1(10);
105 const SIMINT_DBLTYPE const_11 = SIMINT_DBLSET1(11);
106 const SIMINT_DBLTYPE const_12 = SIMINT_DBLSET1(12);
107 const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
108 const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
109 const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
110 const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
111 const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
112 const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
113 const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
114 const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
115 const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
116
117
118 ////////////////////////////////////////
119 // Loop over shells and primitives
120 ////////////////////////////////////////
121
122 real_abcd = 0;
123 istart = 0;
124 for(ab = 0; ab < P.nshell12_clip; ++ab)
125 {
126 const int iend = istart + P.nprim12[ab];
127
128 cd = 0;
129 jstart = 0;
130
131 for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
132 {
133 const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
134 int jend = jstart;
135 for(i = 0; i < nshellbatch; i++)
136 jend += Q.nprim12[cd+i];
137
138 // Clear the beginning of the workspace (where we are accumulating integrals)
139 memset(work, 0, SIMINT_NSHELL_SIMD * 3591 * sizeof(double));
140 abcd = 0;
141
142
143 for(i = istart; i < iend; ++i)
144 {
145 SIMINT_DBLTYPE bra_screen_max; // only used if check_screen
146
147 if(check_screen)
148 {
149 // Skip this whole thing if always insignificant
150 if((P.screen[i] * Q.screen_max) < screen_tol)
151 continue;
152 bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
153 }
154
155 icd = 0;
156 iprimcd = 0;
157 nprim_icd = Q.nprim12[cd];
158 double * restrict PRIM_PTR_INT__p_s_i_s = INT__p_s_i_s + abcd * 84;
159 double * restrict PRIM_PTR_INT__p_s_k_s = INT__p_s_k_s + abcd * 108;
160 double * restrict PRIM_PTR_INT__p_s_l_s = INT__p_s_l_s + abcd * 135;
161 double * restrict PRIM_PTR_INT__p_s_m_s = INT__p_s_m_s + abcd * 165;
162 double * restrict PRIM_PTR_INT__p_s_n_s = INT__p_s_n_s + abcd * 198;
163 double * restrict PRIM_PTR_INT__p_s_o_s = INT__p_s_o_s + abcd * 234;
164 double * restrict PRIM_PTR_INT__p_s_q_s = INT__p_s_q_s + abcd * 273;
165 double * restrict PRIM_PTR_INT__d_s_i_s = INT__d_s_i_s + abcd * 168;
166 double * restrict PRIM_PTR_INT__d_s_k_s = INT__d_s_k_s + abcd * 216;
167 double * restrict PRIM_PTR_INT__d_s_l_s = INT__d_s_l_s + abcd * 270;
168 double * restrict PRIM_PTR_INT__d_s_m_s = INT__d_s_m_s + abcd * 330;
169 double * restrict PRIM_PTR_INT__d_s_n_s = INT__d_s_n_s + abcd * 396;
170 double * restrict PRIM_PTR_INT__d_s_o_s = INT__d_s_o_s + abcd * 468;
171 double * restrict PRIM_PTR_INT__d_s_q_s = INT__d_s_q_s + abcd * 546;
172
173
174
175 // Load these one per loop over i
176 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
177 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
178 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
179
180 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
181
182 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
183 {
184 // calculate the shell offsets
185 // these are the offset from the shell pointed to by cd
186 // for each element
187 int shelloffsets[SIMINT_SIMD_LEN] = {0};
188 int lastoffset = 0;
189 const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
190
191 if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
192 {
193 // Handle if the first element of the vector is a new shell
194 if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
195 {
196 nprim_icd += Q.nprim12[cd + (++icd)];
197 PRIM_PTR_INT__p_s_i_s += 84;
198 PRIM_PTR_INT__p_s_k_s += 108;
199 PRIM_PTR_INT__p_s_l_s += 135;
200 PRIM_PTR_INT__p_s_m_s += 165;
201 PRIM_PTR_INT__p_s_n_s += 198;
202 PRIM_PTR_INT__p_s_o_s += 234;
203 PRIM_PTR_INT__p_s_q_s += 273;
204 PRIM_PTR_INT__d_s_i_s += 168;
205 PRIM_PTR_INT__d_s_k_s += 216;
206 PRIM_PTR_INT__d_s_l_s += 270;
207 PRIM_PTR_INT__d_s_m_s += 330;
208 PRIM_PTR_INT__d_s_n_s += 396;
209 PRIM_PTR_INT__d_s_o_s += 468;
210 PRIM_PTR_INT__d_s_q_s += 546;
211 }
212 iprimcd++;
213 for(n = 1; n < SIMINT_SIMD_LEN; ++n)
214 {
215 if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
216 {
217 shelloffsets[n] = shelloffsets[n-1] + 1;
218 lastoffset++;
219 nprim_icd += Q.nprim12[cd + (++icd)];
220 }
221 else
222 shelloffsets[n] = shelloffsets[n-1];
223 iprimcd++;
224 }
225 }
226 else
227 iprimcd += SIMINT_SIMD_LEN;
228
229 // Do we have to compute this vector (or has it been screened out)?
230 // (not_screened != 0 means we have to do this vector)
231 if(check_screen)
232 {
233 const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
234 if(vmax < screen_tol)
235 {
236 PRIM_PTR_INT__p_s_i_s += lastoffset*84;
237 PRIM_PTR_INT__p_s_k_s += lastoffset*108;
238 PRIM_PTR_INT__p_s_l_s += lastoffset*135;
239 PRIM_PTR_INT__p_s_m_s += lastoffset*165;
240 PRIM_PTR_INT__p_s_n_s += lastoffset*198;
241 PRIM_PTR_INT__p_s_o_s += lastoffset*234;
242 PRIM_PTR_INT__p_s_q_s += lastoffset*273;
243 PRIM_PTR_INT__d_s_i_s += lastoffset*168;
244 PRIM_PTR_INT__d_s_k_s += lastoffset*216;
245 PRIM_PTR_INT__d_s_l_s += lastoffset*270;
246 PRIM_PTR_INT__d_s_m_s += lastoffset*330;
247 PRIM_PTR_INT__d_s_n_s += lastoffset*396;
248 PRIM_PTR_INT__d_s_o_s += lastoffset*468;
249 PRIM_PTR_INT__d_s_q_s += lastoffset*546;
250 continue;
251 }
252 }
253
254 const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
255 const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
256 const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
257 const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
258
259
260 /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
261 SIMINT_DBLTYPE PQ[3];
262 PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
263 PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
264 PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
265 SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
266 R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
267 R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
268
269 const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
270 const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
271 const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
272 const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
273 const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
274 const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
275 const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
276
277 // NOTE: Minus sign!
278 const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
279 SIMINT_DBLTYPE aop_PQ[3];
280 aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
281 aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
282 aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
283
284 SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
285 SIMINT_DBLTYPE aoq_PQ[3];
286 aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
287 aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
288 aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
289 // Put a minus sign here so we don't have to in RR routines
290 a_over_q = SIMINT_NEG(a_over_q);
291
292
293 //////////////////////////////////////////////
294 // Fjt function section
295 // Maximum v value: 14
296 //////////////////////////////////////////////
297 // The parameter to the Fjt function
298 const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
299
300
301 const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
302
303
304 boys_F_split(PRIM_INT__s_s_s_s, F_x, 14);
305 SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
306 prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
307 for(n = 0; n <= 14; n++)
308 PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
309
310 //////////////////////////////////////////////
311 // Primitive integrals: Vertical recurrance
312 //////////////////////////////////////////////
313
314 const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
315 const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
316 const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
317 const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
318 const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
319 const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
320 const SIMINT_DBLTYPE vrr_const_6_over_2q = SIMINT_MUL(const_6, one_over_2q);
321 const SIMINT_DBLTYPE vrr_const_7_over_2q = SIMINT_MUL(const_7, one_over_2q);
322 const SIMINT_DBLTYPE vrr_const_8_over_2q = SIMINT_MUL(const_8, one_over_2q);
323 const SIMINT_DBLTYPE vrr_const_9_over_2q = SIMINT_MUL(const_9, one_over_2q);
324 const SIMINT_DBLTYPE vrr_const_10_over_2q = SIMINT_MUL(const_10, one_over_2q);
325 const SIMINT_DBLTYPE vrr_const_11_over_2q = SIMINT_MUL(const_11, one_over_2q);
326 const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
327 const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
328 const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
329 const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
330 const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
331 const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
332 const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
333 const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
334 const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
335 const SIMINT_DBLTYPE vrr_const_10_over_2pq = SIMINT_MUL(const_10, one_over_2pq);
336 const SIMINT_DBLTYPE vrr_const_11_over_2pq = SIMINT_MUL(const_11, one_over_2pq);
337 const SIMINT_DBLTYPE vrr_const_12_over_2pq = SIMINT_MUL(const_12, one_over_2pq);
338
339
340
341 // Forming PRIM_INT__s_s_p_s[14 * 3];
342 for(n = 0; n < 14; ++n) // loop over orders of auxiliary function
343 {
344
345 PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
346 PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]);
347
348 PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
349 PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 1]);
350
351 PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
352 PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 2]);
353
354 }
355
356
357
358 // Forming PRIM_INT__s_s_d_s[13 * 6];
359 for(n = 0; n < 13; ++n) // loop over orders of auxiliary function
360 {
361
362 PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
363 PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 0]);
364 PRIM_INT__s_s_d_s[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_d_s[n * 6 + 0]);
365
366 PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
367 PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 3]);
368 PRIM_INT__s_s_d_s[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_d_s[n * 6 + 3]);
369
370 PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
371 PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_d_s[n * 6 + 5]);
372 PRIM_INT__s_s_d_s[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_d_s[n * 6 + 5]);
373
374 }
375
376
377
378 // Forming PRIM_INT__s_s_f_s[12 * 10];
379 for(n = 0; n < 12; ++n) // loop over orders of auxiliary function
380 {
381
382 PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
383 PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 0]);
384 PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]), PRIM_INT__s_s_f_s[n * 10 + 0]);
385
386 PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
387 PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 1]);
388
389 PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
390 PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 2]);
391
392 PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
393 PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 3]);
394
395 PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
396 PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 5]);
397
398 PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
399 PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 6]);
400 PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_p_s[n * 3 + 1]), PRIM_INT__s_s_f_s[n * 10 + 6]);
401
402 PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
403 PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 7]);
404
405 PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
406 PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 9]);
407 PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_p_s[n * 3 + 2]), PRIM_INT__s_s_f_s[n * 10 + 9]);
408
409 }
410
411
412 VRR_K_s_s_g_s(
413 PRIM_INT__s_s_g_s,
414 PRIM_INT__s_s_f_s,
415 PRIM_INT__s_s_d_s,
416 Q_PA,
417 a_over_q,
418 aoq_PQ,
419 one_over_2q,
420 11);
421
422
423 VRR_K_s_s_h_s(
424 PRIM_INT__s_s_h_s,
425 PRIM_INT__s_s_g_s,
426 PRIM_INT__s_s_f_s,
427 Q_PA,
428 a_over_q,
429 aoq_PQ,
430 one_over_2q,
431 10);
432
433
434 ostei_general_vrr1_K(6, 9,
435 one_over_2q, a_over_q, aoq_PQ, Q_PA,
436 PRIM_INT__s_s_h_s, PRIM_INT__s_s_g_s, PRIM_INT__s_s_i_s);
437
438
439 ostei_general_vrr_I(1, 0, 6, 0, 2,
440 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
441 PRIM_INT__s_s_i_s, NULL, NULL, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_i_s);
442
443
444 ostei_general_vrr1_K(7, 8,
445 one_over_2q, a_over_q, aoq_PQ, Q_PA,
446 PRIM_INT__s_s_i_s, PRIM_INT__s_s_h_s, PRIM_INT__s_s_k_s);
447
448
449 ostei_general_vrr_I(1, 0, 7, 0, 2,
450 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
451 PRIM_INT__s_s_k_s, NULL, NULL, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_k_s);
452
453
454 ostei_general_vrr_I(1, 0, 5, 0, 2,
455 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
456 PRIM_INT__s_s_h_s, NULL, NULL, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_h_s);
457
458
459 ostei_general_vrr_I(2, 0, 6, 0, 1,
460 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
461 PRIM_INT__p_s_i_s, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_i_s);
462
463
464 ostei_general_vrr1_K(8, 7,
465 one_over_2q, a_over_q, aoq_PQ, Q_PA,
466 PRIM_INT__s_s_k_s, PRIM_INT__s_s_i_s, PRIM_INT__s_s_l_s);
467
468
469 ostei_general_vrr_I(1, 0, 8, 0, 2,
470 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
471 PRIM_INT__s_s_l_s, NULL, NULL, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_l_s);
472
473
474 ostei_general_vrr_I(2, 0, 7, 0, 1,
475 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
476 PRIM_INT__p_s_k_s, PRIM_INT__s_s_k_s, NULL, PRIM_INT__p_s_i_s, NULL, PRIM_INT__d_s_k_s);
477
478
479 ostei_general_vrr1_K(9, 6,
480 one_over_2q, a_over_q, aoq_PQ, Q_PA,
481 PRIM_INT__s_s_l_s, PRIM_INT__s_s_k_s, PRIM_INT__s_s_m_s);
482
483
484 ostei_general_vrr_I(1, 0, 9, 0, 2,
485 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
486 PRIM_INT__s_s_m_s, NULL, NULL, PRIM_INT__s_s_l_s, NULL, PRIM_INT__p_s_m_s);
487
488
489 ostei_general_vrr_I(2, 0, 8, 0, 1,
490 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
491 PRIM_INT__p_s_l_s, PRIM_INT__s_s_l_s, NULL, PRIM_INT__p_s_k_s, NULL, PRIM_INT__d_s_l_s);
492
493
494 ostei_general_vrr1_K(10, 5,
495 one_over_2q, a_over_q, aoq_PQ, Q_PA,
496 PRIM_INT__s_s_m_s, PRIM_INT__s_s_l_s, PRIM_INT__s_s_n_s);
497
498
499 ostei_general_vrr_I(1, 0, 10, 0, 2,
500 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
501 PRIM_INT__s_s_n_s, NULL, NULL, PRIM_INT__s_s_m_s, NULL, PRIM_INT__p_s_n_s);
502
503
504 ostei_general_vrr_I(2, 0, 9, 0, 1,
505 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
506 PRIM_INT__p_s_m_s, PRIM_INT__s_s_m_s, NULL, PRIM_INT__p_s_l_s, NULL, PRIM_INT__d_s_m_s);
507
508
509 ostei_general_vrr1_K(11, 4,
510 one_over_2q, a_over_q, aoq_PQ, Q_PA,
511 PRIM_INT__s_s_n_s, PRIM_INT__s_s_m_s, PRIM_INT__s_s_o_s);
512
513
514 ostei_general_vrr_I(1, 0, 11, 0, 2,
515 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
516 PRIM_INT__s_s_o_s, NULL, NULL, PRIM_INT__s_s_n_s, NULL, PRIM_INT__p_s_o_s);
517
518
519 ostei_general_vrr_I(2, 0, 10, 0, 1,
520 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
521 PRIM_INT__p_s_n_s, PRIM_INT__s_s_n_s, NULL, PRIM_INT__p_s_m_s, NULL, PRIM_INT__d_s_n_s);
522
523
524 ostei_general_vrr1_K(12, 3,
525 one_over_2q, a_over_q, aoq_PQ, Q_PA,
526 PRIM_INT__s_s_o_s, PRIM_INT__s_s_n_s, PRIM_INT__s_s_q_s);
527
528
529 ostei_general_vrr_I(1, 0, 12, 0, 2,
530 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
531 PRIM_INT__s_s_q_s, NULL, NULL, PRIM_INT__s_s_o_s, NULL, PRIM_INT__p_s_q_s);
532
533
534 ostei_general_vrr_I(2, 0, 11, 0, 1,
535 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
536 PRIM_INT__p_s_o_s, PRIM_INT__s_s_o_s, NULL, PRIM_INT__p_s_n_s, NULL, PRIM_INT__d_s_o_s);
537
538
539 ostei_general_vrr_I(2, 0, 12, 0, 1,
540 one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
541 PRIM_INT__p_s_q_s, PRIM_INT__s_s_q_s, NULL, PRIM_INT__p_s_o_s, NULL, PRIM_INT__d_s_q_s);
542
543
544
545
546 ////////////////////////////////////
547 // Accumulate contracted integrals
548 ////////////////////////////////////
549 if(lastoffset == 0)
550 {
551 contract_all(84, PRIM_INT__p_s_i_s, PRIM_PTR_INT__p_s_i_s);
552 contract_all(108, PRIM_INT__p_s_k_s, PRIM_PTR_INT__p_s_k_s);
553 contract_all(135, PRIM_INT__p_s_l_s, PRIM_PTR_INT__p_s_l_s);
554 contract_all(165, PRIM_INT__p_s_m_s, PRIM_PTR_INT__p_s_m_s);
555 contract_all(198, PRIM_INT__p_s_n_s, PRIM_PTR_INT__p_s_n_s);
556 contract_all(234, PRIM_INT__p_s_o_s, PRIM_PTR_INT__p_s_o_s);
557 contract_all(273, PRIM_INT__p_s_q_s, PRIM_PTR_INT__p_s_q_s);
558 contract_all(168, PRIM_INT__d_s_i_s, PRIM_PTR_INT__d_s_i_s);
559 contract_all(216, PRIM_INT__d_s_k_s, PRIM_PTR_INT__d_s_k_s);
560 contract_all(270, PRIM_INT__d_s_l_s, PRIM_PTR_INT__d_s_l_s);
561 contract_all(330, PRIM_INT__d_s_m_s, PRIM_PTR_INT__d_s_m_s);
562 contract_all(396, PRIM_INT__d_s_n_s, PRIM_PTR_INT__d_s_n_s);
563 contract_all(468, PRIM_INT__d_s_o_s, PRIM_PTR_INT__d_s_o_s);
564 contract_all(546, PRIM_INT__d_s_q_s, PRIM_PTR_INT__d_s_q_s);
565 }
566 else
567 {
568 contract(84, shelloffsets, PRIM_INT__p_s_i_s, PRIM_PTR_INT__p_s_i_s);
569 contract(108, shelloffsets, PRIM_INT__p_s_k_s, PRIM_PTR_INT__p_s_k_s);
570 contract(135, shelloffsets, PRIM_INT__p_s_l_s, PRIM_PTR_INT__p_s_l_s);
571 contract(165, shelloffsets, PRIM_INT__p_s_m_s, PRIM_PTR_INT__p_s_m_s);
572 contract(198, shelloffsets, PRIM_INT__p_s_n_s, PRIM_PTR_INT__p_s_n_s);
573 contract(234, shelloffsets, PRIM_INT__p_s_o_s, PRIM_PTR_INT__p_s_o_s);
574 contract(273, shelloffsets, PRIM_INT__p_s_q_s, PRIM_PTR_INT__p_s_q_s);
575 contract(168, shelloffsets, PRIM_INT__d_s_i_s, PRIM_PTR_INT__d_s_i_s);
576 contract(216, shelloffsets, PRIM_INT__d_s_k_s, PRIM_PTR_INT__d_s_k_s);
577 contract(270, shelloffsets, PRIM_INT__d_s_l_s, PRIM_PTR_INT__d_s_l_s);
578 contract(330, shelloffsets, PRIM_INT__d_s_m_s, PRIM_PTR_INT__d_s_m_s);
579 contract(396, shelloffsets, PRIM_INT__d_s_n_s, PRIM_PTR_INT__d_s_n_s);
580 contract(468, shelloffsets, PRIM_INT__d_s_o_s, PRIM_PTR_INT__d_s_o_s);
581 contract(546, shelloffsets, PRIM_INT__d_s_q_s, PRIM_PTR_INT__d_s_q_s);
582 PRIM_PTR_INT__p_s_i_s += lastoffset*84;
583 PRIM_PTR_INT__p_s_k_s += lastoffset*108;
584 PRIM_PTR_INT__p_s_l_s += lastoffset*135;
585 PRIM_PTR_INT__p_s_m_s += lastoffset*165;
586 PRIM_PTR_INT__p_s_n_s += lastoffset*198;
587 PRIM_PTR_INT__p_s_o_s += lastoffset*234;
588 PRIM_PTR_INT__p_s_q_s += lastoffset*273;
589 PRIM_PTR_INT__d_s_i_s += lastoffset*168;
590 PRIM_PTR_INT__d_s_k_s += lastoffset*216;
591 PRIM_PTR_INT__d_s_l_s += lastoffset*270;
592 PRIM_PTR_INT__d_s_m_s += lastoffset*330;
593 PRIM_PTR_INT__d_s_n_s += lastoffset*396;
594 PRIM_PTR_INT__d_s_o_s += lastoffset*468;
595 PRIM_PTR_INT__d_s_q_s += lastoffset*546;
596 }
597
598 } // close loop over j
599 } // close loop over i
600
601 //Advance to the next batch
602 jstart = SIMINT_SIMD_ROUND(jend);
603
604 //////////////////////////////////////////////
605 // Contracted integrals: Horizontal recurrance
606 //////////////////////////////////////////////
607
608
609 const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
610
611
612 for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
613 {
614 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
615
616 // set up HRR pointers
617 double const * restrict HRR_INT__p_s_i_s = INT__p_s_i_s + abcd * 84;
618 double const * restrict HRR_INT__p_s_k_s = INT__p_s_k_s + abcd * 108;
619 double const * restrict HRR_INT__p_s_l_s = INT__p_s_l_s + abcd * 135;
620 double const * restrict HRR_INT__p_s_m_s = INT__p_s_m_s + abcd * 165;
621 double const * restrict HRR_INT__p_s_n_s = INT__p_s_n_s + abcd * 198;
622 double const * restrict HRR_INT__p_s_o_s = INT__p_s_o_s + abcd * 234;
623 double const * restrict HRR_INT__p_s_q_s = INT__p_s_q_s + abcd * 273;
624 double const * restrict HRR_INT__d_s_i_s = INT__d_s_i_s + abcd * 168;
625 double const * restrict HRR_INT__d_s_k_s = INT__d_s_k_s + abcd * 216;
626 double const * restrict HRR_INT__d_s_l_s = INT__d_s_l_s + abcd * 270;
627 double const * restrict HRR_INT__d_s_m_s = INT__d_s_m_s + abcd * 330;
628 double const * restrict HRR_INT__d_s_n_s = INT__d_s_n_s + abcd * 396;
629 double const * restrict HRR_INT__d_s_o_s = INT__d_s_o_s + abcd * 468;
630 double const * restrict HRR_INT__d_s_q_s = INT__d_s_q_s + abcd * 546;
631 double * restrict HRR_INT__p_p_i_i = INT__p_p_i_i + real_abcd * 7056;
632
633 // form INT__p_p_i_s
634 for(iket = 0; iket < 28; ++iket)
635 {
636 HRR_INT__p_p_i_s[0 * 28 + iket] = HRR_INT__d_s_i_s[0 * 28 + iket] + ( hAB[0] * HRR_INT__p_s_i_s[0 * 28 + iket] );
637
638 HRR_INT__p_p_i_s[1 * 28 + iket] = HRR_INT__d_s_i_s[1 * 28 + iket] + ( hAB[1] * HRR_INT__p_s_i_s[0 * 28 + iket] );
639
640 HRR_INT__p_p_i_s[2 * 28 + iket] = HRR_INT__d_s_i_s[2 * 28 + iket] + ( hAB[2] * HRR_INT__p_s_i_s[0 * 28 + iket] );
641
642 HRR_INT__p_p_i_s[3 * 28 + iket] = HRR_INT__d_s_i_s[1 * 28 + iket] + ( hAB[0] * HRR_INT__p_s_i_s[1 * 28 + iket] );
643
644 HRR_INT__p_p_i_s[4 * 28 + iket] = HRR_INT__d_s_i_s[3 * 28 + iket] + ( hAB[1] * HRR_INT__p_s_i_s[1 * 28 + iket] );
645
646 HRR_INT__p_p_i_s[5 * 28 + iket] = HRR_INT__d_s_i_s[4 * 28 + iket] + ( hAB[2] * HRR_INT__p_s_i_s[1 * 28 + iket] );
647
648 HRR_INT__p_p_i_s[6 * 28 + iket] = HRR_INT__d_s_i_s[2 * 28 + iket] + ( hAB[0] * HRR_INT__p_s_i_s[2 * 28 + iket] );
649
650 HRR_INT__p_p_i_s[7 * 28 + iket] = HRR_INT__d_s_i_s[4 * 28 + iket] + ( hAB[1] * HRR_INT__p_s_i_s[2 * 28 + iket] );
651
652 HRR_INT__p_p_i_s[8 * 28 + iket] = HRR_INT__d_s_i_s[5 * 28 + iket] + ( hAB[2] * HRR_INT__p_s_i_s[2 * 28 + iket] );
653
654 }
655
656
657 // form INT__p_p_k_s
658 for(iket = 0; iket < 36; ++iket)
659 {
660 HRR_INT__p_p_k_s[0 * 36 + iket] = HRR_INT__d_s_k_s[0 * 36 + iket] + ( hAB[0] * HRR_INT__p_s_k_s[0 * 36 + iket] );
661
662 HRR_INT__p_p_k_s[1 * 36 + iket] = HRR_INT__d_s_k_s[1 * 36 + iket] + ( hAB[1] * HRR_INT__p_s_k_s[0 * 36 + iket] );
663
664 HRR_INT__p_p_k_s[2 * 36 + iket] = HRR_INT__d_s_k_s[2 * 36 + iket] + ( hAB[2] * HRR_INT__p_s_k_s[0 * 36 + iket] );
665
666 HRR_INT__p_p_k_s[3 * 36 + iket] = HRR_INT__d_s_k_s[1 * 36 + iket] + ( hAB[0] * HRR_INT__p_s_k_s[1 * 36 + iket] );
667
668 HRR_INT__p_p_k_s[4 * 36 + iket] = HRR_INT__d_s_k_s[3 * 36 + iket] + ( hAB[1] * HRR_INT__p_s_k_s[1 * 36 + iket] );
669
670 HRR_INT__p_p_k_s[5 * 36 + iket] = HRR_INT__d_s_k_s[4 * 36 + iket] + ( hAB[2] * HRR_INT__p_s_k_s[1 * 36 + iket] );
671
672 HRR_INT__p_p_k_s[6 * 36 + iket] = HRR_INT__d_s_k_s[2 * 36 + iket] + ( hAB[0] * HRR_INT__p_s_k_s[2 * 36 + iket] );
673
674 HRR_INT__p_p_k_s[7 * 36 + iket] = HRR_INT__d_s_k_s[4 * 36 + iket] + ( hAB[1] * HRR_INT__p_s_k_s[2 * 36 + iket] );
675
676 HRR_INT__p_p_k_s[8 * 36 + iket] = HRR_INT__d_s_k_s[5 * 36 + iket] + ( hAB[2] * HRR_INT__p_s_k_s[2 * 36 + iket] );
677
678 }
679
680
681 // form INT__p_p_l_s
682 for(iket = 0; iket < 45; ++iket)
683 {
684 HRR_INT__p_p_l_s[0 * 45 + iket] = HRR_INT__d_s_l_s[0 * 45 + iket] + ( hAB[0] * HRR_INT__p_s_l_s[0 * 45 + iket] );
685
686 HRR_INT__p_p_l_s[1 * 45 + iket] = HRR_INT__d_s_l_s[1 * 45 + iket] + ( hAB[1] * HRR_INT__p_s_l_s[0 * 45 + iket] );
687
688 HRR_INT__p_p_l_s[2 * 45 + iket] = HRR_INT__d_s_l_s[2 * 45 + iket] + ( hAB[2] * HRR_INT__p_s_l_s[0 * 45 + iket] );
689
690 HRR_INT__p_p_l_s[3 * 45 + iket] = HRR_INT__d_s_l_s[1 * 45 + iket] + ( hAB[0] * HRR_INT__p_s_l_s[1 * 45 + iket] );
691
692 HRR_INT__p_p_l_s[4 * 45 + iket] = HRR_INT__d_s_l_s[3 * 45 + iket] + ( hAB[1] * HRR_INT__p_s_l_s[1 * 45 + iket] );
693
694 HRR_INT__p_p_l_s[5 * 45 + iket] = HRR_INT__d_s_l_s[4 * 45 + iket] + ( hAB[2] * HRR_INT__p_s_l_s[1 * 45 + iket] );
695
696 HRR_INT__p_p_l_s[6 * 45 + iket] = HRR_INT__d_s_l_s[2 * 45 + iket] + ( hAB[0] * HRR_INT__p_s_l_s[2 * 45 + iket] );
697
698 HRR_INT__p_p_l_s[7 * 45 + iket] = HRR_INT__d_s_l_s[4 * 45 + iket] + ( hAB[1] * HRR_INT__p_s_l_s[2 * 45 + iket] );
699
700 HRR_INT__p_p_l_s[8 * 45 + iket] = HRR_INT__d_s_l_s[5 * 45 + iket] + ( hAB[2] * HRR_INT__p_s_l_s[2 * 45 + iket] );
701
702 }
703
704
705 // form INT__p_p_m_s
706 for(iket = 0; iket < 55; ++iket)
707 {
708 HRR_INT__p_p_m_s[0 * 55 + iket] = HRR_INT__d_s_m_s[0 * 55 + iket] + ( hAB[0] * HRR_INT__p_s_m_s[0 * 55 + iket] );
709
710 HRR_INT__p_p_m_s[1 * 55 + iket] = HRR_INT__d_s_m_s[1 * 55 + iket] + ( hAB[1] * HRR_INT__p_s_m_s[0 * 55 + iket] );
711
712 HRR_INT__p_p_m_s[2 * 55 + iket] = HRR_INT__d_s_m_s[2 * 55 + iket] + ( hAB[2] * HRR_INT__p_s_m_s[0 * 55 + iket] );
713
714 HRR_INT__p_p_m_s[3 * 55 + iket] = HRR_INT__d_s_m_s[1 * 55 + iket] + ( hAB[0] * HRR_INT__p_s_m_s[1 * 55 + iket] );
715
716 HRR_INT__p_p_m_s[4 * 55 + iket] = HRR_INT__d_s_m_s[3 * 55 + iket] + ( hAB[1] * HRR_INT__p_s_m_s[1 * 55 + iket] );
717
718 HRR_INT__p_p_m_s[5 * 55 + iket] = HRR_INT__d_s_m_s[4 * 55 + iket] + ( hAB[2] * HRR_INT__p_s_m_s[1 * 55 + iket] );
719
720 HRR_INT__p_p_m_s[6 * 55 + iket] = HRR_INT__d_s_m_s[2 * 55 + iket] + ( hAB[0] * HRR_INT__p_s_m_s[2 * 55 + iket] );
721
722 HRR_INT__p_p_m_s[7 * 55 + iket] = HRR_INT__d_s_m_s[4 * 55 + iket] + ( hAB[1] * HRR_INT__p_s_m_s[2 * 55 + iket] );
723
724 HRR_INT__p_p_m_s[8 * 55 + iket] = HRR_INT__d_s_m_s[5 * 55 + iket] + ( hAB[2] * HRR_INT__p_s_m_s[2 * 55 + iket] );
725
726 }
727
728
729 // form INT__p_p_n_s
730 for(iket = 0; iket < 66; ++iket)
731 {
732 HRR_INT__p_p_n_s[0 * 66 + iket] = HRR_INT__d_s_n_s[0 * 66 + iket] + ( hAB[0] * HRR_INT__p_s_n_s[0 * 66 + iket] );
733
734 HRR_INT__p_p_n_s[1 * 66 + iket] = HRR_INT__d_s_n_s[1 * 66 + iket] + ( hAB[1] * HRR_INT__p_s_n_s[0 * 66 + iket] );
735
736 HRR_INT__p_p_n_s[2 * 66 + iket] = HRR_INT__d_s_n_s[2 * 66 + iket] + ( hAB[2] * HRR_INT__p_s_n_s[0 * 66 + iket] );
737
738 HRR_INT__p_p_n_s[3 * 66 + iket] = HRR_INT__d_s_n_s[1 * 66 + iket] + ( hAB[0] * HRR_INT__p_s_n_s[1 * 66 + iket] );
739
740 HRR_INT__p_p_n_s[4 * 66 + iket] = HRR_INT__d_s_n_s[3 * 66 + iket] + ( hAB[1] * HRR_INT__p_s_n_s[1 * 66 + iket] );
741
742 HRR_INT__p_p_n_s[5 * 66 + iket] = HRR_INT__d_s_n_s[4 * 66 + iket] + ( hAB[2] * HRR_INT__p_s_n_s[1 * 66 + iket] );
743
744 HRR_INT__p_p_n_s[6 * 66 + iket] = HRR_INT__d_s_n_s[2 * 66 + iket] + ( hAB[0] * HRR_INT__p_s_n_s[2 * 66 + iket] );
745
746 HRR_INT__p_p_n_s[7 * 66 + iket] = HRR_INT__d_s_n_s[4 * 66 + iket] + ( hAB[1] * HRR_INT__p_s_n_s[2 * 66 + iket] );
747
748 HRR_INT__p_p_n_s[8 * 66 + iket] = HRR_INT__d_s_n_s[5 * 66 + iket] + ( hAB[2] * HRR_INT__p_s_n_s[2 * 66 + iket] );
749
750 }
751
752
753 // form INT__p_p_o_s
754 for(iket = 0; iket < 78; ++iket)
755 {
756 HRR_INT__p_p_o_s[0 * 78 + iket] = HRR_INT__d_s_o_s[0 * 78 + iket] + ( hAB[0] * HRR_INT__p_s_o_s[0 * 78 + iket] );
757
758 HRR_INT__p_p_o_s[1 * 78 + iket] = HRR_INT__d_s_o_s[1 * 78 + iket] + ( hAB[1] * HRR_INT__p_s_o_s[0 * 78 + iket] );
759
760 HRR_INT__p_p_o_s[2 * 78 + iket] = HRR_INT__d_s_o_s[2 * 78 + iket] + ( hAB[2] * HRR_INT__p_s_o_s[0 * 78 + iket] );
761
762 HRR_INT__p_p_o_s[3 * 78 + iket] = HRR_INT__d_s_o_s[1 * 78 + iket] + ( hAB[0] * HRR_INT__p_s_o_s[1 * 78 + iket] );
763
764 HRR_INT__p_p_o_s[4 * 78 + iket] = HRR_INT__d_s_o_s[3 * 78 + iket] + ( hAB[1] * HRR_INT__p_s_o_s[1 * 78 + iket] );
765
766 HRR_INT__p_p_o_s[5 * 78 + iket] = HRR_INT__d_s_o_s[4 * 78 + iket] + ( hAB[2] * HRR_INT__p_s_o_s[1 * 78 + iket] );
767
768 HRR_INT__p_p_o_s[6 * 78 + iket] = HRR_INT__d_s_o_s[2 * 78 + iket] + ( hAB[0] * HRR_INT__p_s_o_s[2 * 78 + iket] );
769
770 HRR_INT__p_p_o_s[7 * 78 + iket] = HRR_INT__d_s_o_s[4 * 78 + iket] + ( hAB[1] * HRR_INT__p_s_o_s[2 * 78 + iket] );
771
772 HRR_INT__p_p_o_s[8 * 78 + iket] = HRR_INT__d_s_o_s[5 * 78 + iket] + ( hAB[2] * HRR_INT__p_s_o_s[2 * 78 + iket] );
773
774 }
775
776
777 // form INT__p_p_q_s
778 for(iket = 0; iket < 91; ++iket)
779 {
780 HRR_INT__p_p_q_s[0 * 91 + iket] = HRR_INT__d_s_q_s[0 * 91 + iket] + ( hAB[0] * HRR_INT__p_s_q_s[0 * 91 + iket] );
781
782 HRR_INT__p_p_q_s[1 * 91 + iket] = HRR_INT__d_s_q_s[1 * 91 + iket] + ( hAB[1] * HRR_INT__p_s_q_s[0 * 91 + iket] );
783
784 HRR_INT__p_p_q_s[2 * 91 + iket] = HRR_INT__d_s_q_s[2 * 91 + iket] + ( hAB[2] * HRR_INT__p_s_q_s[0 * 91 + iket] );
785
786 HRR_INT__p_p_q_s[3 * 91 + iket] = HRR_INT__d_s_q_s[1 * 91 + iket] + ( hAB[0] * HRR_INT__p_s_q_s[1 * 91 + iket] );
787
788 HRR_INT__p_p_q_s[4 * 91 + iket] = HRR_INT__d_s_q_s[3 * 91 + iket] + ( hAB[1] * HRR_INT__p_s_q_s[1 * 91 + iket] );
789
790 HRR_INT__p_p_q_s[5 * 91 + iket] = HRR_INT__d_s_q_s[4 * 91 + iket] + ( hAB[2] * HRR_INT__p_s_q_s[1 * 91 + iket] );
791
792 HRR_INT__p_p_q_s[6 * 91 + iket] = HRR_INT__d_s_q_s[2 * 91 + iket] + ( hAB[0] * HRR_INT__p_s_q_s[2 * 91 + iket] );
793
794 HRR_INT__p_p_q_s[7 * 91 + iket] = HRR_INT__d_s_q_s[4 * 91 + iket] + ( hAB[1] * HRR_INT__p_s_q_s[2 * 91 + iket] );
795
796 HRR_INT__p_p_q_s[8 * 91 + iket] = HRR_INT__d_s_q_s[5 * 91 + iket] + ( hAB[2] * HRR_INT__p_s_q_s[2 * 91 + iket] );
797
798 }
799
800
801 // form INT__p_p_i_p
802 ostei_general_hrr_L(1, 1, 6, 1, hCD, HRR_INT__p_p_k_s, HRR_INT__p_p_i_s, HRR_INT__p_p_i_p);
803
804 // form INT__p_p_k_p
805 ostei_general_hrr_L(1, 1, 7, 1, hCD, HRR_INT__p_p_l_s, HRR_INT__p_p_k_s, HRR_INT__p_p_k_p);
806
807 // form INT__p_p_l_p
808 ostei_general_hrr_L(1, 1, 8, 1, hCD, HRR_INT__p_p_m_s, HRR_INT__p_p_l_s, HRR_INT__p_p_l_p);
809
810 // form INT__p_p_m_p
811 ostei_general_hrr_L(1, 1, 9, 1, hCD, HRR_INT__p_p_n_s, HRR_INT__p_p_m_s, HRR_INT__p_p_m_p);
812
813 // form INT__p_p_n_p
814 ostei_general_hrr_L(1, 1, 10, 1, hCD, HRR_INT__p_p_o_s, HRR_INT__p_p_n_s, HRR_INT__p_p_n_p);
815
816 // form INT__p_p_o_p
817 ostei_general_hrr_L(1, 1, 11, 1, hCD, HRR_INT__p_p_q_s, HRR_INT__p_p_o_s, HRR_INT__p_p_o_p);
818
819 // form INT__p_p_i_d
820 ostei_general_hrr_L(1, 1, 6, 2, hCD, HRR_INT__p_p_k_p, HRR_INT__p_p_i_p, HRR_INT__p_p_i_d);
821
822 // form INT__p_p_k_d
823 ostei_general_hrr_L(1, 1, 7, 2, hCD, HRR_INT__p_p_l_p, HRR_INT__p_p_k_p, HRR_INT__p_p_k_d);
824
825 // form INT__p_p_l_d
826 ostei_general_hrr_L(1, 1, 8, 2, hCD, HRR_INT__p_p_m_p, HRR_INT__p_p_l_p, HRR_INT__p_p_l_d);
827
828 // form INT__p_p_m_d
829 ostei_general_hrr_L(1, 1, 9, 2, hCD, HRR_INT__p_p_n_p, HRR_INT__p_p_m_p, HRR_INT__p_p_m_d);
830
831 // form INT__p_p_n_d
832 ostei_general_hrr_L(1, 1, 10, 2, hCD, HRR_INT__p_p_o_p, HRR_INT__p_p_n_p, HRR_INT__p_p_n_d);
833
834 // form INT__p_p_i_f
835 ostei_general_hrr_L(1, 1, 6, 3, hCD, HRR_INT__p_p_k_d, HRR_INT__p_p_i_d, HRR_INT__p_p_i_f);
836
837 // form INT__p_p_k_f
838 ostei_general_hrr_L(1, 1, 7, 3, hCD, HRR_INT__p_p_l_d, HRR_INT__p_p_k_d, HRR_INT__p_p_k_f);
839
840 // form INT__p_p_l_f
841 ostei_general_hrr_L(1, 1, 8, 3, hCD, HRR_INT__p_p_m_d, HRR_INT__p_p_l_d, HRR_INT__p_p_l_f);
842
843 // form INT__p_p_m_f
844 ostei_general_hrr_L(1, 1, 9, 3, hCD, HRR_INT__p_p_n_d, HRR_INT__p_p_m_d, HRR_INT__p_p_m_f);
845
846 // form INT__p_p_i_g
847 ostei_general_hrr_L(1, 1, 6, 4, hCD, HRR_INT__p_p_k_f, HRR_INT__p_p_i_f, HRR_INT__p_p_i_g);
848
849 // form INT__p_p_k_g
850 ostei_general_hrr_L(1, 1, 7, 4, hCD, HRR_INT__p_p_l_f, HRR_INT__p_p_k_f, HRR_INT__p_p_k_g);
851
852 // form INT__p_p_l_g
853 ostei_general_hrr_L(1, 1, 8, 4, hCD, HRR_INT__p_p_m_f, HRR_INT__p_p_l_f, HRR_INT__p_p_l_g);
854
855 // form INT__p_p_i_h
856 ostei_general_hrr_L(1, 1, 6, 5, hCD, HRR_INT__p_p_k_g, HRR_INT__p_p_i_g, HRR_INT__p_p_i_h);
857
858 // form INT__p_p_k_h
859 ostei_general_hrr_L(1, 1, 7, 5, hCD, HRR_INT__p_p_l_g, HRR_INT__p_p_k_g, HRR_INT__p_p_k_h);
860
861 // form INT__p_p_i_i
862 ostei_general_hrr_L(1, 1, 6, 6, hCD, HRR_INT__p_p_k_h, HRR_INT__p_p_i_h, HRR_INT__p_p_i_i);
863
864
865 } // close HRR loop
866
867
868 } // close loop cdbatch
869
870 istart = iend;
871 } // close loop over ab
872
873 return P.nshell12_clip * Q.nshell12_clip;
874 }
875
876