From a519a357c73757620ecaf3f6e931e16da5a5dce0 Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Sun, 15 Dec 2024 17:06:17 -0800 Subject: [PATCH] Improve ft_ao and sr-eri upper bound estimator --- pyscf/lib/vhf/nr_sr_vhf.c | 59 ++++++++++++++++++++++++++++++++++-- pyscf/pbc/df/ft_ao.py | 23 ++++++++------ pyscf/pbc/df/rsdf_builder.py | 31 ++++++++++++------- 3 files changed, 90 insertions(+), 23 deletions(-) diff --git a/pyscf/lib/vhf/nr_sr_vhf.c b/pyscf/lib/vhf/nr_sr_vhf.c index bd4a93fda1..221b8325af 100644 --- a/pyscf/lib/vhf/nr_sr_vhf.c +++ b/pyscf/lib/vhf/nr_sr_vhf.c @@ -890,7 +890,8 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond, int ij, i, j, di, dj, dij, di2, dj2; float ai, aj, aij, ai_aij, a1, ci, cj; float xi, yi, zi, xj, yj, zj, xij, yij, zij; - float dx, dy, dz, r2, v, log_fac, r_guess, theta, theta_r; + float dx, dy, dz, r2, r, v, log_fac, r_guess, theta, theta_r; + float u, ti, tj, ti_fac, tj_fac; double qtmp, tmp; float log_qmax; int shls[4]; @@ -924,20 +925,72 @@ void CVHFnr_sr_int2e_q_cond(int (*intor)(), CINTOpt *cintopt, float *q_cond, dy = yj - yi; dz = zj - zi; aij = ai + aj; - ai_aij = ai / aij; + fi = ai / aij; + fj = aj / aij; + ai_aij = fi; a1 = ai_aij * aj; xij = xi + ai_aij * dx; yij = yi + ai_aij * dy; zij = zi + ai_aij * dz; theta = omega2/(omega2+aij); + // r_guess is a guess for a radius, inside which the potential + // from the estimator is always larger than actual potential r_guess = R_GUESS_FAC / sqrtf(aij * theta); theta_r = theta * r_guess; // log(ci*cj * ((2*li+1)*(2*lj+1))**.5/(4*pi) * (pi/aij)**1.5) log_fac = logf(ci*cj * sqrtf((2*li+1.f)*(2*lj+1.f))/(4*M_PI)) + 1.5f*logf(M_PI/aij) + fac_guess; r2 = dx * dx + dy * dy + dz * dz; - v = (li+lj)*logf(MAX(theta_r, 1.f)) - a1*r2 + log_fac; + 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 +// 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 +// 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] +// g[4] = (t^2 + 3*u) * g[2] + 2*u*t * g[1] +// = (t^2 + 3*u) * (t^2+u) * g[0] + 2*u*t^2 * g[0] +// < (t^2 + 3*u) * (t^2+u) * g[0] + 2*u*(t^2+3*u) * g[0] +// = (t^2 + 3*u)^2 * g[0] +// 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] +// ... +// +// Higher order terms have the following inequality relations +// g[i+1] < (t^2 + (i+1)*u)^(i/2) * g[1] if i is even +// g[i+1] < (t^2 + i*u)^((i+1)/2) * g[0] if i is odd +// +// It can be proven as follows: +// When i is even +// g[i+1] = t * g[i] + i*u * g[i-1] +// < 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] +// 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] +// 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); + } + + v = ti_fac + tj_fac - a1*r2 + log_fac; s_index[ish*Nbas+jsh] = v; s_index[jsh*Nbas+ish] = v; xij_cond[ish*Nbas+jsh] = xij; diff --git a/pyscf/pbc/df/ft_ao.py b/pyscf/pbc/df/ft_ao.py index b85524abb7..5e7e0f7f2c 100644 --- a/pyscf/pbc/df/ft_ao.py +++ b/pyscf/pbc/df/ft_ao.py @@ -690,15 +690,20 @@ def get_ovlp_mask(self, cutoff=None): theta = cell_exps[:,None] * supmol_exps / aij dr = np.linalg.norm(cell_bas_coords[:,None,:] - supmol_bas_coords, axis=2) - aij1 = 1./aij - aij2 = aij**-.5 - dri = supmol_exps*aij1 * dr + aij2 - drj = cell_exps[:,None]*aij1 * dr + aij2 - norm_i = cell_cs * ((2*cell_l+1)/(4*np.pi))**.5 - norm_j = supmol_cs * ((2*supmol_l+1)/(4*np.pi))**.5 - fl = 2*np.pi/rs_cell.vol*dr/theta + 1. - ovlp = (np.pi**1.5 * norm_i[:,None]*norm_j * np.exp(-theta*dr**2) * - dri**cell_l[:,None] * drj**supmol_l * aij1**1.5 * fl) + dri = supmol_exps[None,:]/aij * dr + drj = cell_exps[:,None]/aij * dr + 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 + 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) return ovlp > cutoff @staticmethod diff --git a/pyscf/pbc/df/rsdf_builder.py b/pyscf/pbc/df/rsdf_builder.py index 89e9807a47..8edb519f1f 100644 --- a/pyscf/pbc/df/rsdf_builder.py +++ b/pyscf/pbc/df/rsdf_builder.py @@ -1530,13 +1530,18 @@ def estimate_ft_rcut(rs_cell, precision=None, exclude_dd_block=False): fac /= precision r0 = rs_cell.rcut - dri = aj*aij1 * r0 + 1. - drj = ai*aij1 * r0 + 1. + # 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. - r0 = (np.log(fac * dri**li * drj**lj * fl + 1.) / theta)**.5 + r0 = (np.log(fac * fac_dri * fac_drj * fl + 1.) / theta)**.5 - dri = aj*aij1 * r0 + 1. - drj = ai*aij1 * r0 + 1. + 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 r0 = (np.log(fac * dri**li * drj**lj * fl + 1.) / theta)**.5 rcut = r0 @@ -1563,15 +1568,19 @@ def estimate_ft_rcut(rs_cell, precision=None, exclude_dd_block=False): fac /= precision r0 = rs_cell.rcut - dri = aj*aij1 * r0 + 1. - drj = ai*aij1 * r0 + 1. + 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 - r0 = (np.log(fac * dri**li * drj**lj * fl + 1.) / theta)**.5 + r0 = (np.log(fac * fac_dri * fac_drj * fl + 1.) / theta)**.5 - dri = aj*aij1 * r0 + 1. - drj = ai*aij1 * r0 + 1. + 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. - 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[smooth_mask] = r0 return rcut