gaussian charge model for int1e_grids

sunqm · Sep 17, 2023 · cbebbab · cbebbab
1 parent 9dd786d
commit cbebbab
Show file tree

Hide file tree

Showing 2 changed files with 78 additions and 48 deletions.
diff --git a/src/g1e_grids.c b/src/g1e_grids.c
@@ -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;

diff --git a/testsuite/test_int1e_grids.py b/testsuite/test_int1e_grids.py
@@ -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)