1 #include "simint/boys/boys.h"
2 #include "simint/ostei/gen/ostei_generated.h"
3 #include "simint/vectorization/vectorization.h"
4 #include <math.h>
5 #include <string.h>
6 
7 
ostei_p_p_s_d(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__p_p_s_d)8 int ostei_p_p_s_d(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__p_p_s_d)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__p_p_s_d);
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__p_s_s_d = work + (SIMINT_NSHELL_SIMD * 0);
29     double * const INT__d_s_s_d = work + (SIMINT_NSHELL_SIMD * 18);
30     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*54);
31     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
32     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_p = primwork + 5;
33     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 11;
34     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_p = primwork + 23;
35     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_d = primwork + 41;
36     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 59;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_p = primwork + 77;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_d = primwork + 113;
39     double * const hrrwork = (double *)(primwork + 149);
40 
41 
42     // Create constants
43     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
44     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
45     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
46 
47 
48     ////////////////////////////////////////
49     // Loop over shells and primitives
50     ////////////////////////////////////////
51 
52     real_abcd = 0;
53     istart = 0;
54     for(ab = 0; ab < P.nshell12_clip; ++ab)
55     {
56         const int iend = istart + P.nprim12[ab];
57 
58         cd = 0;
59         jstart = 0;
60 
61         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
62         {
63             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
64             int jend = jstart;
65             for(i = 0; i < nshellbatch; i++)
66                 jend += Q.nprim12[cd+i];
67 
68             // Clear the beginning of the workspace (where we are accumulating integrals)
69             memset(work, 0, SIMINT_NSHELL_SIMD * 54 * sizeof(double));
70             abcd = 0;
71 
72 
73             for(i = istart; i < iend; ++i)
74             {
75                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
76 
77                 if(check_screen)
78                 {
79                     // Skip this whole thing if always insignificant
80                     if((P.screen[i] * Q.screen_max) < screen_tol)
81                         continue;
82                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
83                 }
84 
85                 icd = 0;
86                 iprimcd = 0;
87                 nprim_icd = Q.nprim12[cd];
88                 double * restrict PRIM_PTR_INT__p_s_s_d = INT__p_s_s_d + abcd * 18;
89                 double * restrict PRIM_PTR_INT__d_s_s_d = INT__d_s_s_d + abcd * 36;
90 
91 
92 
93                 // Load these one per loop over i
94                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
95                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
96                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
97 
98                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
99 
100                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
101                 {
102                     // calculate the shell offsets
103                     // these are the offset from the shell pointed to by cd
104                     // for each element
105                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
106                     int lastoffset = 0;
107                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
108 
109                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
110                     {
111                         // Handle if the first element of the vector is a new shell
112                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
113                         {
114                             nprim_icd += Q.nprim12[cd + (++icd)];
115                             PRIM_PTR_INT__p_s_s_d += 18;
116                             PRIM_PTR_INT__d_s_s_d += 36;
117                         }
118                         iprimcd++;
119                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
120                         {
121                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
122                             {
123                                 shelloffsets[n] = shelloffsets[n-1] + 1;
124                                 lastoffset++;
125                                 nprim_icd += Q.nprim12[cd + (++icd)];
126                             }
127                             else
128                                 shelloffsets[n] = shelloffsets[n-1];
129                             iprimcd++;
130                         }
131                     }
132                     else
133                         iprimcd += SIMINT_SIMD_LEN;
134 
135                     // Do we have to compute this vector (or has it been screened out)?
136                     // (not_screened != 0 means we have to do this vector)
137                     if(check_screen)
138                     {
139                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
140                         if(vmax < screen_tol)
141                         {
142                             PRIM_PTR_INT__p_s_s_d += lastoffset*18;
143                             PRIM_PTR_INT__d_s_s_d += lastoffset*36;
144                             continue;
145                         }
146                     }
147 
148                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
149                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
150                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
151                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
152 
153 
154                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
155                     SIMINT_DBLTYPE PQ[3];
156                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
157                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
158                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
159                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
160                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
161                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
162 
163                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
164                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
165                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
166                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
167                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
168                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
169                     const SIMINT_DBLTYPE Q_PB[3] = { SIMINT_DBLLOAD(Q.PB_x, j), SIMINT_DBLLOAD(Q.PB_y, j), SIMINT_DBLLOAD(Q.PB_z, j) };
170 
171                     // NOTE: Minus sign!
172                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
173                     SIMINT_DBLTYPE aop_PQ[3];
174                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
175                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
176                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
177 
178                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
179                     SIMINT_DBLTYPE aoq_PQ[3];
180                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
181                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
182                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
183                     // Put a minus sign here so we don't have to in RR routines
184                     a_over_q = SIMINT_NEG(a_over_q);
185 
186 
187                     //////////////////////////////////////////////
188                     // Fjt function section
189                     // Maximum v value: 4
190                     //////////////////////////////////////////////
191                     // The parameter to the Fjt function
192                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
193 
194 
195                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
196 
197 
198                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 4);
199                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
200                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
201                     for(n = 0; n <= 4; n++)
202                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
203 
204                     //////////////////////////////////////////////
205                     // Primitive integrals: Vertical recurrance
206                     //////////////////////////////////////////////
207 
208                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
209                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
210                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
211                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
212 
213 
214 
215                     // Forming PRIM_INT__p_s_s_s[4 * 3];
216                     for(n = 0; n < 4; ++n)  // loop over orders of auxiliary function
217                     {
218 
219                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
220                         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]);
221 
222                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
223                         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]);
224 
225                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
226                         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]);
227 
228                     }
229 
230 
231 
232                     // Forming PRIM_INT__p_s_s_p[2 * 9];
233                     for(n = 0; n < 2; ++n)  // loop over orders of auxiliary function
234                     {
235 
236                         PRIM_INT__p_s_s_p[n * 9 + 0] = SIMINT_MUL(Q_PB[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
237                         PRIM_INT__p_s_s_p[n * 9 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_p[n * 9 + 0]);
238                         PRIM_INT__p_s_s_p[n * 9 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_p[n * 9 + 0]);
239 
240                         PRIM_INT__p_s_s_p[n * 9 + 1] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
241                         PRIM_INT__p_s_s_p[n * 9 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_p[n * 9 + 1]);
242 
243                         PRIM_INT__p_s_s_p[n * 9 + 2] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
244                         PRIM_INT__p_s_s_p[n * 9 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_p[n * 9 + 2]);
245 
246                         PRIM_INT__p_s_s_p[n * 9 + 3] = SIMINT_MUL(Q_PB[0], PRIM_INT__p_s_s_s[n * 3 + 1]);
247                         PRIM_INT__p_s_s_p[n * 9 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_p[n * 9 + 3]);
248 
249                         PRIM_INT__p_s_s_p[n * 9 + 4] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
250                         PRIM_INT__p_s_s_p[n * 9 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_p[n * 9 + 4]);
251                         PRIM_INT__p_s_s_p[n * 9 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_p[n * 9 + 4]);
252 
253                         PRIM_INT__p_s_s_p[n * 9 + 5] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
254                         PRIM_INT__p_s_s_p[n * 9 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_p[n * 9 + 5]);
255 
256                         PRIM_INT__p_s_s_p[n * 9 + 6] = SIMINT_MUL(Q_PB[0], PRIM_INT__p_s_s_s[n * 3 + 2]);
257                         PRIM_INT__p_s_s_p[n * 9 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_p[n * 9 + 6]);
258 
259                         PRIM_INT__p_s_s_p[n * 9 + 7] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_s[n * 3 + 2]);
260                         PRIM_INT__p_s_s_p[n * 9 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_p[n * 9 + 7]);
261 
262                         PRIM_INT__p_s_s_p[n * 9 + 8] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
263                         PRIM_INT__p_s_s_p[n * 9 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_p[n * 9 + 8]);
264                         PRIM_INT__p_s_s_p[n * 9 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__p_s_s_p[n * 9 + 8]);
265 
266                     }
267 
268 
269 
270                     // Forming PRIM_INT__s_s_s_p[2 * 3];
271                     for(n = 0; n < 2; ++n)  // loop over orders of auxiliary function
272                     {
273 
274                         PRIM_INT__s_s_s_p[n * 3 + 0] = SIMINT_MUL(Q_PB[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
275                         PRIM_INT__s_s_s_p[n * 3 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_p[n * 3 + 0]);
276 
277                         PRIM_INT__s_s_s_p[n * 3 + 1] = SIMINT_MUL(Q_PB[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
278                         PRIM_INT__s_s_s_p[n * 3 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_p[n * 3 + 1]);
279 
280                         PRIM_INT__s_s_s_p[n * 3 + 2] = SIMINT_MUL(Q_PB[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
281                         PRIM_INT__s_s_s_p[n * 3 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_p[n * 3 + 2]);
282 
283                     }
284 
285 
286 
287                     // Forming PRIM_INT__p_s_s_d[1 * 18];
288                     for(n = 0; n < 1; ++n)  // loop over orders of auxiliary function
289                     {
290 
291                         PRIM_INT__p_s_s_d[n * 18 + 0] = SIMINT_MUL(Q_PB[0], PRIM_INT__p_s_s_p[n * 9 + 0]);
292                         PRIM_INT__p_s_s_d[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_p[(n+1) * 9 + 0], PRIM_INT__p_s_s_d[n * 18 + 0]);
293                         PRIM_INT__p_s_s_d[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_s_d[n * 18 + 0]);
294                         PRIM_INT__p_s_s_d[n * 18 + 0] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_d[n * 18 + 0]);
295 
296                         PRIM_INT__p_s_s_d[n * 18 + 1] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_p[n * 9 + 0]);
297                         PRIM_INT__p_s_s_d[n * 18 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_p[(n+1) * 9 + 0], PRIM_INT__p_s_s_d[n * 18 + 1]);
298 
299                         PRIM_INT__p_s_s_d[n * 18 + 2] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 0]);
300                         PRIM_INT__p_s_s_d[n * 18 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 0], PRIM_INT__p_s_s_d[n * 18 + 2]);
301 
302                         PRIM_INT__p_s_s_d[n * 18 + 3] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_p[n * 9 + 1]);
303                         PRIM_INT__p_s_s_d[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_p[(n+1) * 9 + 1], PRIM_INT__p_s_s_d[n * 18 + 3]);
304                         PRIM_INT__p_s_s_d[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_s_d[n * 18 + 3]);
305 
306                         PRIM_INT__p_s_s_d[n * 18 + 4] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 1]);
307                         PRIM_INT__p_s_s_d[n * 18 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 1], PRIM_INT__p_s_s_d[n * 18 + 4]);
308 
309                         PRIM_INT__p_s_s_d[n * 18 + 5] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 2]);
310                         PRIM_INT__p_s_s_d[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 2], PRIM_INT__p_s_s_d[n * 18 + 5]);
311                         PRIM_INT__p_s_s_d[n * 18 + 5] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__p_s_s_d[n * 18 + 5]);
312 
313                         PRIM_INT__p_s_s_d[n * 18 + 6] = SIMINT_MUL(Q_PB[0], PRIM_INT__p_s_s_p[n * 9 + 3]);
314                         PRIM_INT__p_s_s_d[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_p[(n+1) * 9 + 3], PRIM_INT__p_s_s_d[n * 18 + 6]);
315                         PRIM_INT__p_s_s_d[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_s_d[n * 18 + 6]);
316 
317                         PRIM_INT__p_s_s_d[n * 18 + 7] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_p[n * 9 + 3]);
318                         PRIM_INT__p_s_s_d[n * 18 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_p[(n+1) * 9 + 3], PRIM_INT__p_s_s_d[n * 18 + 7]);
319                         PRIM_INT__p_s_s_d[n * 18 + 7] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_d[n * 18 + 7]);
320 
321                         PRIM_INT__p_s_s_d[n * 18 + 8] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 3]);
322                         PRIM_INT__p_s_s_d[n * 18 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 3], PRIM_INT__p_s_s_d[n * 18 + 8]);
323 
324                         PRIM_INT__p_s_s_d[n * 18 + 9] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_p[n * 9 + 4]);
325                         PRIM_INT__p_s_s_d[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_p[(n+1) * 9 + 4], PRIM_INT__p_s_s_d[n * 18 + 9]);
326                         PRIM_INT__p_s_s_d[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_s_d[n * 18 + 9]);
327                         PRIM_INT__p_s_s_d[n * 18 + 9] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__p_s_s_d[n * 18 + 9]);
328 
329                         PRIM_INT__p_s_s_d[n * 18 + 10] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 4]);
330                         PRIM_INT__p_s_s_d[n * 18 + 10] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 4], PRIM_INT__p_s_s_d[n * 18 + 10]);
331 
332                         PRIM_INT__p_s_s_d[n * 18 + 11] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 5]);
333                         PRIM_INT__p_s_s_d[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 5], PRIM_INT__p_s_s_d[n * 18 + 11]);
334                         PRIM_INT__p_s_s_d[n * 18 + 11] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__p_s_s_d[n * 18 + 11]);
335 
336                         PRIM_INT__p_s_s_d[n * 18 + 12] = SIMINT_MUL(Q_PB[0], PRIM_INT__p_s_s_p[n * 9 + 6]);
337                         PRIM_INT__p_s_s_d[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__p_s_s_p[(n+1) * 9 + 6], PRIM_INT__p_s_s_d[n * 18 + 12]);
338                         PRIM_INT__p_s_s_d[n * 18 + 12] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_s_d[n * 18 + 12]);
339 
340                         PRIM_INT__p_s_s_d[n * 18 + 13] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_p[n * 9 + 6]);
341                         PRIM_INT__p_s_s_d[n * 18 + 13] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_p[(n+1) * 9 + 6], PRIM_INT__p_s_s_d[n * 18 + 13]);
342 
343                         PRIM_INT__p_s_s_d[n * 18 + 14] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 6]);
344                         PRIM_INT__p_s_s_d[n * 18 + 14] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 6], PRIM_INT__p_s_s_d[n * 18 + 14]);
345                         PRIM_INT__p_s_s_d[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 0], PRIM_INT__p_s_s_d[n * 18 + 14]);
346 
347                         PRIM_INT__p_s_s_d[n * 18 + 15] = SIMINT_MUL(Q_PB[1], PRIM_INT__p_s_s_p[n * 9 + 7]);
348                         PRIM_INT__p_s_s_d[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__p_s_s_p[(n+1) * 9 + 7], PRIM_INT__p_s_s_d[n * 18 + 15]);
349                         PRIM_INT__p_s_s_d[n * 18 + 15] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_s_d[n * 18 + 15]);
350 
351                         PRIM_INT__p_s_s_d[n * 18 + 16] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 7]);
352                         PRIM_INT__p_s_s_d[n * 18 + 16] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 7], PRIM_INT__p_s_s_d[n * 18 + 16]);
353                         PRIM_INT__p_s_s_d[n * 18 + 16] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 1], PRIM_INT__p_s_s_d[n * 18 + 16]);
354 
355                         PRIM_INT__p_s_s_d[n * 18 + 17] = SIMINT_MUL(Q_PB[2], PRIM_INT__p_s_s_p[n * 9 + 8]);
356                         PRIM_INT__p_s_s_d[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__p_s_s_p[(n+1) * 9 + 8], PRIM_INT__p_s_s_d[n * 18 + 17]);
357                         PRIM_INT__p_s_s_d[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__p_s_s_d[n * 18 + 17]);
358                         PRIM_INT__p_s_s_d[n * 18 + 17] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__s_s_s_p[(n+1) * 3 + 2], PRIM_INT__p_s_s_d[n * 18 + 17]);
359 
360                     }
361 
362 
363 
364                     // Forming PRIM_INT__d_s_s_s[3 * 6];
365                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
366                     {
367 
368                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
369                         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]);
370                         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]);
371 
372                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
373                         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]);
374 
375                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
376                         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]);
377 
378                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
379                         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]);
380                         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]);
381 
382                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
383                         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]);
384 
385                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
386                         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]);
387                         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]);
388 
389                     }
390 
391 
392 
393                     // Forming PRIM_INT__d_s_s_p[2 * 18];
394                     for(n = 0; n < 2; ++n)  // loop over orders of auxiliary function
395                     {
396 
397                         PRIM_INT__d_s_s_p[n * 18 + 0] = SIMINT_MUL(Q_PB[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
398                         PRIM_INT__d_s_s_p[n * 18 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_s_p[n * 18 + 0]);
399                         PRIM_INT__d_s_s_p[n * 18 + 0] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_p[n * 18 + 0]);
400 
401                         PRIM_INT__d_s_s_p[n * 18 + 1] = SIMINT_MUL(Q_PB[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
402                         PRIM_INT__d_s_s_p[n * 18 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_s_p[n * 18 + 1]);
403 
404                         PRIM_INT__d_s_s_p[n * 18 + 2] = SIMINT_MUL(Q_PB[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
405                         PRIM_INT__d_s_s_p[n * 18 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__d_s_s_p[n * 18 + 2]);
406 
407                         PRIM_INT__d_s_s_p[n * 18 + 3] = SIMINT_MUL(Q_PB[0], PRIM_INT__d_s_s_s[n * 6 + 1]);
408                         PRIM_INT__d_s_s_p[n * 18 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_s_p[n * 18 + 3]);
409                         PRIM_INT__d_s_s_p[n * 18 + 3] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_p[n * 18 + 3]);
410 
411                         PRIM_INT__d_s_s_p[n * 18 + 4] = SIMINT_MUL(Q_PB[1], PRIM_INT__d_s_s_s[n * 6 + 1]);
412                         PRIM_INT__d_s_s_p[n * 18 + 4] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_s_p[n * 18 + 4]);
413                         PRIM_INT__d_s_s_p[n * 18 + 4] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_p[n * 18 + 4]);
414 
415                         PRIM_INT__d_s_s_p[n * 18 + 5] = SIMINT_MUL(Q_PB[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
416                         PRIM_INT__d_s_s_p[n * 18 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__d_s_s_p[n * 18 + 5]);
417 
418                         PRIM_INT__d_s_s_p[n * 18 + 6] = SIMINT_MUL(Q_PB[0], PRIM_INT__d_s_s_s[n * 6 + 2]);
419                         PRIM_INT__d_s_s_p[n * 18 + 6] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_s_p[n * 18 + 6]);
420                         PRIM_INT__d_s_s_p[n * 18 + 6] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_s_p[n * 18 + 6]);
421 
422                         PRIM_INT__d_s_s_p[n * 18 + 7] = SIMINT_MUL(Q_PB[1], PRIM_INT__d_s_s_s[n * 6 + 2]);
423                         PRIM_INT__d_s_s_p[n * 18 + 7] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_s_p[n * 18 + 7]);
424 
425                         PRIM_INT__d_s_s_p[n * 18 + 8] = SIMINT_MUL(Q_PB[2], PRIM_INT__d_s_s_s[n * 6 + 2]);
426                         PRIM_INT__d_s_s_p[n * 18 + 8] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 2], PRIM_INT__d_s_s_p[n * 18 + 8]);
427                         PRIM_INT__d_s_s_p[n * 18 + 8] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__d_s_s_p[n * 18 + 8]);
428 
429                         PRIM_INT__d_s_s_p[n * 18 + 9] = SIMINT_MUL(Q_PB[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
430                         PRIM_INT__d_s_s_p[n * 18 + 9] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_s_p[n * 18 + 9]);
431 
432                         PRIM_INT__d_s_s_p[n * 18 + 10] = SIMINT_MUL(Q_PB[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
433                         PRIM_INT__d_s_s_p[n * 18 + 10] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_s_p[n * 18 + 10]);
434                         PRIM_INT__d_s_s_p[n * 18 + 10] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_p[n * 18 + 10]);
435 
436                         PRIM_INT__d_s_s_p[n * 18 + 11] = SIMINT_MUL(Q_PB[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
437                         PRIM_INT__d_s_s_p[n * 18 + 11] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__d_s_s_p[n * 18 + 11]);
438 
439                         PRIM_INT__d_s_s_p[n * 18 + 12] = SIMINT_MUL(Q_PB[0], PRIM_INT__d_s_s_s[n * 6 + 4]);
440                         PRIM_INT__d_s_s_p[n * 18 + 12] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_s_p[n * 18 + 12]);
441 
442                         PRIM_INT__d_s_s_p[n * 18 + 13] = SIMINT_MUL(Q_PB[1], PRIM_INT__d_s_s_s[n * 6 + 4]);
443                         PRIM_INT__d_s_s_p[n * 18 + 13] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_s_p[n * 18 + 13]);
444                         PRIM_INT__d_s_s_p[n * 18 + 13] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_s_p[n * 18 + 13]);
445 
446                         PRIM_INT__d_s_s_p[n * 18 + 14] = SIMINT_MUL(Q_PB[2], PRIM_INT__d_s_s_s[n * 6 + 4]);
447                         PRIM_INT__d_s_s_p[n * 18 + 14] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 4], PRIM_INT__d_s_s_p[n * 18 + 14]);
448                         PRIM_INT__d_s_s_p[n * 18 + 14] = SIMINT_FMADD( vrr_const_1_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__d_s_s_p[n * 18 + 14]);
449 
450                         PRIM_INT__d_s_s_p[n * 18 + 15] = SIMINT_MUL(Q_PB[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
451                         PRIM_INT__d_s_s_p[n * 18 + 15] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_s_p[n * 18 + 15]);
452 
453                         PRIM_INT__d_s_s_p[n * 18 + 16] = SIMINT_MUL(Q_PB[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
454                         PRIM_INT__d_s_s_p[n * 18 + 16] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_s_p[n * 18 + 16]);
455 
456                         PRIM_INT__d_s_s_p[n * 18 + 17] = SIMINT_MUL(Q_PB[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
457                         PRIM_INT__d_s_s_p[n * 18 + 17] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__d_s_s_p[n * 18 + 17]);
458                         PRIM_INT__d_s_s_p[n * 18 + 17] = SIMINT_FMADD( vrr_const_2_over_2pq, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__d_s_s_p[n * 18 + 17]);
459 
460                     }
461 
462 
463                     VRR_L_d_s_s_d(
464                             PRIM_INT__d_s_s_d,
465                             PRIM_INT__d_s_s_p,
466                             PRIM_INT__d_s_s_s,
467                             PRIM_INT__p_s_s_p,
468                             Q_PB,
469                             a_over_q,
470                             aoq_PQ,
471                             one_over_2pq,
472                             one_over_2q,
473                             1);
474 
475 
476 
477 
478                     ////////////////////////////////////
479                     // Accumulate contracted integrals
480                     ////////////////////////////////////
481                     if(lastoffset == 0)
482                     {
483                         contract_all(18, PRIM_INT__p_s_s_d, PRIM_PTR_INT__p_s_s_d);
484                         contract_all(36, PRIM_INT__d_s_s_d, PRIM_PTR_INT__d_s_s_d);
485                     }
486                     else
487                     {
488                         contract(18, shelloffsets, PRIM_INT__p_s_s_d, PRIM_PTR_INT__p_s_s_d);
489                         contract(36, shelloffsets, PRIM_INT__d_s_s_d, PRIM_PTR_INT__d_s_s_d);
490                         PRIM_PTR_INT__p_s_s_d += lastoffset*18;
491                         PRIM_PTR_INT__d_s_s_d += lastoffset*36;
492                     }
493 
494                 }  // close loop over j
495             }  // close loop over i
496 
497             //Advance to the next batch
498             jstart = SIMINT_SIMD_ROUND(jend);
499 
500             //////////////////////////////////////////////
501             // Contracted integrals: Horizontal recurrance
502             //////////////////////////////////////////////
503 
504 
505             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
506 
507 
508             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
509             {
510 
511                 // set up HRR pointers
512                 double const * restrict HRR_INT__p_s_s_d = INT__p_s_s_d + abcd * 18;
513                 double const * restrict HRR_INT__d_s_s_d = INT__d_s_s_d + abcd * 36;
514                 double * restrict HRR_INT__p_p_s_d = INT__p_p_s_d + real_abcd * 54;
515 
516                 // form INT__p_p_s_d
517                 for(iket = 0; iket < 6; ++iket)
518                 {
519                     HRR_INT__p_p_s_d[0 * 6 + iket] = HRR_INT__d_s_s_d[0 * 6 + iket] + ( hAB[0] * HRR_INT__p_s_s_d[0 * 6 + iket] );
520 
521                     HRR_INT__p_p_s_d[1 * 6 + iket] = HRR_INT__d_s_s_d[1 * 6 + iket] + ( hAB[1] * HRR_INT__p_s_s_d[0 * 6 + iket] );
522 
523                     HRR_INT__p_p_s_d[2 * 6 + iket] = HRR_INT__d_s_s_d[2 * 6 + iket] + ( hAB[2] * HRR_INT__p_s_s_d[0 * 6 + iket] );
524 
525                     HRR_INT__p_p_s_d[3 * 6 + iket] = HRR_INT__d_s_s_d[1 * 6 + iket] + ( hAB[0] * HRR_INT__p_s_s_d[1 * 6 + iket] );
526 
527                     HRR_INT__p_p_s_d[4 * 6 + iket] = HRR_INT__d_s_s_d[3 * 6 + iket] + ( hAB[1] * HRR_INT__p_s_s_d[1 * 6 + iket] );
528 
529                     HRR_INT__p_p_s_d[5 * 6 + iket] = HRR_INT__d_s_s_d[4 * 6 + iket] + ( hAB[2] * HRR_INT__p_s_s_d[1 * 6 + iket] );
530 
531                     HRR_INT__p_p_s_d[6 * 6 + iket] = HRR_INT__d_s_s_d[2 * 6 + iket] + ( hAB[0] * HRR_INT__p_s_s_d[2 * 6 + iket] );
532 
533                     HRR_INT__p_p_s_d[7 * 6 + iket] = HRR_INT__d_s_s_d[4 * 6 + iket] + ( hAB[1] * HRR_INT__p_s_s_d[2 * 6 + iket] );
534 
535                     HRR_INT__p_p_s_d[8 * 6 + iket] = HRR_INT__d_s_s_d[5 * 6 + iket] + ( hAB[2] * HRR_INT__p_s_s_d[2 * 6 + iket] );
536 
537                 }
538 
539 
540 
541             }  // close HRR loop
542 
543 
544         }   // close loop cdbatch
545 
546         istart = iend;
547     }  // close loop over ab
548 
549     return P.nshell12_clip * Q.nshell12_clip;
550 }
551 
552