Skip to content

Commit

Permalink
bugfix: error in int2c2e (pyscf#212)
Browse files Browse the repository at this point in the history
  • Loading branch information
wxj6000 authored Sep 22, 2024
1 parent cb5a8b6 commit aa533f6
Show file tree
Hide file tree
Showing 5 changed files with 221 additions and 114 deletions.
112 changes: 0 additions & 112 deletions gpu4pyscf/df/int3c2e.py
Original file line number Diff line number Diff line change
Expand Up @@ -1453,118 +1453,6 @@ def get_int3c2e(mol, auxmol=None, auxbasis='weigend+etb', direct_scf_tol=1e-13,

return int3c.transpose([2,1,0])

def get_int2c2e_sorted(mol, auxmol, intopt=None, direct_scf_tol=1e-13, aosym=None, omega=None, stream=None):
'''
Generated int2c2e consistent with pyscf
'''
if omega is None: omega = 0.0
if stream is None: stream = cupy.cuda.get_current_stream()
if intopt is None:
intopt = VHFOpt(mol, auxmol, 'int2e')
intopt.build(direct_scf_tol, diag_block_with_triu=True, aosym=False)
naux = auxmol.nao
rows, cols = np.tril_indices(naux)

nbins = 1

nao_cart = intopt.mol.nao
naux_cart = intopt.auxmol.nao
norb_cart = nao_cart + naux_cart + 1

int2c = cupy.zeros([naux_cart, naux_cart], order='F')
ao_offsets = np.array([nao_cart+1, nao_cart, nao_cart+1, nao_cart], dtype=np.int32)
strides = np.array([1, naux_cart, naux_cart, naux_cart*naux_cart], dtype=np.int32)
for k_id, log_q_k in enumerate(intopt.aux_log_qs):
bins_locs_k = _make_s_index_offsets(log_q_k, nbins)
cp_k_id = k_id + len(intopt.log_qs)
for l_id, log_q_l in enumerate(intopt.aux_log_qs):
if k_id > l_id: continue
bins_locs_l = _make_s_index_offsets(log_q_l, nbins)
cp_l_id = l_id + len(intopt.log_qs)
err = libgint.GINTfill_int2e(
ctypes.cast(stream.ptr, ctypes.c_void_p),
intopt.bpcache,
ctypes.cast(int2c.data.ptr, ctypes.c_void_p),
ctypes.c_int(norb_cart),
strides.ctypes.data_as(ctypes.c_void_p),
ao_offsets.ctypes.data_as(ctypes.c_void_p),
bins_locs_k.ctypes.data_as(ctypes.c_void_p),
bins_locs_l.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(nbins),
ctypes.c_int(cp_k_id),
ctypes.c_int(cp_l_id),
ctypes.c_double(omega))

if err != 0:
raise RuntimeError("int2c2e failed\n")

int2c[rows, cols] = int2c[cols, rows]
if not mol.cart:
coeff = intopt.aux_cart2sph
int2c = coeff.T @ int2c @ coeff

return int2c

def get_int2c2e_ip_sorted(mol, auxmol, intopt=None, direct_scf_tol=1e-13, intor=None, aosym=None, stream=None):
'''
TODO: WIP
'''
if stream is None: stream = cupy.cuda.get_current_stream()
if intopt is None:
intopt = VHFOpt(mol, auxmol, 'int2e')
intopt.build(direct_scf_tol, diag_block_with_triu=True, aosym=False)

nbins = 1

nao_cart = intopt.mol.nao
naux_cart = intopt.auxmol.nao
norb_cart = nao_cart + naux_cart + 1
rows, cols = np.tril_indices(naux_cart)

int2c = cupy.zeros([naux_cart, naux_cart], order='F')
ao_offsets = np.array([nao_cart+1, nao_cart, nao_cart+1, nao_cart], dtype=np.int32)
strides = np.array([1, naux_cart, naux_cart, naux_cart*naux_cart], dtype=np.int32)
for k_id, log_q_k in enumerate(intopt.aux_log_qs):
bins_locs_k = _make_s_index_offsets(log_q_k, nbins)
cp_k_id = k_id + len(intopt.log_qs)
for l_id, log_q_l in enumerate(intopt.aux_log_qs):
if k_id > l_id: continue
bins_locs_l = _make_s_index_offsets(log_q_l, nbins)
cp_l_id = l_id + len(intopt.log_qs)
err = libgint.GINTfill_int2e(
ctypes.cast(stream.ptr, ctypes.c_void_p),
intopt.bpcache,
ctypes.cast(int2c.data.ptr, ctypes.c_void_p),
ctypes.c_int(norb_cart),
strides.ctypes.data_as(ctypes.c_void_p),
ao_offsets.ctypes.data_as(ctypes.c_void_p),
bins_locs_k.ctypes.data_as(ctypes.c_void_p),
bins_locs_l.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(nbins),
ctypes.c_int(cp_k_id),
ctypes.c_int(cp_l_id))

if err != 0:
raise RuntimeError("int2c2e failed\n")

int2c[rows, cols] = int2c[cols, rows]
if not auxmol.cart:
coeff = intopt.aux_cart2sph
int2c = coeff.T @ int2c @ coeff

return int2c

def get_int2c2e(mol, auxmol, direct_scf_tol=1e-13):
'''
Generate int2c2e on GPU
'''
intopt = VHFOpt(mol, auxmol, 'int2e')
intopt.build(direct_scf_tol, diag_block_with_triu=True, aosym=True)
int2c = get_int2c2e_sorted(mol, auxmol, intopt=intopt)
aux_idx = np.argsort(intopt.aux_ao_idx)
int2c = int2c[np.ix_(aux_idx, aux_idx)]
return int2c

def sort_mol(mol0, cart=True, log=None):
'''
# Sort basis according to angular momentum and contraction patterns so
Expand Down
1 change: 0 additions & 1 deletion gpu4pyscf/lib/gint/nr_fill_ao_ints.cu
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,6 @@ int GINTfill_int2e(cudaStream_t stream, BasisProdCache *bpcache, double *eri, in
ContractionProdType *cp_ij = bpcache->cptype + cp_ij_id;
ContractionProdType *cp_kl = bpcache->cptype + cp_kl_id;
GINTEnvVars envs;

int ng[4] = {0,0,0,0};
GINTinit_EnvVars(&envs, cp_ij, cp_kl, ng);
envs.omega = omega;
Expand Down
143 changes: 143 additions & 0 deletions gpu4pyscf/scf/int2c2e.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
# Copyright 2024 The GPU4PySCF Authors. All Rights Reserved.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.

import ctypes
import cupy
import numpy as np
from gpu4pyscf.scf.hf import _make_s_index_offsets
from gpu4pyscf.lib.cupy_helper import load_library, take_last2d
from gpu4pyscf.df.int3c2e import VHFOpt, make_fake_mol

libgint = load_library('libgint')

def get_int2c2e_sorted(mol, intopt=None, direct_scf_tol=1e-13, aosym=None, omega=None, stream=None):
'''
Generated int2c2e consistent with pyscf
'''
if omega is None: omega = 0.0
if stream is None: stream = cupy.cuda.get_current_stream()
if intopt is None:
# Reuse int3c2e
intopt = VHFOpt(mol, mol, 'int2e')
intopt.build(direct_scf_tol, diag_block_with_triu=True, aosym=True)
nao = mol.nao
rows, cols = np.tril_indices(nao)

nao_cart = intopt.mol.nao
norb_cart = nao_cart + 1

int2c = cupy.zeros([nao_cart, nao_cart], order='F')
ao_offsets = np.array([nao_cart+1, nao_cart, nao_cart+1, nao_cart], dtype=np.int32)
strides = np.array([1, nao_cart, nao_cart, nao_cart*nao_cart], dtype=np.int32)
log_cutoff = np.log(direct_scf_tol)
for k_id, log_q_k in enumerate(intopt.aux_log_qs):
#bins_locs_k = _make_s_index_offsets(log_q_k, nbins)
bins_locs_k = np.array([0, len(log_q_k)], dtype=np.int32)
bins_floor_k = np.array([100], dtype=np.double)
cp_k_id = k_id + len(intopt.log_qs)
for l_id, log_q_l in enumerate(intopt.aux_log_qs):
if k_id > l_id: continue
#bins_locs_l = _make_s_index_offsets(log_q_l, nbins)
bins_locs_l = np.array([0, len(log_q_l)], dtype=np.int32)
bins_floor_l = np.array([100], dtype=np.double)
cp_l_id = l_id + len(intopt.log_qs)
nbins_locs_k = len(bins_locs_k) - 1
nbins_locs_l = len(bins_locs_l) - 1
err = libgint.GINTfill_int2e(
ctypes.cast(stream.ptr, ctypes.c_void_p),
intopt.bpcache,
ctypes.cast(int2c.data.ptr, ctypes.c_void_p),
ctypes.c_int(norb_cart),
strides.ctypes.data_as(ctypes.c_void_p),
ao_offsets.ctypes.data_as(ctypes.c_void_p),
bins_locs_k.ctypes.data_as(ctypes.c_void_p),
bins_locs_l.ctypes.data_as(ctypes.c_void_p),
bins_floor_k.ctypes.data_as(ctypes.c_void_p),
bins_floor_l.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(nbins_locs_k),
ctypes.c_int(nbins_locs_l),
ctypes.c_int(cp_k_id),
ctypes.c_int(cp_l_id),
ctypes.c_double(log_cutoff),
ctypes.c_double(omega))

if err != 0:
raise RuntimeError("int2c2e failed\n")

int2c[rows, cols] = int2c[cols, rows]
if not mol.cart:
coeff = intopt.cart2sph
int2c = coeff.T @ int2c @ coeff

return int2c

def get_int2c2e_ip_sorted(mol, auxmol, intopt=None, direct_scf_tol=1e-13, intor=None, aosym=None, stream=None):
'''
TODO: WIP
'''
if stream is None: stream = cupy.cuda.get_current_stream()
if intopt is None:
intopt = VHFOpt(mol, auxmol, 'int2e')
intopt.build(direct_scf_tol, diag_block_with_triu=True, aosym=False)

nbins = 1

nao_cart = intopt.mol.nao
naux_cart = intopt.auxmol.nao
norb_cart = nao_cart + naux_cart + 1
rows, cols = np.tril_indices(naux_cart)

int2c = cupy.zeros([naux_cart, naux_cart], order='F')
ao_offsets = np.array([nao_cart+1, nao_cart, nao_cart+1, nao_cart], dtype=np.int32)
strides = np.array([1, naux_cart, naux_cart, naux_cart*naux_cart], dtype=np.int32)
for k_id, log_q_k in enumerate(intopt.aux_log_qs):
bins_locs_k = _make_s_index_offsets(log_q_k, nbins)
cp_k_id = k_id + len(intopt.log_qs)
for l_id, log_q_l in enumerate(intopt.aux_log_qs):
if k_id > l_id: continue
bins_locs_l = _make_s_index_offsets(log_q_l, nbins)
cp_l_id = l_id + len(intopt.log_qs)
err = libgint.GINTfill_int2e(
ctypes.cast(stream.ptr, ctypes.c_void_p),
intopt.bpcache,
ctypes.cast(int2c.data.ptr, ctypes.c_void_p),
ctypes.c_int(norb_cart),
strides.ctypes.data_as(ctypes.c_void_p),
ao_offsets.ctypes.data_as(ctypes.c_void_p),
bins_locs_k.ctypes.data_as(ctypes.c_void_p),
bins_locs_l.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(nbins),
ctypes.c_int(cp_k_id),
ctypes.c_int(cp_l_id))

if err != 0:
raise RuntimeError("int2c2e failed\n")

int2c[rows, cols] = int2c[cols, rows]
if not auxmol.cart:
coeff = intopt.aux_cart2sph
int2c = coeff.T @ int2c @ coeff

return int2c

def get_int2c2e(mol, direct_scf_tol=1e-13):
'''
Generate int2c2e on GPU
'''
intopt = VHFOpt(mol, mol, 'int2e')
intopt.build(direct_scf_tol, diag_block_with_triu=True, aosym=True)
int2c = get_int2c2e_sorted(mol, intopt=intopt)
int2c = take_last2d(int2c, intopt.rev_ao_idx)
return int2c
3 changes: 2 additions & 1 deletion gpu4pyscf/scf/int4c2e.py
Original file line number Diff line number Diff line change
Expand Up @@ -438,4 +438,5 @@ def loop_int4c2e_general(intopt, ip_type='', direct_scf_tol=1e-13, omega=None, s
raise RuntimeError(f'GINT_fill_int4c2e general failed, err={err}')
if cp_ij_id == cp_kl_id:
int4c *= 0.5
yield i0,i1,j0,j1,k0,k1,l0,l1,int4c
yield i0,i1,j0,j1,k0,k1,l0,l1,int4c

76 changes: 76 additions & 0 deletions gpu4pyscf/scf/tests/test_int2c2e.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
# Copyright 2023 The GPU4PySCF Authors. All Rights Reserved.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.

import pyscf
from pyscf import lib, df
from pyscf.scf import _vhf
from pyscf.gto.moleintor import getints, make_cintopt
from pyscf.df.grad.rhf import _int3c_wrapper
import numpy as np
import cupy
import unittest

from gpu4pyscf.scf import int2c2e
from gpu4pyscf.lib.cupy_helper import load_library
libgint = load_library('libgint')

'''
compare int2c2e by pyscf and gpu4pyscf
'''

def setUpModule():
global mol_cart, mol_sph
mol_cart = pyscf.M(atom='''
H -0.7570000000 -0.0000000000 -0.4696000000
H 0.7570000000 0.0000000000 -0.4696000000
''',
basis= 'ccpvdz',
verbose=1,
output = '/dev/null')
mol_cart.build()

mol_sph = pyscf.M(atom='''
H -0.7570000000 -0.0000000000 -0.4696000000
H 0.7570000000 0.0000000000 -0.4696000000
''',
basis= 'ccpvdz',
verbose=1,
output = '/dev/null')
mol_sph.build()

omega = 0.2

def tearDownModule():
global mol_sph, mol_cart
mol_sph.stdout.close()
mol_cart.stdout.close()
del mol_sph, mol_cart

class KnownValues(unittest.TestCase):
def test_int2c2e_sph(self):
nao = mol_sph.nao
eri = mol_sph.intor('int2c2e_sph').reshape((nao,)*2)
int2c = int2c2e.get_int2c2e(mol_sph)
assert np.linalg.norm(eri - int2c.get()) < 1e-10

def test_int2c2e_cart(self):
nao = mol_cart.nao
eri = mol_cart.intor('int2c2e_cart').reshape((nao,)*2)
int2c = int2c2e.get_int2c2e(mol_cart)
assert np.linalg.norm(eri - int2c.get()) < 1e-10

if __name__ == "__main__":
print("Full Tests for int2c")
unittest.main()

0 comments on commit aa533f6

Please sign in to comment.