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