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_k_i_p_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__k_i_p_p)8 int ostei_k_i_p_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__k_i_p_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__k_i_p_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__k_s_p_s = work + (SIMINT_NSHELL_SIMD * 0);
30     double * const INT__k_s_d_s = work + (SIMINT_NSHELL_SIMD * 108);
31     double * const INT__l_s_p_s = work + (SIMINT_NSHELL_SIMD * 324);
32     double * const INT__l_s_d_s = work + (SIMINT_NSHELL_SIMD * 459);
33     double * const INT__m_s_p_s = work + (SIMINT_NSHELL_SIMD * 729);
34     double * const INT__m_s_d_s = work + (SIMINT_NSHELL_SIMD * 894);
35     double * const INT__n_s_p_s = work + (SIMINT_NSHELL_SIMD * 1224);
36     double * const INT__n_s_d_s = work + (SIMINT_NSHELL_SIMD * 1422);
37     double * const INT__o_s_p_s = work + (SIMINT_NSHELL_SIMD * 1818);
38     double * const INT__o_s_d_s = work + (SIMINT_NSHELL_SIMD * 2052);
39     double * const INT__q_s_p_s = work + (SIMINT_NSHELL_SIMD * 2520);
40     double * const INT__q_s_d_s = work + (SIMINT_NSHELL_SIMD * 2793);
41     double * const INT__r_s_p_s = work + (SIMINT_NSHELL_SIMD * 3339);
42     double * const INT__r_s_d_s = work + (SIMINT_NSHELL_SIMD * 3654);
43     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*4284);
44     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
45     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 16;
46     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 61;
47     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 145;
48     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 275;
49     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 455;
50     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 686;
51     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_p_s = primwork + 966;
52     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 1134;
53     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_p_s = primwork + 1458;
54     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_d_s = primwork + 1674;
55     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_s_s = primwork + 1890;
56     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_p_s = primwork + 2250;
57     SIMINT_DBLTYPE * const restrict PRIM_INT__l_s_d_s = primwork + 2520;
58     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_s_s = primwork + 2790;
59     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_p_s = primwork + 3175;
60     SIMINT_DBLTYPE * const restrict PRIM_INT__m_s_d_s = primwork + 3505;
61     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_s_s = primwork + 3835;
62     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_p_s = primwork + 4231;
63     SIMINT_DBLTYPE * const restrict PRIM_INT__n_s_d_s = primwork + 4627;
64     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_s_s = primwork + 5023;
65     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_p_s = primwork + 5413;
66     SIMINT_DBLTYPE * const restrict PRIM_INT__o_s_d_s = primwork + 5881;
67     SIMINT_DBLTYPE * const restrict PRIM_INT__q_s_s_s = primwork + 6349;
68     SIMINT_DBLTYPE * const restrict PRIM_INT__q_s_p_s = primwork + 6713;
69     SIMINT_DBLTYPE * const restrict PRIM_INT__q_s_d_s = primwork + 7259;
70     SIMINT_DBLTYPE * const restrict PRIM_INT__r_s_s_s = primwork + 7805;
71     SIMINT_DBLTYPE * const restrict PRIM_INT__r_s_p_s = primwork + 8120;
72     SIMINT_DBLTYPE * const restrict PRIM_INT__r_s_d_s = primwork + 8750;
73     double * const hrrwork = (double *)(primwork + 9380);
74     double * const HRR_INT__k_p_p_s = hrrwork + 0;
75     double * const HRR_INT__k_p_d_s = hrrwork + 324;
76     double * const HRR_INT__k_d_p_s = hrrwork + 972;
77     double * const HRR_INT__k_d_d_s = hrrwork + 1620;
78     double * const HRR_INT__k_f_p_s = hrrwork + 2916;
79     double * const HRR_INT__k_f_d_s = hrrwork + 3996;
80     double * const HRR_INT__k_g_p_s = hrrwork + 6156;
81     double * const HRR_INT__k_g_d_s = hrrwork + 7776;
82     double * const HRR_INT__k_h_p_s = hrrwork + 11016;
83     double * const HRR_INT__k_h_d_s = hrrwork + 13284;
84     double * const HRR_INT__k_i_p_s = hrrwork + 17820;
85     double * const HRR_INT__k_i_d_s = hrrwork + 20844;
86     double * const HRR_INT__l_p_p_s = hrrwork + 26892;
87     double * const HRR_INT__l_p_d_s = hrrwork + 27297;
88     double * const HRR_INT__l_d_p_s = hrrwork + 28107;
89     double * const HRR_INT__l_d_d_s = hrrwork + 28917;
90     double * const HRR_INT__l_f_p_s = hrrwork + 30537;
91     double * const HRR_INT__l_f_d_s = hrrwork + 31887;
92     double * const HRR_INT__l_g_p_s = hrrwork + 34587;
93     double * const HRR_INT__l_g_d_s = hrrwork + 36612;
94     double * const HRR_INT__l_h_p_s = hrrwork + 40662;
95     double * const HRR_INT__l_h_d_s = hrrwork + 43497;
96     double * const HRR_INT__m_p_p_s = hrrwork + 49167;
97     double * const HRR_INT__m_p_d_s = hrrwork + 49662;
98     double * const HRR_INT__m_d_p_s = hrrwork + 50652;
99     double * const HRR_INT__m_d_d_s = hrrwork + 51642;
100     double * const HRR_INT__m_f_p_s = hrrwork + 53622;
101     double * const HRR_INT__m_f_d_s = hrrwork + 55272;
102     double * const HRR_INT__m_g_p_s = hrrwork + 58572;
103     double * const HRR_INT__m_g_d_s = hrrwork + 61047;
104     double * const HRR_INT__n_p_p_s = hrrwork + 65997;
105     double * const HRR_INT__n_p_d_s = hrrwork + 66591;
106     double * const HRR_INT__n_d_p_s = hrrwork + 67779;
107     double * const HRR_INT__n_d_d_s = hrrwork + 68967;
108     double * const HRR_INT__n_f_p_s = hrrwork + 71343;
109     double * const HRR_INT__n_f_d_s = hrrwork + 73323;
110     double * const HRR_INT__o_p_p_s = hrrwork + 77283;
111     double * const HRR_INT__o_p_d_s = hrrwork + 77985;
112     double * const HRR_INT__o_d_p_s = hrrwork + 79389;
113     double * const HRR_INT__o_d_d_s = hrrwork + 80793;
114     double * const HRR_INT__q_p_p_s = hrrwork + 83601;
115     double * const HRR_INT__q_p_d_s = hrrwork + 84420;
116 
117 
118     // Create constants
119     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
120     const SIMINT_DBLTYPE const_10 = SIMINT_DBLSET1(10);
121     const SIMINT_DBLTYPE const_11 = SIMINT_DBLSET1(11);
122     const SIMINT_DBLTYPE const_12 = SIMINT_DBLSET1(12);
123     const SIMINT_DBLTYPE const_13 = SIMINT_DBLSET1(13);
124     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
125     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
126     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
127     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
128     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
129     const SIMINT_DBLTYPE const_7 = SIMINT_DBLSET1(7);
130     const SIMINT_DBLTYPE const_8 = SIMINT_DBLSET1(8);
131     const SIMINT_DBLTYPE const_9 = SIMINT_DBLSET1(9);
132     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
133 
134 
135     ////////////////////////////////////////
136     // Loop over shells and primitives
137     ////////////////////////////////////////
138 
139     real_abcd = 0;
140     istart = 0;
141     for(ab = 0; ab < P.nshell12_clip; ++ab)
142     {
143         const int iend = istart + P.nprim12[ab];
144 
145         cd = 0;
146         jstart = 0;
147 
148         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
149         {
150             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
151             int jend = jstart;
152             for(i = 0; i < nshellbatch; i++)
153                 jend += Q.nprim12[cd+i];
154 
155             // Clear the beginning of the workspace (where we are accumulating integrals)
156             memset(work, 0, SIMINT_NSHELL_SIMD * 4284 * sizeof(double));
157             abcd = 0;
158 
159 
160             for(i = istart; i < iend; ++i)
161             {
162                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
163 
164                 if(check_screen)
165                 {
166                     // Skip this whole thing if always insignificant
167                     if((P.screen[i] * Q.screen_max) < screen_tol)
168                         continue;
169                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
170                 }
171 
172                 icd = 0;
173                 iprimcd = 0;
174                 nprim_icd = Q.nprim12[cd];
175                 double * restrict PRIM_PTR_INT__k_s_p_s = INT__k_s_p_s + abcd * 108;
176                 double * restrict PRIM_PTR_INT__k_s_d_s = INT__k_s_d_s + abcd * 216;
177                 double * restrict PRIM_PTR_INT__l_s_p_s = INT__l_s_p_s + abcd * 135;
178                 double * restrict PRIM_PTR_INT__l_s_d_s = INT__l_s_d_s + abcd * 270;
179                 double * restrict PRIM_PTR_INT__m_s_p_s = INT__m_s_p_s + abcd * 165;
180                 double * restrict PRIM_PTR_INT__m_s_d_s = INT__m_s_d_s + abcd * 330;
181                 double * restrict PRIM_PTR_INT__n_s_p_s = INT__n_s_p_s + abcd * 198;
182                 double * restrict PRIM_PTR_INT__n_s_d_s = INT__n_s_d_s + abcd * 396;
183                 double * restrict PRIM_PTR_INT__o_s_p_s = INT__o_s_p_s + abcd * 234;
184                 double * restrict PRIM_PTR_INT__o_s_d_s = INT__o_s_d_s + abcd * 468;
185                 double * restrict PRIM_PTR_INT__q_s_p_s = INT__q_s_p_s + abcd * 273;
186                 double * restrict PRIM_PTR_INT__q_s_d_s = INT__q_s_d_s + abcd * 546;
187                 double * restrict PRIM_PTR_INT__r_s_p_s = INT__r_s_p_s + abcd * 315;
188                 double * restrict PRIM_PTR_INT__r_s_d_s = INT__r_s_d_s + abcd * 630;
189 
190 
191 
192                 // Load these one per loop over i
193                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
194                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
195                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
196 
197                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
198 
199                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
200                 {
201                     // calculate the shell offsets
202                     // these are the offset from the shell pointed to by cd
203                     // for each element
204                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
205                     int lastoffset = 0;
206                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
207 
208                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
209                     {
210                         // Handle if the first element of the vector is a new shell
211                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
212                         {
213                             nprim_icd += Q.nprim12[cd + (++icd)];
214                             PRIM_PTR_INT__k_s_p_s += 108;
215                             PRIM_PTR_INT__k_s_d_s += 216;
216                             PRIM_PTR_INT__l_s_p_s += 135;
217                             PRIM_PTR_INT__l_s_d_s += 270;
218                             PRIM_PTR_INT__m_s_p_s += 165;
219                             PRIM_PTR_INT__m_s_d_s += 330;
220                             PRIM_PTR_INT__n_s_p_s += 198;
221                             PRIM_PTR_INT__n_s_d_s += 396;
222                             PRIM_PTR_INT__o_s_p_s += 234;
223                             PRIM_PTR_INT__o_s_d_s += 468;
224                             PRIM_PTR_INT__q_s_p_s += 273;
225                             PRIM_PTR_INT__q_s_d_s += 546;
226                             PRIM_PTR_INT__r_s_p_s += 315;
227                             PRIM_PTR_INT__r_s_d_s += 630;
228                         }
229                         iprimcd++;
230                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
231                         {
232                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
233                             {
234                                 shelloffsets[n] = shelloffsets[n-1] + 1;
235                                 lastoffset++;
236                                 nprim_icd += Q.nprim12[cd + (++icd)];
237                             }
238                             else
239                                 shelloffsets[n] = shelloffsets[n-1];
240                             iprimcd++;
241                         }
242                     }
243                     else
244                         iprimcd += SIMINT_SIMD_LEN;
245 
246                     // Do we have to compute this vector (or has it been screened out)?
247                     // (not_screened != 0 means we have to do this vector)
248                     if(check_screen)
249                     {
250                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
251                         if(vmax < screen_tol)
252                         {
253                             PRIM_PTR_INT__k_s_p_s += lastoffset*108;
254                             PRIM_PTR_INT__k_s_d_s += lastoffset*216;
255                             PRIM_PTR_INT__l_s_p_s += lastoffset*135;
256                             PRIM_PTR_INT__l_s_d_s += lastoffset*270;
257                             PRIM_PTR_INT__m_s_p_s += lastoffset*165;
258                             PRIM_PTR_INT__m_s_d_s += lastoffset*330;
259                             PRIM_PTR_INT__n_s_p_s += lastoffset*198;
260                             PRIM_PTR_INT__n_s_d_s += lastoffset*396;
261                             PRIM_PTR_INT__o_s_p_s += lastoffset*234;
262                             PRIM_PTR_INT__o_s_d_s += lastoffset*468;
263                             PRIM_PTR_INT__q_s_p_s += lastoffset*273;
264                             PRIM_PTR_INT__q_s_d_s += lastoffset*546;
265                             PRIM_PTR_INT__r_s_p_s += lastoffset*315;
266                             PRIM_PTR_INT__r_s_d_s += lastoffset*630;
267                             continue;
268                         }
269                     }
270 
271                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
272                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
273                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
274                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
275 
276 
277                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
278                     SIMINT_DBLTYPE PQ[3];
279                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
280                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
281                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
282                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
283                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
284                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
285 
286                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
287                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
288                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
289                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
290                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
291                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
292                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
293 
294                     // NOTE: Minus sign!
295                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
296                     SIMINT_DBLTYPE aop_PQ[3];
297                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
298                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
299                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
300 
301                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
302                     SIMINT_DBLTYPE aoq_PQ[3];
303                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
304                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
305                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
306                     // Put a minus sign here so we don't have to in RR routines
307                     a_over_q = SIMINT_NEG(a_over_q);
308 
309 
310                     //////////////////////////////////////////////
311                     // Fjt function section
312                     // Maximum v value: 15
313                     //////////////////////////////////////////////
314                     // The parameter to the Fjt function
315                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
316 
317 
318                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
319 
320 
321                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 15);
322                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
323                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
324                     for(n = 0; n <= 15; n++)
325                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
326 
327                     //////////////////////////////////////////////
328                     // Primitive integrals: Vertical recurrance
329                     //////////////////////////////////////////////
330 
331                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
332                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
333                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
334                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
335                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
336                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
337                     const SIMINT_DBLTYPE vrr_const_7_over_2p = SIMINT_MUL(const_7, one_over_2p);
338                     const SIMINT_DBLTYPE vrr_const_8_over_2p = SIMINT_MUL(const_8, one_over_2p);
339                     const SIMINT_DBLTYPE vrr_const_9_over_2p = SIMINT_MUL(const_9, one_over_2p);
340                     const SIMINT_DBLTYPE vrr_const_10_over_2p = SIMINT_MUL(const_10, one_over_2p);
341                     const SIMINT_DBLTYPE vrr_const_11_over_2p = SIMINT_MUL(const_11, one_over_2p);
342                     const SIMINT_DBLTYPE vrr_const_12_over_2p = SIMINT_MUL(const_12, one_over_2p);
343                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
344                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
345                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
346                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
347                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
348                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
349                     const SIMINT_DBLTYPE vrr_const_6_over_2pq = SIMINT_MUL(const_6, one_over_2pq);
350                     const SIMINT_DBLTYPE vrr_const_7_over_2pq = SIMINT_MUL(const_7, one_over_2pq);
351                     const SIMINT_DBLTYPE vrr_const_8_over_2pq = SIMINT_MUL(const_8, one_over_2pq);
352                     const SIMINT_DBLTYPE vrr_const_9_over_2pq = SIMINT_MUL(const_9, one_over_2pq);
353                     const SIMINT_DBLTYPE vrr_const_10_over_2pq = SIMINT_MUL(const_10, one_over_2pq);
354                     const SIMINT_DBLTYPE vrr_const_11_over_2pq = SIMINT_MUL(const_11, one_over_2pq);
355                     const SIMINT_DBLTYPE vrr_const_12_over_2pq = SIMINT_MUL(const_12, one_over_2pq);
356                     const SIMINT_DBLTYPE vrr_const_13_over_2pq = SIMINT_MUL(const_13, one_over_2pq);
357 
358 
359 
360                     // Forming PRIM_INT__p_s_s_s[15 * 3];
361                     for(n = 0; n < 15; ++n)  // loop over orders of auxiliary function
362                     {
363 
364                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
365                         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]);
366 
367                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
368                         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]);
369 
370                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
371                         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]);
372 
373                     }
374 
375 
376 
377                     // Forming PRIM_INT__d_s_s_s[14 * 6];
378                     for(n = 0; n < 14; ++n)  // loop over orders of auxiliary function
379                     {
380 
381                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
382                         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]);
383                         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]);
384 
385                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
386                         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]);
387                         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]);
388 
389                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
390                         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]);
391                         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]);
392 
393                     }
394 
395 
396 
397                     // Forming PRIM_INT__f_s_s_s[13 * 10];
398                     for(n = 0; n < 13; ++n)  // loop over orders of auxiliary function
399                     {
400 
401                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
402                         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]);
403                         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]);
404 
405                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
406                         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]);
407 
408                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
409                         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]);
410 
411                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
412                         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]);
413                         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]);
414 
415                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
416                         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]);
417 
418                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
419                         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]);
420                         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]);
421 
422                     }
423 
424 
425                     VRR_I_g_s_s_s(
426                             PRIM_INT__g_s_s_s,
427                             PRIM_INT__f_s_s_s,
428                             PRIM_INT__d_s_s_s,
429                             P_PA,
430                             a_over_p,
431                             aop_PQ,
432                             one_over_2p,
433                             12);
434 
435 
436                     VRR_I_h_s_s_s(
437                             PRIM_INT__h_s_s_s,
438                             PRIM_INT__g_s_s_s,
439                             PRIM_INT__f_s_s_s,
440                             P_PA,
441                             a_over_p,
442                             aop_PQ,
443                             one_over_2p,
444                             11);
445 
446 
447                     ostei_general_vrr1_I(6, 10,
448                             one_over_2p, a_over_p, aop_PQ, P_PA,
449                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
450 
451 
452                     ostei_general_vrr1_I(7, 9,
453                             one_over_2p, a_over_p, aop_PQ, P_PA,
454                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
455 
456 
457                     ostei_general_vrr_K(7, 0, 1, 0, 2,
458                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
459                             PRIM_INT__k_s_s_s, NULL, NULL, PRIM_INT__i_s_s_s, NULL, PRIM_INT__k_s_p_s);
460 
461 
462                     ostei_general_vrr_K(6, 0, 1, 0, 2,
463                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
464                             PRIM_INT__i_s_s_s, NULL, NULL, PRIM_INT__h_s_s_s, NULL, PRIM_INT__i_s_p_s);
465 
466 
467                     ostei_general_vrr_K(7, 0, 2, 0, 1,
468                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
469                             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);
470 
471 
472                     ostei_general_vrr1_I(8, 8,
473                             one_over_2p, a_over_p, aop_PQ, P_PA,
474                             PRIM_INT__k_s_s_s, PRIM_INT__i_s_s_s, PRIM_INT__l_s_s_s);
475 
476 
477                     ostei_general_vrr_K(8, 0, 1, 0, 2,
478                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
479                             PRIM_INT__l_s_s_s, NULL, NULL, PRIM_INT__k_s_s_s, NULL, PRIM_INT__l_s_p_s);
480 
481 
482                     ostei_general_vrr_K(8, 0, 2, 0, 1,
483                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
484                             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);
485 
486 
487                     ostei_general_vrr1_I(9, 7,
488                             one_over_2p, a_over_p, aop_PQ, P_PA,
489                             PRIM_INT__l_s_s_s, PRIM_INT__k_s_s_s, PRIM_INT__m_s_s_s);
490 
491 
492                     ostei_general_vrr_K(9, 0, 1, 0, 2,
493                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
494                             PRIM_INT__m_s_s_s, NULL, NULL, PRIM_INT__l_s_s_s, NULL, PRIM_INT__m_s_p_s);
495 
496 
497                     ostei_general_vrr_K(9, 0, 2, 0, 1,
498                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
499                             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);
500 
501 
502                     ostei_general_vrr1_I(10, 6,
503                             one_over_2p, a_over_p, aop_PQ, P_PA,
504                             PRIM_INT__m_s_s_s, PRIM_INT__l_s_s_s, PRIM_INT__n_s_s_s);
505 
506 
507                     ostei_general_vrr_K(10, 0, 1, 0, 2,
508                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
509                             PRIM_INT__n_s_s_s, NULL, NULL, PRIM_INT__m_s_s_s, NULL, PRIM_INT__n_s_p_s);
510 
511 
512                     ostei_general_vrr_K(10, 0, 2, 0, 1,
513                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
514                             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);
515 
516 
517                     ostei_general_vrr1_I(11, 5,
518                             one_over_2p, a_over_p, aop_PQ, P_PA,
519                             PRIM_INT__n_s_s_s, PRIM_INT__m_s_s_s, PRIM_INT__o_s_s_s);
520 
521 
522                     ostei_general_vrr_K(11, 0, 1, 0, 2,
523                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
524                             PRIM_INT__o_s_s_s, NULL, NULL, PRIM_INT__n_s_s_s, NULL, PRIM_INT__o_s_p_s);
525 
526 
527                     ostei_general_vrr_K(11, 0, 2, 0, 1,
528                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
529                             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);
530 
531 
532                     ostei_general_vrr1_I(12, 4,
533                             one_over_2p, a_over_p, aop_PQ, P_PA,
534                             PRIM_INT__o_s_s_s, PRIM_INT__n_s_s_s, PRIM_INT__q_s_s_s);
535 
536 
537                     ostei_general_vrr_K(12, 0, 1, 0, 2,
538                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
539                             PRIM_INT__q_s_s_s, NULL, NULL, PRIM_INT__o_s_s_s, NULL, PRIM_INT__q_s_p_s);
540 
541 
542                     ostei_general_vrr_K(12, 0, 2, 0, 1,
543                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
544                             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);
545 
546 
547                     ostei_general_vrr1_I(13, 3,
548                             one_over_2p, a_over_p, aop_PQ, P_PA,
549                             PRIM_INT__q_s_s_s, PRIM_INT__o_s_s_s, PRIM_INT__r_s_s_s);
550 
551 
552                     ostei_general_vrr_K(13, 0, 1, 0, 2,
553                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
554                             PRIM_INT__r_s_s_s, NULL, NULL, PRIM_INT__q_s_s_s, NULL, PRIM_INT__r_s_p_s);
555 
556 
557                     ostei_general_vrr_K(13, 0, 2, 0, 1,
558                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
559                             PRIM_INT__r_s_p_s, PRIM_INT__r_s_s_s, NULL, PRIM_INT__q_s_p_s, NULL, PRIM_INT__r_s_d_s);
560 
561 
562 
563 
564                     ////////////////////////////////////
565                     // Accumulate contracted integrals
566                     ////////////////////////////////////
567                     if(lastoffset == 0)
568                     {
569                         contract_all(108, PRIM_INT__k_s_p_s, PRIM_PTR_INT__k_s_p_s);
570                         contract_all(216, PRIM_INT__k_s_d_s, PRIM_PTR_INT__k_s_d_s);
571                         contract_all(135, PRIM_INT__l_s_p_s, PRIM_PTR_INT__l_s_p_s);
572                         contract_all(270, PRIM_INT__l_s_d_s, PRIM_PTR_INT__l_s_d_s);
573                         contract_all(165, PRIM_INT__m_s_p_s, PRIM_PTR_INT__m_s_p_s);
574                         contract_all(330, PRIM_INT__m_s_d_s, PRIM_PTR_INT__m_s_d_s);
575                         contract_all(198, PRIM_INT__n_s_p_s, PRIM_PTR_INT__n_s_p_s);
576                         contract_all(396, PRIM_INT__n_s_d_s, PRIM_PTR_INT__n_s_d_s);
577                         contract_all(234, PRIM_INT__o_s_p_s, PRIM_PTR_INT__o_s_p_s);
578                         contract_all(468, PRIM_INT__o_s_d_s, PRIM_PTR_INT__o_s_d_s);
579                         contract_all(273, PRIM_INT__q_s_p_s, PRIM_PTR_INT__q_s_p_s);
580                         contract_all(546, PRIM_INT__q_s_d_s, PRIM_PTR_INT__q_s_d_s);
581                         contract_all(315, PRIM_INT__r_s_p_s, PRIM_PTR_INT__r_s_p_s);
582                         contract_all(630, PRIM_INT__r_s_d_s, PRIM_PTR_INT__r_s_d_s);
583                     }
584                     else
585                     {
586                         contract(108, shelloffsets, PRIM_INT__k_s_p_s, PRIM_PTR_INT__k_s_p_s);
587                         contract(216, shelloffsets, PRIM_INT__k_s_d_s, PRIM_PTR_INT__k_s_d_s);
588                         contract(135, shelloffsets, PRIM_INT__l_s_p_s, PRIM_PTR_INT__l_s_p_s);
589                         contract(270, shelloffsets, PRIM_INT__l_s_d_s, PRIM_PTR_INT__l_s_d_s);
590                         contract(165, shelloffsets, PRIM_INT__m_s_p_s, PRIM_PTR_INT__m_s_p_s);
591                         contract(330, shelloffsets, PRIM_INT__m_s_d_s, PRIM_PTR_INT__m_s_d_s);
592                         contract(198, shelloffsets, PRIM_INT__n_s_p_s, PRIM_PTR_INT__n_s_p_s);
593                         contract(396, shelloffsets, PRIM_INT__n_s_d_s, PRIM_PTR_INT__n_s_d_s);
594                         contract(234, shelloffsets, PRIM_INT__o_s_p_s, PRIM_PTR_INT__o_s_p_s);
595                         contract(468, shelloffsets, PRIM_INT__o_s_d_s, PRIM_PTR_INT__o_s_d_s);
596                         contract(273, shelloffsets, PRIM_INT__q_s_p_s, PRIM_PTR_INT__q_s_p_s);
597                         contract(546, shelloffsets, PRIM_INT__q_s_d_s, PRIM_PTR_INT__q_s_d_s);
598                         contract(315, shelloffsets, PRIM_INT__r_s_p_s, PRIM_PTR_INT__r_s_p_s);
599                         contract(630, shelloffsets, PRIM_INT__r_s_d_s, PRIM_PTR_INT__r_s_d_s);
600                         PRIM_PTR_INT__k_s_p_s += lastoffset*108;
601                         PRIM_PTR_INT__k_s_d_s += lastoffset*216;
602                         PRIM_PTR_INT__l_s_p_s += lastoffset*135;
603                         PRIM_PTR_INT__l_s_d_s += lastoffset*270;
604                         PRIM_PTR_INT__m_s_p_s += lastoffset*165;
605                         PRIM_PTR_INT__m_s_d_s += lastoffset*330;
606                         PRIM_PTR_INT__n_s_p_s += lastoffset*198;
607                         PRIM_PTR_INT__n_s_d_s += lastoffset*396;
608                         PRIM_PTR_INT__o_s_p_s += lastoffset*234;
609                         PRIM_PTR_INT__o_s_d_s += lastoffset*468;
610                         PRIM_PTR_INT__q_s_p_s += lastoffset*273;
611                         PRIM_PTR_INT__q_s_d_s += lastoffset*546;
612                         PRIM_PTR_INT__r_s_p_s += lastoffset*315;
613                         PRIM_PTR_INT__r_s_d_s += lastoffset*630;
614                     }
615 
616                 }  // close loop over j
617             }  // close loop over i
618 
619             //Advance to the next batch
620             jstart = SIMINT_SIMD_ROUND(jend);
621 
622             //////////////////////////////////////////////
623             // Contracted integrals: Horizontal recurrance
624             //////////////////////////////////////////////
625 
626 
627             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
628 
629 
630             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
631             {
632                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
633 
634                 // set up HRR pointers
635                 double const * restrict HRR_INT__k_s_p_s = INT__k_s_p_s + abcd * 108;
636                 double const * restrict HRR_INT__k_s_d_s = INT__k_s_d_s + abcd * 216;
637                 double const * restrict HRR_INT__l_s_p_s = INT__l_s_p_s + abcd * 135;
638                 double const * restrict HRR_INT__l_s_d_s = INT__l_s_d_s + abcd * 270;
639                 double const * restrict HRR_INT__m_s_p_s = INT__m_s_p_s + abcd * 165;
640                 double const * restrict HRR_INT__m_s_d_s = INT__m_s_d_s + abcd * 330;
641                 double const * restrict HRR_INT__n_s_p_s = INT__n_s_p_s + abcd * 198;
642                 double const * restrict HRR_INT__n_s_d_s = INT__n_s_d_s + abcd * 396;
643                 double const * restrict HRR_INT__o_s_p_s = INT__o_s_p_s + abcd * 234;
644                 double const * restrict HRR_INT__o_s_d_s = INT__o_s_d_s + abcd * 468;
645                 double const * restrict HRR_INT__q_s_p_s = INT__q_s_p_s + abcd * 273;
646                 double const * restrict HRR_INT__q_s_d_s = INT__q_s_d_s + abcd * 546;
647                 double const * restrict HRR_INT__r_s_p_s = INT__r_s_p_s + abcd * 315;
648                 double const * restrict HRR_INT__r_s_d_s = INT__r_s_d_s + abcd * 630;
649                 double * restrict HRR_INT__k_i_p_p = INT__k_i_p_p + real_abcd * 9072;
650 
651                 // form INT__k_p_p_s
652                 ostei_general_hrr_J(7, 1, 1, 0, hAB, HRR_INT__l_s_p_s, HRR_INT__k_s_p_s, HRR_INT__k_p_p_s);
653 
654                 // form INT__k_p_d_s
655                 ostei_general_hrr_J(7, 1, 2, 0, hAB, HRR_INT__l_s_d_s, HRR_INT__k_s_d_s, HRR_INT__k_p_d_s);
656 
657                 // form INT__l_p_p_s
658                 ostei_general_hrr_J(8, 1, 1, 0, hAB, HRR_INT__m_s_p_s, HRR_INT__l_s_p_s, HRR_INT__l_p_p_s);
659 
660                 // form INT__l_p_d_s
661                 ostei_general_hrr_J(8, 1, 2, 0, hAB, HRR_INT__m_s_d_s, HRR_INT__l_s_d_s, HRR_INT__l_p_d_s);
662 
663                 // form INT__m_p_p_s
664                 ostei_general_hrr_J(9, 1, 1, 0, hAB, HRR_INT__n_s_p_s, HRR_INT__m_s_p_s, HRR_INT__m_p_p_s);
665 
666                 // form INT__m_p_d_s
667                 ostei_general_hrr_J(9, 1, 2, 0, hAB, HRR_INT__n_s_d_s, HRR_INT__m_s_d_s, HRR_INT__m_p_d_s);
668 
669                 // form INT__n_p_p_s
670                 ostei_general_hrr_J(10, 1, 1, 0, hAB, HRR_INT__o_s_p_s, HRR_INT__n_s_p_s, HRR_INT__n_p_p_s);
671 
672                 // form INT__n_p_d_s
673                 ostei_general_hrr_J(10, 1, 2, 0, hAB, HRR_INT__o_s_d_s, HRR_INT__n_s_d_s, HRR_INT__n_p_d_s);
674 
675                 // form INT__o_p_p_s
676                 ostei_general_hrr_J(11, 1, 1, 0, hAB, HRR_INT__q_s_p_s, HRR_INT__o_s_p_s, HRR_INT__o_p_p_s);
677 
678                 // form INT__o_p_d_s
679                 ostei_general_hrr_J(11, 1, 2, 0, hAB, HRR_INT__q_s_d_s, HRR_INT__o_s_d_s, HRR_INT__o_p_d_s);
680 
681                 // form INT__q_p_p_s
682                 ostei_general_hrr_J(12, 1, 1, 0, hAB, HRR_INT__r_s_p_s, HRR_INT__q_s_p_s, HRR_INT__q_p_p_s);
683 
684                 // form INT__q_p_d_s
685                 ostei_general_hrr_J(12, 1, 2, 0, hAB, HRR_INT__r_s_d_s, HRR_INT__q_s_d_s, HRR_INT__q_p_d_s);
686 
687                 // form INT__k_d_p_s
688                 ostei_general_hrr_J(7, 2, 1, 0, hAB, HRR_INT__l_p_p_s, HRR_INT__k_p_p_s, HRR_INT__k_d_p_s);
689 
690                 // form INT__k_d_d_s
691                 ostei_general_hrr_J(7, 2, 2, 0, hAB, HRR_INT__l_p_d_s, HRR_INT__k_p_d_s, HRR_INT__k_d_d_s);
692 
693                 // form INT__l_d_p_s
694                 ostei_general_hrr_J(8, 2, 1, 0, hAB, HRR_INT__m_p_p_s, HRR_INT__l_p_p_s, HRR_INT__l_d_p_s);
695 
696                 // form INT__l_d_d_s
697                 ostei_general_hrr_J(8, 2, 2, 0, hAB, HRR_INT__m_p_d_s, HRR_INT__l_p_d_s, HRR_INT__l_d_d_s);
698 
699                 // form INT__m_d_p_s
700                 ostei_general_hrr_J(9, 2, 1, 0, hAB, HRR_INT__n_p_p_s, HRR_INT__m_p_p_s, HRR_INT__m_d_p_s);
701 
702                 // form INT__m_d_d_s
703                 ostei_general_hrr_J(9, 2, 2, 0, hAB, HRR_INT__n_p_d_s, HRR_INT__m_p_d_s, HRR_INT__m_d_d_s);
704 
705                 // form INT__n_d_p_s
706                 ostei_general_hrr_J(10, 2, 1, 0, hAB, HRR_INT__o_p_p_s, HRR_INT__n_p_p_s, HRR_INT__n_d_p_s);
707 
708                 // form INT__n_d_d_s
709                 ostei_general_hrr_J(10, 2, 2, 0, hAB, HRR_INT__o_p_d_s, HRR_INT__n_p_d_s, HRR_INT__n_d_d_s);
710 
711                 // form INT__o_d_p_s
712                 ostei_general_hrr_J(11, 2, 1, 0, hAB, HRR_INT__q_p_p_s, HRR_INT__o_p_p_s, HRR_INT__o_d_p_s);
713 
714                 // form INT__o_d_d_s
715                 ostei_general_hrr_J(11, 2, 2, 0, hAB, HRR_INT__q_p_d_s, HRR_INT__o_p_d_s, HRR_INT__o_d_d_s);
716 
717                 // form INT__k_f_p_s
718                 ostei_general_hrr_J(7, 3, 1, 0, hAB, HRR_INT__l_d_p_s, HRR_INT__k_d_p_s, HRR_INT__k_f_p_s);
719 
720                 // form INT__k_f_d_s
721                 ostei_general_hrr_J(7, 3, 2, 0, hAB, HRR_INT__l_d_d_s, HRR_INT__k_d_d_s, HRR_INT__k_f_d_s);
722 
723                 // form INT__l_f_p_s
724                 ostei_general_hrr_J(8, 3, 1, 0, hAB, HRR_INT__m_d_p_s, HRR_INT__l_d_p_s, HRR_INT__l_f_p_s);
725 
726                 // form INT__l_f_d_s
727                 ostei_general_hrr_J(8, 3, 2, 0, hAB, HRR_INT__m_d_d_s, HRR_INT__l_d_d_s, HRR_INT__l_f_d_s);
728 
729                 // form INT__m_f_p_s
730                 ostei_general_hrr_J(9, 3, 1, 0, hAB, HRR_INT__n_d_p_s, HRR_INT__m_d_p_s, HRR_INT__m_f_p_s);
731 
732                 // form INT__m_f_d_s
733                 ostei_general_hrr_J(9, 3, 2, 0, hAB, HRR_INT__n_d_d_s, HRR_INT__m_d_d_s, HRR_INT__m_f_d_s);
734 
735                 // form INT__n_f_p_s
736                 ostei_general_hrr_J(10, 3, 1, 0, hAB, HRR_INT__o_d_p_s, HRR_INT__n_d_p_s, HRR_INT__n_f_p_s);
737 
738                 // form INT__n_f_d_s
739                 ostei_general_hrr_J(10, 3, 2, 0, hAB, HRR_INT__o_d_d_s, HRR_INT__n_d_d_s, HRR_INT__n_f_d_s);
740 
741                 // form INT__k_g_p_s
742                 ostei_general_hrr_J(7, 4, 1, 0, hAB, HRR_INT__l_f_p_s, HRR_INT__k_f_p_s, HRR_INT__k_g_p_s);
743 
744                 // form INT__k_g_d_s
745                 ostei_general_hrr_J(7, 4, 2, 0, hAB, HRR_INT__l_f_d_s, HRR_INT__k_f_d_s, HRR_INT__k_g_d_s);
746 
747                 // form INT__l_g_p_s
748                 ostei_general_hrr_J(8, 4, 1, 0, hAB, HRR_INT__m_f_p_s, HRR_INT__l_f_p_s, HRR_INT__l_g_p_s);
749 
750                 // form INT__l_g_d_s
751                 ostei_general_hrr_J(8, 4, 2, 0, hAB, HRR_INT__m_f_d_s, HRR_INT__l_f_d_s, HRR_INT__l_g_d_s);
752 
753                 // form INT__m_g_p_s
754                 ostei_general_hrr_J(9, 4, 1, 0, hAB, HRR_INT__n_f_p_s, HRR_INT__m_f_p_s, HRR_INT__m_g_p_s);
755 
756                 // form INT__m_g_d_s
757                 ostei_general_hrr_J(9, 4, 2, 0, hAB, HRR_INT__n_f_d_s, HRR_INT__m_f_d_s, HRR_INT__m_g_d_s);
758 
759                 // form INT__k_h_p_s
760                 ostei_general_hrr_J(7, 5, 1, 0, hAB, HRR_INT__l_g_p_s, HRR_INT__k_g_p_s, HRR_INT__k_h_p_s);
761 
762                 // form INT__k_h_d_s
763                 ostei_general_hrr_J(7, 5, 2, 0, hAB, HRR_INT__l_g_d_s, HRR_INT__k_g_d_s, HRR_INT__k_h_d_s);
764 
765                 // form INT__l_h_p_s
766                 ostei_general_hrr_J(8, 5, 1, 0, hAB, HRR_INT__m_g_p_s, HRR_INT__l_g_p_s, HRR_INT__l_h_p_s);
767 
768                 // form INT__l_h_d_s
769                 ostei_general_hrr_J(8, 5, 2, 0, hAB, HRR_INT__m_g_d_s, HRR_INT__l_g_d_s, HRR_INT__l_h_d_s);
770 
771                 // form INT__k_i_p_s
772                 ostei_general_hrr_J(7, 6, 1, 0, hAB, HRR_INT__l_h_p_s, HRR_INT__k_h_p_s, HRR_INT__k_i_p_s);
773 
774                 // form INT__k_i_d_s
775                 ostei_general_hrr_J(7, 6, 2, 0, hAB, HRR_INT__l_h_d_s, HRR_INT__k_h_d_s, HRR_INT__k_i_d_s);
776 
777                 // form INT__k_i_p_p
778                 for(ibra = 0; ibra < 1008; ++ibra)
779                 {
780                     HRR_INT__k_i_p_p[ibra * 9 + 0] = HRR_INT__k_i_d_s[ibra * 6 + 0] + ( hCD[0] * HRR_INT__k_i_p_s[ibra * 3 + 0] );
781 
782                     HRR_INT__k_i_p_p[ibra * 9 + 1] = HRR_INT__k_i_d_s[ibra * 6 + 1] + ( hCD[1] * HRR_INT__k_i_p_s[ibra * 3 + 0] );
783 
784                     HRR_INT__k_i_p_p[ibra * 9 + 2] = HRR_INT__k_i_d_s[ibra * 6 + 2] + ( hCD[2] * HRR_INT__k_i_p_s[ibra * 3 + 0] );
785 
786                     HRR_INT__k_i_p_p[ibra * 9 + 3] = HRR_INT__k_i_d_s[ibra * 6 + 1] + ( hCD[0] * HRR_INT__k_i_p_s[ibra * 3 + 1] );
787 
788                     HRR_INT__k_i_p_p[ibra * 9 + 4] = HRR_INT__k_i_d_s[ibra * 6 + 3] + ( hCD[1] * HRR_INT__k_i_p_s[ibra * 3 + 1] );
789 
790                     HRR_INT__k_i_p_p[ibra * 9 + 5] = HRR_INT__k_i_d_s[ibra * 6 + 4] + ( hCD[2] * HRR_INT__k_i_p_s[ibra * 3 + 1] );
791 
792                     HRR_INT__k_i_p_p[ibra * 9 + 6] = HRR_INT__k_i_d_s[ibra * 6 + 2] + ( hCD[0] * HRR_INT__k_i_p_s[ibra * 3 + 2] );
793 
794                     HRR_INT__k_i_p_p[ibra * 9 + 7] = HRR_INT__k_i_d_s[ibra * 6 + 4] + ( hCD[1] * HRR_INT__k_i_p_s[ibra * 3 + 2] );
795 
796                     HRR_INT__k_i_p_p[ibra * 9 + 8] = HRR_INT__k_i_d_s[ibra * 6 + 5] + ( hCD[2] * HRR_INT__k_i_p_s[ibra * 3 + 2] );
797 
798                 }
799 
800 
801             }  // close HRR loop
802 
803 
804         }   // close loop cdbatch
805 
806         istart = iend;
807     }  // close loop over ab
808 
809     return P.nshell12_clip * Q.nshell12_clip;
810 }
811 
ostei_i_k_p_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__i_k_p_p)812 int ostei_i_k_p_p(struct simint_multi_shellpair const P,
813                   struct simint_multi_shellpair const Q,
814                   double screen_tol,
815                   double * const restrict work,
816                   double * const restrict INT__i_k_p_p)
817 {
818     double P_AB[3*P.nshell12];
819     struct simint_multi_shellpair P_tmp = P;
820     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
821     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
822     P_tmp.AB_x = P_AB;
823     P_tmp.AB_y = P_AB + P.nshell12;
824     P_tmp.AB_z = P_AB + 2*P.nshell12;
825 
826     for(int i = 0; i < P.nshell12; i++)
827     {
828         P_tmp.AB_x[i] = -P.AB_x[i];
829         P_tmp.AB_y[i] = -P.AB_y[i];
830         P_tmp.AB_z[i] = -P.AB_z[i];
831     }
832 
833     int ret = ostei_k_i_p_p(P_tmp, Q, screen_tol, work, INT__i_k_p_p);
834     double buffer[9072] SIMINT_ALIGN_ARRAY_DBL;
835 
836     for(int q = 0; q < ret; q++)
837     {
838         int idx = 0;
839         for(int a = 0; a < 28; ++a)
840         for(int b = 0; b < 36; ++b)
841         for(int c = 0; c < 3; ++c)
842         for(int d = 0; d < 3; ++d)
843             buffer[idx++] = INT__i_k_p_p[q*9072+b*252+a*9+c*3+d];
844 
845         memcpy(INT__i_k_p_p+q*9072, buffer, 9072*sizeof(double));
846     }
847 
848     return ret;
849 }
850 
851