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_h_d_s_s(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__h_d_s_s)8 int ostei_h_d_s_s(struct simint_multi_shellpair const P,
9                   struct simint_multi_shellpair const Q,
10                   double screen_tol,
11                   double * const restrict work,
12                   double * const restrict INT__h_d_s_s)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__h_d_s_s);
17     int ab, cd, abcd;
18     int istart, jstart;
19     int iprimcd, nprim_icd, icd;
20     const int check_screen = (screen_tol > 0.0);
21     int i, j;
22     int n;
23     int not_screened;
24     int real_abcd;
25     int iket;
26 
27     // partition workspace
28     double * const INT__h_s_s_s = work + (SIMINT_NSHELL_SIMD * 0);
29     double * const INT__i_s_s_s = work + (SIMINT_NSHELL_SIMD * 21);
30     double * const INT__k_s_s_s = work + (SIMINT_NSHELL_SIMD * 49);
31     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*85);
32     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
33     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 8;
34     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 29;
35     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 65;
36     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 115;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 175;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__i_s_s_s = primwork + 238;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__k_s_s_s = primwork + 294;
40     double * const hrrwork = (double *)(primwork + 330);
41     double * const HRR_INT__h_p_s_s = hrrwork + 0;
42     double * const HRR_INT__i_p_s_s = hrrwork + 63;
43 
44 
45     // Create constants
46     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
47     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
48     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
49     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
50     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
51     const SIMINT_DBLTYPE const_6 = SIMINT_DBLSET1(6);
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 * 85 * 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__h_s_s_s = INT__h_s_s_s + abcd * 21;
96                 double * restrict PRIM_PTR_INT__i_s_s_s = INT__i_s_s_s + abcd * 28;
97                 double * restrict PRIM_PTR_INT__k_s_s_s = INT__k_s_s_s + abcd * 36;
98 
99 
100 
101                 // Load these one per loop over i
102                 const SIMINT_DBLTYPE P_alpha = SIMINT_DBLSET1(P.alpha[i]);
103                 const SIMINT_DBLTYPE P_prefac = SIMINT_DBLSET1(P.prefac[i]);
104                 const SIMINT_DBLTYPE Pxyz[3] = { SIMINT_DBLSET1(P.x[i]), SIMINT_DBLSET1(P.y[i]), SIMINT_DBLSET1(P.z[i]) };
105 
106                 const SIMINT_DBLTYPE P_PA[3] = { SIMINT_DBLSET1(P.PA_x[i]), SIMINT_DBLSET1(P.PA_y[i]), SIMINT_DBLSET1(P.PA_z[i]) };
107 
108                 for(j = jstart; j < jend; j += SIMINT_SIMD_LEN)
109                 {
110                     // calculate the shell offsets
111                     // these are the offset from the shell pointed to by cd
112                     // for each element
113                     int shelloffsets[SIMINT_SIMD_LEN] = {0};
114                     int lastoffset = 0;
115                     const int nlane = ( ((j + SIMINT_SIMD_LEN) < jend) ? SIMINT_SIMD_LEN : (jend - j));
116 
117                     if((iprimcd + SIMINT_SIMD_LEN) >= nprim_icd)
118                     {
119                         // Handle if the first element of the vector is a new shell
120                         if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
121                         {
122                             nprim_icd += Q.nprim12[cd + (++icd)];
123                             PRIM_PTR_INT__h_s_s_s += 21;
124                             PRIM_PTR_INT__i_s_s_s += 28;
125                             PRIM_PTR_INT__k_s_s_s += 36;
126                         }
127                         iprimcd++;
128                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
129                         {
130                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
131                             {
132                                 shelloffsets[n] = shelloffsets[n-1] + 1;
133                                 lastoffset++;
134                                 nprim_icd += Q.nprim12[cd + (++icd)];
135                             }
136                             else
137                                 shelloffsets[n] = shelloffsets[n-1];
138                             iprimcd++;
139                         }
140                     }
141                     else
142                         iprimcd += SIMINT_SIMD_LEN;
143 
144                     // Do we have to compute this vector (or has it been screened out)?
145                     // (not_screened != 0 means we have to do this vector)
146                     if(check_screen)
147                     {
148                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
149                         if(vmax < screen_tol)
150                         {
151                             PRIM_PTR_INT__h_s_s_s += lastoffset*21;
152                             PRIM_PTR_INT__i_s_s_s += lastoffset*28;
153                             PRIM_PTR_INT__k_s_s_s += lastoffset*36;
154                             continue;
155                         }
156                     }
157 
158                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
159                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
160                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
161                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
162 
163 
164                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
165                     SIMINT_DBLTYPE PQ[3];
166                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
167                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
168                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
169                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
170                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
171                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
172 
173                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
174                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
175                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
176                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
177                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
178                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
179 
180                     // NOTE: Minus sign!
181                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
182                     SIMINT_DBLTYPE aop_PQ[3];
183                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
184                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
185                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
186 
187 
188                     //////////////////////////////////////////////
189                     // Fjt function section
190                     // Maximum v value: 7
191                     //////////////////////////////////////////////
192                     // The parameter to the Fjt function
193                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
194 
195 
196                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
197 
198 
199                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 7);
200                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
201                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
202                     for(n = 0; n <= 7; n++)
203                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
204 
205                     //////////////////////////////////////////////
206                     // Primitive integrals: Vertical recurrance
207                     //////////////////////////////////////////////
208 
209                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
210                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
211                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
212                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
213                     const SIMINT_DBLTYPE vrr_const_5_over_2p = SIMINT_MUL(const_5, one_over_2p);
214                     const SIMINT_DBLTYPE vrr_const_6_over_2p = SIMINT_MUL(const_6, one_over_2p);
215 
216 
217 
218                     // Forming PRIM_INT__p_s_s_s[7 * 3];
219                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
220                     {
221 
222                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
223                         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]);
224 
225                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
226                         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]);
227 
228                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
229                         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]);
230 
231                     }
232 
233 
234 
235                     // Forming PRIM_INT__d_s_s_s[6 * 6];
236                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
237                     {
238 
239                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
240                         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]);
241                         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]);
242 
243                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
244                         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]);
245                         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]);
246 
247                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
248                         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]);
249                         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]);
250 
251                     }
252 
253 
254 
255                     // Forming PRIM_INT__f_s_s_s[5 * 10];
256                     for(n = 0; n < 5; ++n)  // loop over orders of auxiliary function
257                     {
258 
259                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
260                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 0]);
261                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 0], PRIM_INT__p_s_s_s[n * 3 + 0]), PRIM_INT__f_s_s_s[n * 10 + 0]);
262 
263                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
264                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 1]);
265 
266                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
267                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 0], PRIM_INT__f_s_s_s[n * 10 + 2]);
268 
269                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
270                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 6]);
271                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 1], PRIM_INT__p_s_s_s[n * 3 + 1]), PRIM_INT__f_s_s_s[n * 10 + 6]);
272 
273                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
274                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 7]);
275 
276                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
277                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 9]);
278                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_FMADD( vrr_const_2_over_2p, SIMINT_FMADD(a_over_p, PRIM_INT__p_s_s_s[(n+1) * 3 + 2], PRIM_INT__p_s_s_s[n * 3 + 2]), PRIM_INT__f_s_s_s[n * 10 + 9]);
279 
280                     }
281 
282 
283                     VRR_I_g_s_s_s(
284                             PRIM_INT__g_s_s_s,
285                             PRIM_INT__f_s_s_s,
286                             PRIM_INT__d_s_s_s,
287                             P_PA,
288                             a_over_p,
289                             aop_PQ,
290                             one_over_2p,
291                             4);
292 
293 
294                     VRR_I_h_s_s_s(
295                             PRIM_INT__h_s_s_s,
296                             PRIM_INT__g_s_s_s,
297                             PRIM_INT__f_s_s_s,
298                             P_PA,
299                             a_over_p,
300                             aop_PQ,
301                             one_over_2p,
302                             3);
303 
304 
305                     ostei_general_vrr1_I(6, 2,
306                             one_over_2p, a_over_p, aop_PQ, P_PA,
307                             PRIM_INT__h_s_s_s, PRIM_INT__g_s_s_s, PRIM_INT__i_s_s_s);
308 
309 
310                     ostei_general_vrr1_I(7, 1,
311                             one_over_2p, a_over_p, aop_PQ, P_PA,
312                             PRIM_INT__i_s_s_s, PRIM_INT__h_s_s_s, PRIM_INT__k_s_s_s);
313 
314 
315 
316 
317                     ////////////////////////////////////
318                     // Accumulate contracted integrals
319                     ////////////////////////////////////
320                     if(lastoffset == 0)
321                     {
322                         contract_all(21, PRIM_INT__h_s_s_s, PRIM_PTR_INT__h_s_s_s);
323                         contract_all(28, PRIM_INT__i_s_s_s, PRIM_PTR_INT__i_s_s_s);
324                         contract_all(36, PRIM_INT__k_s_s_s, PRIM_PTR_INT__k_s_s_s);
325                     }
326                     else
327                     {
328                         contract(21, shelloffsets, PRIM_INT__h_s_s_s, PRIM_PTR_INT__h_s_s_s);
329                         contract(28, shelloffsets, PRIM_INT__i_s_s_s, PRIM_PTR_INT__i_s_s_s);
330                         contract(36, shelloffsets, PRIM_INT__k_s_s_s, PRIM_PTR_INT__k_s_s_s);
331                         PRIM_PTR_INT__h_s_s_s += lastoffset*21;
332                         PRIM_PTR_INT__i_s_s_s += lastoffset*28;
333                         PRIM_PTR_INT__k_s_s_s += lastoffset*36;
334                     }
335 
336                 }  // close loop over j
337             }  // close loop over i
338 
339             //Advance to the next batch
340             jstart = SIMINT_SIMD_ROUND(jend);
341 
342             //////////////////////////////////////////////
343             // Contracted integrals: Horizontal recurrance
344             //////////////////////////////////////////////
345 
346 
347             const double hAB[3] = { P.AB_x[ab], P.AB_y[ab], P.AB_z[ab] };
348 
349 
350             for(abcd = 0; abcd < nshellbatch; ++abcd, ++real_abcd)
351             {
352 
353                 // set up HRR pointers
354                 double const * restrict HRR_INT__h_s_s_s = INT__h_s_s_s + abcd * 21;
355                 double const * restrict HRR_INT__i_s_s_s = INT__i_s_s_s + abcd * 28;
356                 double const * restrict HRR_INT__k_s_s_s = INT__k_s_s_s + abcd * 36;
357                 double * restrict HRR_INT__h_d_s_s = INT__h_d_s_s + real_abcd * 126;
358 
359                 // form INT__h_p_s_s
360                 ostei_general_hrr_J(5, 1, 0, 0, hAB, HRR_INT__i_s_s_s, HRR_INT__h_s_s_s, HRR_INT__h_p_s_s);
361 
362                 // form INT__i_p_s_s
363                 ostei_general_hrr_J(6, 1, 0, 0, hAB, HRR_INT__k_s_s_s, HRR_INT__i_s_s_s, HRR_INT__i_p_s_s);
364 
365                 // form INT__h_d_s_s
366                 ostei_general_hrr_J(5, 2, 0, 0, hAB, HRR_INT__i_p_s_s, HRR_INT__h_p_s_s, HRR_INT__h_d_s_s);
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 
ostei_d_h_s_s(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__d_h_s_s)380 int ostei_d_h_s_s(struct simint_multi_shellpair const P,
381                   struct simint_multi_shellpair const Q,
382                   double screen_tol,
383                   double * const restrict work,
384                   double * const restrict INT__d_h_s_s)
385 {
386     double P_AB[3*P.nshell12];
387     struct simint_multi_shellpair P_tmp = P;
388     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
389     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
390     P_tmp.AB_x = P_AB;
391     P_tmp.AB_y = P_AB + P.nshell12;
392     P_tmp.AB_z = P_AB + 2*P.nshell12;
393 
394     for(int i = 0; i < P.nshell12; i++)
395     {
396         P_tmp.AB_x[i] = -P.AB_x[i];
397         P_tmp.AB_y[i] = -P.AB_y[i];
398         P_tmp.AB_z[i] = -P.AB_z[i];
399     }
400 
401     int ret = ostei_h_d_s_s(P_tmp, Q, screen_tol, work, INT__d_h_s_s);
402     double buffer[126] SIMINT_ALIGN_ARRAY_DBL;
403 
404     for(int q = 0; q < ret; q++)
405     {
406         int idx = 0;
407         for(int a = 0; a < 6; ++a)
408         for(int b = 0; b < 21; ++b)
409         for(int c = 0; c < 1; ++c)
410         for(int d = 0; d < 1; ++d)
411             buffer[idx++] = INT__d_h_s_s[q*126+b*6+a*1+c*1+d];
412 
413         memcpy(INT__d_h_s_s+q*126, buffer, 126*sizeof(double));
414     }
415 
416     return ret;
417 }
418 
419