Updates for ft_ao rcut estimation

pyscf/lib/vhf/nr_sr_vhf.c

@@ -167,7 +167,7 @@ void CVHFdot_sr_nrs2ij(int (*intor)(), JKOperator **jkop, JKArray **vjk,
         int ish0 = ishls[0];
         int ish1 = ishls[1];
         int jsh0 = jshls[0];
-        int jsh1 = jshls[1];
+        //int jsh1 = jshls[1];
         int ksh0 = kshls[0];
         int ksh1 = kshls[1];
         int lsh0 = lshls[0];
@@ -304,7 +304,7 @@ void CVHFdot_sr_nrs2kl(int (*intor)(), JKOperator **jkop, JKArray **vjk,
         int ksh0 = kshls[0];
         int ksh1 = kshls[1];
         int lsh0 = lshls[0];
-        int lsh1 = lshls[1];
+        //int lsh1 = lshls[1];
         size_t Nbas = nbas;
         size_t Nbas2 = Nbas * Nbas;
         float *q_ijij = (float *)vhfopt->logq_cond;
@@ -704,15 +704,15 @@ void CVHFnr_sr_direct_drv(int (*intor)(), void (*fdot)(), JKOperator **jkop,
         int ish0 = shls_slice[0];
         int ish1 = shls_slice[1];
         int jsh0 = shls_slice[2];
-        int jsh1 = shls_slice[3];
+        //int jsh1 = shls_slice[3];
         int ksh0 = shls_slice[4];
-        int ksh1 = shls_slice[5];
+        //int ksh1 = shls_slice[5];
         int lsh0 = shls_slice[6];
-        int lsh1 = shls_slice[7];
+        //int lsh1 = shls_slice[7];
         int nish = ish1 - ish0;
-        int njsh = jsh1 - jsh0;
-        int nksh = ksh1 - ksh0;
-        int nlsh = lsh1 - lsh0;
+        //int njsh = jsh1 - jsh0;
+        //int nksh = ksh1 - ksh0;
+        //int nlsh = lsh1 - lsh0;
         assert(njsh == nish);
         assert(nksh == nish);
         assert(nlsh == nish);
@@ -981,11 +981,17 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond,
 //          = (t^2 + i*u)^((i+1)/2) * g[0]
 //          < (t^2 + (i+1)*u)^((i+1)/2) * g[0]
 // The upper bound of the (ij|ij) can be estimated using these inequalities
-                        ti = fj * r;
-                        tj = fi * r;
+//
+// When treating the orbital product in \int P_i(x) P_j(x) Gaussian(x) dx,
+// let F(x) = int P_i(x) Gaussian(x) dx
+// \int P_i(x) P_j(x) Gaussian(x) dx = \int F(x) d P_j(x)
+//     (mostly correct, needs proof) < (t^2+i*u)^(i/2) \int Gaussian(x) d P_j(x)
+//                                   < (t^2+i*u)^(i/2) (t^2+j*u)^(j/2) \int Gaussian(x) dx
+                        ti = fj * r + theta_r;
+                        tj = fi * r + theta_r;
                         u = .5f / aij;
-                        ti_fac = .5f*li * logf(ti*ti + li*u + theta_r);
-                        tj_fac = .5f*lj * logf(tj*tj + lj*u + theta_r);
+                        ti_fac = .5f*li * logf(ti*ti + li*u);
+                        tj_fac = .5f*lj * logf(tj*tj + lj*u);
 
                         v = ti_fac + tj_fac - a1*r2 + log_fac;
                         s_index[ish*Nbas+jsh] = v;

pyscf/lib/vhf/test/test_nr_direct.py

@@ -308,16 +308,16 @@ def test_sr_vhf_q_cond(self):
 
                     s_index = q_cond[2]
                     si_0 = s_index[0,0]
-                    si_others = s_index.diagonal()[1:]
+                    si_others = s_index.diagonal()
                     with mol.with_short_range_coulomb(omega):
                         ints = [abs(mol.intor_by_shell('int2e', (0,0,i,i))).max()
-                                for i in range(1, mol.nbas)]
+                                for i in range(mol.nbas)]
 
                         aij = akl = a * 2
                         omega2 = mol.omega**2
                         theta = 1/(2/aij+1/omega2)
-                        rr = rs**2
-                        estimator = rr * numpy.exp(si_0 + si_others - theta*rr)
+                        rr = numpy.append(0, rs)**2
+                        estimator = numpy.exp(si_0 + si_others - theta*rr)
                         assert all(estimator / ints > 1)
 
 

pyscf/pbc/df/ft_ao.py

@@ -772,7 +772,6 @@ def estimate_rcut(cell, precision=None):
     lj = ls
     cj = cs
     aij = ai + aj
-    lij = li + lj
     norm_ang = ((2*li+1)*(2*lj+1))**.5/(4*np.pi)
     c1 = ci * cj * norm_ang
     theta = ai * aj / aij

pyscf/pbc/df/rsdf_builder.py

@@ -1521,7 +1521,6 @@ def estimate_ft_rcut(rs_cell, precision=None, exclude_dd_block=False):
     lj = ls
     cj = cs
     aij = ai + aj
-    lij = li + lj
     norm_ang = ((2*li+1)*(2*lj+1))**.5/(4*np.pi)
     c1 = ci * cj * norm_ang
     theta = ai * aj / aij
@@ -1538,7 +1537,7 @@ def estimate_ft_rcut(rs_cell, precision=None, exclude_dd_block=False):
     fac_dri = (li * .5/aij + (aj/aij * r0)**2) ** (li/2)
     fac_drj = (lj * .5/aij + (ai/aij * r0)**2) ** (lj/2)
     fl = 2*np.pi/rs_cell.vol*r0/theta + 1.
-    r0 = (np.log(fac * dri**li * drj**lj * fl + 1.) / theta)**.5
+    r0 = (np.log(fac * fac_dri * fac_drj * fl + 1.) / theta)**.5
     rcut = r0
 
     if exclude_dd_block:
@@ -1554,7 +1553,6 @@ def estimate_ft_rcut(rs_cell, precision=None, exclude_dd_block=False):
             lj = ls[smooth_mask]
             cj = cs[smooth_mask]
             aij = ai + aj
-            lij = li + lj
             norm_ang = ((2*li+1)*(2*lj+1))**.5/(4*np.pi)
             c1 = ci * cj * norm_ang
             theta = ai * aj / aij

pyscf/pbc/gto/cell.py

@@ -402,8 +402,10 @@ def _estimate_rcut(alpha, l, c, precision=INTEGRAL_PRECISION):
     # For kinetic operator, basis becomes 2*a*r*|orig-basis>.
     # A penalty term 4*a^2*r^2 is included in the estimation.
     fac *= 4*alpha**2
-    r0 = (np.log(fac * r0 * (r0*.5+a1)**(2*l+2) + 1.) / theta)**.5
-    r0 = (np.log(fac * r0 * (r0*.5+a1)**(2*l+2) + 1.) / theta)**.5
+    # To encounter the impact of r^n in integrals, see the comments in
+    # lib/vhf/nr_sr_vhf.c for the upper bound estimation
+    r0 = (np.log(fac * r0 * ((r0*.5)**2+a1)**(l+1) + 1.) / theta)**.5
+    r0 = (np.log(fac * r0 * ((r0*.5)**2+a1)**(l+1) + 1.) / theta)**.5
     return r0
 
 def bas_rcut(cell, bas_id, precision=None):