From aae58feab0f97b2acbe9c8401c78ee86c91c1a49 Mon Sep 17 00:00:00 2001 From: Qiming Sun Date: Thu, 2 Jan 2025 07:35:10 -0800 Subject: [PATCH] Add unittest pipeline --- .github/workflows/unittest.yml | 19 ++ pyscf/pbc/prop/efg/rhf.py | 166 +++++++++--------- ...t_integral.py => test_pbc_efg_integral.py} | 8 +- ..._rhf.py => test_pbc_rhf_polarizability.py} | 36 ++-- pyscf/prop/gtensor/test/test_uks.py | 76 +++++--- pyscf/prop/gtensor/uks.py | 3 +- .../test/test_magnetizability_integral.py | 123 +++++++++++++ pyscf/prop/nmr/test/test_rks_nmr.py | 3 +- pyscf/prop/nmr/test/test_uks_nmr.py | 7 +- pyscf/prop/polarizability/rhf.py | 62 ------- ...test_rhf.py => test_rhf_polarizability.py} | 56 +++++- ..._rhf.py => test_rhf_rotational_gtensor.py} | 0 pyscf/prop/ssc/test/test_scf_integral.py | 40 +++-- 13 files changed, 384 insertions(+), 215 deletions(-) create mode 100644 .github/workflows/unittest.yml rename pyscf/pbc/prop/efg/test/{test_integral.py => test_pbc_efg_integral.py} (95%) rename pyscf/pbc/prop/polarizability/test/{test_rhf.py => test_pbc_rhf_polarizability.py} (76%) create mode 100644 pyscf/prop/magnetizability/test/test_magnetizability_integral.py rename pyscf/prop/polarizability/test/{test_rhf.py => test_rhf_polarizability.py} (62%) rename pyscf/prop/rotational_gtensor/test/{test_rhf.py => test_rhf_rotational_gtensor.py} (100%) diff --git a/.github/workflows/unittest.yml b/.github/workflows/unittest.yml new file mode 100644 index 0000000..11a649b --- /dev/null +++ b/.github/workflows/unittest.yml @@ -0,0 +1,19 @@ +name: unittest + +on: [push, pull_request] + +jobs: + linux-build: + runs-on: ubuntu-22.04 + strategy: + fail-fast: false + steps: + - uses: actions/checkout@v3 + - name: Set up Python 3.10 + uses: actions/setup-python@v4 + with: + python-version: 3.10 + - name: Install + run: pip install pyscf packaging pytest + - name: Test + run: pytest -v pyscf diff --git a/pyscf/pbc/prop/efg/rhf.py b/pyscf/pbc/prop/efg/rhf.py index a67beb7..bd1f6c9 100644 --- a/pyscf/pbc/prop/efg/rhf.py +++ b/pyscf/pbc/prop/efg/rhf.py @@ -29,14 +29,16 @@ ''' -import numpy +import numpy as np import scipy.special from pyscf import lib -from pyscf.pbc import gto +from pyscf.gto.mole import gaussian_int +from pyscf.pbc import gto as pbcgto from pyscf.pbc import scf -from pyscf.pbc import tools +from pyscf.pbc import tools as pbctools from pyscf.pbc import df from pyscf.pbc.df import aft +from pyscf.pbc.df.gdf_builder import estimate_eta_for_ke_cutoff from pyscf.prop.efg import rhf as rhf_efg def kernel(method, efg_nuc=None): @@ -47,14 +49,14 @@ def kernel(method, efg_nuc=None): dm = method.make_rdm1() if isinstance(method, scf.khf.KSCF): - if isinstance(dm[0][0], numpy.ndarray) and dm[0][0].ndim == 2: + if isinstance(dm[0][0], np.ndarray) and dm[0][0].ndim == 2: dm = dm[0] + dm[1] # KUHF density matrix elif isinstance(method, scf.hf.SCF): - if isinstance(dm[0], numpy.ndarray) and dm[0].ndim == 2: + if isinstance(dm[0], np.ndarray) and dm[0].ndim == 2: dm = dm[0] + dm[1] # UHF density matrix else: mo = method.mo_coeff - if isinstance(dm[0][0], numpy.ndarray) and dm[0][0].ndim == 2: + if isinstance(dm[0][0], np.ndarray) and dm[0][0].ndim == 2: dm_a = [lib.einsum('pi,ij,qj->pq', c, dm[0][k], c.conj()) for k, c in enumerate(mo)] dm_b = [lib.einsum('pi,ij,qj->pq', c, dm[1][k], c.conj()) @@ -86,7 +88,7 @@ def kernel(method, efg_nuc=None): rhf_efg._analyze(cell, atm_id, v, log) - return numpy.asarray(efg) + return np.asarray(efg) EFG = kernel @@ -98,45 +100,47 @@ def _get_quad_nuc(cell, atm_id): Lall = cell.get_lattice_Ls(rcut=ew_cut) rLij = coords[atm_id,:] - coords + Lall[:,None,:] - rr = numpy.einsum('Ljx,Ljy->Ljxy', rLij, rLij) - r = numpy.sqrt(numpy.einsum('Ljxx->Lj', rr)) + rr = np.einsum('Ljx,Ljy->Ljxy', rLij, rLij) + r = np.sqrt(np.einsum('Ljxx->Lj', rr)) r[r<1e-16] = 1e60 - idx = numpy.arange(3) + idx = np.arange(3) erfc_part = scipy.special.erfc(ew_eta * r) / r**5 - ewovrl = 3 * chargs[atm_id] * numpy.einsum('Ljxy,j,Lj->xy', rr, chargs, erfc_part) + ewovrl = 3 * chargs[atm_id] * np.einsum('Ljxy,j,Lj->xy', rr, chargs, erfc_part) ewovrl[idx,idx] -= ewovrl.trace() / 3 - exp_part = numpy.exp(-ew_eta**2 * r**2) / r**4 - ewovrl_part = (2./numpy.sqrt(numpy.pi) * ew_eta * chargs[atm_id] * - numpy.einsum('Ljxy,j,Lj->xy', rr, chargs, exp_part)) + exp_part = np.exp(-ew_eta**2 * r**2) / r**4 + ewovrl_part = (2./np.sqrt(np.pi) * ew_eta * chargs[atm_id] * + np.einsum('Ljxy,j,Lj->xy', rr, chargs, exp_part)) ewovrl += ewovrl_part ewovrl += 2*ewovrl_part ewovrl[idx,idx] -= ewovrl_part.trace() - exp_part = numpy.exp(-ew_eta**2 * r**2) / r**2 - ewovrl += (4./numpy.sqrt(numpy.pi) * ew_eta**3 * chargs[atm_id] * - numpy.einsum('Ljxy,j,Lj->xy', rr, chargs, exp_part)) + exp_part = np.exp(-ew_eta**2 * r**2) / r**2 + ewovrl += (4./np.sqrt(np.pi) * ew_eta**3 * chargs[atm_id] * + np.einsum('Ljxy,j,Lj->xy', rr, chargs, exp_part)) # Fermi contact term - ewovrl_fc = 4./3*ew_eta**3/numpy.pi**.5 * numpy.exp(-ew_eta**2*r**2) - ewovrl[idx,idx] -= numpy.einsum('j,Lj->', chargs, ewovrl_fc) + ewovrl_fc = 4./3*ew_eta**3/np.pi**.5 * np.exp(-ew_eta**2*r**2) + ewovrl[idx,idx] -= np.einsum('j,Lj->', chargs, ewovrl_fc) - mesh = gto.cell._cut_mesh_for_ewald(cell, cell.mesh) + log_precision = np.log(cell.precision / (chargs.sum()*16*np.pi**2)) + ke_cutoff = -2*ew_eta**2*log_precision + mesh = cell.cutoff_to_mesh(ke_cutoff) Gv, Gvbase, weights = cell.get_Gv_weights(mesh) - GG = numpy.einsum('gx,gy->gxy', Gv, Gv) - absG2 = numpy.einsum('gxx->g', GG) + GG = np.einsum('gx,gy->gxy', Gv, Gv) + absG2 = np.einsum('gxx->g', GG) # Corresponding to FC term, that makes the tensor zero trace - idx = numpy.arange(3) + idx = np.arange(3) GG[:,idx,idx] -= 1./3 * absG2[:,None] absG2[absG2==0] = 1e200 - coulG = 4*numpy.pi / absG2 + coulG = 4*np.pi / absG2 coulG *= weights - expG2 = numpy.exp(-absG2/(4*ew_eta**2)) + expG2 = np.exp(-absG2/(4*ew_eta**2)) coulG *= expG2 SI = cell.get_SI(Gv) - ZSI = numpy.einsum("i,ij->j", chargs, SI) - ewg =-chargs[atm_id] * numpy.einsum('ixy,i,i,i->xy', + ZSI = np.einsum("i,ij->j", chargs, SI) + ewg =-chargs[atm_id] * np.einsum('ixy,i,i,i->xy', GG, SI[atm_id].conj(), ZSI, coulG).real return ewovrl + ewg @@ -150,7 +154,7 @@ def _fft_quad_integrals(mydf, dm, efg_nuc): mesh = mydf.mesh kpts = mydf.kpts - kpts_lst = numpy.reshape(kpts, (-1,3)) + kpts_lst = np.reshape(kpts, (-1,3)) nkpts = len(kpts_lst) nao = cell.nao_nr() dm_kpts = dm.reshape((nkpts,nao,nao), order='C') @@ -158,24 +162,24 @@ def _fft_quad_integrals(mydf, dm, efg_nuc): ni = mydf._numint hermi = 1 make_rho, nset, nao = ni._gen_rho_evaluator(cell, dm_kpts, hermi) - ngrids = numpy.prod(mesh) - rhoR = numpy.zeros(ngrids) + ngrids = np.prod(mesh) + rhoR = np.zeros(ngrids) for ao_ks_etc, p0, p1 in mydf.aoR_loop(mydf.grids, kpts_lst): ao_ks, mask = ao_ks_etc[0], ao_ks_etc[2] rhoR[p0:p1] += make_rho(0, ao_ks, mask, 'LDA') ao_ks = None - rhoG = tools.fft(rhoR, mesh) + rhoG = pbctools.fft(rhoR, mesh) Gv = cell.get_Gv(mesh) - coulG = tools.get_coulG(cell, mesh=mesh, Gv=Gv) - GG = numpy.einsum('gx,gy->gxy', Gv, Gv) - absG2 = numpy.einsum('gxx->g', GG) + coulG = pbctools.get_coulG(cell, mesh=mesh, Gv=Gv) + GG = np.einsum('gx,gy->gxy', Gv, Gv) + absG2 = np.einsum('gxx->g', GG) # Corresponding to FC term, that makes the tensor traceless - idx = numpy.arange(3) + idx = np.arange(3) GG[:,idx,idx] -= 1./3 * absG2[:,None] - vG = 1./ngrids * numpy.einsum('g,g,gxy->gxy', rhoG, coulG, GG) + vG = 1./ngrids * np.einsum('g,g,gxy->gxy', rhoG, coulG, GG) SI = cell.get_SI(Gv) efg_e = lib.einsum('zg,gxy->zxy', SI[efg_nuc], vG.conj()).real return efg_e @@ -192,47 +196,60 @@ def _aft_quad_integrals(mydf, dm, efg_nuc): mesh = mydf.mesh kpts = mydf.kpts - kpts_lst = numpy.reshape(kpts, (-1,3)) + kpts_lst = np.reshape(kpts, (-1,3)) nkpts = len(kpts_lst) nao = cell.nao_nr() dm_kpts = dm.reshape((nkpts,nao,nao), order='C') - ngrids = numpy.prod(mesh) - rhoG = numpy.zeros(ngrids, dtype=numpy.complex128) - kpt_allow = numpy.zeros(3) + ngrids = np.prod(mesh) + rhoG = np.zeros(ngrids, dtype=np.complex128) + kpt_allow = np.zeros(3) max_memory = max(2000, mydf.max_memory-lib.current_memory()[0]) for aoaoks, p0, p1 in mydf.ft_loop(mesh, kpt_allow, kpts_lst, max_memory=max_memory, aosym='s1'): - rhoG[p0:p1] = numpy.einsum('kgpq,kqp->g', aoaoks.reshape(nkpts,p1-p0,nao,nao), + rhoG[p0:p1] = np.einsum('kgpq,kqp->g', aoaoks.reshape(nkpts,p1-p0,nao,nao), dm_kpts) t1 = log.timer_debug1('contracting Vnuc [%s:%s]'%(p0, p1), *t1) t0 = log.timer_debug1('contracting Vnuc', *t0) rhoG *= 1./nkpts Gv, Gvbase, kws = cell.get_Gv_weights(mesh) - coulG = tools.get_coulG(cell, mesh=mesh, Gv=Gv) - GG = numpy.einsum('gx,gy->gxy', Gv, Gv) - absG2 = numpy.einsum('gxx->g', GG) + coulG = pbctools.get_coulG(cell, mesh=mesh, Gv=Gv) + GG = np.einsum('gx,gy->gxy', Gv, Gv) + absG2 = np.einsum('gxx->g', GG) # Corresponding to FC term, that makes the tensor traceless - idx = numpy.arange(3) + idx = np.arange(3) GG[:,idx,idx] -= 1./3 * absG2[:,None] - vG = 1./cell.vol * numpy.einsum('g,g,gxy->gxy', rhoG, coulG, GG) + vG = 1./cell.vol * np.einsum('g,g,gxy->gxy', rhoG, coulG, GG) - if mydf.eta == 0: - SI = cell.get_SI(Gv) - efg_e = numpy.einsum('zg,gxy->zxy', SI[efg_nuc], vG.conj()).real + assert cell.dimension == 3 + a = cell.lattice_vectors() + ke_cutoff = min(pbctools.mesh_to_cutoff(a, mesh)) + eta = estimate_eta_for_ke_cutoff(cell, ke_cutoff, cell.precision) - else: - nuccell = aft._compensate_nuccell(mydf) - # PP-loc part1 is handled by fakenuc in _int_nuc_vloc - efg_e = _int_nuc_vloc(mydf, nuccell, kpts_lst, dm_kpts) - t0 = log.timer_debug1('vnuc pass1: analytic int', *t0) + nuccell = _compensate_nuccell(cell, eta) + # PP-loc part1 is handled by fakenuc in _int_nuc_vloc + efg_e = _int_nuc_vloc(mydf, nuccell, kpts_lst, dm_kpts) + t0 = log.timer_debug1('vnuc pass1: analytic int', *t0) - aoaux = df.ft_ao.ft_ao(nuccell, Gv) - efg_e += numpy.einsum('gz,gxy->zxy', aoaux[:,efg_nuc], vG.conj()).real + aoaux = df.ft_ao.ft_ao(nuccell, Gv) + efg_e += np.einsum('gz,gxy->zxy', aoaux[:,efg_nuc], vG.conj()).real return efg_e +def _compensate_nuccell(cell, eta): + '''A cell of the compensated Gaussian charges for nucleus''' + modchg_cell = cell.copy(deep=False) + half_sph_norm = .5/np.sqrt(np.pi) + norm = half_sph_norm/gaussian_int(2, eta) + chg_env = [eta, norm] + ptr_eta = cell._env.size + ptr_norm = ptr_eta + 1 + chg_bas = [[ia, 0, 1, 1, 0, ptr_eta, ptr_norm, 0] for ia in range(cell.natm)] + modchg_cell._atm = cell._atm + modchg_cell._bas = np.asarray(chg_bas, dtype=np.int32) + modchg_cell._env = np.hstack((cell._env, chg_env)) + return modchg_cell def _int_nuc_vloc(mydf, nuccell, kpts, dm_kpts): '''Vnuc - Vloc''' @@ -242,10 +259,10 @@ def _int_nuc_vloc(mydf, nuccell, kpts, dm_kpts): # Use the 3c2e code with steep s gaussians to mimic nuclear density fakenuc = aft._fake_nuc(cell) fakenuc._atm, fakenuc._bas, fakenuc._env = \ - gto.conc_env(nuccell._atm, nuccell._bas, nuccell._env, - fakenuc._atm, fakenuc._bas, fakenuc._env) + pbcgto.conc_env(nuccell._atm, nuccell._bas, nuccell._env, + fakenuc._atm, fakenuc._bas, fakenuc._env) - kptij_lst = numpy.hstack((kpts,kpts)).reshape(-1,2,3) + kptij_lst = np.hstack((kpts,kpts)).reshape(-1,2,3) v3c = df.incore.aux_e2(cell, fakenuc, 'int3c2e_ipip1', aosym='s1', comp=9, kptij_lst=kptij_lst) v3c += df.incore.aux_e2(cell, fakenuc, 'int3c2e_ipvip1', aosym='s1', comp=9, @@ -254,46 +271,33 @@ def _int_nuc_vloc(mydf, nuccell, kpts, dm_kpts): nao = cell.nao_nr() natm = cell.natm v3c = v3c.reshape(nkpts,3,3,nao,nao,natm*2) - efg_loc = 1./nkpts * numpy.einsum('kxypqz,kqp->zxy', v3c, dm_kpts) - efg_loc+= 1./nkpts * numpy.einsum('kyxqpz,kqp->zxy', v3c, dm_kpts) + efg_loc = 1./nkpts * np.einsum('kxypqz,kqp->zxy', v3c, dm_kpts) + efg_loc+= 1./nkpts * np.einsum('kyxqpz,kqp->zxy', v3c, dm_kpts) v3c = None # Fermi contact fc = df.incore.aux_e2(cell, fakenuc, 'int3c1e', aosym='s1', kptij_lst=kptij_lst) fc = fc.reshape(nkpts,nao,nao,natm*2) - vfc = 1./nkpts * numpy.einsum('kpqz,kqp->z', fc, dm_kpts) + vfc = 1./nkpts * np.einsum('kpqz,kqp->z', fc, dm_kpts) for i in range(3): - efg_loc[:,i,i] += 4./3*numpy.pi * vfc + efg_loc[:,i,i] += 4./3*np.pi * vfc nuc_part = efg_loc[:natm] modchg_part = efg_loc[natm:] efg_loc = nuc_part - modchg_part if cell.dimension == 3: - fac = numpy.pi/cell.vol - nucbar = numpy.array([fac/nuccell.bas_exp(i)[0] for i in range(natm)]) + fac = np.pi/cell.vol + nucbar = np.array([fac/nuccell.bas_exp(i)[0] for i in range(natm)]) ipip = cell.pbc_intor('int1e_ipipovlp', 9, lib.HERMITIAN, kpts) ipvip = cell.pbc_intor('int1e_ipovlpip', 9, lib.HERMITIAN, kpts) d2rho_bar = 0 for k in range(nkpts): v = (ipip[k] + ipvip[k]).reshape(3,3,nao,nao) - d2rho_bar += numpy.einsum('xypq,qp->xy', v, dm_kpts[k]) - d2rho_bar += numpy.einsum('yxqp,qp->xy', v, dm_kpts[k]) + d2rho_bar += np.einsum('xypq,qp->xy', v, dm_kpts[k]) + d2rho_bar += np.einsum('yxqp,qp->xy', v, dm_kpts[k]) d2rho_bar *= 1./nkpts - efg_loc -= numpy.einsum('z,xy->zxy', nucbar, d2rho_bar) + efg_loc -= np.einsum('z,xy->zxy', nucbar, d2rho_bar) return efg_loc.real - - -if __name__ == '__main__': - cell = gto.Cell() - cell.verbose = 4 - cell.atom = 'H 1 0.8 0; H 0. 1.3 0' - cell.a = numpy.eye(3) * 3 - cell.basis = 'gth-szv' - cell.pseudo = 'gth-pade' - #cell.mesh = [24]*3 - cell.build() - mf = scf.RHF(cell).run() - v = EFG(mf) diff --git a/pyscf/pbc/prop/efg/test/test_integral.py b/pyscf/pbc/prop/efg/test/test_pbc_efg_integral.py similarity index 95% rename from pyscf/pbc/prop/efg/test/test_integral.py rename to pyscf/pbc/prop/efg/test/test_pbc_efg_integral.py index 9494f86..db3a347 100644 --- a/pyscf/pbc/prop/efg/test/test_integral.py +++ b/pyscf/pbc/prop/efg/test/test_pbc_efg_integral.py @@ -39,7 +39,9 @@ def ewald_deriv1(cell, atm_id): np.einsum('Ljx,j,Lj->x', rLij, chargs, np.exp(-ew_eta**2 * r**2) / r**2)) - mesh = gto.cell._cut_mesh_for_ewald(cell, cell.mesh) + log_precision = np.log(cell.precision / (chargs.sum()*16*np.pi**2)) + ke_cutoff = -2*ew_eta**2*log_precision + mesh = cell.cutoff_to_mesh(ke_cutoff) Gv, Gvbase, weights = cell.get_Gv_weights(mesh) absG2 = np.einsum('gi,gi->g', Gv, Gv) absG2[absG2==0] = 1e200 @@ -80,7 +82,9 @@ def ewald_deriv2(cell, atm_id): ewovrl += (4./np.sqrt(np.pi) * ew_eta**3 * chargs[atm_id] * np.einsum('Ljxy,j,Lj->xy', rr, chargs, exp_part)) - mesh = gto.cell._cut_mesh_for_ewald(cell, cell.mesh) + log_precision = np.log(cell.precision / (chargs.sum()*16*np.pi**2)) + ke_cutoff = -2*ew_eta**2*log_precision + mesh = cell.cutoff_to_mesh(ke_cutoff) Gv, Gvbase, weights = cell.get_Gv_weights(mesh) GG = np.einsum('gi,gj->gij', Gv, Gv) absG2 = np.einsum('gi,gi->g', Gv, Gv) diff --git a/pyscf/pbc/prop/polarizability/test/test_rhf.py b/pyscf/pbc/prop/polarizability/test/test_pbc_rhf_polarizability.py similarity index 76% rename from pyscf/pbc/prop/polarizability/test/test_rhf.py rename to pyscf/pbc/prop/polarizability/test/test_pbc_rhf_polarizability.py index 14dc3c2..c161057 100644 --- a/pyscf/pbc/prop/polarizability/test/test_rhf.py +++ b/pyscf/pbc/prop/polarizability/test/test_pbc_rhf_polarizability.py @@ -19,25 +19,27 @@ from pyscf.pbc import gto, scf from pyscf.pbc.prop.polarizability import rhf -cell = gto.Cell() -cell.atom = """H 0.0 0.0 0.0 - F 0.9 0.0 0.0 - """ -cell.basis = 'sto-3g' -cell.a = [[2.82, 0, 0], [0, 2.82, 0], [0, 0, 2.82]] -cell.dimension = 1 -cell.precision = 1e-10 -cell.output = '/dev/null' -cell.build() +def setUpModule(): + global cell, kmf, polar + cell = gto.Cell() + cell.atom = """H 0.0 0.0 0.0 + F 0.9 0.0 0.0 + """ + cell.basis = 'sto-3g' + cell.a = [[2.82, 0, 0], [0, 2.82, 0], [0, 0, 2.82]] + cell.dimension = 1 + cell.precision = 1e-10 + cell.output = '/dev/null' + cell.build() -kpts = cell.make_kpts([16,1,1]) -kmf = scf.KRHF(cell, kpts=kpts, exxdiv="ewald").density_fit() -kmf.kernel() -polar = rhf.Polarizability(kmf, kpts) + kpts = cell.make_kpts([16,1,1]) + kmf = scf.KRHF(cell, kpts=kpts, exxdiv="ewald").density_fit() + kmf.kernel() + polar = rhf.Polarizability(kmf, kpts) def tearDownModule(): - global cell, kmf - del cell, kmf + global cell, kmf, polar + del cell, kmf, polar class KnownValues(unittest.TestCase): def test_dip_moment(self): @@ -63,7 +65,7 @@ def test_polarizability_with_freq(self): def test_hyper_polarizability(self): e3 = polar.hyper_polarizability() - self.assertAlmostEqual(e3[0,0,0], 3.02691297, 6) + self.assertAlmostEqual(e3[0,0,0], 3.02691297, 5) self.assertAlmostEqual(e3[0,1,1], -1.97889278e-02, 6) self.assertAlmostEqual(e3[0,2,2], e3[0,1,1], 9) diff --git a/pyscf/prop/gtensor/test/test_uks.py b/pyscf/prop/gtensor/test/test_uks.py index b025f06..2ad89ea 100644 --- a/pyscf/prop/gtensor/test/test_uks.py +++ b/pyscf/prop/gtensor/test/test_uks.py @@ -19,40 +19,61 @@ from pyscf import gto, lib, scf, dft from pyscf.prop import gtensor from pyscf.data import nist -nist.ALPHA = 1./137.03599967994 -mol = gto.Mole() -mol.verbose = 7 -mol.output = '/dev/null' -mol.atom = ''' - H 0. , 0. , .917 - F 0. , 0. , 0.''' -mol.basis = 'ccpvdz' -mol.spin = 1 -mol.charge = 1 -mol.build() +def setUpModule(): + global mol, mf, dm0, dm1, ALPHA_backup + ALPHA_backup = nist.ALPHA + nist.ALPHA = 1./137.03599967994 -mf = dft.UKS(mol) -mf.xc = 'b3lyp' -mf.conv_tol_grad = 1e-6 -mf.conv_tol = 1e-12 -mf.kernel() + mol = gto.Mole() + mol.verbose = 7 + mol.output = '/dev/null' + mol.atom = ''' + H 0. , 0. , .917 + F 0. , 0. , 0.''' + mol.basis = 'ccpvdz' + mol.spin = 1 + mol.charge = 1 + mol.build() -nao = mol.nao_nr() -numpy.random.seed(1) -dm0 = numpy.random.random((2,nao,nao)) -dm0 = dm0 + dm0.transpose(0,2,1) -dm1 = numpy.random.random((2,3,nao,nao)) -dm1 = dm1 - dm1.transpose(0,1,3,2) + with lib.temporary_env(dft.radi, ATOM_SPECIFIC_TREUTLER_GRIDS=False): + mf = dft.UKS(mol) + mf.xc = 'b3lyp' + mf.conv_tol_grad = 1e-6 + mf.conv_tol = 1e-12 + mf.kernel() + + nao = mol.nao_nr() + numpy.random.seed(1) + dm0 = numpy.random.random((2,nao,nao)) + dm0 = dm0 + dm0.transpose(0,2,1) + dm0 = numpy.einsum('xpi,xqi->xpq', dm0, dm0, optimize=True) + dm1 = numpy.random.random((2,3,nao,nao)) + dm1 = dm1 - dm1.transpose(0,1,3,2) + +def tearDownModule(): + global mol, mf, dm0, dm1, ALPHA_backup + mol.stdout.close() + del mol, mf, dm0, dm1 + nist.ALPHA = ALPHA_backup class KnowValues(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.original_grids = dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS + dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = False + + @classmethod + def tearDownClass(cls): + dft.radi.ATOM_SPECIFIC_TREUTLER_GRIDS = cls.original_grids + def test_nr_lda_para_soc2e(self): mf1 = copy.copy(mf) mf1.xc = 'lda,vwn' g = gtensor.uks.GTensor(mf1) g.para_soc2e = 'SSO' dat = g.make_para_soc2e(dm0, dm1, 1) - self.assertAlmostEqual(lib.finger(dat), -0.008807083583654644, 9) + self.assertAlmostEqual(lib.fp(dat), -0.201859494, 9) def test_nr_bp86_para_soc2e(self): mf1 = copy.copy(mf) @@ -60,16 +81,15 @@ def test_nr_bp86_para_soc2e(self): g = gtensor.uks.GTensor(mf1) g.para_soc2e = 'SSO' dat = g.make_para_soc2e(dm0, dm1, 1) - self.assertAlmostEqual(lib.finger(dat), -0.0088539747015796387, 9) + self.assertAlmostEqual(lib.fp(dat), -0.201794075, 7) def test_nr_b3lyp_para_soc2e(self): mf1 = copy.copy(mf) mf1.xc = 'b3lyp' g = gtensor.uks.GTensor(mf1) g.para_soc2e = 'SSO' - dm = numpy.einsum('xpi,xqi->xpq', dm0, dm0) - dat = g.make_para_soc2e(dm, dm1, 1) - self.assertAlmostEqual(lib.finger(dat), -0.20729647343641752, 9) + dat = g.make_para_soc2e(dm0, dm1, 1) + self.assertAlmostEqual(lib.fp(dat), -0.207296515, 7) def test_nr_uks(self): g = gtensor.uhf.GTensor(mf) @@ -78,7 +98,7 @@ def test_nr_uks(self): g.so_eff_charge = True g.cphf = False dat = g.kernel() - self.assertAlmostEqual(numpy.linalg.norm(dat), 3.47479197036, 6) + self.assertAlmostEqual(numpy.linalg.norm(dat), 3.4748005, 5) if __name__ == "__main__": diff --git a/pyscf/prop/gtensor/uks.py b/pyscf/prop/gtensor/uks.py index 910bcf7..f1ee195 100644 --- a/pyscf/prop/gtensor/uks.py +++ b/pyscf/prop/gtensor/uks.py @@ -157,7 +157,8 @@ def get_vxc_soc(ni, mol, grids, xc_code, dms, max_memory=2000, verbose=None): make_rho, nset, nao = ni._gen_rho_evaluator(mol, dms, hermi=1) ngrids = len(grids.weights) BLKSIZE = numint.BLKSIZE - blksize = min(int(max_memory/12*1e6/8/nao/BLKSIZE)*BLKSIZE, ngrids) + blksize = int(max_memory*1e6/(12*nao*8*BLKSIZE)) + blksize = max(4, min(blksize, ngrids//BLKSIZE+1, 1200)) * BLKSIZE shls_slice = (0, mol.nbas) ao_loc = mol.ao_loc_nr() diff --git a/pyscf/prop/magnetizability/test/test_magnetizability_integral.py b/pyscf/prop/magnetizability/test/test_magnetizability_integral.py new file mode 100644 index 0000000..fb3a00f --- /dev/null +++ b/pyscf/prop/magnetizability/test/test_magnetizability_integral.py @@ -0,0 +1,123 @@ +import numpy +from pyscf import lib +from pyscf import gto, dft +from pyscf.scf import jk +from pyscf.dft import numint + +def test_integral(): + mol = gto.M(atom='''O 0. 0. 0. + H 0. -0.757 0.587 + H 0. 0.757 0.587''', + basis='ccpvdz') + + nao = mol.nao + dm0 = numpy.random.random((nao,nao)) + dm0 = dm0 + dm0.T + + v1 = mol.intor('int2e_gg1', comp=9).reshape(3,3,nao,nao,nao,nao) + v2 = mol.intor('int2e_g1g2', comp=9).reshape(3,3,nao,nao,nao,nao) + v = v1 + v1.transpose(0,1,4,5,2,3) + assert abs(v1 - v1.transpose(0,1,3,2,4,5)).max() < 1e-8 + assert abs(v1 - v1.transpose(0,1,2,3,5,4)).max() < 1e-8 + assert abs(v2 - v2.transpose(1,0,4,5,2,3)).max() < 1e-8 + + assert abs(v - v.transpose(0,1,3,2,4,5)).max() < 1e-8 + assert abs(v - v.transpose(0,1,2,3,5,4)).max() < 1e-8 + assert abs(v - v.transpose(1,0,4,5,2,3)).max() < 1e-8 + assert abs(v - v.transpose(0,1,3,2,5,4)).max() < 1e-8 + jref = numpy.einsum('xyijkl,ji->xykl', v, dm0) + kref = numpy.einsum('xyijkl,jk->xyil', v, dm0) + + vs = jk.get_jk(mol, [dm0]*4, ['ijkl,ji->s2kl', + 'ijkl,lk->s2ij', + 'ijkl,jk->s1il', + 'ijkl,li->s1kj'], + 'int2e_gg1', 's4', 9, hermi=1) + vj = vs[0] + vs[1] + vk = vs[2] + vs[3] + vj = vj.reshape(3,3,nao,nao) + vk = vk.reshape(3,3,nao,nao) + assert abs(vj - jref).max() < 1e-8 + assert abs(vk - kref).max() < 1e-8 + assert abs(jref - jref.transpose(1,0,3,2)).max() < 1e-8 + assert abs(kref - kref.transpose(1,0,3,2)).max() < 1e-8 + + v += 2 * v2 + jref = numpy.einsum('xyijkl,ji->xykl', v, dm0) + kref = numpy.einsum('xyijkl,jk->xyil', v, dm0) + assert abs(jref - jref.transpose(1,0,3,2)).max() < 1e-8 + assert abs(kref - kref.transpose(1,0,3,2)).max() < 1e-8 + + vs = jk.get_jk(mol, [dm0]*2, ['ijkl,ji->s2kl', + 'ijkl,jk->s1il'], + 'int2e_g1g2', 'aa4', 9, hermi=0) + kk = vs[1].reshape(3,3,nao,nao) + assert abs(kk-kk.transpose(1,0,3,2)).max() < 1e-8 + + vj += vs[0].reshape(3,3,nao,nao) * 2 + vk += vs[1].reshape(3,3,nao,nao) * 2 + assert abs(vj - jref).max() < 1e-8 + assert abs(vk - kref).max() < 1e-8 + + +######################### +# DFT + mol = gto.M(atom=''' + H 0. , 0. , 0. + H 0. , .7 , .8 + H .6 , 0. , .5 + H .9 , .4 , 0. + ''') + mf = dft.RKS(mol).run(conv_tol=1e-12) + grids = mf.grids + ni = mf._numint + shls_slice = (0, mol.nbas) + ao_loc = mol.ao_loc_nr() + + dm = mf.make_rdm1() + make_rho, nset, nao = ni._gen_rho_evaluator(mol, dm) + vmat = numpy.zeros((3,nao,nao)) + ao_deriv = 0 + for ao, mask, weight, coords \ + in ni.block_loop(mol, grids, nao, ao_deriv): + rho = make_rho(0, ao, mask, 'LDA') + vxc = ni.eval_xc('LDA,', rho, 0, deriv=1)[1] + vrho = vxc[0] + aow = numpy.einsum('pi,p->pi', ao, weight*vrho) + giao = mol.eval_gto('GTOval_ig', coords, comp=3, non0tab=mask) + vmat[0] += numint._dot_ao_ao(mol, aow, giao[0], mask, shls_slice, ao_loc) + vmat[1] += numint._dot_ao_ao(mol, aow, giao[1], mask, shls_slice, ao_loc) + vmat[2] += numint._dot_ao_ao(mol, aow, giao[2], mask, shls_slice, ao_loc) + rho = vxc = vrho = aow = None + vref = vmat - vmat.transpose(0,2,1) + assert abs(lib.fp(vref) - -0.236996049404) < 1e-7 + + ao = mol.eval_gto('GTOval', grids.coords) + rho = make_rho(0, ao, grids.non0tab, 'LDA') + vxc = ni.eval_xc('LDA,', rho, 0, deriv=1)[1] + vrho = vxc[0] + + aow = numpy.einsum('pi,p->pi', ao, grids.weights*vrho) + giao = mol.eval_gto('GTOval_ig', grids.coords, comp=3) + vmat = (numint._dot_ao_ao(mol, aow, giao[0], grids.non0tab, shls_slice, ao_loc), + numint._dot_ao_ao(mol, aow, giao[1], grids.non0tab, shls_slice, ao_loc), + numint._dot_ao_ao(mol, aow, giao[2], grids.non0tab, shls_slice, ao_loc)) + vmat = numpy.array(vmat) + vmat = vmat - vmat.transpose(0,2,1) + assert abs(vref - vmat).max() < 1e-9 + + aow = numpy.einsum('pi,p,p,px->pxi', ao, grids.weights, vrho, grids.coords) + vmat = lib.einsum('pxi,pj->xij', aow, ao) + atom_coords = mol.atom_coords() + Rx = atom_coords[:,0] + Ry = atom_coords[:,1] + Rz = atom_coords[:,2] + v1 = numpy.zeros_like(vmat) + v1[0] += (Ry[:,None]-Ry) * vmat[2] + v1[0] -= (Rz[:,None]-Rz) * vmat[1] + v1[1] += (Rz[:,None]-Rz) * vmat[0] + v1[1] -= (Rx[:,None]-Rx) * vmat[2] + v1[2] += (Rx[:,None]-Rx) * vmat[1] + v1[2] -= (Ry[:,None]-Ry) * vmat[0] + v1 *= -.5 + assert abs(vref - v1).max() < 1e-9 diff --git a/pyscf/prop/nmr/test/test_rks_nmr.py b/pyscf/prop/nmr/test/test_rks_nmr.py index 4270be5..3220687 100644 --- a/pyscf/prop/nmr/test/test_rks_nmr.py +++ b/pyscf/prop/nmr/test/test_rks_nmr.py @@ -15,7 +15,8 @@ def setUpModule(): nucmod = {'F': 2}, # gaussian nuclear model basis = '6-31g', ) - mf = dft.RKS(mol).run() + with lib.temporary_env(dft.radi, ATOM_SPECIFIC_TREUTLER_GRIDS=False): + mf = dft.RKS(mol).run() def tearDownModule(): global mol, mf diff --git a/pyscf/prop/nmr/test/test_uks_nmr.py b/pyscf/prop/nmr/test/test_uks_nmr.py index 3dc8359..a0e5913 100644 --- a/pyscf/prop/nmr/test/test_uks_nmr.py +++ b/pyscf/prop/nmr/test/test_uks_nmr.py @@ -15,7 +15,8 @@ def setUpModule(): nucmod = {'F': 2}, # gaussian nuclear model basis = '6-31g', ) - mf = dft.UKS(mol).run() + with lib.temporary_env(dft.radi, ATOM_SPECIFIC_TREUTLER_GRIDS=False): + mf = dft.UKS(mol).run() def tearDownModule(): global mol, mf @@ -42,8 +43,8 @@ def test_he(self): def test_nr_giao_cpscf(self): nmr = mf.NMR() msc = nmr.kernel() - self.assertAlmostEqual(msc[1][0,0], 368.881240, 5) - self.assertAlmostEqual(msc[1][1,1], 368.881240, 5) + self.assertAlmostEqual(msc[1][0,0], 368.881232, 5) + self.assertAlmostEqual(msc[1][1,1], 368.881232, 5) self.assertAlmostEqual(msc[1][2,2], 482.413298, 5) self.assertAlmostEqual(lib.fp(msc), -131.708525548098, 5) diff --git a/pyscf/prop/polarizability/rhf.py b/pyscf/prop/polarizability/rhf.py index 36c038a..632329b 100644 --- a/pyscf/prop/polarizability/rhf.py +++ b/pyscf/prop/polarizability/rhf.py @@ -337,65 +337,3 @@ def vind(mo1): from pyscf import scf scf.hf.RHF.Polarizability = lib.class_as_method(Polarizability) - -if __name__ == '__main__': - from pyscf import gto - from pyscf import scf - mol = gto.Mole() - mol.atom = '''h , 0. 0. 0. - F , 0. 0. .917''' - mol.basis = '631g' - mol.build() - - mf = scf.RHF(mol).run(conv_tol=1e-14) - polar = mf.Polarizability().polarizability() - hpol = mf.Polarizability().hyper_polarizability() - print(polar) - - mf.verbose = 0 - charges = mol.atom_charges() - coords = mol.atom_coords() - charge_center = numpy.einsum('i,ix->x', charges, coords) / charges.sum() - with mol.with_common_orig(charge_center): - ao_dip = mol.intor_symmetric('int1e_r', comp=3) - h1 = mf.get_hcore() - def apply_E(E): - mf.get_hcore = lambda *args, **kwargs: h1 + numpy.einsum('x,xij->ij', E, ao_dip) - mf.run(conv_tol=1e-14) - return mf.dip_moment(mol, mf.make_rdm1(), unit='AU', verbose=0) - e1 = apply_E([ 0.0001, 0, 0]) - e2 = apply_E([-0.0001, 0, 0]) - print((e1 - e2) / 0.0002) - e1 = apply_E([0, 0.0001, 0]) - e2 = apply_E([0,-0.0001, 0]) - print((e1 - e2) / 0.0002) - e1 = apply_E([0, 0, 0.0001]) - e2 = apply_E([0, 0,-0.0001]) - print((e1 - e2) / 0.0002) - - print(hpol) - def apply_E(E): - mf.get_hcore = lambda *args, **kwargs: h1 + numpy.einsum('x,xij->ij', E, ao_dip) - mf.run(conv_tol=1e-14) - return Polarizability(mf).polarizability() - e1 = apply_E([ 0.0001, 0, 0]) - e2 = apply_E([-0.0001, 0, 0]) - print((e1 - e2) / 0.0002) - e1 = apply_E([0, 0.0001, 0]) - e2 = apply_E([0,-0.0001, 0]) - print((e1 - e2) / 0.0002) - e1 = apply_E([0, 0, 0.0001]) - e2 = apply_E([0, 0,-0.0001]) - print((e1 - e2) / 0.0002) - - mol = gto.M(atom='''O 0. 0. 0. - H 0. -0.757 0.587 - H 0. 0.757 0.587''', - basis='6-31g') - mf = scf.RHF(mol).run(conv_tol=1e-14) - print(Polarizability(mf).polarizability()) - print(Polarizability(mf).polarizability_with_freq(freq= 0.)) - - print(Polarizability(mf).polarizability_with_freq(freq= 0.1)) - print(Polarizability(mf).polarizability_with_freq(freq=-0.1)) - diff --git a/pyscf/prop/polarizability/test/test_rhf.py b/pyscf/prop/polarizability/test/test_rhf_polarizability.py similarity index 62% rename from pyscf/prop/polarizability/test/test_rhf.py rename to pyscf/prop/polarizability/test/test_rhf_polarizability.py index 8efc32f..c77cba7 100644 --- a/pyscf/prop/polarizability/test/test_rhf.py +++ b/pyscf/prop/polarizability/test/test_rhf_polarizability.py @@ -59,16 +59,19 @@ def test_polarizability_with_freq(self): ref = numpy.einsum('px,py->xy', v, u1)*2 val = rhf.Polarizability(mf).polarizability_with_freq(freq) # errors mainly due to krylov solver - self.assertAlmostEqual(abs(ref-val).max(), 0, 7) + self.assertAlmostEqual(abs(ref-val).max(), 0, 6) # static property ref = numpy.einsum('px,py->xy', v1, numpy.linalg.solve(a+b, v1))*4 val = rhf.Polarizability(mf).polarizability() - self.assertAlmostEqual(abs(ref-val).max(), 0, 7) + self.assertAlmostEqual(abs(ref-val).max(), 0, 6) val = rhf.Polarizability(mf).polarizability_with_freq(freq=0) # errors mainly due to krylov solver - self.assertAlmostEqual(abs(ref-val).max(), 0, 7) + self.assertAlmostEqual(abs(ref-val).max(), 0, 6) + + val = mf.Polarizability().polarizability() + self.assertAlmostEqual(abs(ref-val).max(), 0, 6) def test_polarizability_with_freq_against_tdhf(self): @@ -98,6 +101,53 @@ def test_polarizability_with_freq_against_tdhf(self): val = rhf.Polarizability(mf).polarizability_with_freq(freq= 0.1) self.assertAlmostEqual(abs(ref - val).max(), 0, 8) + def test_finite_difference(self): + mol = gto.Mole() + mol.atom = '''h , 0. 0. 0. + F , 0. 0. .917''' + mol.basis = '631g' + mol.build() + + mf = scf.RHF(mol).run(conv_tol=1e-14) + polar = mf.Polarizability().polarizability() + hpol = mf.Polarizability().hyper_polarizability() + + mf.verbose = 0 + charges = mol.atom_charges() + coords = mol.atom_coords() + charge_center = numpy.einsum('i,ix->x', charges, coords) / charges.sum() + with mol.with_common_orig(charge_center): + ao_dip = mol.intor_symmetric('int1e_r', comp=3) + h1 = mf.get_hcore() + def apply_E(E): + mf.get_hcore = lambda *args, **kwargs: h1 + numpy.einsum('x,xij->ij', E, ao_dip) + mf.run(conv_tol=1e-14) + return mf.dip_moment(mol, mf.make_rdm1(), unit='AU', verbose=0) + e1 = apply_E([ 0.0001, 0, 0]) + e2 = apply_E([-0.0001, 0, 0]) + self.assertAlmostEqual(abs(polar[0] - (e1 - e2) / 0.0002).max(), 0, 5) + e1 = apply_E([0, 0.0001, 0]) + e2 = apply_E([0,-0.0001, 0]) + self.assertAlmostEqual(abs(polar[1] - (e1 - e2) / 0.0002).max(), 0, 5) + e1 = apply_E([0, 0, 0.0001]) + e2 = apply_E([0, 0,-0.0001]) + self.assertAlmostEqual(abs(polar[2] - (e1 - e2) / 0.0002).max(), 0, 5) + + def apply_E(E): + mf.get_hcore = lambda *args, **kwargs: h1 + numpy.einsum('x,xij->ij', E, ao_dip) + mf.run(conv_tol=1e-14) + return mf.Polarizability().polarizability() + e1 = apply_E([ 0.0001, 0, 0]) + e2 = apply_E([-0.0001, 0, 0]) + self.assertAlmostEqual(abs(hpol[0] - (e1 - e2) / 0.0002).max(), 0, 4) + e1 = apply_E([0, 0.0001, 0]) + e2 = apply_E([0,-0.0001, 0]) + self.assertAlmostEqual(abs(hpol[1] - (e1 - e2) / 0.0002).max(), 0, 4) + e1 = apply_E([0, 0, 0.0001]) + e2 = apply_E([0, 0,-0.0001]) + self.assertAlmostEqual(abs(hpol[2] - (e1 - e2) / 0.0002).max(), 0, 4) + + if __name__ == "__main__": print("Tests for polarizability") unittest.main() diff --git a/pyscf/prop/rotational_gtensor/test/test_rhf.py b/pyscf/prop/rotational_gtensor/test/test_rhf_rotational_gtensor.py similarity index 100% rename from pyscf/prop/rotational_gtensor/test/test_rhf.py rename to pyscf/prop/rotational_gtensor/test/test_rhf_rotational_gtensor.py diff --git a/pyscf/prop/ssc/test/test_scf_integral.py b/pyscf/prop/ssc/test/test_scf_integral.py index 81b26d0..27947b2 100644 --- a/pyscf/prop/ssc/test/test_scf_integral.py +++ b/pyscf/prop/ssc/test/test_scf_integral.py @@ -19,6 +19,9 @@ from pyscf import gto from pyscf import dft from pyscf.prop.ssc import rhf as rhf_ssc +import pyscf +from packaging.version import Version +pyscf_28 = Version(pyscf.__version__) >= Version('2.8.0') # Test numerical integration scheme JCP 73, 5718 (1980); DOI:10.1063/1.440051 # Test field gradients integrals @@ -77,11 +80,12 @@ def test_3c_r1(self): fakemol._env = numpy.hstack((r0, a**2, w*2/numpy.pi**.5/s_cart2sph_factor)) fakemol._built = True - pmol = mol + fakemol - i3c = pmol.intor('int4c1e_sph', shls_slice=(mol.nbas,pmol.nbas, 0,mol.nbas,0,mol.nbas,0,mol.nbas)) - nao = mol.nao_nr() - i3c = i3c.reshape(nao,nao,nao) - self.assertAlmostEqual(abs(i3c - nablav).max(), 0, 9) + if pyscf_28: + pmol = mol + fakemol + i3c = pmol.intor('int4c1e_sph', shls_slice=(mol.nbas,pmol.nbas, 0,mol.nbas,0,mol.nbas,0,mol.nbas)) + nao = mol.nao_nr() + i3c = i3c.reshape(nao,nao,nao) + self.assertAlmostEqual(abs(i3c - nablav).max(), 0, 9) def test1_3c_r3(self): '''3-center vec{r}/r^3 type integral''' @@ -103,11 +107,12 @@ def test1_3c_r3(self): fakemol._env = numpy.hstack((r0, a**2, a**2*w*4/numpy.pi**.5/p_cart2sph_factor)) fakemol._built = True - pmol = mol + fakemol - i3c = pmol.intor('int4c1e_sph', shls_slice=(mol.nbas,pmol.nbas, 0,mol.nbas,0,mol.nbas,0,mol.nbas)) - nao = mol.nao_nr() - i3c = i3c.reshape(3,nao,nao,nao) - self.assertAlmostEqual(abs(i3c - nablav).max(), 0, 9) + if pyscf_28: + pmol = mol + fakemol + i3c = pmol.intor('int4c1e_sph', shls_slice=(mol.nbas,pmol.nbas, 0,mol.nbas,0,mol.nbas,0,mol.nbas)) + nao = mol.nao_nr() + i3c = i3c.reshape(3,nao,nao,nao) + self.assertAlmostEqual(abs(i3c - nablav).max(), 0, 9) def test2_3c_r3(self): @@ -150,13 +155,14 @@ def test2_3c_r3(self): fakemol._env = numpy.hstack((orig1, orig2, a**2, a**2*w*4/numpy.pi**.5/p_cart2sph_factor)) fakemol._built = True - pmol = mol + fakemol - mat = pmol.intor('int4c1e_sph', - shls_slice=(mol.nbas, pmol.nbas, mol.nbas, pmol.nbas, 0, mol.nbas, 0, mol.nbas)) - nao = mol.nao_nr() - mat = mat.reshape(6,6,nao,nao) - self.assertAlmostEqual(mat[0,3,0,3], 0.01931792982, 9) - self.assertAlmostEqual(mat[0,4,0,3], 0.01431132964, 9) + if pyscf_28: + pmol = mol + fakemol + mat = pmol.intor('int4c1e_sph', + shls_slice=(mol.nbas, pmol.nbas, mol.nbas, pmol.nbas, 0, mol.nbas, 0, mol.nbas)) + nao = mol.nao_nr() + mat = mat.reshape(6,6,nao,nao) + self.assertAlmostEqual(mat[0,3,0,3], 0.01931792982, 9) + self.assertAlmostEqual(mat[0,4,0,3], 0.01431132964, 9) # With numerical integration grids = dft.gen_grid.Grids(mol)