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_g_p_i_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_p_i_p)8 int ostei_g_p_i_p(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__g_p_i_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__g_p_i_p);
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__g_s_i_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__g_s_k_s = work + (SIMINT_NSHELL_SIMD * 420);
31     double * const INT__h_s_i_s = work + (SIMINT_NSHELL_SIMD * 960);
32     double * const INT__h_s_k_s = work + (SIMINT_NSHELL_SIMD * 1548);
33     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*2304);
34     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
35     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 13;
36     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 49;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_f_s = primwork + 115;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_g_s = primwork + 215;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_h_s = primwork + 350;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_i_s = primwork + 518;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_k_s = primwork + 714;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_d_s = primwork + 930;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_f_s = primwork + 1020;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_g_s = primwork + 1170;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_h_s = primwork + 1395;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_i_s = primwork + 1710;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_k_s = primwork + 2130;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_f_s = primwork + 2670;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_g_s = primwork + 2910;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_h_s = primwork + 3270;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_i_s = primwork + 3774;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_k_s = primwork + 4446;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_g_s = primwork + 5310;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_h_s = primwork + 5760;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_i_s = primwork + 6390;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_k_s = primwork + 7230;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_h_s = primwork + 8310;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_i_s = primwork + 8940;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_k_s = primwork + 9780;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_i_s = primwork + 10860;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_k_s = primwork + 11448;
62     double * const hrrwork = (double *)(primwork + 12204);
63     double * const HRR_INT__g_p_i_s = hrrwork + 0;
64     double * const HRR_INT__g_p_k_s = hrrwork + 1260;
65 
66 
67     // Create constants
68     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
69     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
70     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
71     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
72     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
73     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
74     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
75     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
76 
77 
78     ////////////////////////////////////////
79     // Loop over shells and primitives
80     ////////////////////////////////////////
81 
82     real_abcd = 0;
83     istart = 0;
84     for(ab = 0; ab < P.nshell12_clip; ++ab)
85     {
86         const int iend = istart + P.nprim12[ab];
87 
88         cd = 0;
89         jstart = 0;
90 
91         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
92         {
93             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
94             int jend = jstart;
95             for(i = 0; i < nshellbatch; i++)
96                 jend += Q.nprim12[cd+i];
97 
98             // Clear the beginning of the workspace (where we are accumulating integrals)
99             memset(work, 0, SIMINT_NSHELL_SIMD * 2304 * sizeof(double));
100             abcd = 0;
101 
102 
103             for(i = istart; i < iend; ++i)
104             {
105                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
106 
107                 if(check_screen)
108                 {
109                     // Skip this whole thing if always insignificant
110                     if((P.screen[i] * Q.screen_max) < screen_tol)
111                         continue;
112                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
113                 }
114 
115                 icd = 0;
116                 iprimcd = 0;
117                 nprim_icd = Q.nprim12[cd];
118                 double * restrict PRIM_PTR_INT__g_s_i_s = INT__g_s_i_s + abcd * 420;
119                 double * restrict PRIM_PTR_INT__g_s_k_s = INT__g_s_k_s + abcd * 540;
120                 double * restrict PRIM_PTR_INT__h_s_i_s = INT__h_s_i_s + abcd * 588;
121                 double * restrict PRIM_PTR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
122 
123 
124 
125                 // Load these one per loop over i
126                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
127                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
128                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
129 
130                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
131 
132                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
133                 {
134                     // calculate the shell offsets
135                     // these are the offset from the shell pointed to by cd
136                     // for each element
137                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
138                     int lastoffset = 0;
139                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
140 
141                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
142                     {
143                         // Handle if the first element of the vector is a new shell
144                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
145                         {
146                             nprim_icd += Q.nprim12[cd + (++icd)];
147                             PRIM_PTR_INT__g_s_i_s += 420;
148                             PRIM_PTR_INT__g_s_k_s += 540;
149                             PRIM_PTR_INT__h_s_i_s += 588;
150                             PRIM_PTR_INT__h_s_k_s += 756;
151                         }
152                         iprimcd++;
153                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
154                         {
155                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
156                             {
157                                 shelloffsets[n] = shelloffsets[n-1] + 1;
158                                 lastoffset++;
159                                 nprim_icd += Q.nprim12[cd + (++icd)];
160                             }
161                             else
162                                 shelloffsets[n] = shelloffsets[n-1];
163                             iprimcd++;
164                         }
165                     }
166                     else
167                         iprimcd += SIMINT_SIMD_LEN;
168 
169                     // Do we have to compute this vector (or has it been screened out)?
170                     // (not_screened != 0 means we have to do this vector)
171                     if(check_screen)
172                     {
173                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
174                         if(vmax < screen_tol)
175                         {
176                             PRIM_PTR_INT__g_s_i_s += lastoffset*420;
177                             PRIM_PTR_INT__g_s_k_s += lastoffset*540;
178                             PRIM_PTR_INT__h_s_i_s += lastoffset*588;
179                             PRIM_PTR_INT__h_s_k_s += lastoffset*756;
180                             continue;
181                         }
182                     }
183 
184                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
185                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
186                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
187                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
188 
189 
190                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
191                     SIMINT_DBLTYPE PQ[3];
192                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
193                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
194                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
195                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
196                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
197                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
198 
199                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
200                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
201                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
202                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
203                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
204                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
205                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
206 
207                     // NOTE: Minus sign!
208                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
209                     SIMINT_DBLTYPE aop_PQ[3];
210                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
211                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
212                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
213 
214                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
215                     SIMINT_DBLTYPE aoq_PQ[3];
216                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
217                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
218                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
219                     // Put a minus sign here so we don't have to in RR routines
220                     a_over_q = SIMINT_NEG(a_over_q);
221 
222 
223                     //////////////////////////////////////////////
224                     // Fjt function section
225                     // Maximum v value: 12
226                     //////////////////////////////////////////////
227                     // The parameter to the Fjt function
228                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
229 
230 
231                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
232 
233 
234                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 12);
235                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
236                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
237                     for(n = 0; n <= 12; n++)
238                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
239 
240                     //////////////////////////////////////////////
241                     // Primitive integrals: Vertical recurrance
242                     //////////////////////////////////////////////
243 
244                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
245                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
246                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
247                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
248                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
249                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
250                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
251                     const SIMINT_DBLTYPE vrr_const_4_over_2q = SIMINT_MUL(const_4, one_over_2q);
252                     const SIMINT_DBLTYPE vrr_const_5_over_2q = SIMINT_MUL(const_5, one_over_2q);
253                     const SIMINT_DBLTYPE vrr_const_6_over_2q = SIMINT_MUL(const_6, one_over_2q);
254                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
255                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
256                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
257                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
258                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
259                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
260                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
261 
262 
263 
264                     // Forming PRIM_INT__s_s_p_s[12 * 3];
265                     for(n = 0; n < 12; ++n)  // loop over orders of auxiliary function
266                     {
267 
268                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
269                         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]);
270 
271                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
272                         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]);
273 
274                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
275                         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]);
276 
277                     }
278 
279 
280 
281                     // Forming PRIM_INT__s_s_d_s[11 * 6];
282                     for(n = 0; n < 11; ++n)  // loop over orders of auxiliary function
283                     {
284 
285                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
286                         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]);
287                         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]);
288 
289                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
290                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 1]);
291 
292                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
293                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 2]);
294 
295                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
296                         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]);
297                         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]);
298 
299                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
300                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 4]);
301 
302                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
303                         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]);
304                         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]);
305 
306                     }
307 
308 
309 
310                     // Forming PRIM_INT__s_s_f_s[10 * 10];
311                     for(n = 0; n < 10; ++n)  // loop over orders of auxiliary function
312                     {
313 
314                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
315                         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]);
316                         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]);
317 
318                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
319                         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]);
320 
321                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
322                         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]);
323 
324                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
325                         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]);
326 
327                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
328                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__s_s_f_s[n * 10 + 4]);
329 
330                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
331                         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]);
332 
333                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
334                         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]);
335                         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]);
336 
337                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
338                         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]);
339 
340                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
341                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 8]);
342 
343                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
344                         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]);
345                         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]);
346 
347                     }
348 
349 
350                     VRR_K_s_s_g_s(
351                             PRIM_INT__s_s_g_s,
352                             PRIM_INT__s_s_f_s,
353                             PRIM_INT__s_s_d_s,
354                             Q_PA,
355                             a_over_q,
356                             aoq_PQ,
357                             one_over_2q,
358                             9);
359 
360 
361                     VRR_K_s_s_h_s(
362                             PRIM_INT__s_s_h_s,
363                             PRIM_INT__s_s_g_s,
364                             PRIM_INT__s_s_f_s,
365                             Q_PA,
366                             a_over_q,
367                             aoq_PQ,
368                             one_over_2q,
369                             8);
370 
371 
372                     ostei_general_vrr1_K(6, 7,
373                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
374                             PRIM_INT__s_s_h_s, PRIM_INT__s_s_g_s, PRIM_INT__s_s_i_s);
375 
376 
377                     ostei_general_vrr_I(1, 0, 6, 0, 5,
378                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
379                             PRIM_INT__s_s_i_s, NULL, NULL, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_i_s);
380 
381 
382                     ostei_general_vrr_I(1, 0, 5, 0, 5,
383                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
384                             PRIM_INT__s_s_h_s, NULL, NULL, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_h_s);
385 
386 
387                     ostei_general_vrr_I(2, 0, 6, 0, 4,
388                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
389                             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);
390 
391 
392                     VRR_I_p_s_g_s(
393                             PRIM_INT__p_s_g_s,
394                             PRIM_INT__s_s_g_s,
395                             PRIM_INT__s_s_f_s,
396                             P_PA,
397                             aop_PQ,
398                             one_over_2pq,
399                             5);
400 
401 
402                     ostei_general_vrr_I(2, 0, 5, 0, 4,
403                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
404                             PRIM_INT__p_s_h_s, PRIM_INT__s_s_h_s, NULL, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_h_s);
405 
406 
407                     ostei_general_vrr_I(3, 0, 6, 0, 3,
408                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
409                             PRIM_INT__d_s_i_s, PRIM_INT__p_s_i_s, NULL, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_i_s);
410 
411 
412                     VRR_I_p_s_f_s(
413                             PRIM_INT__p_s_f_s,
414                             PRIM_INT__s_s_f_s,
415                             PRIM_INT__s_s_d_s,
416                             P_PA,
417                             aop_PQ,
418                             one_over_2pq,
419                             5);
420 
421 
422                     ostei_general_vrr_I(2, 0, 4, 0, 4,
423                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
424                             PRIM_INT__p_s_g_s, PRIM_INT__s_s_g_s, NULL, PRIM_INT__p_s_f_s, NULL, PRIM_INT__d_s_g_s);
425 
426 
427                     ostei_general_vrr_I(3, 0, 5, 0, 3,
428                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
429                             PRIM_INT__d_s_h_s, PRIM_INT__p_s_h_s, NULL, PRIM_INT__d_s_g_s, NULL, PRIM_INT__f_s_h_s);
430 
431 
432                     ostei_general_vrr_I(4, 0, 6, 0, 2,
433                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
434                             PRIM_INT__f_s_i_s, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_h_s, NULL, PRIM_INT__g_s_i_s);
435 
436 
437                     ostei_general_vrr1_K(7, 6,
438                             one_over_2q, a_over_q, aoq_PQ, Q_PA,
439                             PRIM_INT__s_s_i_s, PRIM_INT__s_s_h_s, PRIM_INT__s_s_k_s);
440 
441 
442                     ostei_general_vrr_I(1, 0, 7, 0, 5,
443                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
444                             PRIM_INT__s_s_k_s, NULL, NULL, PRIM_INT__s_s_i_s, NULL, PRIM_INT__p_s_k_s);
445 
446 
447                     ostei_general_vrr_I(2, 0, 7, 0, 4,
448                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
449                             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);
450 
451 
452                     ostei_general_vrr_I(3, 0, 7, 0, 3,
453                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
454                             PRIM_INT__d_s_k_s, PRIM_INT__p_s_k_s, NULL, PRIM_INT__d_s_i_s, NULL, PRIM_INT__f_s_k_s);
455 
456 
457                     ostei_general_vrr_I(4, 0, 7, 0, 2,
458                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
459                             PRIM_INT__f_s_k_s, PRIM_INT__d_s_k_s, NULL, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_k_s);
460 
461 
462 
463                     // Forming PRIM_INT__p_s_d_s[5 * 18];
464                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
465                     {
466 
467                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
468                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
469                         PRIM_INT__p_s_d_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 0]);
470 
471                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 1]);
472                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 1]);
473                         PRIM_INT__p_s_d_s[n * 18 + 1] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 1]);
474 
475                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 2]);
476                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 2]);
477                         PRIM_INT__p_s_d_s[n * 18 + 2] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 2]);
478 
479                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
480                         PRIM_INT__p_s_d_s[n * 18 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 3]);
481 
482                         PRIM_INT__p_s_d_s[n * 18 + 4] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 4]);
483                         PRIM_INT__p_s_d_s[n * 18 + 4] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 4]);
484 
485                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
486                         PRIM_INT__p_s_d_s[n * 18 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 5]);
487 
488                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
489                         PRIM_INT__p_s_d_s[n * 18 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 6]);
490 
491                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 1]);
492                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 7]);
493                         PRIM_INT__p_s_d_s[n * 18 + 7] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 7]);
494 
495                         PRIM_INT__p_s_d_s[n * 18 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 2]);
496                         PRIM_INT__p_s_d_s[n * 18 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 8]);
497 
498                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
499                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 9]);
500                         PRIM_INT__p_s_d_s[n * 18 + 9] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 9]);
501 
502                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 4]);
503                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 10]);
504                         PRIM_INT__p_s_d_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 10]);
505 
506                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
507                         PRIM_INT__p_s_d_s[n * 18 + 11] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 11]);
508 
509                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
510                         PRIM_INT__p_s_d_s[n * 18 + 12] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__p_s_d_s[n * 18 + 12]);
511 
512                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
513                         PRIM_INT__p_s_d_s[n * 18 + 13] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__p_s_d_s[n * 18 + 13]);
514 
515                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 2]);
516                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 2], PRIM_INT__p_s_d_s[n * 18 + 14]);
517                         PRIM_INT__p_s_d_s[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__p_s_d_s[n * 18 + 14]);
518 
519                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
520                         PRIM_INT__p_s_d_s[n * 18 + 15] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__p_s_d_s[n * 18 + 15]);
521 
522                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 4]);
523                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 4], PRIM_INT__p_s_d_s[n * 18 + 16]);
524                         PRIM_INT__p_s_d_s[n * 18 + 16] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__p_s_d_s[n * 18 + 16]);
525 
526                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
527                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__p_s_d_s[n * 18 + 17]);
528                         PRIM_INT__p_s_d_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__p_s_d_s[n * 18 + 17]);
529 
530                     }
531 
532 
533                     VRR_I_d_s_f_s(
534                             PRIM_INT__d_s_f_s,
535                             PRIM_INT__p_s_f_s,
536                             PRIM_INT__s_s_f_s,
537                             PRIM_INT__p_s_d_s,
538                             P_PA,
539                             a_over_p,
540                             aop_PQ,
541                             one_over_2p,
542                             one_over_2pq,
543                             4);
544 
545 
546                     ostei_general_vrr_I(3, 0, 4, 0, 3,
547                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
548                             PRIM_INT__d_s_g_s, PRIM_INT__p_s_g_s, NULL, PRIM_INT__d_s_f_s, NULL, PRIM_INT__f_s_g_s);
549 
550 
551                     ostei_general_vrr_I(4, 0, 5, 0, 2,
552                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
553                             PRIM_INT__f_s_h_s, PRIM_INT__d_s_h_s, NULL, PRIM_INT__f_s_g_s, NULL, PRIM_INT__g_s_h_s);
554 
555 
556                     ostei_general_vrr_I(5, 0, 6, 0, 1,
557                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
558                             PRIM_INT__g_s_i_s, PRIM_INT__f_s_i_s, NULL, PRIM_INT__g_s_h_s, NULL, PRIM_INT__h_s_i_s);
559 
560 
561                     ostei_general_vrr_I(5, 0, 7, 0, 1,
562                             one_over_2p, a_over_p, one_over_2pq, aop_PQ, P_PA,
563                             PRIM_INT__g_s_k_s, PRIM_INT__f_s_k_s, NULL, PRIM_INT__g_s_i_s, NULL, PRIM_INT__h_s_k_s);
564 
565 
566 
567 
568                     ////////////////////////////////////
569                     // Accumulate contracted integrals
570                     ////////////////////////////////////
571                     if(lastoffset == 0)
572                     {
573                         contract_all(420, PRIM_INT__g_s_i_s, PRIM_PTR_INT__g_s_i_s);
574                         contract_all(540, PRIM_INT__g_s_k_s, PRIM_PTR_INT__g_s_k_s);
575                         contract_all(588, PRIM_INT__h_s_i_s, PRIM_PTR_INT__h_s_i_s);
576                         contract_all(756, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
577                     }
578                     else
579                     {
580                         contract(420, shelloffsets, PRIM_INT__g_s_i_s, PRIM_PTR_INT__g_s_i_s);
581                         contract(540, shelloffsets, PRIM_INT__g_s_k_s, PRIM_PTR_INT__g_s_k_s);
582                         contract(588, shelloffsets, PRIM_INT__h_s_i_s, PRIM_PTR_INT__h_s_i_s);
583                         contract(756, shelloffsets, PRIM_INT__h_s_k_s, PRIM_PTR_INT__h_s_k_s);
584                         PRIM_PTR_INT__g_s_i_s += lastoffset*420;
585                         PRIM_PTR_INT__g_s_k_s += lastoffset*540;
586                         PRIM_PTR_INT__h_s_i_s += lastoffset*588;
587                         PRIM_PTR_INT__h_s_k_s += lastoffset*756;
588                     }
589 
590                 }  // close loop over j
591             }  // close loop over i
592 
593             //Advance to the next batch
594             jstart = SIMINT_SIMD_ROUND(jend);
595 
596             //////////////////////////////////////////////
597             // Contracted integrals: Horizontal recurrance
598             //////////////////////////////////////////////
599 
600 
601             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
602 
603 
604             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
605             {
606                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
607 
608                 // set up HRR pointers
609                 double const * restrict HRR_INT__g_s_i_s = INT__g_s_i_s + abcd * 420;
610                 double const * restrict HRR_INT__g_s_k_s = INT__g_s_k_s + abcd * 540;
611                 double const * restrict HRR_INT__h_s_i_s = INT__h_s_i_s + abcd * 588;
612                 double const * restrict HRR_INT__h_s_k_s = INT__h_s_k_s + abcd * 756;
613                 double * restrict HRR_INT__g_p_i_p = INT__g_p_i_p + real_abcd * 3780;
614 
615                 // form INT__g_p_i_s
616                 HRR_J_g_p(
617                     HRR_INT__g_p_i_s,
618                     HRR_INT__g_s_i_s,
619                     HRR_INT__h_s_i_s,
620                     hAB, 28);
621 
622                 // form INT__g_p_k_s
623                 HRR_J_g_p(
624                     HRR_INT__g_p_k_s,
625                     HRR_INT__g_s_k_s,
626                     HRR_INT__h_s_k_s,
627                     hAB, 36);
628 
629                 // form INT__g_p_i_p
630                 ostei_general_hrr_L(4, 1, 6, 1, hCD, HRR_INT__g_p_k_s, HRR_INT__g_p_i_s, HRR_INT__g_p_i_p);
631 
632 
633             }  // close HRR loop
634 
635 
636         }   // close loop cdbatch
637 
638         istart = iend;
639     }  // close loop over ab
640 
641     return P.nshell12_clip * Q.nshell12_clip;
642 }
643 
ostei_p_g_i_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__p_g_i_p)644 int ostei_p_g_i_p(struct simint_multi_shellpair const P,
645                   struct simint_multi_shellpair const Q,
646                   double screen_tol,
647                   double * const restrict work,
648                   double * const restrict INT__p_g_i_p)
649 {
650     double P_AB[3*P.nshell12];
651     struct simint_multi_shellpair P_tmp = P;
652     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
653     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
654     P_tmp.AB_x = P_AB;
655     P_tmp.AB_y = P_AB + P.nshell12;
656     P_tmp.AB_z = P_AB + 2*P.nshell12;
657 
658     for(int i = 0; i < P.nshell12; i++)
659     {
660         P_tmp.AB_x[i] = -P.AB_x[i];
661         P_tmp.AB_y[i] = -P.AB_y[i];
662         P_tmp.AB_z[i] = -P.AB_z[i];
663     }
664 
665     int ret = ostei_g_p_i_p(P_tmp, Q, screen_tol, work, INT__p_g_i_p);
666     double buffer[3780] SIMINT_ALIGN_ARRAY_DBL;
667 
668     for(int q = 0; q < ret; q++)
669     {
670         int idx = 0;
671         for(int a = 0; a < 3; ++a)
672         for(int b = 0; b < 15; ++b)
673         for(int c = 0; c < 28; ++c)
674         for(int d = 0; d < 3; ++d)
675             buffer[idx++] = INT__p_g_i_p[q*3780+b*252+a*84+c*3+d];
676 
677         memcpy(INT__p_g_i_p+q*3780, buffer, 3780*sizeof(double));
678     }
679 
680     return ret;
681 }
682 
ostei_g_p_p_i(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_p_p_i)683 int ostei_g_p_p_i(struct simint_multi_shellpair const P,
684                   struct simint_multi_shellpair const Q,
685                   double screen_tol,
686                   double * const restrict work,
687                   double * const restrict INT__g_p_p_i)
688 {
689     double Q_AB[3*Q.nshell12];
690     struct simint_multi_shellpair Q_tmp = Q;
691     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
692     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
693     Q_tmp.AB_x = Q_AB;
694     Q_tmp.AB_y = Q_AB + Q.nshell12;
695     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
696 
697     for(int i = 0; i < Q.nshell12; i++)
698     {
699         Q_tmp.AB_x[i] = -Q.AB_x[i];
700         Q_tmp.AB_y[i] = -Q.AB_y[i];
701         Q_tmp.AB_z[i] = -Q.AB_z[i];
702     }
703 
704     int ret = ostei_g_p_i_p(P, Q_tmp, screen_tol, work, INT__g_p_p_i);
705     double buffer[3780] SIMINT_ALIGN_ARRAY_DBL;
706 
707     for(int q = 0; q < ret; q++)
708     {
709         int idx = 0;
710         for(int a = 0; a < 15; ++a)
711         for(int b = 0; b < 3; ++b)
712         for(int c = 0; c < 3; ++c)
713         for(int d = 0; d < 28; ++d)
714             buffer[idx++] = INT__g_p_p_i[q*3780+a*252+b*84+d*3+c];
715 
716         memcpy(INT__g_p_p_i+q*3780, buffer, 3780*sizeof(double));
717     }
718 
719     return ret;
720 }
721 
ostei_p_g_p_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_g_p_i)722 int ostei_p_g_p_i(struct simint_multi_shellpair const P,
723                   struct simint_multi_shellpair const Q,
724                   double screen_tol,
725                   double * const restrict work,
726                   double * const restrict INT__p_g_p_i)
727 {
728     double P_AB[3*P.nshell12];
729     struct simint_multi_shellpair P_tmp = P;
730     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
731     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
732     P_tmp.AB_x = P_AB;
733     P_tmp.AB_y = P_AB + P.nshell12;
734     P_tmp.AB_z = P_AB + 2*P.nshell12;
735 
736     for(int i = 0; i < P.nshell12; i++)
737     {
738         P_tmp.AB_x[i] = -P.AB_x[i];
739         P_tmp.AB_y[i] = -P.AB_y[i];
740         P_tmp.AB_z[i] = -P.AB_z[i];
741     }
742 
743     double Q_AB[3*Q.nshell12];
744     struct simint_multi_shellpair Q_tmp = Q;
745     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
746     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
747     Q_tmp.AB_x = Q_AB;
748     Q_tmp.AB_y = Q_AB + Q.nshell12;
749     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
750 
751     for(int i = 0; i < Q.nshell12; i++)
752     {
753         Q_tmp.AB_x[i] = -Q.AB_x[i];
754         Q_tmp.AB_y[i] = -Q.AB_y[i];
755         Q_tmp.AB_z[i] = -Q.AB_z[i];
756     }
757 
758     int ret = ostei_g_p_i_p(P_tmp, Q_tmp, screen_tol, work, INT__p_g_p_i);
759     double buffer[3780] SIMINT_ALIGN_ARRAY_DBL;
760 
761     for(int q = 0; q < ret; q++)
762     {
763         int idx = 0;
764         for(int a = 0; a < 3; ++a)
765         for(int b = 0; b < 15; ++b)
766         for(int c = 0; c < 3; ++c)
767         for(int d = 0; d < 28; ++d)
768             buffer[idx++] = INT__p_g_p_i[q*3780+b*252+a*84+d*3+c];
769 
770         memcpy(INT__p_g_p_i+q*3780, buffer, 3780*sizeof(double));
771     }
772 
773     return ret;
774 }
775 
776