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_s_s_d_p(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__s_s_d_p)8 int ostei_s_s_d_p(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__s_s_d_p)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__s_s_d_p);
17     int ab, cd, abcd;
18     int istart, jstart;
19     int iprimcd, nprim_icd, icd;
20     const int check_screen = (screen_tol > 0.0);
21     int i, j;
22     int n;
23     int not_screened;
24     int real_abcd;
25     int ibra;
26 
27     // partition workspace
28     double * const INT__s_s_d_s = work + (SIMINT_NSHELL_SIMD * 0);
29     double * const INT__s_s_f_s = work + (SIMINT_NSHELL_SIMD * 6);
30     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*16);
31     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
32     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_p_s = primwork + 4;
33     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_d_s = primwork + 13;
34     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_f_s = primwork + 25;
35     double * const hrrwork = (double *)(primwork + 35);
36 
37 
38     // Create constants
39     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
40     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
41     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
42 
43 
44     ////////////////////////////////////////
45     // Loop over shells and primitives
46     ////////////////////////////////////////
47 
48     real_abcd = 0;
49     istart = 0;
50     for(ab = 0; ab < P.nshell12_clip; ++ab)
51     {
52         const int iend = istart + P.nprim12[ab];
53 
54         cd = 0;
55         jstart = 0;
56 
57         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
58         {
59             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
60             int jend = jstart;
61             for(i = 0; i < nshellbatch; i++)
62                 jend += Q.nprim12[cd+i];
63 
64             // Clear the beginning of the workspace (where we are accumulating integrals)
65             memset(work, 0, SIMINT_NSHELL_SIMD * 16 * sizeof(double));
66             abcd = 0;
67 
68 
69             for(i = istart; i < iend; ++i)
70             {
71                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
72 
73                 if(check_screen)
74                 {
75                     // Skip this whole thing if always insignificant
76                     if((P.screen[i] * Q.screen_max) < screen_tol)
77                         continue;
78                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
79                 }
80 
81                 icd = 0;
82                 iprimcd = 0;
83                 nprim_icd = Q.nprim12[cd];
84                 double * restrict PRIM_PTR_INT__s_s_d_s = INT__s_s_d_s + abcd * 6;
85                 double * restrict PRIM_PTR_INT__s_s_f_s = INT__s_s_f_s + abcd * 10;
86 
87 
88 
89                 // Load these one per loop over i
90                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
91                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
92                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
93 
94 
95                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
96                 {
97                     // calculate the shell offsets
98                     // these are the offset from the shell pointed to by cd
99                     // for each element
100                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
101                     int lastoffset = 0;
102                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
103 
104                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
105                     {
106                         // Handle if the first element of the vector is a new shell
107                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
108                         {
109                             nprim_icd += Q.nprim12[cd + (++icd)];
110                             PRIM_PTR_INT__s_s_d_s += 6;
111                             PRIM_PTR_INT__s_s_f_s += 10;
112                         }
113                         iprimcd++;
114                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
115                         {
116                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
117                             {
118                                 shelloffsets[n] = shelloffsets[n-1] + 1;
119                                 lastoffset++;
120                                 nprim_icd += Q.nprim12[cd + (++icd)];
121                             }
122                             else
123                                 shelloffsets[n] = shelloffsets[n-1];
124                             iprimcd++;
125                         }
126                     }
127                     else
128                         iprimcd += SIMINT_SIMD_LEN;
129 
130                     // Do we have to compute this vector (or has it been screened out)?
131                     // (not_screened != 0 means we have to do this vector)
132                     if(check_screen)
133                     {
134                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
135                         if(vmax < screen_tol)
136                         {
137                             PRIM_PTR_INT__s_s_d_s += lastoffset*6;
138                             PRIM_PTR_INT__s_s_f_s += lastoffset*10;
139                             continue;
140                         }
141                     }
142 
143                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
144                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
145                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
146                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
147 
148 
149                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
150                     SIMINT_DBLTYPE PQ[3];
151                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
152                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
153                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
154                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
155                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
156                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
157 
158                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
159                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
160                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
161                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
162                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
163                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
164                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
165 
166                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
167                     SIMINT_DBLTYPE aoq_PQ[3];
168                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
169                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
170                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
171                     // Put a minus sign here so we don't have to in RR routines
172                     a_over_q = SIMINT_NEG(a_over_q);
173 
174 
175                     //////////////////////////////////////////////
176                     // Fjt function section
177                     // Maximum v value: 3
178                     //////////////////////////////////////////////
179                     // The parameter to the Fjt function
180                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
181 
182 
183                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
184 
185 
186                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 3);
187                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
188                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
189                     for(n = 0; n <= 3; n++)
190                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
191 
192                     //////////////////////////////////////////////
193                     // Primitive integrals: Vertical recurrance
194                     //////////////////////////////////////////////
195 
196                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
197                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
198 
199 
200 
201                     // Forming PRIM_INT__s_s_p_s[3 * 3];
202                     for(n = 0; n < 3; ++n)  // loop over orders of auxiliary function
203                     {
204 
205                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
206                         PRIM_INT__s_s_p_s[n * 3 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]);
207 
208                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
209                         PRIM_INT__s_s_p_s[n * 3 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 1]);
210 
211                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
212                         PRIM_INT__s_s_p_s[n * 3 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_p_s[n * 3 + 2]);
213 
214                     }
215 
216 
217 
218                     // Forming PRIM_INT__s_s_d_s[2 * 6];
219                     for(n = 0; n < 2; ++n)  // loop over orders of auxiliary function
220                     {
221 
222                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_p_s[n * 3 + 0]);
223                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 0]);
224                         PRIM_INT__s_s_d_s[n * 6 + 0] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 0]);
225 
226                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 0]);
227                         PRIM_INT__s_s_d_s[n * 6 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 1]);
228 
229                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 0]);
230                         PRIM_INT__s_s_d_s[n * 6 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_d_s[n * 6 + 2]);
231 
232                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_p_s[n * 3 + 1]);
233                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 3]);
234                         PRIM_INT__s_s_d_s[n * 6 + 3] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 3]);
235 
236                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 1]);
237                         PRIM_INT__s_s_d_s[n * 6 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_d_s[n * 6 + 4]);
238 
239                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_p_s[n * 3 + 2]);
240                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_d_s[n * 6 + 5]);
241                         PRIM_INT__s_s_d_s[n * 6 + 5] = SIMINT_FMADD( vrr_const_1_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_s_s[(n+1) * 1 + 0], PRIM_INT__s_s_s_s[n * 1 + 0]), PRIM_INT__s_s_d_s[n * 6 + 5]);
242 
243                     }
244 
245 
246 
247                     // Forming PRIM_INT__s_s_f_s[1 * 10];
248                     for(n = 0; n < 1; ++n)  // loop over orders of auxiliary function
249                     {
250 
251                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 0]);
252                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 0]);
253                         PRIM_INT__s_s_f_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 0], PRIM_INT__s_s_p_s[n * 3 + 0]), PRIM_INT__s_s_f_s[n * 10 + 0]);
254 
255                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 0]);
256                         PRIM_INT__s_s_f_s[n * 10 + 1] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 1]);
257 
258                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 0]);
259                         PRIM_INT__s_s_f_s[n * 10 + 2] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 0], PRIM_INT__s_s_f_s[n * 10 + 2]);
260 
261                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 3]);
262                         PRIM_INT__s_s_f_s[n * 10 + 3] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 3]);
263 
264                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 1]);
265                         PRIM_INT__s_s_f_s[n * 10 + 4] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 1], PRIM_INT__s_s_f_s[n * 10 + 4]);
266 
267                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_MUL(Q_PA[0], PRIM_INT__s_s_d_s[n * 6 + 5]);
268                         PRIM_INT__s_s_f_s[n * 10 + 5] = SIMINT_FMADD( aoq_PQ[0], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 5]);
269 
270                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 3]);
271                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 6]);
272                         PRIM_INT__s_s_f_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 1], PRIM_INT__s_s_p_s[n * 3 + 1]), PRIM_INT__s_s_f_s[n * 10 + 6]);
273 
274                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 3]);
275                         PRIM_INT__s_s_f_s[n * 10 + 7] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 3], PRIM_INT__s_s_f_s[n * 10 + 7]);
276 
277                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_MUL(Q_PA[1], PRIM_INT__s_s_d_s[n * 6 + 5]);
278                         PRIM_INT__s_s_f_s[n * 10 + 8] = SIMINT_FMADD( aoq_PQ[1], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 8]);
279 
280                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_MUL(Q_PA[2], PRIM_INT__s_s_d_s[n * 6 + 5]);
281                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( aoq_PQ[2], PRIM_INT__s_s_d_s[(n+1) * 6 + 5], PRIM_INT__s_s_f_s[n * 10 + 9]);
282                         PRIM_INT__s_s_f_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2q, SIMINT_FMADD(a_over_q, PRIM_INT__s_s_p_s[(n+1) * 3 + 2], PRIM_INT__s_s_p_s[n * 3 + 2]), PRIM_INT__s_s_f_s[n * 10 + 9]);
283 
284                     }
285 
286 
287 
288 
289                     ////////////////////////////////////
290                     // Accumulate contracted integrals
291                     ////////////////////////////////////
292                     if(lastoffset == 0)
293                     {
294                         contract_all(6, PRIM_INT__s_s_d_s, PRIM_PTR_INT__s_s_d_s);
295                         contract_all(10, PRIM_INT__s_s_f_s, PRIM_PTR_INT__s_s_f_s);
296                     }
297                     else
298                     {
299                         contract(6, shelloffsets, PRIM_INT__s_s_d_s, PRIM_PTR_INT__s_s_d_s);
300                         contract(10, shelloffsets, PRIM_INT__s_s_f_s, PRIM_PTR_INT__s_s_f_s);
301                         PRIM_PTR_INT__s_s_d_s += lastoffset*6;
302                         PRIM_PTR_INT__s_s_f_s += lastoffset*10;
303                     }
304 
305                 }  // close loop over j
306             }  // close loop over i
307 
308             //Advance to the next batch
309             jstart = SIMINT_SIMD_ROUND(jend);
310 
311             //////////////////////////////////////////////
312             // Contracted integrals: Horizontal recurrance
313             //////////////////////////////////////////////
314 
315 
316 
317 
318             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
319             {
320                 const double hCD[3] = { Q.AB_x[cd+abcd], Q.AB_y[cd+abcd], Q.AB_z[cd+abcd] };
321 
322                 // set up HRR pointers
323                 double const * restrict HRR_INT__s_s_d_s = INT__s_s_d_s + abcd * 6;
324                 double const * restrict HRR_INT__s_s_f_s = INT__s_s_f_s + abcd * 10;
325                 double * restrict HRR_INT__s_s_d_p = INT__s_s_d_p + real_abcd * 18;
326 
327                 // form INT__s_s_d_p
328                 for(ibra = 0; ibra < 1; ++ibra)
329                 {
330                     HRR_INT__s_s_d_p[ibra * 18 + 0] = HRR_INT__s_s_f_s[ibra * 10 + 0] + ( hCD[0] * HRR_INT__s_s_d_s[ibra * 6 + 0] );
331 
332                     HRR_INT__s_s_d_p[ibra * 18 + 1] = HRR_INT__s_s_f_s[ibra * 10 + 1] + ( hCD[1] * HRR_INT__s_s_d_s[ibra * 6 + 0] );
333 
334                     HRR_INT__s_s_d_p[ibra * 18 + 2] = HRR_INT__s_s_f_s[ibra * 10 + 2] + ( hCD[2] * HRR_INT__s_s_d_s[ibra * 6 + 0] );
335 
336                     HRR_INT__s_s_d_p[ibra * 18 + 3] = HRR_INT__s_s_f_s[ibra * 10 + 1] + ( hCD[0] * HRR_INT__s_s_d_s[ibra * 6 + 1] );
337 
338                     HRR_INT__s_s_d_p[ibra * 18 + 4] = HRR_INT__s_s_f_s[ibra * 10 + 3] + ( hCD[1] * HRR_INT__s_s_d_s[ibra * 6 + 1] );
339 
340                     HRR_INT__s_s_d_p[ibra * 18 + 5] = HRR_INT__s_s_f_s[ibra * 10 + 4] + ( hCD[2] * HRR_INT__s_s_d_s[ibra * 6 + 1] );
341 
342                     HRR_INT__s_s_d_p[ibra * 18 + 6] = HRR_INT__s_s_f_s[ibra * 10 + 2] + ( hCD[0] * HRR_INT__s_s_d_s[ibra * 6 + 2] );
343 
344                     HRR_INT__s_s_d_p[ibra * 18 + 7] = HRR_INT__s_s_f_s[ibra * 10 + 4] + ( hCD[1] * HRR_INT__s_s_d_s[ibra * 6 + 2] );
345 
346                     HRR_INT__s_s_d_p[ibra * 18 + 8] = HRR_INT__s_s_f_s[ibra * 10 + 5] + ( hCD[2] * HRR_INT__s_s_d_s[ibra * 6 + 2] );
347 
348                     HRR_INT__s_s_d_p[ibra * 18 + 9] = HRR_INT__s_s_f_s[ibra * 10 + 3] + ( hCD[0] * HRR_INT__s_s_d_s[ibra * 6 + 3] );
349 
350                     HRR_INT__s_s_d_p[ibra * 18 + 10] = HRR_INT__s_s_f_s[ibra * 10 + 6] + ( hCD[1] * HRR_INT__s_s_d_s[ibra * 6 + 3] );
351 
352                     HRR_INT__s_s_d_p[ibra * 18 + 11] = HRR_INT__s_s_f_s[ibra * 10 + 7] + ( hCD[2] * HRR_INT__s_s_d_s[ibra * 6 + 3] );
353 
354                     HRR_INT__s_s_d_p[ibra * 18 + 12] = HRR_INT__s_s_f_s[ibra * 10 + 4] + ( hCD[0] * HRR_INT__s_s_d_s[ibra * 6 + 4] );
355 
356                     HRR_INT__s_s_d_p[ibra * 18 + 13] = HRR_INT__s_s_f_s[ibra * 10 + 7] + ( hCD[1] * HRR_INT__s_s_d_s[ibra * 6 + 4] );
357 
358                     HRR_INT__s_s_d_p[ibra * 18 + 14] = HRR_INT__s_s_f_s[ibra * 10 + 8] + ( hCD[2] * HRR_INT__s_s_d_s[ibra * 6 + 4] );
359 
360                     HRR_INT__s_s_d_p[ibra * 18 + 15] = HRR_INT__s_s_f_s[ibra * 10 + 5] + ( hCD[0] * HRR_INT__s_s_d_s[ibra * 6 + 5] );
361 
362                     HRR_INT__s_s_d_p[ibra * 18 + 16] = HRR_INT__s_s_f_s[ibra * 10 + 8] + ( hCD[1] * HRR_INT__s_s_d_s[ibra * 6 + 5] );
363 
364                     HRR_INT__s_s_d_p[ibra * 18 + 17] = HRR_INT__s_s_f_s[ibra * 10 + 9] + ( hCD[2] * HRR_INT__s_s_d_s[ibra * 6 + 5] );
365 
366                 }
367 
368 
369             }  // close HRR loop
370 
371 
372         }   // close loop cdbatch
373 
374         istart = iend;
375     }  // close loop over ab
376 
377     return P.nshell12_clip * Q.nshell12_clip;
378 }
379 
380