Skip to content

Commit

Permalink
Improve ft_ao and sr-eri upper bound estimator
Browse files Browse the repository at this point in the history
  • Loading branch information
@sunqm
sunqm committed Dec 16, 2024
1 parent b159818 commit a519a35
Show file tree
Hide file tree
Showing 3 changed files with 90 additions and 23 deletions.
59 changes: 56 additions & 3 deletions pyscf/lib/vhf/nr_sr_vhf.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down Expand Up @@ -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;
Expand Down
23 changes: 14 additions & 9 deletions pyscf/pbc/df/ft_ao.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
31 changes: 20 additions & 11 deletions pyscf/pbc/df/rsdf_builder.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand Down

0 comments on commit a519a35

Please sign in to comment.