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_s_f_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_s_f_s)8 int ostei_h_s_f_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_s_f_s)
13 {
14 
15     SIMINT_ASSUME_ALIGN_DBL(work);
16     SIMINT_ASSUME_ALIGN_DBL(INT__h_s_f_s);
17     memset(INT__h_s_f_s, 0, P.nshell12_clip * Q.nshell12_clip * 210 * sizeof(double));
18 
19     int ab, cd, abcd;
20     int istart, jstart;
21     int iprimcd, nprim_icd, icd;
22     const int check_screen = (screen_tol > 0.0);
23     int i, j;
24     int n;
25     int not_screened;
26 
27     // partition workspace
28     SIMINT_DBLTYPE * const primwork = (SIMINT_DBLTYPE *)(work + SIMINT_NSHELL_SIMD*0);
29     SIMINT_DBLTYPE * const restrict PRIM_INT__s_s_s_s = primwork + 0;
30     SIMINT_DBLTYPE * const restrict PRIM_INT__p_s_s_s = primwork + 9;
31     SIMINT_DBLTYPE * const restrict PRIM_INT__d_s_s_s = primwork + 33;
32     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_s_s = primwork + 75;
33     SIMINT_DBLTYPE * const restrict PRIM_INT__f_s_p_s = primwork + 135;
34     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_s_s = primwork + 225;
35     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_p_s = primwork + 300;
36     SIMINT_DBLTYPE * const restrict PRIM_INT__g_s_d_s = primwork + 435;
37     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_s_s = primwork + 615;
38     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_p_s = primwork + 699;
39     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_d_s = primwork + 888;
40     SIMINT_DBLTYPE * const restrict PRIM_INT__h_s_f_s = primwork + 1140;
41     double * const hrrwork = (double *)(primwork + 1350);
42 
43 
44     // Create constants
45     const SIMINT_DBLTYPE const_1 = SIMINT_DBLSET1(1);
46     const SIMINT_DBLTYPE const_2 = SIMINT_DBLSET1(2);
47     const SIMINT_DBLTYPE const_3 = SIMINT_DBLSET1(3);
48     const SIMINT_DBLTYPE const_4 = SIMINT_DBLSET1(4);
49     const SIMINT_DBLTYPE const_5 = SIMINT_DBLSET1(5);
50     const SIMINT_DBLTYPE one_half = SIMINT_DBLSET1(0.5);
51 
52 
53     ////////////////////////////////////////
54     // Loop over shells and primitives
55     ////////////////////////////////////////
56 
57     abcd = 0;
58     istart = 0;
59     for(ab = 0; ab < P.nshell12_clip; ++ab)
60     {
61         const int iend = istart + P.nprim12[ab];
62 
63         cd = 0;
64         jstart = 0;
65 
66         for(cd = 0; cd < Q.nshell12_clip; cd += SIMINT_NSHELL_SIMD)
67         {
68             const int nshellbatch = ((cd + SIMINT_NSHELL_SIMD) > Q.nshell12_clip) ? Q.nshell12_clip - cd : SIMINT_NSHELL_SIMD;
69             int jend = jstart;
70             for(i = 0; i < nshellbatch; i++)
71                 jend += Q.nprim12[cd+i];
72 
73 
74             for(i = istart; i < iend; ++i)
75             {
76                 SIMINT_DBLTYPE bra_screen_max;  // only used if check_screen
77 
78                 if(check_screen)
79                 {
80                     // Skip this whole thing if always insignificant
81                     if((P.screen[i] * Q.screen_max) < screen_tol)
82                         continue;
83                     bra_screen_max = SIMINT_DBLSET1(P.screen[i]);
84                 }
85 
86                 icd = 0;
87                 iprimcd = 0;
88                 nprim_icd = Q.nprim12[cd];
89                 double * restrict PRIM_PTR_INT__h_s_f_s = INT__h_s_f_s + abcd * 210;
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__h_s_f_s += 210;
116                         }
117                         iprimcd++;
118                         for(n = 1; n < SIMINT_SIMD_LEN; ++n)
119                         {
120                             if(iprimcd >= nprim_icd && ((icd+1) < nshellbatch))
121                             {
122                                 shelloffsets[n] = shelloffsets[n-1] + 1;
123                                 lastoffset++;
124                                 nprim_icd += Q.nprim12[cd + (++icd)];
125                             }
126                             else
127                                 shelloffsets[n] = shelloffsets[n-1];
128                             iprimcd++;
129                         }
130                     }
131                     else
132                         iprimcd += SIMINT_SIMD_LEN;
133 
134                     // Do we have to compute this vector (or has it been screened out)?
135                     // (not_screened != 0 means we have to do this vector)
136                     if(check_screen)
137                     {
138                         const double vmax = vector_max(SIMINT_MUL(bra_screen_max, SIMINT_DBLLOAD(Q.screen, j)));
139                         if(vmax < screen_tol)
140                         {
141                             PRIM_PTR_INT__h_s_f_s += lastoffset*210;
142                             continue;
143                         }
144                     }
145 
146                     const SIMINT_DBLTYPE Q_alpha = SIMINT_DBLLOAD(Q.alpha, j);
147                     const SIMINT_DBLTYPE PQalpha_mul = SIMINT_MUL(P_alpha, Q_alpha);
148                     const SIMINT_DBLTYPE PQalpha_sum = SIMINT_ADD(P_alpha, Q_alpha);
149                     const SIMINT_DBLTYPE one_over_PQalpha_sum = SIMINT_DIV(const_1, PQalpha_sum);
150 
151 
152                     /* construct R2 = (Px - Qx)**2 + (Py - Qy)**2 + (Pz -Qz)**2 */
153                     SIMINT_DBLTYPE PQ[3];
154                     PQ[0] = SIMINT_SUB(Pxyz[0], SIMINT_DBLLOAD(Q.x, j));
155                     PQ[1] = SIMINT_SUB(Pxyz[1], SIMINT_DBLLOAD(Q.y, j));
156                     PQ[2] = SIMINT_SUB(Pxyz[2], SIMINT_DBLLOAD(Q.z, j));
157                     SIMINT_DBLTYPE R2 = SIMINT_MUL(PQ[0], PQ[0]);
158                     R2 = SIMINT_FMADD(PQ[1], PQ[1], R2);
159                     R2 = SIMINT_FMADD(PQ[2], PQ[2], R2);
160 
161                     const SIMINT_DBLTYPE alpha = SIMINT_MUL(PQalpha_mul, one_over_PQalpha_sum); // alpha from MEST
162                     const SIMINT_DBLTYPE one_over_p = SIMINT_DIV(const_1, P_alpha);
163                     const SIMINT_DBLTYPE one_over_q = SIMINT_DIV(const_1, Q_alpha);
164                     const SIMINT_DBLTYPE one_over_2p = SIMINT_MUL(one_half, one_over_p);
165                     const SIMINT_DBLTYPE one_over_2q = SIMINT_MUL(one_half, one_over_q);
166                     const SIMINT_DBLTYPE one_over_2pq = SIMINT_MUL(one_half, one_over_PQalpha_sum);
167                     const SIMINT_DBLTYPE Q_PA[3] = { SIMINT_DBLLOAD(Q.PA_x, j), SIMINT_DBLLOAD(Q.PA_y, j), SIMINT_DBLLOAD(Q.PA_z, j) };
168 
169                     // NOTE: Minus sign!
170                     const SIMINT_DBLTYPE a_over_p = SIMINT_MUL(SIMINT_NEG(alpha), one_over_p);
171                     SIMINT_DBLTYPE aop_PQ[3];
172                     aop_PQ[0] = SIMINT_MUL(a_over_p, PQ[0]);
173                     aop_PQ[1] = SIMINT_MUL(a_over_p, PQ[1]);
174                     aop_PQ[2] = SIMINT_MUL(a_over_p, PQ[2]);
175 
176                     SIMINT_DBLTYPE a_over_q = SIMINT_MUL(alpha, one_over_q);
177                     SIMINT_DBLTYPE aoq_PQ[3];
178                     aoq_PQ[0] = SIMINT_MUL(a_over_q, PQ[0]);
179                     aoq_PQ[1] = SIMINT_MUL(a_over_q, PQ[1]);
180                     aoq_PQ[2] = SIMINT_MUL(a_over_q, PQ[2]);
181                     // Put a minus sign here so we don't have to in RR routines
182                     a_over_q = SIMINT_NEG(a_over_q);
183 
184 
185                     //////////////////////////////////////////////
186                     // Fjt function section
187                     // Maximum v value: 8
188                     //////////////////////////////////////////////
189                     // The parameter to the Fjt function
190                     const SIMINT_DBLTYPE F_x = SIMINT_MUL(R2, alpha);
191 
192 
193                     const SIMINT_DBLTYPE Q_prefac = mask_load(nlane, Q.prefac + j);
194 
195 
196                     boys_F_split(PRIM_INT__s_s_s_s, F_x, 8);
197                     SIMINT_DBLTYPE prefac = SIMINT_SQRT(one_over_PQalpha_sum);
198                     prefac = SIMINT_MUL(SIMINT_MUL(P_prefac, Q_prefac), prefac);
199                     for(n = 0; n <= 8; n++)
200                         PRIM_INT__s_s_s_s[n] = SIMINT_MUL(PRIM_INT__s_s_s_s[n], prefac);
201 
202                     //////////////////////////////////////////////
203                     // Primitive integrals: Vertical recurrance
204                     //////////////////////////////////////////////
205 
206                     const SIMINT_DBLTYPE vrr_const_1_over_2p = one_over_2p;
207                     const SIMINT_DBLTYPE vrr_const_2_over_2p = SIMINT_MUL(const_2, one_over_2p);
208                     const SIMINT_DBLTYPE vrr_const_3_over_2p = SIMINT_MUL(const_3, one_over_2p);
209                     const SIMINT_DBLTYPE vrr_const_4_over_2p = SIMINT_MUL(const_4, one_over_2p);
210                     const SIMINT_DBLTYPE vrr_const_1_over_2q = one_over_2q;
211                     const SIMINT_DBLTYPE vrr_const_2_over_2q = SIMINT_MUL(const_2, one_over_2q);
212                     const SIMINT_DBLTYPE vrr_const_1_over_2pq = one_over_2pq;
213                     const SIMINT_DBLTYPE vrr_const_2_over_2pq = SIMINT_MUL(const_2, one_over_2pq);
214                     const SIMINT_DBLTYPE vrr_const_3_over_2pq = SIMINT_MUL(const_3, one_over_2pq);
215                     const SIMINT_DBLTYPE vrr_const_4_over_2pq = SIMINT_MUL(const_4, one_over_2pq);
216                     const SIMINT_DBLTYPE vrr_const_5_over_2pq = SIMINT_MUL(const_5, one_over_2pq);
217 
218 
219 
220                     // Forming PRIM_INT__p_s_s_s[8 * 3];
221                     for(n = 0; n < 8; ++n)  // loop over orders of auxiliary function
222                     {
223 
224                         PRIM_INT__p_s_s_s[n * 3 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__s_s_s_s[n * 1 + 0]);
225                         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]);
226 
227                         PRIM_INT__p_s_s_s[n * 3 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__s_s_s_s[n * 1 + 0]);
228                         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]);
229 
230                         PRIM_INT__p_s_s_s[n * 3 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__s_s_s_s[n * 1 + 0]);
231                         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]);
232 
233                     }
234 
235 
236 
237                     // Forming PRIM_INT__d_s_s_s[7 * 6];
238                     for(n = 0; n < 7; ++n)  // loop over orders of auxiliary function
239                     {
240 
241                         PRIM_INT__d_s_s_s[n * 6 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__p_s_s_s[n * 3 + 0]);
242                         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]);
243                         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]);
244 
245                         PRIM_INT__d_s_s_s[n * 6 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 0]);
246                         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]);
247 
248                         PRIM_INT__d_s_s_s[n * 6 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 0]);
249                         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]);
250 
251                         PRIM_INT__d_s_s_s[n * 6 + 3] = SIMINT_MUL(P_PA[1], PRIM_INT__p_s_s_s[n * 3 + 1]);
252                         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]);
253                         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]);
254 
255                         PRIM_INT__d_s_s_s[n * 6 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 1]);
256                         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]);
257 
258                         PRIM_INT__d_s_s_s[n * 6 + 5] = SIMINT_MUL(P_PA[2], PRIM_INT__p_s_s_s[n * 3 + 2]);
259                         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]);
260                         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]);
261 
262                     }
263 
264 
265 
266                     // Forming PRIM_INT__f_s_s_s[6 * 10];
267                     for(n = 0; n < 6; ++n)  // loop over orders of auxiliary function
268                     {
269 
270                         PRIM_INT__f_s_s_s[n * 10 + 0] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 0]);
271                         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]);
272                         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]);
273 
274                         PRIM_INT__f_s_s_s[n * 10 + 1] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 0]);
275                         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]);
276 
277                         PRIM_INT__f_s_s_s[n * 10 + 2] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 0]);
278                         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]);
279 
280                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 3]);
281                         PRIM_INT__f_s_s_s[n * 10 + 3] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 3], PRIM_INT__f_s_s_s[n * 10 + 3]);
282 
283                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 1]);
284                         PRIM_INT__f_s_s_s[n * 10 + 4] = SIMINT_FMADD( aop_PQ[2], PRIM_INT__d_s_s_s[(n+1) * 6 + 1], PRIM_INT__f_s_s_s[n * 10 + 4]);
285 
286                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_MUL(P_PA[0], PRIM_INT__d_s_s_s[n * 6 + 5]);
287                         PRIM_INT__f_s_s_s[n * 10 + 5] = SIMINT_FMADD( aop_PQ[0], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 5]);
288 
289                         PRIM_INT__f_s_s_s[n * 10 + 6] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 3]);
290                         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]);
291                         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]);
292 
293                         PRIM_INT__f_s_s_s[n * 10 + 7] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 3]);
294                         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]);
295 
296                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_MUL(P_PA[1], PRIM_INT__d_s_s_s[n * 6 + 5]);
297                         PRIM_INT__f_s_s_s[n * 10 + 8] = SIMINT_FMADD( aop_PQ[1], PRIM_INT__d_s_s_s[(n+1) * 6 + 5], PRIM_INT__f_s_s_s[n * 10 + 8]);
298 
299                         PRIM_INT__f_s_s_s[n * 10 + 9] = SIMINT_MUL(P_PA[2], PRIM_INT__d_s_s_s[n * 6 + 5]);
300                         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]);
301                         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]);
302 
303                     }
304 
305 
306                     VRR_I_g_s_s_s(
307                             PRIM_INT__g_s_s_s,
308                             PRIM_INT__f_s_s_s,
309                             PRIM_INT__d_s_s_s,
310                             P_PA,
311                             a_over_p,
312                             aop_PQ,
313                             one_over_2p,
314                             5);
315 
316 
317                     VRR_I_h_s_s_s(
318                             PRIM_INT__h_s_s_s,
319                             PRIM_INT__g_s_s_s,
320                             PRIM_INT__f_s_s_s,
321                             P_PA,
322                             a_over_p,
323                             aop_PQ,
324                             one_over_2p,
325                             4);
326 
327 
328                     ostei_general_vrr_K(5, 0, 1, 0, 3,
329                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
330                             PRIM_INT__h_s_s_s, NULL, NULL, PRIM_INT__g_s_s_s, NULL, PRIM_INT__h_s_p_s);
331 
332 
333                     VRR_K_g_s_p_s(
334                             PRIM_INT__g_s_p_s,
335                             PRIM_INT__g_s_s_s,
336                             PRIM_INT__f_s_s_s,
337                             Q_PA,
338                             aoq_PQ,
339                             one_over_2pq,
340                             3);
341 
342 
343                     ostei_general_vrr_K(5, 0, 2, 0, 2,
344                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
345                             PRIM_INT__h_s_p_s, PRIM_INT__h_s_s_s, NULL, PRIM_INT__g_s_p_s, NULL, PRIM_INT__h_s_d_s);
346 
347 
348                     VRR_K_f_s_p_s(
349                             PRIM_INT__f_s_p_s,
350                             PRIM_INT__f_s_s_s,
351                             PRIM_INT__d_s_s_s,
352                             Q_PA,
353                             aoq_PQ,
354                             one_over_2pq,
355                             3);
356 
357 
358                     ostei_general_vrr_K(4, 0, 2, 0, 2,
359                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
360                             PRIM_INT__g_s_p_s, PRIM_INT__g_s_s_s, NULL, PRIM_INT__f_s_p_s, NULL, PRIM_INT__g_s_d_s);
361 
362 
363                     ostei_general_vrr_K(5, 0, 3, 0, 1,
364                             one_over_2q, a_over_q, one_over_2pq, aoq_PQ, Q_PA,
365                             PRIM_INT__h_s_d_s, PRIM_INT__h_s_p_s, NULL, PRIM_INT__g_s_d_s, NULL, PRIM_INT__h_s_f_s);
366 
367 
368 
369 
370                     ////////////////////////////////////
371                     // Accumulate contracted integrals
372                     ////////////////////////////////////
373                     if(lastoffset == 0)
374                     {
375                         contract_all(210, PRIM_INT__h_s_f_s, PRIM_PTR_INT__h_s_f_s);
376                     }
377                     else
378                     {
379                         contract(210, shelloffsets, PRIM_INT__h_s_f_s, PRIM_PTR_INT__h_s_f_s);
380                         PRIM_PTR_INT__h_s_f_s += lastoffset*210;
381                     }
382 
383                 }  // close loop over j
384             }  // close loop over i
385 
386             //Advance to the next batch
387             jstart = SIMINT_SIMD_ROUND(jend);
388             abcd += nshellbatch;
389 
390         }   // close loop cdbatch
391 
392         istart = iend;
393     }  // close loop over ab
394 
395     return P.nshell12_clip * Q.nshell12_clip;
396 }
397 
ostei_s_h_f_s(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__s_h_f_s)398 int ostei_s_h_f_s(struct simint_multi_shellpair const P,
399                   struct simint_multi_shellpair const Q,
400                   double screen_tol,
401                   double * const restrict work,
402                   double * const restrict INT__s_h_f_s)
403 {
404     double P_AB[3*P.nshell12];
405     struct simint_multi_shellpair P_tmp = P;
406     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
407     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
408     P_tmp.AB_x = P_AB;
409     P_tmp.AB_y = P_AB + P.nshell12;
410     P_tmp.AB_z = P_AB + 2*P.nshell12;
411 
412     for(int i = 0; i < P.nshell12; i++)
413     {
414         P_tmp.AB_x[i] = -P.AB_x[i];
415         P_tmp.AB_y[i] = -P.AB_y[i];
416         P_tmp.AB_z[i] = -P.AB_z[i];
417     }
418 
419     int ret = ostei_h_s_f_s(P_tmp, Q, screen_tol, work, INT__s_h_f_s);
420 
421     return ret;
422 }
423 
ostei_h_s_s_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__h_s_s_f)424 int ostei_h_s_s_f(struct simint_multi_shellpair const P,
425                   struct simint_multi_shellpair const Q,
426                   double screen_tol,
427                   double * const restrict work,
428                   double * const restrict INT__h_s_s_f)
429 {
430     double Q_AB[3*Q.nshell12];
431     struct simint_multi_shellpair Q_tmp = Q;
432     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
433     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
434     Q_tmp.AB_x = Q_AB;
435     Q_tmp.AB_y = Q_AB + Q.nshell12;
436     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
437 
438     for(int i = 0; i < Q.nshell12; i++)
439     {
440         Q_tmp.AB_x[i] = -Q.AB_x[i];
441         Q_tmp.AB_y[i] = -Q.AB_y[i];
442         Q_tmp.AB_z[i] = -Q.AB_z[i];
443     }
444 
445     int ret = ostei_h_s_f_s(P, Q_tmp, screen_tol, work, INT__h_s_s_f);
446 
447     return ret;
448 }
449 
ostei_s_h_s_f(struct simint_multi_shellpair const P,struct simint_multi_shellpair const Q,double screen_tol,double * const restrict work,double * const restrict INT__s_h_s_f)450 int ostei_s_h_s_f(struct simint_multi_shellpair const P,
451                   struct simint_multi_shellpair const Q,
452                   double screen_tol,
453                   double * const restrict work,
454                   double * const restrict INT__s_h_s_f)
455 {
456     double P_AB[3*P.nshell12];
457     struct simint_multi_shellpair P_tmp = P;
458     P_tmp.PA_x = P.PB_x;  P_tmp.PA_y = P.PB_y;  P_tmp.PA_z = P.PB_z;
459     P_tmp.PB_x = P.PA_x;  P_tmp.PB_y = P.PA_y;  P_tmp.PB_z = P.PA_z;
460     P_tmp.AB_x = P_AB;
461     P_tmp.AB_y = P_AB + P.nshell12;
462     P_tmp.AB_z = P_AB + 2*P.nshell12;
463 
464     for(int i = 0; i < P.nshell12; i++)
465     {
466         P_tmp.AB_x[i] = -P.AB_x[i];
467         P_tmp.AB_y[i] = -P.AB_y[i];
468         P_tmp.AB_z[i] = -P.AB_z[i];
469     }
470 
471     double Q_AB[3*Q.nshell12];
472     struct simint_multi_shellpair Q_tmp = Q;
473     Q_tmp.PA_x = Q.PB_x;  Q_tmp.PA_y = Q.PB_y;  Q_tmp.PA_z = Q.PB_z;
474     Q_tmp.PB_x = Q.PA_x;  Q_tmp.PB_y = Q.PA_y;  Q_tmp.PB_z = Q.PA_z;
475     Q_tmp.AB_x = Q_AB;
476     Q_tmp.AB_y = Q_AB + Q.nshell12;
477     Q_tmp.AB_z = Q_AB + 2*Q.nshell12;
478 
479     for(int i = 0; i < Q.nshell12; i++)
480     {
481         Q_tmp.AB_x[i] = -Q.AB_x[i];
482         Q_tmp.AB_y[i] = -Q.AB_y[i];
483         Q_tmp.AB_z[i] = -Q.AB_z[i];
484     }
485 
486     int ret = ostei_h_s_f_s(P_tmp, Q_tmp, screen_tol, work, INT__s_h_s_f);
487 
488     return ret;
489 }
490 
491