1 #include <math.h>
2 #include <string.h>
3 #include "simint/recur_lookup.h"
4 #include "simint/osoei/osoei.h"
5
6 #define S_IJ(i,j) (s_ij[((i)*(nam2) + j)])
7
simint_compute_osoei_overlap(struct simint_shell const * sh1,struct simint_shell const * sh2,double * restrict integrals)8 int simint_compute_osoei_overlap(struct simint_shell const * sh1,
9 struct simint_shell const * sh2,
10 double * restrict integrals)
11 {
12 const int am1 = sh1->am;
13 const int am2 = sh2->am;
14
15 const int nam1 = am1 + 1;
16 const int nam2 = am2 + 1;
17 const int nam12 = nam1*nam2;
18
19 // Workspace for calculating terms
20 double work[3*nam12];
21
22 const double xyz1[3] = { sh1->x, sh1->y, sh1->z };
23 const double xyz2[3] = { sh2->x, sh2->y, sh2->z };
24
25 const int ncart1 = ((am1+1)*(am1+2))/2;
26 const int ncart2 = ((am2+1)*(am2+2))/2;
27
28 // Mostly needed just for the ordering
29 const int arrstart1 = am_recur_map[sh1->am];
30 const int arrstart2 = am_recur_map[sh2->am];
31 struct RecurInfo const * aminfo1 = &recurinfo_array[arrstart1];
32 struct RecurInfo const * aminfo2 = &recurinfo_array[arrstart2];
33
34 memset(integrals, 0, ncart1*ncart2*sizeof(double));
35
36 for(int i = 0; i < sh1->nprim; i++)
37 {
38 for(int j = 0; j < sh2->nprim; j++)
39 {
40 const double prefac = sh1->coef[i] * sh2->coef[j];
41 simint_osoei_overlap_terms(sh1->alpha[i], xyz1, sh2->alpha[j], xyz2,
42 sh1->am+1, sh2->am+1, work);
43
44 size_t outidx = 0;
45
46 for(int n = 0; n < ncart1; n++)
47 for(int m = 0; m < ncart2; m++)
48 {
49 const int8_t * ijk1 = aminfo1[n].ijk;
50 const int8_t * ijk2 = aminfo2[m].ijk;
51 const int xidx = ijk1[0]*nam2 + ijk2[0];
52 const int yidx = ijk1[1]*nam2 + ijk2[1];
53 const int zidx = ijk1[2]*nam2 + ijk2[2];
54
55 const double val = work[0*nam12 + xidx] *
56 work[1*nam12 + yidx] *
57 work[2*nam12 + zidx];
58
59 // remember: a and b are indices of primitives
60 integrals[outidx++] += prefac * val;
61 }
62 }
63 }
64
65 return 1;
66 }
67
68