Fix bugs in int1e_grids for Coulomb attenuated operators

sunqm · Jun 4, 2021 · 8d74777 · 8d74777
1 parent a621bcd
commit 8d74777
Showing 1 changed file with 6 additions and 5 deletions.
diff --git a/src/g1e_grids.c b/src/g1e_grids.c
@@ -100,14 +100,15 @@ FINT CINTg0_1e_grids(double *g, const double fac, CINTEnvVars *envs,
                         x = aij * RGSQUARE(rijrg, ig);
                         CINTrys_roots(nroots, x, ubuf, wbuf);
                         for (i = 0; i < nroots; i++) {
+                                // transform to t^2
                                 u[ig+GRID_BLKSIZE*i] = ubuf[i] / (ubuf[i] + 1);
                                 w[ig+GRID_BLKSIZE*i] = wbuf[i] * fac1;
                         }
                 }
         } else if (omega < 0) { // short-range part of range-separated Coulomb
                 theta = omega * omega / (omega * omega + aij);
                 sqrt_theta = sqrt(theta);
-                fac1 = fac * sqrt_theta / aij;
+                fac1 = fac / aij;
                 for (ig = 0; ig < bgrids; ig++) {
                         x = aij * RGSQUARE(rijrg, ig);
                         if (theta * x > envs->expcutoff) {
@@ -126,13 +127,13 @@ FINT CINTg0_1e_grids(double *g, const double fac, CINTEnvVars *envs,
                 }
         } else {  // long-range part of range-separated Coulomb
                 theta = omega * omega / (omega * omega + aij);
-                aij *= theta;
                 fac1 = fac * sqrt(theta) / aij;
                 for (ig = 0; ig < bgrids; ig++) {
-                        x = aij * RGSQUARE(rijrg, ig);
+                        x = aij * theta * RGSQUARE(rijrg, ig);
                         CINTrys_roots(nroots, x, ubuf, wbuf);
                         for (i = 0; i < nroots; i++) {
-                                u[ig+GRID_BLKSIZE*i] = ubuf[i] / (ubuf[i] + 1 - ubuf[i] * theta);
+                                // u stores t^2 = tau^2 * theta
+                                u[ig+GRID_BLKSIZE*i] = ubuf[i] / (ubuf[i] + 1) * theta;
                                 w[ig+GRID_BLKSIZE*i] = wbuf[i] * fac1;
                         }
                 }
@@ -180,7 +181,7 @@ FINT CINTg0_1e_grids(double *g, const double fac, CINTEnvVars *envs,
         double *rirg;
         MALLOC_ALIGN8_INSTACK(rirg, GRID_BLKSIZE*2);
         double *t2 = rirg + GRID_BLKSIZE;
-        double aij2 = 0.5 / aij;
+        double aij2 = 0.5 / envs->aij;
 
         for (n = 0; n < nroots; n++) {
                 for (ig = 0; ig < bgrids; ig++) {