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_i_i_g_s(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_i_g_s)8 int ostei_i_i_g_s(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__i_i_g_s)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__i_i_g_s);
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 
27     // partition workspace
28     double * const INT__i_s_g_s = work + (SIMINT_NSHELL_SIMD * 0);
29     double * const INT__k_s_g_s = work + (SIMINT_NSHELL_SIMD * 420);
30     double * const INT__l_s_g_s = work + (SIMINT_NSHELL_SIMD * 960);
31     double * const INT__m_s_g_s = work + (SIMINT_NSHELL_SIMD * 1635);
32     double * const INT__n_s_g_s = work + (SIMINT_NSHELL_SIMD * 2460);
33     double * const INT__o_s_g_s = work + (SIMINT_NSHELL_SIMD * 3450);
34     double * const INT__q_s_g_s = work + (SIMINT_NSHELL_SIMD * 4620);
35     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*5985);
36     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 17;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 65;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 155;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 295;
41     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 415;
42     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 610;
43     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 790;
44     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 1060;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 1312;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 1564;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 1942;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 2362;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 2670;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_d_s = primwork + 3006;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_f_s = primwork + 3510;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_g_s = primwork + 4070;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 4490;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 4850;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 5282;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_f_s = primwork + 5930;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_g_s = primwork + 6650;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_s_s = primwork + 7190;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_p_s = primwork + 7595;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_d_s = primwork + 8135;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_f_s = primwork + 8945;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_g_s = primwork + 9845;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_s_s = primwork + 10520;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_p_s = primwork + 10960;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_d_s = primwork + 11620;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_f_s = primwork + 12610;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_g_s = primwork + 13710;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_s_s = primwork + 14535;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_p_s = primwork + 14997;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_d_s = primwork + 15789;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_f_s = primwork + 16977;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_g_s = primwork + 18297;
73     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_s_s = primwork + 19287;
74     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_p_s = primwork + 19755;
75     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_d_s = primwork + 20691;
76     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_f_s = primwork + 22095;
77     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_g_s = primwork + 23655;
78     SIMINT_DBLTYPE * const restrict PRIM_INT__q_s_s_s = primwork + 24825;
79     SIMINT_DBLTYPE * const restrict PRIM_INT__q_s_p_s = primwork + 25280;
80     SIMINT_DBLTYPE * const restrict PRIM_INT__q_s_d_s = primwork + 26372;
81     SIMINT_DBLTYPE * const restrict PRIM_INT__q_s_f_s = primwork + 28010;
82     SIMINT_DBLTYPE * const restrict PRIM_INT__q_s_g_s = primwork + 29830;
83     double * const hrrwork = (double *)(primwork + 31195);
84     double * const HRR_INT__i_p_g_s = hrrwork + 0;
85     double * const HRR_INT__i_d_g_s = hrrwork + 1260;
86     double * const HRR_INT__i_f_g_s = hrrwork + 3780;
87     double * const HRR_INT__i_g_g_s = hrrwork + 7980;
88     double * const HRR_INT__i_h_g_s = hrrwork + 14280;
89     double * const HRR_INT__k_p_g_s = hrrwork + 23100;
90     double * const HRR_INT__k_d_g_s = hrrwork + 24720;
91     double * const HRR_INT__k_f_g_s = hrrwork + 27960;
92     double * const HRR_INT__k_g_g_s = hrrwork + 33360;
93     double * const HRR_INT__k_h_g_s = hrrwork + 41460;
94     double * const HRR_INT__l_p_g_s = hrrwork + 52800;
95     double * const HRR_INT__l_d_g_s = hrrwork + 54825;
96     double * const HRR_INT__l_f_g_s = hrrwork + 58875;
97     double * const HRR_INT__l_g_g_s = hrrwork + 65625;
98     double * const HRR_INT__m_p_g_s = hrrwork + 75750;
99     double * const HRR_INT__m_d_g_s = hrrwork + 78225;
100     double * const HRR_INT__m_f_g_s = hrrwork + 83175;
101     double * const HRR_INT__n_p_g_s = hrrwork + 91425;
102     double * const HRR_INT__n_d_g_s = hrrwork + 94395;
103     double * const HRR_INT__o_p_g_s = hrrwork + 100335;
104 
105 
106     // Create constants
107     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
108     const SIMINT_DBLTYPE const_10 = SIMINT_DBLSET1(10);
109     const SIMINT_DBLTYPE const_11 = SIMINT_DBLSET1(11);
110     const SIMINT_DBLTYPE const_12 = SIMINT_DBLSET1(12);
111     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
112     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
113     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
114     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
115     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
116     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
117     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
118     const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
119     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
120 
121 
122     ////////////////////////////////////////
123     // Loop over shells and primitives
124     ////////////////////////////////////////
125 
126     real_abcd = 0;
127     istart = 0;
128     for(ab = 0; ab < P.nshell12_clip; ++ab)
129     {
130         const int iend = istart + P.nprim12[ab];
131 
132         cd = 0;
133         jstart = 0;
134 
135         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
136         {
137             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
138             int jend = jstart;
139             for(i = 0; i < nshellbatch; i++)
140                 jend += Q.nprim12[cd+i];
141 
142             // Clear the beginning of the workspace (where we are accumulating integrals)
143             memset(work, 0, SIMINT_NSHELL_SIMD * 5985 * sizeof(double));
144             abcd = 0;
145 
146 
147             for(i = istart; i < iend; ++i)
148             {
149                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
150 
151                 if(check_screen)
152                 {
153                     // Skip this whole thing if always insignificant
154                     if((P.screen[i] * Q.screen_max) < screen_tol)
155                         continue;
156                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
157                 }
158 
159                 icd = 0;
160                 iprimcd = 0;
161                 nprim_icd = Q.nprim12[cd];
162                 double * restrict PRIM_PTR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
163                 double * restrict PRIM_PTR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
164                 double * restrict PRIM_PTR_INT__l_s_g_s = INT__l_s_g_s + abcd * 675;
165                 double * restrict PRIM_PTR_INT__m_s_g_s = INT__m_s_g_s + abcd * 825;
166                 double * restrict PRIM_PTR_INT__n_s_g_s = INT__n_s_g_s + abcd * 990;
167                 double * restrict PRIM_PTR_INT__o_s_g_s = INT__o_s_g_s + abcd * 1170;
168                 double * restrict PRIM_PTR_INT__q_s_g_s = INT__q_s_g_s + abcd * 1365;
169 
170 
171 
172                 // Load these one per loop over i
173                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
174                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
175                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
176 
177                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
178 
179                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
180                 {
181                     // calculate the shell offsets
182                     // these are the offset from the shell pointed to by cd
183                     // for each element
184                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
185                     int lastoffset = 0;
186                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
187 
188                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
189                     {
190                         // Handle if the first element of the vector is a new shell
191                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
192                         {
193                             nprim_icd += Q.nprim12[cd + (++icd)];
194                             PRIM_PTR_INT__i_s_g_s += 420;
195                             PRIM_PTR_INT__k_s_g_s += 540;
196                             PRIM_PTR_INT__l_s_g_s += 675;
197                             PRIM_PTR_INT__m_s_g_s += 825;
198                             PRIM_PTR_INT__n_s_g_s += 990;
199                             PRIM_PTR_INT__o_s_g_s += 1170;
200                             PRIM_PTR_INT__q_s_g_s += 1365;
201                         }
202                         iprimcd++;
203                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
204                         {
205                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
206                             {
207                                 shelloffsets[n] = shelloffsets[n-1] + 1;
208                                 lastoffset++;
209                                 nprim_icd += Q.nprim12[cd + (++icd)];
210                             }
211                             else
212                                 shelloffsets[n] = shelloffsets[n-1];
213                             iprimcd++;
214                         }
215                     }
216                     else
217                         iprimcd += SIMINT_SIMD_LEN;
218 
219                     // Do we have to compute this vector (or has it been screened out)?
220                     // (not_screened != 0 means we have to do this vector)
221                     if(check_screen)
222                     {
223                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
224                         if(vmax < screen_tol)
225                         {
226                             PRIM_PTR_INT__i_s_g_s += lastoffset*420;
227                             PRIM_PTR_INT__k_s_g_s += lastoffset*540;
228                             PRIM_PTR_INT__l_s_g_s += lastoffset*675;
229                             PRIM_PTR_INT__m_s_g_s += lastoffset*825;
230                             PRIM_PTR_INT__n_s_g_s += lastoffset*990;
231                             PRIM_PTR_INT__o_s_g_s += lastoffset*1170;
232                             PRIM_PTR_INT__q_s_g_s += lastoffset*1365;
233                             continue;
234                         }
235                     }
236 
237                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
238                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
239                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
240                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
241 
242 
243                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
244                     SIMINT_DBLTYPE PQ[3];
245                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
246                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
247                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
248                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
249                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
250                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
251 
252                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
253                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
254                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
255                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
256                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
257                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
258                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
259 
260                     // NOTE: Minus sign!
261                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
262                     SIMINT_DBLTYPE aop_PQ[3];
263                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
264                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
265                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
266 
267                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
268                     SIMINT_DBLTYPE aoq_PQ[3];
269                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
270                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
271                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
272                     // Put a minus sign here so we don't have to in RR routines
273                     a_over_q = SIMINT_NEG(a_over_q);
274 
275 
276                     //////////////////////////////////////////////
277                     // Fjt function section
278                     // Maximum v value: 16
279                     //////////////////////////////////////////////
280                     // The parameter to the Fjt function
281                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
282 
283 
284                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
285 
286 
287                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 16);
288                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
289                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
290                     for(n = 0; n <= 16; n++)
291                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
292 
293                     //////////////////////////////////////////////
294                     // Primitive integrals: Vertical recurrance
295                     //////////////////////////////////////////////
296 
297                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
298                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
299                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
300                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
301                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
302                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
303                     const SIMINT_DBLTYPE vrr_const_7_over_2p = SIMINT_MUL(const_7, one_over_2p);
304                     const SIMINT_DBLTYPE vrr_const_8_over_2p = SIMINT_MUL(const_8, one_over_2p);
305                     const SIMINT_DBLTYPE vrr_const_9_over_2p = SIMINT_MUL(const_9, one_over_2p);
306                     const SIMINT_DBLTYPE vrr_const_10_over_2p = SIMINT_MUL(const_10, one_over_2p);
307                     const SIMINT_DBLTYPE vrr_const_11_over_2p = SIMINT_MUL(const_11, one_over_2p);
308                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
309                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
310                     const SIMINT_DBLTYPE vrr_const_3_over_2q = SIMINT_MUL(const_3, one_over_2q);
311                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
312                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
313                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
314                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
315                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
316                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
317                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
318                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
319                     const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
320                     const SIMINT_DBLTYPE vrr_const_10_over_2pq = SIMINT_MUL(const_10, one_over_2pq);
321                     const SIMINT_DBLTYPE vrr_const_11_over_2pq = SIMINT_MUL(const_11, one_over_2pq);
322                     const SIMINT_DBLTYPE vrr_const_12_over_2pq = SIMINT_MUL(const_12, one_over_2pq);
323 
324 
325 
326                     // Forming PRIM_INT__p_s_s_s[16 * 3];
327                     for(n = 0; n < 16; ++n)  // loop over orders of auxiliary function
328                     {
329 
330                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
331                         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]);
332 
333                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
334                         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]);
335 
336                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
337                         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]);
338 
339                     }
340 
341 
342 
343                     // Forming PRIM_INT__d_s_s_s[15 * 6];
344                     for(n = 0; n < 15; ++n)  // loop over orders of auxiliary function
345                     {
346 
347                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
348                         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]);
349                         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]);
350 
351                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
352                         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]);
353 
354                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
355                         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]);
356 
357                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
358                         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]);
359                         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]);
360 
361                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
362                         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]);
363 
364                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
365                         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]);
366                         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]);
367 
368                     }
369 
370 
371 
372                     // Forming PRIM_INT__f_s_s_s[14 * 10];
373                     for(n = 0; n < 14; ++n)  // loop over orders of auxiliary function
374                     {
375 
376                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
377                         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]);
378                         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]);
379 
380                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
381                         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]);
382 
383                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
384                         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]);
385 
386                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
387                         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]);
388 
389                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
390                         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]);
391 
392                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
393                         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]);
394 
395                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
396                         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]);
397                         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]);
398 
399                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
400                         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]);
401 
402                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
403                         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]);
404 
405                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
406                         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]);
407                         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]);
408 
409                     }
410 
411 
412                     VRR_I_g_s_s_s(
413                             PRIM_INT__g_s_s_s,
414                             PRIM_INT__f_s_s_s,
415                             PRIM_INT__d_s_s_s,
416                             P_PA,
417                             a_over_p,
418                             aop_PQ,
419                             one_over_2p,
420                             13);
421 
422 
423                     VRR_I_h_s_s_s(
424                             PRIM_INT__h_s_s_s,
425                             PRIM_INT__g_s_s_s,
426                             PRIM_INT__f_s_s_s,
427                             P_PA,
428                             a_over_p,
429                             aop_PQ,
430                             one_over_2p,
431                             12);
432 
433 
434                     ostei_general_vrr1_I(6, 11,
435                             one_over_2p, a_over_p, aop_PQ, P_PA,
436                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
437 
438 
439                     ostei_general_vrr_K(6, 0, 1, 0, 4,
440                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
441                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
442 
443 
444                     ostei_general_vrr_K(5, 0, 1, 0, 4,
445                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
446                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
447 
448 
449                     ostei_general_vrr_K(6, 0, 2, 0, 3,
450                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
451                             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);
452 
453 
454                     VRR_K_g_s_p_s(
455                             PRIM_INT__g_s_p_s,
456                             PRIM_INT__g_s_s_s,
457                             PRIM_INT__f_s_s_s,
458                             Q_PA,
459                             aoq_PQ,
460                             one_over_2pq,
461                             4);
462 
463 
464                     ostei_general_vrr_K(5, 0, 2, 0, 3,
465                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
466                             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);
467 
468 
469                     ostei_general_vrr_K(6, 0, 3, 0, 2,
470                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
471                             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);
472 
473 
474                     VRR_K_f_s_p_s(
475                             PRIM_INT__f_s_p_s,
476                             PRIM_INT__f_s_s_s,
477                             PRIM_INT__d_s_s_s,
478                             Q_PA,
479                             aoq_PQ,
480                             one_over_2pq,
481                             4);
482 
483 
484                     ostei_general_vrr_K(4, 0, 2, 0, 3,
485                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
486                             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);
487 
488 
489                     ostei_general_vrr_K(5, 0, 3, 0, 2,
490                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
491                             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);
492 
493 
494                     ostei_general_vrr_K(6, 0, 4, 0, 1,
495                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
496                             PRIM_INT__i_s_f_s, PRIM_INT__i_s_d_s, NULL, PRIM_INT__h_s_f_s, NULL, PRIM_INT__i_s_g_s);
497 
498 
499                     ostei_general_vrr1_I(7, 10,
500                             one_over_2p, a_over_p, aop_PQ, P_PA,
501                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
502 
503 
504                     ostei_general_vrr_K(7, 0, 1, 0, 4,
505                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
506                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
507 
508 
509                     ostei_general_vrr_K(7, 0, 2, 0, 3,
510                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
511                             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);
512 
513 
514                     ostei_general_vrr_K(7, 0, 3, 0, 2,
515                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
516                             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);
517 
518 
519                     ostei_general_vrr_K(7, 0, 4, 0, 1,
520                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
521                             PRIM_INT__k_s_f_s, PRIM_INT__k_s_d_s, NULL, PRIM_INT__i_s_f_s, NULL, PRIM_INT__k_s_g_s);
522 
523 
524                     ostei_general_vrr1_I(8, 9,
525                             one_over_2p, a_over_p, aop_PQ, P_PA,
526                             PRIM_INT__k_s_s_s, PRIM_INT__i_s_s_s, PRIM_INT__l_s_s_s);
527 
528 
529                     ostei_general_vrr_K(8, 0, 1, 0, 4,
530                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
531                             PRIM_INT__l_s_s_s, NULL, NULL, PRIM_INT__k_s_s_s, NULL, PRIM_INT__l_s_p_s);
532 
533 
534                     ostei_general_vrr_K(8, 0, 2, 0, 3,
535                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
536                             PRIM_INT__l_s_p_s, PRIM_INT__l_s_s_s, NULL, PRIM_INT__k_s_p_s, NULL, PRIM_INT__l_s_d_s);
537 
538 
539                     ostei_general_vrr_K(8, 0, 3, 0, 2,
540                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
541                             PRIM_INT__l_s_d_s, PRIM_INT__l_s_p_s, NULL, PRIM_INT__k_s_d_s, NULL, PRIM_INT__l_s_f_s);
542 
543 
544                     ostei_general_vrr_K(8, 0, 4, 0, 1,
545                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
546                             PRIM_INT__l_s_f_s, PRIM_INT__l_s_d_s, NULL, PRIM_INT__k_s_f_s, NULL, PRIM_INT__l_s_g_s);
547 
548 
549                     ostei_general_vrr1_I(9, 8,
550                             one_over_2p, a_over_p, aop_PQ, P_PA,
551                             PRIM_INT__l_s_s_s, PRIM_INT__k_s_s_s, PRIM_INT__m_s_s_s);
552 
553 
554                     ostei_general_vrr_K(9, 0, 1, 0, 4,
555                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
556                             PRIM_INT__m_s_s_s, NULL, NULL, PRIM_INT__l_s_s_s, NULL, PRIM_INT__m_s_p_s);
557 
558 
559                     ostei_general_vrr_K(9, 0, 2, 0, 3,
560                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
561                             PRIM_INT__m_s_p_s, PRIM_INT__m_s_s_s, NULL, PRIM_INT__l_s_p_s, NULL, PRIM_INT__m_s_d_s);
562 
563 
564                     ostei_general_vrr_K(9, 0, 3, 0, 2,
565                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
566                             PRIM_INT__m_s_d_s, PRIM_INT__m_s_p_s, NULL, PRIM_INT__l_s_d_s, NULL, PRIM_INT__m_s_f_s);
567 
568 
569                     ostei_general_vrr_K(9, 0, 4, 0, 1,
570                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
571                             PRIM_INT__m_s_f_s, PRIM_INT__m_s_d_s, NULL, PRIM_INT__l_s_f_s, NULL, PRIM_INT__m_s_g_s);
572 
573 
574                     ostei_general_vrr1_I(10, 7,
575                             one_over_2p, a_over_p, aop_PQ, P_PA,
576                             PRIM_INT__m_s_s_s, PRIM_INT__l_s_s_s, PRIM_INT__n_s_s_s);
577 
578 
579                     ostei_general_vrr_K(10, 0, 1, 0, 4,
580                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
581                             PRIM_INT__n_s_s_s, NULL, NULL, PRIM_INT__m_s_s_s, NULL, PRIM_INT__n_s_p_s);
582 
583 
584                     ostei_general_vrr_K(10, 0, 2, 0, 3,
585                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
586                             PRIM_INT__n_s_p_s, PRIM_INT__n_s_s_s, NULL, PRIM_INT__m_s_p_s, NULL, PRIM_INT__n_s_d_s);
587 
588 
589                     ostei_general_vrr_K(10, 0, 3, 0, 2,
590                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
591                             PRIM_INT__n_s_d_s, PRIM_INT__n_s_p_s, NULL, PRIM_INT__m_s_d_s, NULL, PRIM_INT__n_s_f_s);
592 
593 
594                     ostei_general_vrr_K(10, 0, 4, 0, 1,
595                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
596                             PRIM_INT__n_s_f_s, PRIM_INT__n_s_d_s, NULL, PRIM_INT__m_s_f_s, NULL, PRIM_INT__n_s_g_s);
597 
598 
599                     ostei_general_vrr1_I(11, 6,
600                             one_over_2p, a_over_p, aop_PQ, P_PA,
601                             PRIM_INT__n_s_s_s, PRIM_INT__m_s_s_s, PRIM_INT__o_s_s_s);
602 
603 
604                     ostei_general_vrr_K(11, 0, 1, 0, 4,
605                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
606                             PRIM_INT__o_s_s_s, NULL, NULL, PRIM_INT__n_s_s_s, NULL, PRIM_INT__o_s_p_s);
607 
608 
609                     ostei_general_vrr_K(11, 0, 2, 0, 3,
610                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
611                             PRIM_INT__o_s_p_s, PRIM_INT__o_s_s_s, NULL, PRIM_INT__n_s_p_s, NULL, PRIM_INT__o_s_d_s);
612 
613 
614                     ostei_general_vrr_K(11, 0, 3, 0, 2,
615                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
616                             PRIM_INT__o_s_d_s, PRIM_INT__o_s_p_s, NULL, PRIM_INT__n_s_d_s, NULL, PRIM_INT__o_s_f_s);
617 
618 
619                     ostei_general_vrr_K(11, 0, 4, 0, 1,
620                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
621                             PRIM_INT__o_s_f_s, PRIM_INT__o_s_d_s, NULL, PRIM_INT__n_s_f_s, NULL, PRIM_INT__o_s_g_s);
622 
623 
624                     ostei_general_vrr1_I(12, 5,
625                             one_over_2p, a_over_p, aop_PQ, P_PA,
626                             PRIM_INT__o_s_s_s, PRIM_INT__n_s_s_s, PRIM_INT__q_s_s_s);
627 
628 
629                     ostei_general_vrr_K(12, 0, 1, 0, 4,
630                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
631                             PRIM_INT__q_s_s_s, NULL, NULL, PRIM_INT__o_s_s_s, NULL, PRIM_INT__q_s_p_s);
632 
633 
634                     ostei_general_vrr_K(12, 0, 2, 0, 3,
635                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
636                             PRIM_INT__q_s_p_s, PRIM_INT__q_s_s_s, NULL, PRIM_INT__o_s_p_s, NULL, PRIM_INT__q_s_d_s);
637 
638 
639                     ostei_general_vrr_K(12, 0, 3, 0, 2,
640                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
641                             PRIM_INT__q_s_d_s, PRIM_INT__q_s_p_s, NULL, PRIM_INT__o_s_d_s, NULL, PRIM_INT__q_s_f_s);
642 
643 
644                     ostei_general_vrr_K(12, 0, 4, 0, 1,
645                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
646                             PRIM_INT__q_s_f_s, PRIM_INT__q_s_d_s, NULL, PRIM_INT__o_s_f_s, NULL, PRIM_INT__q_s_g_s);
647 
648 
649 
650 
651                     ////////////////////////////////////
652                     // Accumulate contracted integrals
653                     ////////////////////////////////////
654                     if(lastoffset == 0)
655                     {
656                         contract_all(420, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
657                         contract_all(540, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
658                         contract_all(675, PRIM_INT__l_s_g_s, PRIM_PTR_INT__l_s_g_s);
659                         contract_all(825, PRIM_INT__m_s_g_s, PRIM_PTR_INT__m_s_g_s);
660                         contract_all(990, PRIM_INT__n_s_g_s, PRIM_PTR_INT__n_s_g_s);
661                         contract_all(1170, PRIM_INT__o_s_g_s, PRIM_PTR_INT__o_s_g_s);
662                         contract_all(1365, PRIM_INT__q_s_g_s, PRIM_PTR_INT__q_s_g_s);
663                     }
664                     else
665                     {
666                         contract(420, shelloffsets, PRIM_INT__i_s_g_s, PRIM_PTR_INT__i_s_g_s);
667                         contract(540, shelloffsets, PRIM_INT__k_s_g_s, PRIM_PTR_INT__k_s_g_s);
668                         contract(675, shelloffsets, PRIM_INT__l_s_g_s, PRIM_PTR_INT__l_s_g_s);
669                         contract(825, shelloffsets, PRIM_INT__m_s_g_s, PRIM_PTR_INT__m_s_g_s);
670                         contract(990, shelloffsets, PRIM_INT__n_s_g_s, PRIM_PTR_INT__n_s_g_s);
671                         contract(1170, shelloffsets, PRIM_INT__o_s_g_s, PRIM_PTR_INT__o_s_g_s);
672                         contract(1365, shelloffsets, PRIM_INT__q_s_g_s, PRIM_PTR_INT__q_s_g_s);
673                         PRIM_PTR_INT__i_s_g_s += lastoffset*420;
674                         PRIM_PTR_INT__k_s_g_s += lastoffset*540;
675                         PRIM_PTR_INT__l_s_g_s += lastoffset*675;
676                         PRIM_PTR_INT__m_s_g_s += lastoffset*825;
677                         PRIM_PTR_INT__n_s_g_s += lastoffset*990;
678                         PRIM_PTR_INT__o_s_g_s += lastoffset*1170;
679                         PRIM_PTR_INT__q_s_g_s += lastoffset*1365;
680                     }
681 
682                 }  // close loop over j
683             }  // close loop over i
684 
685             //Advance to the next batch
686             jstart = SIMINT_SIMD_ROUND(jend);
687 
688             //////////////////////////////////////////////
689             // Contracted integrals: Horizontal recurrance
690             //////////////////////////////////////////////
691 
692 
693             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
694 
695 
696             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
697             {
698 
699                 // set up HRR pointers
700                 double const * restrict HRR_INT__i_s_g_s = INT__i_s_g_s + abcd * 420;
701                 double const * restrict HRR_INT__k_s_g_s = INT__k_s_g_s + abcd * 540;
702                 double const * restrict HRR_INT__l_s_g_s = INT__l_s_g_s + abcd * 675;
703                 double const * restrict HRR_INT__m_s_g_s = INT__m_s_g_s + abcd * 825;
704                 double const * restrict HRR_INT__n_s_g_s = INT__n_s_g_s + abcd * 990;
705                 double const * restrict HRR_INT__o_s_g_s = INT__o_s_g_s + abcd * 1170;
706                 double const * restrict HRR_INT__q_s_g_s = INT__q_s_g_s + abcd * 1365;
707                 double * restrict HRR_INT__i_i_g_s = INT__i_i_g_s + real_abcd * 11760;
708 
709                 // form INT__i_p_g_s
710                 ostei_general_hrr_J(6, 1, 4, 0, hAB, HRR_INT__k_s_g_s, HRR_INT__i_s_g_s, HRR_INT__i_p_g_s);
711 
712                 // form INT__k_p_g_s
713                 ostei_general_hrr_J(7, 1, 4, 0, hAB, HRR_INT__l_s_g_s, HRR_INT__k_s_g_s, HRR_INT__k_p_g_s);
714 
715                 // form INT__l_p_g_s
716                 ostei_general_hrr_J(8, 1, 4, 0, hAB, HRR_INT__m_s_g_s, HRR_INT__l_s_g_s, HRR_INT__l_p_g_s);
717 
718                 // form INT__m_p_g_s
719                 ostei_general_hrr_J(9, 1, 4, 0, hAB, HRR_INT__n_s_g_s, HRR_INT__m_s_g_s, HRR_INT__m_p_g_s);
720 
721                 // form INT__n_p_g_s
722                 ostei_general_hrr_J(10, 1, 4, 0, hAB, HRR_INT__o_s_g_s, HRR_INT__n_s_g_s, HRR_INT__n_p_g_s);
723 
724                 // form INT__o_p_g_s
725                 ostei_general_hrr_J(11, 1, 4, 0, hAB, HRR_INT__q_s_g_s, HRR_INT__o_s_g_s, HRR_INT__o_p_g_s);
726 
727                 // form INT__i_d_g_s
728                 ostei_general_hrr_J(6, 2, 4, 0, hAB, HRR_INT__k_p_g_s, HRR_INT__i_p_g_s, HRR_INT__i_d_g_s);
729 
730                 // form INT__k_d_g_s
731                 ostei_general_hrr_J(7, 2, 4, 0, hAB, HRR_INT__l_p_g_s, HRR_INT__k_p_g_s, HRR_INT__k_d_g_s);
732 
733                 // form INT__l_d_g_s
734                 ostei_general_hrr_J(8, 2, 4, 0, hAB, HRR_INT__m_p_g_s, HRR_INT__l_p_g_s, HRR_INT__l_d_g_s);
735 
736                 // form INT__m_d_g_s
737                 ostei_general_hrr_J(9, 2, 4, 0, hAB, HRR_INT__n_p_g_s, HRR_INT__m_p_g_s, HRR_INT__m_d_g_s);
738 
739                 // form INT__n_d_g_s
740                 ostei_general_hrr_J(10, 2, 4, 0, hAB, HRR_INT__o_p_g_s, HRR_INT__n_p_g_s, HRR_INT__n_d_g_s);
741 
742                 // form INT__i_f_g_s
743                 ostei_general_hrr_J(6, 3, 4, 0, hAB, HRR_INT__k_d_g_s, HRR_INT__i_d_g_s, HRR_INT__i_f_g_s);
744 
745                 // form INT__k_f_g_s
746                 ostei_general_hrr_J(7, 3, 4, 0, hAB, HRR_INT__l_d_g_s, HRR_INT__k_d_g_s, HRR_INT__k_f_g_s);
747 
748                 // form INT__l_f_g_s
749                 ostei_general_hrr_J(8, 3, 4, 0, hAB, HRR_INT__m_d_g_s, HRR_INT__l_d_g_s, HRR_INT__l_f_g_s);
750 
751                 // form INT__m_f_g_s
752                 ostei_general_hrr_J(9, 3, 4, 0, hAB, HRR_INT__n_d_g_s, HRR_INT__m_d_g_s, HRR_INT__m_f_g_s);
753 
754                 // form INT__i_g_g_s
755                 ostei_general_hrr_J(6, 4, 4, 0, hAB, HRR_INT__k_f_g_s, HRR_INT__i_f_g_s, HRR_INT__i_g_g_s);
756 
757                 // form INT__k_g_g_s
758                 ostei_general_hrr_J(7, 4, 4, 0, hAB, HRR_INT__l_f_g_s, HRR_INT__k_f_g_s, HRR_INT__k_g_g_s);
759 
760                 // form INT__l_g_g_s
761                 ostei_general_hrr_J(8, 4, 4, 0, hAB, HRR_INT__m_f_g_s, HRR_INT__l_f_g_s, HRR_INT__l_g_g_s);
762 
763                 // form INT__i_h_g_s
764                 ostei_general_hrr_J(6, 5, 4, 0, hAB, HRR_INT__k_g_g_s, HRR_INT__i_g_g_s, HRR_INT__i_h_g_s);
765 
766                 // form INT__k_h_g_s
767                 ostei_general_hrr_J(7, 5, 4, 0, hAB, HRR_INT__l_g_g_s, HRR_INT__k_g_g_s, HRR_INT__k_h_g_s);
768 
769                 // form INT__i_i_g_s
770                 ostei_general_hrr_J(6, 6, 4, 0, hAB, HRR_INT__k_h_g_s, HRR_INT__i_h_g_s, HRR_INT__i_i_g_s);
771 
772 
773             }  // close HRR loop
774 
775 
776         }   // close loop cdbatch
777 
778         istart = iend;
779     }  // close loop over ab
780 
781     return P.nshell12_clip * Q.nshell12_clip;
782 }
783 
ostei_i_i_s_g(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_i_s_g)784 int ostei_i_i_s_g(struct simint_multi_shellpair const P,
785                   struct simint_multi_shellpair const Q,
786                   double screen_tol,
787                   double * const restrict work,
788                   double * const restrict INT__i_i_s_g)
789 {
790     double Q_AB[3*Q.nshell12];
791     struct simint_multi_shellpair Q_tmp = Q;
792     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
793     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
794     Q_tmp.AB_x = Q_AB;
795     Q_tmp.AB_y = Q_AB + Q.nshell12;
796     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
797 
798     for(int i = 0; i < Q.nshell12; i++)
799     {
800         Q_tmp.AB_x[i] = -Q.AB_x[i];
801         Q_tmp.AB_y[i] = -Q.AB_y[i];
802         Q_tmp.AB_z[i] = -Q.AB_z[i];
803     }
804 
805     int ret = ostei_i_i_g_s(P, Q_tmp, screen_tol, work, INT__i_i_s_g);
806     double buffer[11760] SIMINT_ALIGN_ARRAY_DBL;
807 
808     for(int q = 0; q < ret; q++)
809     {
810         int idx = 0;
811         for(int a = 0; a < 28; ++a)
812         for(int b = 0; b < 28; ++b)
813         for(int c = 0; c < 1; ++c)
814         for(int d = 0; d < 15; ++d)
815             buffer[idx++] = INT__i_i_s_g[q*11760+a*420+b*15+d*1+c];
816 
817         memcpy(INT__i_i_s_g+q*11760, buffer, 11760*sizeof(double));
818     }
819 
820     return ret;
821 }
822 
823