Skip to content

Commit

Permalink
Remove int3c2e kl operation, which causes an edge case bug (#284)
Browse files Browse the repository at this point in the history 
* Add a edge case test that computes int1e with fakemol, that supposes to fail

* remove xkl - xk

---------

Co-authored-by: Xiaojie Wu <[email protected]>
  • Loading branch information
@henryw7 @wxj6000
henryw7 and wxj6000 authored Dec 14, 2024
1 parent 9d168d5 commit bf03d85
Show file tree
Hide file tree
Showing 10 changed files with 80 additions and 92 deletions.
43 changes: 42 additions & 1 deletion gpu4pyscf/df/tests/test_int3c2e.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,11 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.

import pyscf
from pyscf import df
from pyscf import df, gto
from pyscf.gto.moleintor import getints, make_cintopt
from pyscf.df.grad.rhf import _int3c_wrapper
import numpy as np
import cupy as cp
import unittest
from gpu4pyscf.df import int3c2e
from gpu4pyscf.lib.cupy_helper import load_library
Expand Down Expand Up @@ -147,6 +148,46 @@ def test_int1e_iprinvip(self):
h1ao = mol.intor('int1e_iprinvip', comp=9) # <\nabla|1/r|>
assert np.linalg.norm(int3c[:,:,:,i] - h1ao) < 1e-7

def test_int1e_edge_case(self):
mol = gto.M(
atom =
"""
H 0 0 0
H -1 0.0 0
""",
basis =
"""
H S
0.1612777588 1.0000000
""")
mol.build()

xs = np.arange(-3.01, 3.00, 5.0)
ys = np.arange(-3.02, 3.00, 10.0)
zs = np.arange(-3.03, 3.00, 10.0)
grid_points = pyscf.lib.cartesian_prod([xs, ys, zs])

nao = mol.nao
ngrids = grid_points.shape[0]
dm = np.ones((nao, nao))
charges = np.ones(ngrids)

int3c2e_ip2 = mol._add_suffix('int3c2e_ip2')
cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, int3c2e_ip2)
fakemol = gto.fakemol_for_charges(grid_points)
q_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e_ip2, aosym='s1', cintopt=cintopt)
dq_cpu = np.einsum('xijk,ij,k->xk', q_nj, dm, charges)

dm = cp.asarray(dm)
charges = cp.asarray(charges)
intopt = int3c2e.VHFOpt(mol, fakemol, 'int2e')
intopt.build(1e-14, diag_block_with_triu=True, aosym=False)
coeff = intopt.coeff
dm_cart = coeff @ dm @ coeff.T
dq_gpu, _ = int3c2e.get_int3c2e_ip_jk(intopt, 0, 'ip2', charges, None, dm_cart)

cp.testing.assert_allclose(dq_cpu, dq_gpu, atol = 1e-10)

if __name__ == "__main__":
print("Full Tests for int3c")
unittest.main()
12 changes: 3 additions & 9 deletions gpu4pyscf/lib/gint/g2e.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -479,15 +479,9 @@ static void GINTg0_int3c2e(GINTEnvVars envs, double* __restrict__ g,
const double xi = bas_x[ish];
const double yi = bas_y[ish];
const double zi = bas_z[ish];
const double xk = bas_x[ksh];
const double yk = bas_y[ksh];
const double zk = bas_z[ksh];
const double xijxi = xij - xi;
const double yijyi = yij - yi;
const double zijzi = zij - zi;
const double xklxk = xkl - xk;
const double yklyk = ykl - yk;
const double zklzk = zkl - zk;

int nmax = envs.li_ceil + envs.lj_ceil;
int mmax = envs.lk_ceil + envs.ll_ceil;
Expand Down Expand Up @@ -563,9 +557,9 @@ static void GINTg0_int3c2e(GINTEnvVars envs, double* __restrict__ g,
//}
const double tmp3 = tmp1 * aij;
const double b01 = b00 + tmp4 * aij;
const double c0px = xklxk + tmp3 * xijxkl;
const double c0py = yklyk + tmp3 * yijykl;
const double c0pz = zklzk + tmp3 * zijzkl;
const double c0px = tmp3 * xijxkl;
const double c0py = tmp3 * yijykl;
const double c0pz = tmp3 * zijzkl;
double s0x = gx[i];
double s0y = gy[i];
double s0z = gz[i];
Expand Down
15 changes: 4 additions & 11 deletions gpu4pyscf/lib/gint/g3c2e.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,18 +151,11 @@ static void GINTfill_int3c2e_kernel0010(GINTEnvVars envs, ERITensor eri, BasisPr
double* __restrict__ x12 = c_bpcache.x12;
double* __restrict__ y12 = c_bpcache.y12;
double* __restrict__ z12 = c_bpcache.z12;

const int nbas = c_bpcache.nbas;
double* __restrict__ bas_x = c_bpcache.bas_coords;
double* __restrict__ bas_y = bas_x + nbas;
double* __restrict__ bas_z = bas_y + nbas;

double gout0 = 0;
double gout1 = 0;
double gout2 = 0;
const double xk = bas_x[ksh];
const double yk = bas_y[ksh];
const double zk = bas_z[ksh];

const int prim_ij0 = prim_ij;
const int prim_ij1 = prim_ij + nprim_ij;
const int prim_kl0 = prim_kl;
Expand Down Expand Up @@ -208,9 +201,9 @@ static void GINTfill_int3c2e_kernel0010(GINTEnvVars envs, ERITensor eri, BasisPr
//const double b00 = u2 * tmp4;
//const double tmp1 = aij * u2 / (u2 * aijkl + a1);;
const double tmp3 = aij * u2 / (u2 * aijkl + a1);;
const double c0px = xkl - xk + tmp3 * xijxkl;
const double c0py = ykl - yk + tmp3 * yijykl;
const double c0pz = zkl - zk + tmp3 * zijzkl;
const double c0px = tmp3 * xijxkl;
const double c0py = tmp3 * yijykl;
const double c0pz = tmp3 * zijzkl;
const double g_0 = 1;
const double g_1 = c0px;
const double g_2 = 1;
Expand Down
11 changes: 5 additions & 6 deletions gpu4pyscf/lib/gint/g3c2e_ip1ip2.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -165,12 +165,11 @@ static void GINTfill_int3c2e_ip1ip2_kernel000(GINTEnvVars envs, ERITensor eri, B
double gzx = 0;
double gzy = 0;
double gzz = 0;

const double xi = bas_x[ish];
const double yi = bas_y[ish];
const double zi = bas_z[ish];
const double xk = bas_x[ksh];
const double yk = bas_y[ksh];
const double zk = bas_z[ksh];

const int prim_ij0 = prim_ij;
const int prim_ij1 = prim_ij + nprim_ij;
const int prim_kl0 = prim_kl;
Expand Down Expand Up @@ -217,9 +216,9 @@ static void GINTfill_int3c2e_ip1ip2_kernel000(GINTEnvVars envs, ERITensor eri, B
const double c00y = yij - yi - tmp2 * yijykl;
const double c00z = zij - zi - tmp2 * zijzkl;
const double tmp3 = tmp1 * aij;
const double c0px = xkl - xk + tmp3 * xijxkl;
const double c0py = ykl - yk + tmp3 * yijykl;
const double c0pz = zkl - zk + tmp3 * zijzkl;
const double c0px = tmp3 * xijxkl;
const double c0py = tmp3 * yijykl;
const double c0pz = tmp3 * zijzkl;
const double g_0 = 1;
const double g_1 = c00x;
const double g_2 = c0px;
Expand Down
19 changes: 5 additions & 14 deletions gpu4pyscf/lib/gint/g3c2e_ip2.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,17 +138,10 @@ static void GINTfill_int3c2e_ip2_kernel000(GINTEnvVars envs, ERITensor eri, Basi
double* __restrict__ z12 = c_bpcache.z12;
double* __restrict__ a1 = c_bpcache.a1;

const int nbas = c_bpcache.nbas;
double* __restrict__ bas_x = c_bpcache.bas_coords;
double* __restrict__ bas_y = bas_x + nbas;
double* __restrict__ bas_z = bas_y + nbas;

double gout0 = 0;
double gout1 = 0;
double gout2 = 0;
const double xk = bas_x[ksh];
const double yk = bas_y[ksh];
const double zk = bas_z[ksh];

const int prim_ij0 = prim_ij;
const int prim_ij1 = prim_ij + nprim_ij;
const int prim_kl0 = prim_kl;
Expand Down Expand Up @@ -191,13 +184,11 @@ static void GINTfill_int3c2e_ip2_kernel000(GINTEnvVars envs, ERITensor eri, Basi
}
root0 /= root0 + 1 - root0 * theta;
const double u2 = a0 * root0;
const double tmp4 = .5 / (u2 * aijkl + a1);
const double b00 = u2 * tmp4;
const double tmp1 = 2 * b00;
const double tmp1 = u2 / (u2 * aijkl + a1);
const double tmp3 = tmp1 * aij;
const double c0px = xkl - xk + tmp3 * xijxkl;
const double c0py = ykl - yk + tmp3 * yijykl;
const double c0pz = zkl - zk + tmp3 * zijzkl;
const double c0px = tmp3 * xijxkl;
const double c0py = tmp3 * yijykl;
const double c0pz = tmp3 * zijzkl;
const double g_0 = 1;
const double g_1 = c0px;
const double g_2 = 1;
Expand Down
15 changes: 4 additions & 11 deletions gpu4pyscf/lib/gint/g3c2e_ipip2.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -147,11 +147,6 @@ static void GINTfill_int3c2e_ipip2_kernel000(GINTEnvVars envs, ERITensor eri, Ba
double* __restrict__ x12 = c_bpcache.x12;
double* __restrict__ y12 = c_bpcache.y12;
double* __restrict__ z12 = c_bpcache.z12;

const int nbas = c_bpcache.nbas;
double* __restrict__ bas_x = c_bpcache.bas_coords;
double* __restrict__ bas_y = bas_x + nbas;
double* __restrict__ bas_z = bas_y + nbas;

double gxx = 0;
double gxy = 0;
Expand All @@ -162,9 +157,7 @@ static void GINTfill_int3c2e_ipip2_kernel000(GINTEnvVars envs, ERITensor eri, Ba
double gzx = 0;
double gzy = 0;
double gzz = 0;
double xk = bas_x[ksh];
double yk = bas_y[ksh];
double zk = bas_z[ksh];

const int prim_ij0 = prim_ij;
const int prim_ij1 = prim_ij + nprim_ij;
const int prim_kl0 = prim_kl;
Expand Down Expand Up @@ -208,9 +201,9 @@ static void GINTfill_int3c2e_ipip2_kernel000(GINTEnvVars envs, ERITensor eri, Ba
const double tmp1 = 2 * b00;
const double tmp3 = tmp1 * aij;
const double b01 = b00 + tmp4 * aij;
const double c0px = xkl - xk + tmp3 * xijxkl;
const double c0py = ykl - yk + tmp3 * yijykl;
const double c0pz = zkl - zk + tmp3 * zijzkl;
const double c0px = tmp3 * xijxkl;
const double c0py = tmp3 * yijykl;
const double c0pz = tmp3 * zijzkl;
const double g_0 = 1;
const double g_1 = c0px;
const double g_2 = c0px * c0px + b01;
Expand Down
18 changes: 5 additions & 13 deletions gpu4pyscf/lib/gvhf/g3c2e_ip2.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,17 +163,11 @@ static void GINTint3c2e_ip2_jk_kernel001(GINTEnvVars envs, JKMatrix jk, BasisPro
double* __restrict__ a1 = c_bpcache.a1;
int ij, kl;
int prim_ij0, prim_ij1, prim_kl0, prim_kl1;
int nbas = c_bpcache.nbas;
double* __restrict__ bas_x = c_bpcache.bas_coords;
double* __restrict__ bas_y = bas_x + nbas;
double* __restrict__ bas_z = bas_y + nbas;

double gout0 = 0;
double gout1 = 0;
double gout2 = 0;
double xk = bas_x[ksh];
double yk = bas_y[ksh];
double zk = bas_z[ksh];

prim_ij0 = prim_ij;
prim_ij1 = prim_ij + nprim_ij;
prim_kl0 = prim_kl;
Expand Down Expand Up @@ -216,13 +210,11 @@ static void GINTint3c2e_ip2_jk_kernel001(GINTEnvVars envs, JKMatrix jk, BasisPro
}
root0 /= root0 + 1 - root0 * theta;
double u2 = a0 * root0;
double tmp4 = .5 / (u2 * aijkl + a1);
double b00 = u2 * tmp4;
double tmp1 = 2 * b00;
double tmp1 = u2 / (u2 * aijkl + a1);
double tmp3 = tmp1 * aij;
double c0px = xkl - xk + tmp3 * xijxkl;
double c0py = ykl - yk + tmp3 * yijykl;
double c0pz = zkl - zk + tmp3 * zijzkl;
double c0px = tmp3 * xijxkl;
double c0py = tmp3 * yijykl;
double c0pz = tmp3 * zijzkl;
double g_0 = 1;
double g_1 = c0px;
double g_2 = 1;
Expand Down
18 changes: 5 additions & 13 deletions gpu4pyscf/lib/gvhf/g3c2e_pass1.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -155,17 +155,11 @@ static void GINTint3c2e_pass1_j_kernel0010(GINTEnvVars envs, JKMatrix jk, BasisP
}
int ij, kl;
int prim_ij0, prim_ij1, prim_kl0, prim_kl1;
int nbas = c_bpcache.nbas;
double* __restrict__ bas_x = c_bpcache.bas_coords;
double* __restrict__ bas_y = bas_x + nbas;
double* __restrict__ bas_z = bas_y + nbas;

double gout0 = 0;
double gout1 = 0;
double gout2 = 0;
double xk = bas_x[ksh];
double yk = bas_y[ksh];
double zk = bas_z[ksh];

prim_ij0 = prim_ij;
prim_ij1 = prim_ij + nprim_ij;
prim_kl0 = prim_kl;
Expand Down Expand Up @@ -204,13 +198,11 @@ static void GINTint3c2e_pass1_j_kernel0010(GINTEnvVars envs, JKMatrix jk, BasisP
root0 = fmt1 / (fmt0 - fmt1);
}
double u2 = a0 * root0;
double tmp4 = .5 / (u2 * aijkl + a1);
double b00 = u2 * tmp4;
double tmp1 = 2 * b00;
double tmp1 = u2 / (u2 * aijkl + a1);
double tmp3 = tmp1 * aij;
double c0px = xkl - xk + tmp3 * xijxkl;
double c0py = ykl - yk + tmp3 * yijykl;
double c0pz = zkl - zk + tmp3 * zijzkl;
double c0px = tmp3 * xijxkl;
double c0py = tmp3 * yijykl;
double c0pz = tmp3 * zijzkl;
double g_0 = 1;
double g_1 = c0px;
double g_2 = 1;
Expand Down
19 changes: 6 additions & 13 deletions gpu4pyscf/lib/gvhf/g3c2e_pass2.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,19 +152,14 @@ static void GINTint3c2e_pass2_j_kernel0010(GINTEnvVars envs, JKMatrix jk, BasisP
double* __restrict__ z12 = c_bpcache.z12;
int ij, kl;
int prim_ij0, prim_ij1, prim_kl0, prim_kl1;
int nbas = c_bpcache.nbas;
double* __restrict__ bas_x = c_bpcache.bas_coords;
double* __restrict__ bas_y = bas_x + nbas;
double* __restrict__ bas_z = bas_y + nbas;

if (ish == jsh){
norm *= .5;
}
double gout0 = 0;
double gout1 = 0;
double gout2 = 0;
double xk = bas_x[ksh];
double yk = bas_y[ksh];
double zk = bas_z[ksh];

prim_ij0 = prim_ij;
prim_ij1 = prim_ij + nprim_ij;
prim_kl0 = prim_kl;
Expand Down Expand Up @@ -203,13 +198,11 @@ static void GINTint3c2e_pass2_j_kernel0010(GINTEnvVars envs, JKMatrix jk, BasisP
root0 = fmt1 / (fmt0 - fmt1);
}
double u2 = a0 * root0;
double tmp4 = .5 / (u2 * aijkl + a1);
double b00 = u2 * tmp4;
double tmp1 = 2 * b00;
double tmp1 = u2 / (u2 * aijkl + a1);
double tmp3 = tmp1 * aij;
double c0px = xkl - xk + tmp3 * xijxkl;
double c0py = ykl - yk + tmp3 * yijykl;
double c0pz = zkl - zk + tmp3 * zijzkl;
double c0px = tmp3 * xijxkl;
double c0py = tmp3 * yijykl;
double c0pz = tmp3 * zijzkl;
double g_0 = 1;
double g_1 = c0px;
double g_2 = 1;
Expand Down
2 changes: 1 addition & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ def initialize_with_default_plat_name(self):
install_requires=[
'pyscf~=2.7.0',
'pyscf-dispersion',
f'cupy-cuda{CUDA_VERSION}',
f'cupy-cuda{CUDA_VERSION}>=13.0', # Due to expm in cupyx.scipy.linalg and cutensor 2.0
'geometric',
f'gpu4pyscf-libxc-cuda{CUDA_VERSION}>=0.5',
]
Expand Down

0 comments on commit bf03d85

Please sign in to comment.