Improve ft_ao and sr-eri upper bound estimator

sunqm · Dec 17, 2024 · b4e71e4 · b4e71e4
1 parent c18dc14
commit b4e71e4
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 19 deletions.
diff --git a/pyscf/lib/vhf/nr_sr_vhf.c b/pyscf/lib/vhf/nr_sr_vhf.c
@@ -945,9 +945,11 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond,
                         r = sqrtf(r2);
 // u = .5/aij * (1 - root * akl/(aij+akl));
 // t = R_PA - root * akl/(aij+akl) * R_PQ;
-// The VRR iteration for constructing g[i] in libcint is 
+// The VRR iteration for constructing g[i] in libcint is
 //     g[i+1] = t * g[i] + i*u * g[i-1] = (t^2 + i*u) * g[i-1] + (i-1)*u*t * g[i-2]
-// To estimate the upper bound of g[n], t > 0 can be assumed. The first few terms are
+//
+// To estimate the upper bound of g[n] for a relative large separation of two
+// basis functions, t > 0 is assumed. The first few terms are
 // g[1] = t * g[0]
 // g[2] = t * g[1] + u * g[0] = (t^2+u) * g[0]
 // g[3] = t * g[2] + 2*u * g[1] = (t^2 + 3*u) * g[1]
@@ -958,6 +960,7 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond,
 // g[5] = t * g[4] + 4*u * g[3]
 //        < (t^2 + 3*u)^2 * g[1] + 4*u * (t^2 + 3*u) * g[1]
 //        < (t^2 + 5*u)^2 * g[1]
+//        < (t^2 + 5*u)^2.5 * g[0]
 // ...
 //
 // Higher order terms have the following inequality relations
@@ -970,25 +973,19 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond,
 //        < t * (t^2+(i-1)*u)^(i/2) * g[0] + i*u * (t^2+(i-1)*u)^(i/2-1) * g[1]
 //          = (t^2+(2*i-1)*u) * (t^2+(i-1)*u)^(i/2-1) * g[1]
 //          < (t^2 + (i+1)*u)^(i/2) * g[1]
+//          < (t^2 + (i+1)*u)^((i+1)/2) * g[0]
 // When i is odd
 // g[i+1] = t * g[i] + i*u * g[i-1]
 //        < t * (t^2+i*u)^((i-1)/2) * g[1] + i*u * (t^2+(i-2)*u)^((i-1)/2) * g[0]
 //        < t^2 * (t^2+i*u)^((i-1)/2) * g[0] + i*u * (t^2+i*u)^((i-1)/2) * g[0]
 //          = (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;
-                        u = .5 / aij;
-                        if (li % 2 == 1) {
-                                ti_fac = (li-1)/2 * logf(ti*ti + li*u + theta_r);
-                        } else {
-                                ti_fac = li/2 * logf(ti*ti + (li-1)*u + theta_r);
-                        }
-                        if (lj % 2 == 1) {
-                                tj_fac = (lj-1)/2 * logf(tj*tj + lj*u + theta_r);
-                        } else {
-                                tj_fac = lj/2 * logf(tj*tj + (lj-1)*u + 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);
 
                         v = ti_fac + tj_fac - a1*r2 + log_fac;
                         s_index[ish*Nbas+jsh] = v;

diff --git a/pyscf/pbc/df/ft_ao.py b/pyscf/pbc/df/ft_ao.py
@@ -695,12 +695,8 @@ def get_ovlp_mask(self, cutoff=None):
         li = cell_l[:,None]
         lj = supmol_l[None,:]
         odd_l = ls % 2 == 1
-        # li is even: ((li-1) * .5/aij + dri**2) ** (li//2)
-        # li is odd : (li * .5/aij + dri**2) ** ((li+1)//2) * dri
-        fac_dri = ((li-1) * .5/aij + dri**2) ** (li//2)
-        fac_drj = ((lj-1) * .5/aij + drj**2) ** (lj//2)
-        fac_dri[odd_l,:] = (li*.5/aij + dri**2)**((li-1)//2) * dri
-        fac_drj[:,odd_l] = (lj*.5/aij + drj**2)**((lj-1)//2) * drj
+        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.
         ovlp = (norm[:,None]*norm * np.pi**1.5 * aij**-1.5 *
                 cp.exp(-theta*dr**2) * fac_dri * fac_drj * fl)