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_f_d_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_f_d_p)8 int ostei_g_f_d_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_f_d_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__g_f_d_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_d_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__g_s_f_s = work + (SIMINT_NSHELL_SIMD * 90);
31     double * const INT__h_s_d_s = work + (SIMINT_NSHELL_SIMD * 240);
32     double * const INT__h_s_f_s = work + (SIMINT_NSHELL_SIMD * 366);
33     double * const INT__i_s_d_s = work + (SIMINT_NSHELL_SIMD * 576);
34     double * const INT__i_s_f_s = work + (SIMINT_NSHELL_SIMD * 744);
35     double * const INT__k_s_d_s = work + (SIMINT_NSHELL_SIMD * 1024);
36     double * const INT__k_s_f_s = work + (SIMINT_NSHELL_SIMD * 1240);
37     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*1600);
38     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 11;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 41;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_p_s = primwork + 95;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 149;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 229;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_d_s = primwork + 319;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 439;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 544;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 679;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_f_s = primwork + 859;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 1009;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 1135;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 1324;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 1576;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 1786;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 1926;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 2178;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 2514;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 2794;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 2938;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 3262;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_f_s = primwork + 3694;
61     double * const hrrwork = (double *)(primwork + 4054);
62     double * const HRR_INT__g_p_d_s = hrrwork + 0;
63     double * const HRR_INT__g_p_f_s = hrrwork + 270;
64     double * const HRR_INT__g_d_d_s = hrrwork + 720;
65     double * const HRR_INT__g_d_f_s = hrrwork + 1260;
66     double * const HRR_INT__g_f_d_s = hrrwork + 2160;
67     double * const HRR_INT__g_f_f_s = hrrwork + 3060;
68     double * const HRR_INT__h_p_d_s = hrrwork + 4560;
69     double * const HRR_INT__h_p_f_s = hrrwork + 4938;
70     double * const HRR_INT__h_d_d_s = hrrwork + 5568;
71     double * const HRR_INT__h_d_f_s = hrrwork + 6324;
72     double * const HRR_INT__i_p_d_s = hrrwork + 7584;
73     double * const HRR_INT__i_p_f_s = hrrwork + 8088;
74 
75 
76     // Create constants
77     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
78     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
79     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
80     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
81     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
82     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
83     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
84     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
85 
86 
87     ////////////////////////////////////////
88     // Loop over shells and primitives
89     ////////////////////////////////////////
90 
91     real_abcd = 0;
92     istart = 0;
93     for(ab = 0; ab < P.nshell12_clip; ++ab)
94     {
95         const int iend = istart + P.nprim12[ab];
96 
97         cd = 0;
98         jstart = 0;
99 
100         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
101         {
102             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
103             int jend = jstart;
104             for(i = 0; i < nshellbatch; i++)
105                 jend += Q.nprim12[cd+i];
106 
107             // Clear the beginning of the workspace (where we are accumulating integrals)
108             memset(work, 0, SIMINT_NSHELL_SIMD * 1600 * sizeof(double));
109             abcd = 0;
110 
111 
112             for(i = istart; i < iend; ++i)
113             {
114                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
115 
116                 if(check_screen)
117                 {
118                     // Skip this whole thing if always insignificant
119                     if((P.screen[i] * Q.screen_max) < screen_tol)
120                         continue;
121                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
122                 }
123 
124                 icd = 0;
125                 iprimcd = 0;
126                 nprim_icd = Q.nprim12[cd];
127                 double * restrict PRIM_PTR_INT__g_s_d_s = INT__g_s_d_s + abcd * 90;
128                 double * restrict PRIM_PTR_INT__g_s_f_s = INT__g_s_f_s + abcd * 150;
129                 double * restrict PRIM_PTR_INT__h_s_d_s = INT__h_s_d_s + abcd * 126;
130                 double * restrict PRIM_PTR_INT__h_s_f_s = INT__h_s_f_s + abcd * 210;
131                 double * restrict PRIM_PTR_INT__i_s_d_s = INT__i_s_d_s + abcd * 168;
132                 double * restrict PRIM_PTR_INT__i_s_f_s = INT__i_s_f_s + abcd * 280;
133                 double * restrict PRIM_PTR_INT__k_s_d_s = INT__k_s_d_s + abcd * 216;
134                 double * restrict PRIM_PTR_INT__k_s_f_s = INT__k_s_f_s + abcd * 360;
135 
136 
137 
138                 // Load these one per loop over i
139                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
140                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
141                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
142 
143                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
144 
145                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
146                 {
147                     // calculate the shell offsets
148                     // these are the offset from the shell pointed to by cd
149                     // for each element
150                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
151                     int lastoffset = 0;
152                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
153 
154                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
155                     {
156                         // Handle if the first element of the vector is a new shell
157                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
158                         {
159                             nprim_icd += Q.nprim12[cd + (++icd)];
160                             PRIM_PTR_INT__g_s_d_s += 90;
161                             PRIM_PTR_INT__g_s_f_s += 150;
162                             PRIM_PTR_INT__h_s_d_s += 126;
163                             PRIM_PTR_INT__h_s_f_s += 210;
164                             PRIM_PTR_INT__i_s_d_s += 168;
165                             PRIM_PTR_INT__i_s_f_s += 280;
166                             PRIM_PTR_INT__k_s_d_s += 216;
167                             PRIM_PTR_INT__k_s_f_s += 360;
168                         }
169                         iprimcd++;
170                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
171                         {
172                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
173                             {
174                                 shelloffsets[n] = shelloffsets[n-1] + 1;
175                                 lastoffset++;
176                                 nprim_icd += Q.nprim12[cd + (++icd)];
177                             }
178                             else
179                                 shelloffsets[n] = shelloffsets[n-1];
180                             iprimcd++;
181                         }
182                     }
183                     else
184                         iprimcd += SIMINT_SIMD_LEN;
185 
186                     // Do we have to compute this vector (or has it been screened out)?
187                     // (not_screened != 0 means we have to do this vector)
188                     if(check_screen)
189                     {
190                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
191                         if(vmax < screen_tol)
192                         {
193                             PRIM_PTR_INT__g_s_d_s += lastoffset*90;
194                             PRIM_PTR_INT__g_s_f_s += lastoffset*150;
195                             PRIM_PTR_INT__h_s_d_s += lastoffset*126;
196                             PRIM_PTR_INT__h_s_f_s += lastoffset*210;
197                             PRIM_PTR_INT__i_s_d_s += lastoffset*168;
198                             PRIM_PTR_INT__i_s_f_s += lastoffset*280;
199                             PRIM_PTR_INT__k_s_d_s += lastoffset*216;
200                             PRIM_PTR_INT__k_s_f_s += lastoffset*360;
201                             continue;
202                         }
203                     }
204 
205                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
206                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
207                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
208                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
209 
210 
211                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
212                     SIMINT_DBLTYPE PQ[3];
213                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
214                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
215                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
216                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
217                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
218                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
219 
220                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
221                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
222                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
223                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
224                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
225                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
226                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
227 
228                     // NOTE: Minus sign!
229                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
230                     SIMINT_DBLTYPE aop_PQ[3];
231                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
232                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
233                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
234 
235                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
236                     SIMINT_DBLTYPE aoq_PQ[3];
237                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
238                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
239                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
240                     // Put a minus sign here so we don't have to in RR routines
241                     a_over_q = SIMINT_NEG(a_over_q);
242 
243 
244                     //////////////////////////////////////////////
245                     // Fjt function section
246                     // Maximum v value: 10
247                     //////////////////////////////////////////////
248                     // The parameter to the Fjt function
249                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
250 
251 
252                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
253 
254 
255                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 10);
256                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
257                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
258                     for(n = 0; n <= 10; n++)
259                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
260 
261                     //////////////////////////////////////////////
262                     // Primitive integrals: Vertical recurrance
263                     //////////////////////////////////////////////
264 
265                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
266                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
267                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
268                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
269                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
270                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
271                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
272                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
273                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
274                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
275                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
276                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
277                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
278                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
279                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
280 
281 
282 
283                     // Forming PRIM_INT__p_s_s_s[10 * 3];
284                     for(n = 0; n < 10; ++n)  // loop over orders of auxiliary function
285                     {
286 
287                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
288                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]);
289 
290                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
291                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 1]);
292 
293                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
294                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_s[n * 3 + 2]);
295 
296                     }
297 
298 
299 
300                     // Forming PRIM_INT__d_s_s_s[9 * 6];
301                     for(n = 0; n < 9; ++n)  // loop over orders of auxiliary function
302                     {
303 
304                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
305                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 0]);
306                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 0]);
307 
308                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
309                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 1]);
310 
311                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
312                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_s[n * 6 + 2]);
313 
314                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
315                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_s[n * 6 + 3]);
316                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 3]);
317 
318                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
319                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_s[n * 6 + 4]);
320 
321                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
322                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_s_s[n * 6 + 5]);
323                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_FMADD( vrr_const_1_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__d_s_s_s[n * 6 + 5]);
324 
325                     }
326 
327 
328 
329                     // Forming PRIM_INT__f_s_s_s[8 * 10];
330                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
331                     {
332 
333                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
334                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 0]);
335                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__f_s_s_s[n * 10 + 0]);
336 
337                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
338                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 1]);
339 
340                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
341                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 2]);
342 
343                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
344                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 3]);
345 
346                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
347                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__f_s_s_s[n * 10 + 4]);
348 
349                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
350                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 5]);
351 
352                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
353                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 6]);
354                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__f_s_s_s[n * 10 + 6]);
355 
356                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
357                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 7]);
358 
359                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
360                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 8]);
361 
362                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
363                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 9]);
364                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__f_s_s_s[n * 10 + 9]);
365 
366                     }
367 
368 
369                     VRR_I_g_s_s_s(
370                             PRIM_INT__g_s_s_s,
371                             PRIM_INT__f_s_s_s,
372                             PRIM_INT__d_s_s_s,
373                             P_PA,
374                             a_over_p,
375                             aop_PQ,
376                             one_over_2p,
377                             7);
378 
379 
380                     VRR_K_g_s_p_s(
381                             PRIM_INT__g_s_p_s,
382                             PRIM_INT__g_s_s_s,
383                             PRIM_INT__f_s_s_s,
384                             Q_PA,
385                             aoq_PQ,
386                             one_over_2pq,
387                             3);
388 
389 
390                     VRR_K_f_s_p_s(
391                             PRIM_INT__f_s_p_s,
392                             PRIM_INT__f_s_s_s,
393                             PRIM_INT__d_s_s_s,
394                             Q_PA,
395                             aoq_PQ,
396                             one_over_2pq,
397                             3);
398 
399 
400                     ostei_general_vrr_K(4, 0, 2, 0, 2,
401                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
402                             PRIM_INT__g_s_p_s, PRIM_INT__g_s_s_s, NULL, PRIM_INT__f_s_p_s, NULL, PRIM_INT__g_s_d_s);
403 
404 
405 
406                     // Forming PRIM_INT__d_s_p_s[3 * 18];
407                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
408                     {
409 
410                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
411                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
412                         PRIM_INT__d_s_p_s[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 0]);
413 
414                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
415                         PRIM_INT__d_s_p_s[n * 18 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 1]);
416 
417                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
418                         PRIM_INT__d_s_p_s[n * 18 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_p_s[n * 18 + 2]);
419 
420                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
421                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 3]);
422                         PRIM_INT__d_s_p_s[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 3]);
423 
424                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
425                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 4]);
426                         PRIM_INT__d_s_p_s[n * 18 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 4]);
427 
428                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
429                         PRIM_INT__d_s_p_s[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_p_s[n * 18 + 5]);
430 
431                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
432                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 6]);
433                         PRIM_INT__d_s_p_s[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 6]);
434 
435                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
436                         PRIM_INT__d_s_p_s[n * 18 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 7]);
437 
438                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
439                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_p_s[n * 18 + 8]);
440                         PRIM_INT__d_s_p_s[n * 18 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_p_s[n * 18 + 8]);
441 
442                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
443                         PRIM_INT__d_s_p_s[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 9]);
444 
445                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
446                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 10]);
447                         PRIM_INT__d_s_p_s[n * 18 + 10] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 10]);
448 
449                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
450                         PRIM_INT__d_s_p_s[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_p_s[n * 18 + 11]);
451 
452                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
453                         PRIM_INT__d_s_p_s[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 12]);
454 
455                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
456                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 13]);
457                         PRIM_INT__d_s_p_s[n * 18 + 13] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 13]);
458 
459                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
460                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_p_s[n * 18 + 14]);
461                         PRIM_INT__d_s_p_s[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_p_s[n * 18 + 14]);
462 
463                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_MUL(Q_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
464                         PRIM_INT__d_s_p_s[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 15]);
465 
466                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_MUL(Q_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
467                         PRIM_INT__d_s_p_s[n * 18 + 16] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 16]);
468 
469                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_MUL(Q_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
470                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_p_s[n * 18 + 17]);
471                         PRIM_INT__d_s_p_s[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_p_s[n * 18 + 17]);
472 
473                     }
474 
475 
476                     VRR_K_f_s_d_s(
477                             PRIM_INT__f_s_d_s,
478                             PRIM_INT__f_s_p_s,
479                             PRIM_INT__f_s_s_s,
480                             PRIM_INT__d_s_p_s,
481                             Q_PA,
482                             a_over_q,
483                             aoq_PQ,
484                             one_over_2pq,
485                             one_over_2q,
486                             2);
487 
488 
489                     ostei_general_vrr_K(4, 0, 3, 0, 1,
490                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
491                             PRIM_INT__g_s_d_s, PRIM_INT__g_s_p_s, NULL, PRIM_INT__f_s_d_s, NULL, PRIM_INT__g_s_f_s);
492 
493 
494                     VRR_I_h_s_s_s(
495                             PRIM_INT__h_s_s_s,
496                             PRIM_INT__g_s_s_s,
497                             PRIM_INT__f_s_s_s,
498                             P_PA,
499                             a_over_p,
500                             aop_PQ,
501                             one_over_2p,
502                             6);
503 
504 
505                     ostei_general_vrr_K(5, 0, 1, 0, 3,
506                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
507                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
508 
509 
510                     ostei_general_vrr_K(5, 0, 2, 0, 2,
511                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
512                             PRIM_INT__h_s_p_s, PRIM_INT__h_s_s_s, NULL, PRIM_INT__g_s_p_s, NULL, PRIM_INT__h_s_d_s);
513 
514 
515                     ostei_general_vrr_K(5, 0, 3, 0, 1,
516                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
517                             PRIM_INT__h_s_d_s, PRIM_INT__h_s_p_s, NULL, PRIM_INT__g_s_d_s, NULL, PRIM_INT__h_s_f_s);
518 
519 
520                     ostei_general_vrr1_I(6, 5,
521                             one_over_2p, a_over_p, aop_PQ, P_PA,
522                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
523 
524 
525                     ostei_general_vrr_K(6, 0, 1, 0, 3,
526                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
527                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
528 
529 
530                     ostei_general_vrr_K(6, 0, 2, 0, 2,
531                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
532                             PRIM_INT__i_s_p_s, PRIM_INT__i_s_s_s, NULL, PRIM_INT__h_s_p_s, NULL, PRIM_INT__i_s_d_s);
533 
534 
535                     ostei_general_vrr_K(6, 0, 3, 0, 1,
536                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
537                             PRIM_INT__i_s_d_s, PRIM_INT__i_s_p_s, NULL, PRIM_INT__h_s_d_s, NULL, PRIM_INT__i_s_f_s);
538 
539 
540                     ostei_general_vrr1_I(7, 4,
541                             one_over_2p, a_over_p, aop_PQ, P_PA,
542                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
543 
544 
545                     ostei_general_vrr_K(7, 0, 1, 0, 3,
546                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
547                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
548 
549 
550                     ostei_general_vrr_K(7, 0, 2, 0, 2,
551                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
552                             PRIM_INT__k_s_p_s, PRIM_INT__k_s_s_s, NULL, PRIM_INT__i_s_p_s, NULL, PRIM_INT__k_s_d_s);
553 
554 
555                     ostei_general_vrr_K(7, 0, 3, 0, 1,
556                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
557                             PRIM_INT__k_s_d_s, PRIM_INT__k_s_p_s, NULL, PRIM_INT__i_s_d_s, NULL, PRIM_INT__k_s_f_s);
558 
559 
560 
561 
562                     ////////////////////////////////////
563                     // Accumulate contracted integrals
564                     ////////////////////////////////////
565                     if(lastoffset == 0)
566                     {
567                         contract_all(90, PRIM_INT__g_s_d_s, PRIM_PTR_INT__g_s_d_s);
568                         contract_all(150, PRIM_INT__g_s_f_s, PRIM_PTR_INT__g_s_f_s);
569                         contract_all(126, PRIM_INT__h_s_d_s, PRIM_PTR_INT__h_s_d_s);
570                         contract_all(210, PRIM_INT__h_s_f_s, PRIM_PTR_INT__h_s_f_s);
571                         contract_all(168, PRIM_INT__i_s_d_s, PRIM_PTR_INT__i_s_d_s);
572                         contract_all(280, PRIM_INT__i_s_f_s, PRIM_PTR_INT__i_s_f_s);
573                         contract_all(216, PRIM_INT__k_s_d_s, PRIM_PTR_INT__k_s_d_s);
574                         contract_all(360, PRIM_INT__k_s_f_s, PRIM_PTR_INT__k_s_f_s);
575                     }
576                     else
577                     {
578                         contract(90, shelloffsets, PRIM_INT__g_s_d_s, PRIM_PTR_INT__g_s_d_s);
579                         contract(150, shelloffsets, PRIM_INT__g_s_f_s, PRIM_PTR_INT__g_s_f_s);
580                         contract(126, shelloffsets, PRIM_INT__h_s_d_s, PRIM_PTR_INT__h_s_d_s);
581                         contract(210, shelloffsets, PRIM_INT__h_s_f_s, PRIM_PTR_INT__h_s_f_s);
582                         contract(168, shelloffsets, PRIM_INT__i_s_d_s, PRIM_PTR_INT__i_s_d_s);
583                         contract(280, shelloffsets, PRIM_INT__i_s_f_s, PRIM_PTR_INT__i_s_f_s);
584                         contract(216, shelloffsets, PRIM_INT__k_s_d_s, PRIM_PTR_INT__k_s_d_s);
585                         contract(360, shelloffsets, PRIM_INT__k_s_f_s, PRIM_PTR_INT__k_s_f_s);
586                         PRIM_PTR_INT__g_s_d_s += lastoffset*90;
587                         PRIM_PTR_INT__g_s_f_s += lastoffset*150;
588                         PRIM_PTR_INT__h_s_d_s += lastoffset*126;
589                         PRIM_PTR_INT__h_s_f_s += lastoffset*210;
590                         PRIM_PTR_INT__i_s_d_s += lastoffset*168;
591                         PRIM_PTR_INT__i_s_f_s += lastoffset*280;
592                         PRIM_PTR_INT__k_s_d_s += lastoffset*216;
593                         PRIM_PTR_INT__k_s_f_s += lastoffset*360;
594                     }
595 
596                 }  // close loop over j
597             }  // close loop over i
598 
599             //Advance to the next batch
600             jstart = SIMINT_SIMD_ROUND(jend);
601 
602             //////////////////////////////////////////////
603             // Contracted integrals: Horizontal recurrance
604             //////////////////////////////////////////////
605 
606 
607             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
608 
609 
610             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
611             {
612                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
613 
614                 // set up HRR pointers
615                 double const * restrict HRR_INT__g_s_d_s = INT__g_s_d_s + abcd * 90;
616                 double const * restrict HRR_INT__g_s_f_s = INT__g_s_f_s + abcd * 150;
617                 double const * restrict HRR_INT__h_s_d_s = INT__h_s_d_s + abcd * 126;
618                 double const * restrict HRR_INT__h_s_f_s = INT__h_s_f_s + abcd * 210;
619                 double const * restrict HRR_INT__i_s_d_s = INT__i_s_d_s + abcd * 168;
620                 double const * restrict HRR_INT__i_s_f_s = INT__i_s_f_s + abcd * 280;
621                 double const * restrict HRR_INT__k_s_d_s = INT__k_s_d_s + abcd * 216;
622                 double const * restrict HRR_INT__k_s_f_s = INT__k_s_f_s + abcd * 360;
623                 double * restrict HRR_INT__g_f_d_p = INT__g_f_d_p + real_abcd * 2700;
624 
625                 // form INT__g_p_d_s
626                 HRR_J_g_p(
627                     HRR_INT__g_p_d_s,
628                     HRR_INT__g_s_d_s,
629                     HRR_INT__h_s_d_s,
630                     hAB, 6);
631 
632                 // form INT__g_p_f_s
633                 HRR_J_g_p(
634                     HRR_INT__g_p_f_s,
635                     HRR_INT__g_s_f_s,
636                     HRR_INT__h_s_f_s,
637                     hAB, 10);
638 
639                 // form INT__h_p_d_s
640                 ostei_general_hrr_J(5, 1, 2, 0, hAB, HRR_INT__i_s_d_s, HRR_INT__h_s_d_s, HRR_INT__h_p_d_s);
641 
642                 // form INT__h_p_f_s
643                 ostei_general_hrr_J(5, 1, 3, 0, hAB, HRR_INT__i_s_f_s, HRR_INT__h_s_f_s, HRR_INT__h_p_f_s);
644 
645                 // form INT__i_p_d_s
646                 ostei_general_hrr_J(6, 1, 2, 0, hAB, HRR_INT__k_s_d_s, HRR_INT__i_s_d_s, HRR_INT__i_p_d_s);
647 
648                 // form INT__i_p_f_s
649                 ostei_general_hrr_J(6, 1, 3, 0, hAB, HRR_INT__k_s_f_s, HRR_INT__i_s_f_s, HRR_INT__i_p_f_s);
650 
651                 // form INT__g_d_d_s
652                 ostei_general_hrr_J(4, 2, 2, 0, hAB, HRR_INT__h_p_d_s, HRR_INT__g_p_d_s, HRR_INT__g_d_d_s);
653 
654                 // form INT__g_d_f_s
655                 ostei_general_hrr_J(4, 2, 3, 0, hAB, HRR_INT__h_p_f_s, HRR_INT__g_p_f_s, HRR_INT__g_d_f_s);
656 
657                 // form INT__h_d_d_s
658                 ostei_general_hrr_J(5, 2, 2, 0, hAB, HRR_INT__i_p_d_s, HRR_INT__h_p_d_s, HRR_INT__h_d_d_s);
659 
660                 // form INT__h_d_f_s
661                 ostei_general_hrr_J(5, 2, 3, 0, hAB, HRR_INT__i_p_f_s, HRR_INT__h_p_f_s, HRR_INT__h_d_f_s);
662 
663                 // form INT__g_f_d_s
664                 ostei_general_hrr_J(4, 3, 2, 0, hAB, HRR_INT__h_d_d_s, HRR_INT__g_d_d_s, HRR_INT__g_f_d_s);
665 
666                 // form INT__g_f_f_s
667                 ostei_general_hrr_J(4, 3, 3, 0, hAB, HRR_INT__h_d_f_s, HRR_INT__g_d_f_s, HRR_INT__g_f_f_s);
668 
669                 // form INT__g_f_d_p
670                 for(ibra = 0; ibra < 150; ++ibra)
671                 {
672                     HRR_INT__g_f_d_p[ibra * 18 + 0] = HRR_INT__g_f_f_s[ibra * 10 + 0] + ( hCD[0] * HRR_INT__g_f_d_s[ibra * 6 + 0] );
673 
674                     HRR_INT__g_f_d_p[ibra * 18 + 1] = HRR_INT__g_f_f_s[ibra * 10 + 1] + ( hCD[1] * HRR_INT__g_f_d_s[ibra * 6 + 0] );
675 
676                     HRR_INT__g_f_d_p[ibra * 18 + 2] = HRR_INT__g_f_f_s[ibra * 10 + 2] + ( hCD[2] * HRR_INT__g_f_d_s[ibra * 6 + 0] );
677 
678                     HRR_INT__g_f_d_p[ibra * 18 + 3] = HRR_INT__g_f_f_s[ibra * 10 + 1] + ( hCD[0] * HRR_INT__g_f_d_s[ibra * 6 + 1] );
679 
680                     HRR_INT__g_f_d_p[ibra * 18 + 4] = HRR_INT__g_f_f_s[ibra * 10 + 3] + ( hCD[1] * HRR_INT__g_f_d_s[ibra * 6 + 1] );
681 
682                     HRR_INT__g_f_d_p[ibra * 18 + 5] = HRR_INT__g_f_f_s[ibra * 10 + 4] + ( hCD[2] * HRR_INT__g_f_d_s[ibra * 6 + 1] );
683 
684                     HRR_INT__g_f_d_p[ibra * 18 + 6] = HRR_INT__g_f_f_s[ibra * 10 + 2] + ( hCD[0] * HRR_INT__g_f_d_s[ibra * 6 + 2] );
685 
686                     HRR_INT__g_f_d_p[ibra * 18 + 7] = HRR_INT__g_f_f_s[ibra * 10 + 4] + ( hCD[1] * HRR_INT__g_f_d_s[ibra * 6 + 2] );
687 
688                     HRR_INT__g_f_d_p[ibra * 18 + 8] = HRR_INT__g_f_f_s[ibra * 10 + 5] + ( hCD[2] * HRR_INT__g_f_d_s[ibra * 6 + 2] );
689 
690                     HRR_INT__g_f_d_p[ibra * 18 + 9] = HRR_INT__g_f_f_s[ibra * 10 + 3] + ( hCD[0] * HRR_INT__g_f_d_s[ibra * 6 + 3] );
691 
692                     HRR_INT__g_f_d_p[ibra * 18 + 10] = HRR_INT__g_f_f_s[ibra * 10 + 6] + ( hCD[1] * HRR_INT__g_f_d_s[ibra * 6 + 3] );
693 
694                     HRR_INT__g_f_d_p[ibra * 18 + 11] = HRR_INT__g_f_f_s[ibra * 10 + 7] + ( hCD[2] * HRR_INT__g_f_d_s[ibra * 6 + 3] );
695 
696                     HRR_INT__g_f_d_p[ibra * 18 + 12] = HRR_INT__g_f_f_s[ibra * 10 + 4] + ( hCD[0] * HRR_INT__g_f_d_s[ibra * 6 + 4] );
697 
698                     HRR_INT__g_f_d_p[ibra * 18 + 13] = HRR_INT__g_f_f_s[ibra * 10 + 7] + ( hCD[1] * HRR_INT__g_f_d_s[ibra * 6 + 4] );
699 
700                     HRR_INT__g_f_d_p[ibra * 18 + 14] = HRR_INT__g_f_f_s[ibra * 10 + 8] + ( hCD[2] * HRR_INT__g_f_d_s[ibra * 6 + 4] );
701 
702                     HRR_INT__g_f_d_p[ibra * 18 + 15] = HRR_INT__g_f_f_s[ibra * 10 + 5] + ( hCD[0] * HRR_INT__g_f_d_s[ibra * 6 + 5] );
703 
704                     HRR_INT__g_f_d_p[ibra * 18 + 16] = HRR_INT__g_f_f_s[ibra * 10 + 8] + ( hCD[1] * HRR_INT__g_f_d_s[ibra * 6 + 5] );
705 
706                     HRR_INT__g_f_d_p[ibra * 18 + 17] = HRR_INT__g_f_f_s[ibra * 10 + 9] + ( hCD[2] * HRR_INT__g_f_d_s[ibra * 6 + 5] );
707 
708                 }
709 
710 
711             }  // close HRR loop
712 
713 
714         }   // close loop cdbatch
715 
716         istart = iend;
717     }  // close loop over ab
718 
719     return P.nshell12_clip * Q.nshell12_clip;
720 }
721 
ostei_f_g_d_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__f_g_d_p)722 int ostei_f_g_d_p(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__f_g_d_p)
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     int ret = ostei_g_f_d_p(P_tmp, Q, screen_tol, work, INT__f_g_d_p);
744     double buffer[2700] SIMINT_ALIGN_ARRAY_DBL;
745 
746     for(int q = 0; q < ret; q++)
747     {
748         int idx = 0;
749         for(int a = 0; a < 10; ++a)
750         for(int b = 0; b < 15; ++b)
751         for(int c = 0; c < 6; ++c)
752         for(int d = 0; d < 3; ++d)
753             buffer[idx++] = INT__f_g_d_p[q*2700+b*180+a*18+c*3+d];
754 
755         memcpy(INT__f_g_d_p+q*2700, buffer, 2700*sizeof(double));
756     }
757 
758     return ret;
759 }
760 
ostei_g_f_p_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__g_f_p_d)761 int ostei_g_f_p_d(struct simint_multi_shellpair const P,
762                   struct simint_multi_shellpair const Q,
763                   double screen_tol,
764                   double * const restrict work,
765                   double * const restrict INT__g_f_p_d)
766 {
767     double Q_AB[3*Q.nshell12];
768     struct simint_multi_shellpair Q_tmp = Q;
769     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
770     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
771     Q_tmp.AB_x = Q_AB;
772     Q_tmp.AB_y = Q_AB + Q.nshell12;
773     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
774 
775     for(int i = 0; i < Q.nshell12; i++)
776     {
777         Q_tmp.AB_x[i] = -Q.AB_x[i];
778         Q_tmp.AB_y[i] = -Q.AB_y[i];
779         Q_tmp.AB_z[i] = -Q.AB_z[i];
780     }
781 
782     int ret = ostei_g_f_d_p(P, Q_tmp, screen_tol, work, INT__g_f_p_d);
783     double buffer[2700] SIMINT_ALIGN_ARRAY_DBL;
784 
785     for(int q = 0; q < ret; q++)
786     {
787         int idx = 0;
788         for(int a = 0; a < 15; ++a)
789         for(int b = 0; b < 10; ++b)
790         for(int c = 0; c < 3; ++c)
791         for(int d = 0; d < 6; ++d)
792             buffer[idx++] = INT__g_f_p_d[q*2700+a*180+b*18+d*3+c];
793 
794         memcpy(INT__g_f_p_d+q*2700, buffer, 2700*sizeof(double));
795     }
796 
797     return ret;
798 }
799 
ostei_f_g_p_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__f_g_p_d)800 int ostei_f_g_p_d(struct simint_multi_shellpair const P,
801                   struct simint_multi_shellpair const Q,
802                   double screen_tol,
803                   double * const restrict work,
804                   double * const restrict INT__f_g_p_d)
805 {
806     double P_AB[3*P.nshell12];
807     struct simint_multi_shellpair P_tmp = P;
808     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
809     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
810     P_tmp.AB_x = P_AB;
811     P_tmp.AB_y = P_AB + P.nshell12;
812     P_tmp.AB_z = P_AB + 2*P.nshell12;
813 
814     for(int i = 0; i < P.nshell12; i++)
815     {
816         P_tmp.AB_x[i] = -P.AB_x[i];
817         P_tmp.AB_y[i] = -P.AB_y[i];
818         P_tmp.AB_z[i] = -P.AB_z[i];
819     }
820 
821     double Q_AB[3*Q.nshell12];
822     struct simint_multi_shellpair Q_tmp = Q;
823     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
824     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
825     Q_tmp.AB_x = Q_AB;
826     Q_tmp.AB_y = Q_AB + Q.nshell12;
827     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
828 
829     for(int i = 0; i < Q.nshell12; i++)
830     {
831         Q_tmp.AB_x[i] = -Q.AB_x[i];
832         Q_tmp.AB_y[i] = -Q.AB_y[i];
833         Q_tmp.AB_z[i] = -Q.AB_z[i];
834     }
835 
836     int ret = ostei_g_f_d_p(P_tmp, Q_tmp, screen_tol, work, INT__f_g_p_d);
837     double buffer[2700] SIMINT_ALIGN_ARRAY_DBL;
838 
839     for(int q = 0; q < ret; q++)
840     {
841         int idx = 0;
842         for(int a = 0; a < 10; ++a)
843         for(int b = 0; b < 15; ++b)
844         for(int c = 0; c < 3; ++c)
845         for(int d = 0; d < 6; ++d)
846             buffer[idx++] = INT__f_g_p_d[q*2700+b*180+a*18+d*3+c];
847 
848         memcpy(INT__f_g_p_d+q*2700, buffer, 2700*sizeof(double));
849     }
850 
851     return ret;
852 }
853 
854