diff --git a/pyscf/lib/vhf/nr_sr_vhf.c b/pyscf/lib/vhf/nr_sr_vhf.c
index f45a9dd5d8..12979e9832 100644
--- a/pyscf/lib/vhf/nr_sr_vhf.c
+++ b/pyscf/lib/vhf/nr_sr_vhf.c
@@ -891,7 +891,7 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond,
         float ai, aj, aij, ai_aij, a1, ci, cj;
         float xi, yi, zi, xj, yj, zj, xij, yij, zij;
         float dx, dy, dz, r2, r, v, log_fac, r_guess, theta, theta_r;
-        float u, ti, tj, ti_fac, tj_fac;
+        float fi, fj, u, ti, tj, ti_fac, tj_fac;
         double qtmp, tmp;
         float log_qmax;
         int shls[4];
diff --git a/pyscf/pbc/df/ft_ao.py b/pyscf/pbc/df/ft_ao.py
index c13214ed49..e6ae10b3bf 100644
--- a/pyscf/pbc/df/ft_ao.py
+++ b/pyscf/pbc/df/ft_ao.py
@@ -697,6 +697,8 @@ def get_ovlp_mask(self, cutoff=None):
         fac_dri = (li * .5/aij + dri**2) ** (li*.5)
         fac_drj = (lj * .5/aij + drj**2) ** (lj*.5)
         fl = 2*np.pi/vol * (dr/theta[:,None,None,:]) + 1.
+        # See also lib/vhf/nr_sr_vhf.c the comments of the upper bound
+        # estimation for CVHFnr_sr_int2e_q_cond.
         ovlp = ((cell_cs[:,None]*((2*li+1)/(4*np.pi))**.5) *
                 (supmol_cs*((2*lj+1)/(4*np.pi))**.5) * (np.pi/aij)**1.5 *
                 np.exp(-theta*dr**2) * fac_dri * fac_drj * fl)
@@ -774,18 +776,17 @@ def estimate_rcut(cell, precision=None):
     norm_ang = ((2*li+1)*(2*lj+1))**.5/(4*np.pi)
     c1 = ci * cj * norm_ang
     theta = ai * aj / aij
-    aij1 = aij**-.5
-    fac = np.pi**1.5*c1 * aij1**(lij+3) * (2*aij/np.pi)**.25 * aij**lij
+    fac = c1 * (np.pi/aij)**1.5 * (2*aij/np.pi)**.25
     fac /= precision
 
     r0 = cell.rcut
-    dri = aj*aij1 * r0 + 1.
-    drj = ai*aij1 * r0 + 1.
-    fl = 2*np.pi/cell.vol * r0/theta
-    r0 = (np.log(fac * dri**li * drj**lj * fl + 1.) / theta)**.5
-
-    dri = aj*aij1 * r0 + 1.
-    drj = ai*aij1 * r0 + 1.
-    fl = 2*np.pi/cell.vol * r0/theta
-    r0 = (np.log(fac * dri**li * drj**lj * fl + 1.) / theta)**.5
+    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/cell.vol * r0/theta + 1.
+    r0 = (np.log(fac * fac_dri * fac_drj * fl + 1.) / theta)**.5
+
+    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/cell.vol * r0/theta + 1.
+    r0 = (np.log(fac * fac_dri * fac_drj * fl + 1.) / theta)**.5
     return r0
diff --git a/pyscf/pbc/df/rsdf_builder.py b/pyscf/pbc/df/rsdf_builder.py
index bed680b007..e7e93014d2 100644
--- a/pyscf/pbc/df/rsdf_builder.py
+++ b/pyscf/pbc/df/rsdf_builder.py
@@ -1525,24 +1525,19 @@ def estimate_ft_rcut(rs_cell, precision=None, exclude_dd_block=False):
     norm_ang = ((2*li+1)*(2*lj+1))**.5/(4*np.pi)
     c1 = ci * cj * norm_ang
     theta = ai * aj / aij
-    aij1 = aij**-.5
-    fac = np.pi**1.5*c1 * aij1**(lij+3) * (2*aij/np.pi)**.25 * aij**lij
+    fac = c1 * (np.pi/aij)**1.5 * (2*aij/np.pi)**.25
     fac /= precision
 
     r0 = rs_cell.rcut
     # See also the estimator implemented in lib/vhf/nr_sr_vhf.c
-    dri = aj*aij1 * r0
-    drj = ai*aij1 * r0
-    fac_dri = (li * .5/aij + dri**2) ** (li/2)
-    fac_drj = (lj * .5/aij + drj**2) ** (lj/2)
-    fl = 2*np.pi*r0/theta + 1.
+    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 * fac_dri * fac_drj * fl + 1.) / theta)**.5
 
-    dri = aj*aij1 * r0
-    drj = ai*aij1 * r0
-    fac_dri = (li * .5/aij + dri**2) ** (li/2)
-    fac_drj = (lj * .5/aij + drj**2) ** (lj/2)
-    fl = 2*np.pi/rs_cell.vol*r0/theta
+    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
     rcut = r0
 
@@ -1563,23 +1558,18 @@ def estimate_ft_rcut(rs_cell, precision=None, exclude_dd_block=False):
             norm_ang = ((2*li+1)*(2*lj+1))**.5/(4*np.pi)
             c1 = ci * cj * norm_ang
             theta = ai * aj / aij
-            aij1 = aij**-.5
-            fac = np.pi**1.5*c1 * aij1**(lij+3) * (2*aij/np.pi)**.25 * aij**lij
+            fac = c1 * (np.pi/aij)**1.5 * (2*aij/np.pi)**.25
             fac /= precision
 
             r0 = rs_cell.rcut
-            dri = aj*aij1 * r0
-            drj = ai*aij1 * r0
-            fac_dri = (li * .5/aij + dri**2) ** (li/2)
-            fac_drj = (lj * .5/aij + drj**2) ** (lj/2)
-            fl = 2*np.pi/rs_cell.vol*r0/theta
+            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 * fac_dri * fac_drj * fl + 1.) / theta)**.5
 
-            dri = aj*aij1 * r0
-            drj = ai*aij1 * r0
-            fac_dri = (li * .5/aij + dri**2) ** (li/2)
-            fac_drj = (lj * .5/aij + drj**2) ** (lj/2)
-            fl = 2*np.pi*r0/theta + 1.
+            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 * fac_dri * fac_drj * fl + 1.) / theta)**.5
             rcut[smooth_mask] = r0
     return rcut