diff --git a/pyscf/lib/vhf/nr_sr_vhf.c b/pyscf/lib/vhf/nr_sr_vhf.c index 221b8325af..f45a9dd5d8 100644 --- 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 index 5e7e0f7f2c..95f63b969c 100644 --- 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)