Skip to content

Commit

Permalink
gaussian charge model for int1e_grids
Browse files Browse the repository at this point in the history
sunqm committed Sep 17, 2023
1 parent 9dd786d commit cbebbab
Showing 2 changed files with 78 additions and 48 deletions.
73 changes: 47 additions & 26 deletions src/g1e_grids.c
Original file line number Diff line number Diff line change
@@ -88,10 +88,15 @@ FINT CINTg0_1e_grids(double *g, double cutoff,
}

#ifdef WITH_RANGE_COULOMB
const double omega = envs->env[PTR_RANGE_OMEGA];
double theta, sqrt_theta;
double omega = envs->env[PTR_RANGE_OMEGA];
#else
double omega = 0.;
#endif
double zeta = envs->env[PTR_RINV_ZETA];
double omega2, theta, sqrt_theta, a0, tau2;

if (omega == 0) {
assert(zeta >= 0);
if (omega == 0. && zeta == 0.) {
fac1 = envs->fac[0] / aij;
for (ig = 0; ig < bgrids; ig++) {
x = aij * RGSQUARE(rijrg, ig);
@@ -102,17 +107,27 @@ FINT CINTg0_1e_grids(double *g, double cutoff,
w[ig+GRID_BLKSIZE*i] = wbuf[i] * fac1;
}
}
} else if (omega < 0) { // short-range part of range-separated Coulomb
theta = omega * omega / (omega * omega + aij);
sqrt_theta = sqrt(theta);
} else if (omega < 0.) { // short-range part of range-separated Coulomb
a0 = aij;
fac1 = envs->fac[0] / aij;
// FIXME:
if (zeta == 0.) {
tau2 = 1.;
omega2 = omega * omega;
theta = omega2 / (omega2 + aij);
} else { // zeta > 0.
tau2 = zeta / (zeta + aij);
a0 *= tau2;
fac1 *= sqrt(tau2);
omega2 = omega * omega;
theta = omega2 / (omega2 + a0);
}
sqrt_theta = sqrt(theta);
// very small erfc() leads to ~0 weights. They can cause
// numerical issue in sr_rys_roots Use this cutoff as a
// temporary solution to avoid the numerical issue
// numerical issue in sr_rys_roots. Use this cutoff as a
// temporary solution to avoid numerical issues
double temp_cutoff = MIN(cutoff, EXPCUTOFF_SR);
for (ig = 0; ig < bgrids; ig++) {
x = aij * RGSQUARE(rijrg, ig);
x = a0 * RGSQUARE(rijrg, ig);
if (theta * x > temp_cutoff) {
// very small erfc() leads to ~0 weights
for (i = 0; i < nroots; i++) {
@@ -122,16 +137,33 @@ FINT CINTg0_1e_grids(double *g, double cutoff,
} else {
CINTsr_rys_roots(nroots, x, sqrt_theta, ubuf, wbuf);
for (i = 0; i < nroots; i++) {
u[ig+GRID_BLKSIZE*i] = ubuf[i] / (ubuf[i] + 1);
u[ig+GRID_BLKSIZE*i] = ubuf[i] / (ubuf[i] + 1) * tau2;
w[ig+GRID_BLKSIZE*i] = wbuf[i] * fac1;
}
}
}
} else { // long-range part of range-separated Coulomb
theta = omega * omega / (omega * omega + aij);
fac1 = envs->fac[0] * sqrt(theta) / aij;
} else {
// * long-range part of range-separated Coulomb
// * or Gaussian charge model, with rho(r) = Norm exp(-zeta*r^2)
a0 = aij;
fac1 = envs->fac[0] / aij;
if (zeta == 0.) { // omega > 0.
omega2 = omega * omega;
theta = omega2 / (omega2 + aij);
a0 *= theta;
fac1 *= sqrt(theta);
} else if (omega == 0.) { // zeta > 0.
theta = zeta / (zeta + aij);
a0 *= theta;
fac1 *= sqrt(theta);
} else { // omega > 0. && zeta > 0.
omega2 = omega * omega;
theta = omega2*zeta / (omega2*zeta + (zeta+omega2)*aij);
a0 *= theta;
fac1 *= sqrt(theta);
}
for (ig = 0; ig < bgrids; ig++) {
x = aij * theta * RGSQUARE(rijrg, ig);
x = a0 * RGSQUARE(rijrg, ig);
CINTrys_roots(nroots, x, ubuf, wbuf);
for (i = 0; i < nroots; i++) {
// u stores t^2 = tau^2 * theta
@@ -140,17 +172,6 @@ FINT CINTg0_1e_grids(double *g, double cutoff,
}
}
}
#else
fac1 = envs->fac[0] / aij;
for (ig = 0; ig < bgrids; ig++) {
x = aij * RGSQUARE(rijrg, ig);
CINTrys_roots(nroots, x, ubuf, wbuf);
for (i = 0; i < nroots; i++) {
u[ig+GRID_BLKSIZE*i] = ubuf[i] / (ubuf[i] + 1);
w[ig+GRID_BLKSIZE*i] = wbuf[i] * fac1;
}
}
#endif
FINT nmax = envs->li_ceil + envs->lj_ceil;
if (nmax == 0) {
return 1;
53 changes: 31 additions & 22 deletions testsuite/test_int1e_grids.py
Original file line number Diff line number Diff line change
@@ -206,7 +206,7 @@ def test_mol1():
import time
import pyscf
from pyscf import df
mol = pyscf.M(atom='''
mol0 = pyscf.M(atom='''
C 0.16146 -4.47867 0.00000
H -0.89492 5.29077 0.00000
H 0.47278 4.59602 0.88488
@@ -216,27 +216,36 @@ def test_mol1():
numpy.random.seed(12)
ngrids = 201
grids = numpy.random.random((ngrids, 3)) * 12 - 5
fmol = pyscf.gto.fakemol_for_charges(grids)
ref = df.incore.aux_e2(mol, fmol, intor='int3c2e').transpose(2,0,1)
j3c = mol.intor('int1e_grids', grids=grids)
print(abs(j3c - ref).max())

ref = df.incore.aux_e2(mol, fmol, intor='int3c2e_ip1').transpose(0,3,1,2)
j3c = mol.intor('int1e_grids_ip', grids=grids)
print(abs(j3c - ref).max())

ref = df.incore.aux_e2(mol, fmol, intor='int3c2e_ip1_cart').transpose(0,3,1,2)
j3c = mol.intor('int1e_grids_ip_cart', grids=grids)
print(abs(j3c - ref).max())

ref = df.incore.aux_e2(mol, fmol, intor='int3c2e_ip1_spinor').transpose(0,3,1,2)
j3c = mol.intor('int1e_grids_ip_spinor', grids=grids)
print(abs(j3c - ref).max())

ref = df.r_incore.aux_e2(mol, fmol, intor='int3c2e_spsp1_spinor').transpose(2,0,1)
j3c = mol.intor('int1e_grids_spvsp_spinor', grids=grids)
print(abs(j3c - ref).max())

for omega in (0, 0.1, -0.1):
for zeta in (0, 10, 1e16):
print('omega, zeta', omega, zeta)
if zeta == 0:
expnt = 1e16
else:
expnt = zeta
mol = mol0.copy()
mol.omega = omega
mol.set_rinv_zeta(zeta)
fmol = pyscf.gto.fakemol_for_charges(grids, expnt)
ref = df.incore.aux_e2(mol, fmol, intor='int3c2e').transpose(2,0,1)
j3c = mol.intor('int1e_grids', grids=grids)
print(abs(j3c - ref).max())

ref = df.incore.aux_e2(mol, fmol, intor='int3c2e_ip1').transpose(0,3,1,2)
j3c = mol.intor('int1e_grids_ip', grids=grids)
print(abs(j3c - ref).max())

ref = df.incore.aux_e2(mol, fmol, intor='int3c2e_ip1_cart').transpose(0,3,1,2)
j3c = mol.intor('int1e_grids_ip_cart', grids=grids)
print(abs(j3c - ref).max())

ref = df.incore.aux_e2(mol, fmol, intor='int3c2e_ip1_spinor').transpose(0,3,1,2)
j3c = mol.intor('int1e_grids_ip_spinor', grids=grids)
print(abs(j3c - ref).max())

ref = df.r_incore.aux_e2(mol, fmol, intor='int3c2e_spsp1_spinor').transpose(2,0,1)
j3c = mol.intor('int1e_grids_spvsp_spinor', grids=grids)
print(abs(j3c - ref).max())

test_int1e_grids_sph1('cint1e_grids_sph', 0, 1, 9)
test_int1e_grids_sph('cint1e_grids_ip_sph', 0, 1, 9)

0 comments on commit cbebbab

Please sign in to comment.