From 884c93c1c251671a978c160b1e7f1bed9db1e062 Mon Sep 17 00:00:00 2001 From: "xiaojie.wu" Date: Tue, 23 Apr 2024 03:14:42 +0800 Subject: [PATCH] pricise unit test for solvent models --- gpu4pyscf/solvent/grad/smd.py | 3 +- gpu4pyscf/solvent/smd.py | 4 +- gpu4pyscf/solvent/tests/test_pcm.py | 37 +++++++++-------- gpu4pyscf/solvent/tests/test_pcm_grad.py | 53 +++++++++++++----------- gpu4pyscf/solvent/tests/test_smd.py | 2 - gpu4pyscf/solvent/tests/test_smd_grad.py | 2 +- 6 files changed, 52 insertions(+), 49 deletions(-) diff --git a/gpu4pyscf/solvent/grad/smd.py b/gpu4pyscf/solvent/grad/smd.py index da164fc2..aafd7f92 100644 --- a/gpu4pyscf/solvent/grad/smd.py +++ b/gpu4pyscf/solvent/grad/smd.py @@ -264,7 +264,8 @@ def kernel(self, *args, dm=None, atmlst=None, **kwargs): self.de_solvent = pcm_grad.grad_qv(self.base.with_solvent, dm) self.de_solvent+= pcm_grad.grad_solver(self.base.with_solvent, dm) self.de_solvent+= pcm_grad.grad_nuc(self.base.with_solvent, dm) - self.de_cds = get_cds(self.base.with_solvent) + #self.de_cds = get_cds(self.base.with_solvent) + self.de_cds = smd.get_cds_legacy(self.base.with_solvent)[1] self.de = self.de_solute + self.de_solvent + self.de_cds if self.verbose >= logger.NOTE: diff --git a/gpu4pyscf/solvent/smd.py b/gpu4pyscf/solvent/smd.py index 9ddec2e4..8518cd8c 100644 --- a/gpu4pyscf/solvent/smd.py +++ b/gpu4pyscf/solvent/smd.py @@ -539,7 +539,7 @@ def get_cds_legacy(smdobj): sola, solb, solc, solg, solh, soln, icds, ctypes.byref(gcds), ctypes.byref(areacds), dcds) - return gcds.value / 627.509451 + return gcds.value / 627.509451, dcds class SMD(pcm.PCM): _keys = { @@ -606,7 +606,7 @@ def dump_flags(self, verbose=None): return self def get_cds(self): - return get_cds_legacy(self) + return get_cds_legacy(self)[0] def nuc_grad_method(self, grad_method): from gpu4pyscf.solvent.grad import smd as smd_grad diff --git a/gpu4pyscf/solvent/tests/test_pcm.py b/gpu4pyscf/solvent/tests/test_pcm.py index 02fffe64..db9d506c 100644 --- a/gpu4pyscf/solvent/tests/test_pcm.py +++ b/gpu4pyscf/solvent/tests/test_pcm.py @@ -35,6 +35,7 @@ def setUpModule(): mol.basis = 'sto3g' mol.output = '/dev/null' mol.build(verbose=0) + mol.nelectron = mol.nao * 2 epsilon = 35.9 lebedev_order = 3 @@ -56,48 +57,48 @@ def _energy_with_solvent(mf, method): class KnownValues(unittest.TestCase): def test_CPCM(self): e_tot = _energy_with_solvent(scf.RHF(mol), 'C-PCM') - print(f"Energy error in RHF with C-PCM: {numpy.abs(e_tot - -74.9690902442)}") - assert numpy.abs(e_tot - -74.9690902442) < 1e-9 + print(f"Energy error in RHF with C-PCM: {numpy.abs(e_tot - -71.19244927767662)}") + assert numpy.abs(e_tot - -71.19244927767662) < 1e-5 def test_COSMO(self): e_tot = _energy_with_solvent(scf.RHF(mol), 'COSMO') - print(f"Energy error in RHF with COSMO: {numpy.abs(e_tot - -74.96900351922464)}") - assert numpy.abs(e_tot - -74.96900351922464) < 1e-9 + print(f"Energy error in RHF with COSMO: {numpy.abs(e_tot - -71.16259314943571)}") + assert numpy.abs(e_tot - -71.16259314943571) < 1e-5 def test_IEFPCM(self): e_tot = _energy_with_solvent(scf.RHF(mol), 'IEF-PCM') - print(f"Energy error in RHF with IEF-PCM: {numpy.abs(e_tot - -74.9690111344)}") - assert numpy.abs(e_tot - -74.9690111344) < 1e-9 + print(f"Energy error in RHF with IEF-PCM: {numpy.abs(e_tot - -71.19244457024647)}") + assert numpy.abs(e_tot - -71.19244457024647) < 1e-5 def test_SSVPE(self): e_tot = _energy_with_solvent(scf.RHF(mol), 'SS(V)PE') - print(f"Energy error in RHF with SS(V)PE: {numpy.abs(e_tot - -74.9689577454)}") - assert numpy.abs(e_tot - -74.9689577454) < 1e-9 + print(f"Energy error in RHF with SS(V)PE: {numpy.abs(e_tot - -71.13576912425178)}") + assert numpy.abs(e_tot - -71.13576912425178) < 1e-5 def test_uhf(self): e_tot = _energy_with_solvent(scf.UHF(mol), 'IEF-PCM') - print(f"Energy error in UHF with IEF-PCM: {numpy.abs(e_tot - -74.96901113434953)}") - assert numpy.abs(e_tot - -74.96901113434953) < 1e-9 + print(f"Energy error in UHF with IEF-PCM: {numpy.abs(e_tot - -71.19244457024645)}") + assert numpy.abs(e_tot - -71.19244457024645) < 1e-5 def test_rks(self): e_tot = _energy_with_solvent(dft.RKS(mol, xc='b3lyp'), 'IEF-PCM') - print(f"Energy error in RKS with IEF-PCM: {numpy.abs(e_tot - -75.3182692148)}") - assert numpy.abs(e_tot - -75.3182692148) < 1e-6 + print(f"Energy error in RKS with IEF-PCM: {numpy.abs(e_tot - -71.67007402042326)}") + assert numpy.abs(e_tot - -71.67007402042326) < 1e-5 def test_uks(self): e_tot = _energy_with_solvent(dft.UKS(mol, xc='b3lyp'), 'IEF-PCM') - print(f"Energy error in UKS with IEF-PCM: {numpy.abs(e_tot - -75.3182692148)}") - assert numpy.abs(e_tot - -75.3182692148) < 1e-6 + print(f"Energy error in UKS with IEF-PCM: {numpy.abs(e_tot - -71.67007402042326)}") + assert numpy.abs(e_tot - -71.67007402042326) < 1e-5 def test_dfrks(self): e_tot = _energy_with_solvent(dft.RKS(mol, xc='b3lyp').density_fit(), 'IEF-PCM') - print(f"Energy error in DFRKS with IEF-PCM: {numpy.abs(e_tot - -75.31863727142068)}") - assert numpy.abs(e_tot - -75.31863727142068) < 1e-9 + print(f"Energy error in DFRKS with IEF-PCM: {numpy.abs(e_tot - -71.67135250643568)}") + assert numpy.abs(e_tot - -71.67135250643568) < 1e-5 def test_dfuks(self): e_tot = _energy_with_solvent(dft.UKS(mol, xc='b3lyp').density_fit(), 'IEF-PCM') - print(f"Energy error in DFUKS with IEF-PCM: {numpy.abs(e_tot - -75.31863727142068)}") - assert numpy.abs(e_tot - -75.31863727142068) < 1e-9 + print(f"Energy error in DFUKS with IEF-PCM: {numpy.abs(e_tot - -71.67135250643567)}") + assert numpy.abs(e_tot - -71.67135250643567) < 1e-5 def test_to_cpu(self): mf = dft.RKS(mol, xc='b3lyp') diff --git a/gpu4pyscf/solvent/tests/test_pcm_grad.py b/gpu4pyscf/solvent/tests/test_pcm_grad.py index 919e06c6..f7f7cd8f 100644 --- a/gpu4pyscf/solvent/tests/test_pcm_grad.py +++ b/gpu4pyscf/solvent/tests/test_pcm_grad.py @@ -36,6 +36,7 @@ def setUpModule(): mol.basis = 'sto3g' mol.output = '/dev/null' mol.build(verbose=0) + mol.nelectron = mol.nao * 2 epsilon = 35.9 lebedev_order = 3 @@ -126,48 +127,50 @@ def get_FADS(mol): def test_grad_CPCM(self): grad = _grad_with_solvent('C-PCM') - g0 = numpy.asarray([ - [0.49773047433563E-15, -0.12128126037559E-15, -0.58936988992306E-01], - [0.22810111996954E-01, -0.68951901317025E-17, 0.29468494708267E-01], - [-0.22810111996957E-01, 0.12949813945902E-15, 0.29468494708266E-01]]) - + g0 = numpy.asarray( + [[ 4.65578319e-15, 5.62862593e-17, -1.61722589e+00], + [ 1.07512481e+00, 5.66523976e-17, 8.08612943e-01], + [-1.07512481e+00, -7.81228374e-17, 8.08612943e-01]] + ) print(f"Gradient error in RHF with CPCM: {numpy.linalg.norm(g0 - grad)}") - assert numpy.linalg.norm(g0 - grad) < 1e-9 + assert numpy.linalg.norm(g0 - grad) < 1e-6 def test_grad_COSMO(self): grad = _grad_with_solvent('COSMO') g0 = numpy.asarray( - [[-1.33560836e-16, 8.70874355e-17, -5.89638726e-02], - [ 2.28202396e-02, 2.63784344e-17, 2.94819363e-02], - [-2.28202396e-02, -1.08799896e-16, 2.94819363e-02]]) - + [[-8.53959617e-16, -4.87015595e-16, -1.61739114e+00], + [ 1.07538942e+00, 7.78180254e-16, 8.08695569e-01], + [-1.07538942e+00, -1.70254021e-16, 8.08695569e-01]]) print(f"Gradient error in RHF with COSMO: {numpy.linalg.norm(g0 - grad)}") - assert numpy.linalg.norm(g0 - grad) < 1e-9 + assert numpy.linalg.norm(g0 - grad) < 1e-6 def test_grad_IEFPCM(self): grad = _grad_with_solvent('IEF-PCM') - g0 = numpy.asarray([ - [0.18357915015649E-14, 0.14192681822347E-15, -0.58988087999658E-01], - [0.22822709179063E-01, -0.10002010417168E-15, 0.29494044211805E-01], - [-0.22822709179066E-01, -0.31051364515588E-16, 0.29494044211806E-01]]) + g0 = numpy.asarray( + [[-4.41438069e-15, 2.20049192e-16, -1.61732554e+00], + [ 1.07584098e+00, -5.28912700e-16, 8.08662770e-01], + [-1.07584098e+00, 2.81699314e-16, 8.08662770e-01]] + ) print(f"Gradient error in RHF with IEFPCM: {numpy.linalg.norm(g0 - grad)}") - assert numpy.linalg.norm(g0 - grad) < 1e-9 + assert numpy.linalg.norm(g0 - grad) < 1e-6 def test_grad_SSVPE(self): grad = _grad_with_solvent('SS(V)PE') - g0 = numpy.asarray([ - [0.76104817971710E-15, 0.11185701540547E-15, -0.58909172879217E-01], - [0.22862990009767E-01, -0.13861633974903E-15, 0.29454586651678E-01], - [-0.22862990009769E-01, 0.34988765678591E-16, 0.29454586651679E-01]]) + g0 = numpy.asarray( + [[ 3.42479745e-15, -1.00280742e-16, -1.61117735e+00], + [ 1.07135985e+00, -6.97375148e-16, 8.05588676e-01], + [-1.07135985e+00, 7.91425487e-16, 8.05588676e-01]] + ) print(f"Gradient error in RHF with SS(V)PE: {numpy.linalg.norm(g0 - grad)}") - assert numpy.linalg.norm(g0 - grad) < 1e-9 + assert numpy.linalg.norm(g0 - grad) < 1e-6 def test_uhf_grad_IEFPCM(self): grad = _grad_with_solvent('IEF-PCM', unrestricted=True) - g0 = numpy.asarray([ - [0.18357915015649E-14, 0.14192681822347E-15, -0.58988087999658E-01], - [0.22822709179063E-01, -0.10002010417168E-15, 0.29494044211805E-01], - [-0.22822709179066E-01, -0.31051364515588E-16, 0.29494044211806E-01]]) + g0 = numpy.asarray( + [[-5.46822686e-16, -3.41150050e-17, -1.61732554e+00], + [ 1.07584098e+00, -1.52839767e-16, 8.08662770e-01], + [-1.07584098e+00, 1.51295204e-16, 8.08662770e-01]] + ) print(f"Gradient error in UHF with IEFPCM: {numpy.linalg.norm(g0 - grad)}") assert numpy.linalg.norm(g0 - grad) < 1e-6 diff --git a/gpu4pyscf/solvent/tests/test_smd.py b/gpu4pyscf/solvent/tests/test_smd.py index a75caa46..68c65bd4 100644 --- a/gpu4pyscf/solvent/tests/test_smd.py +++ b/gpu4pyscf/solvent/tests/test_smd.py @@ -114,7 +114,6 @@ def test_cds_solvent(self): smdobj.sasa_ng = 590 smdobj.solvent = 'toluene' e_cds = smdobj.get_cds() - print(e_cds) assert numpy.abs(e_cds - -0.0013479524949097355) < 1e-8 def test_cds_water(self): @@ -122,7 +121,6 @@ def test_cds_water(self): smdobj.sasa_ng = 590 smdobj.solvent = 'water' e_cds = smdobj.get_cds() - print(e_cds) assert numpy.abs(e_cds - 0.002298448590009083) < 1e-8 def test_smd_solvent(self): diff --git a/gpu4pyscf/solvent/tests/test_smd_grad.py b/gpu4pyscf/solvent/tests/test_smd_grad.py index 58729434..9ee40636 100644 --- a/gpu4pyscf/solvent/tests/test_smd_grad.py +++ b/gpu4pyscf/solvent/tests/test_smd_grad.py @@ -76,7 +76,7 @@ def _check_grad(atom, solvent='water'): smdobj = smd.SMD(mol) smdobj.solvent = solvent - grad_cds = smd_grad.get_cds(smdobj) + grad_cds = smd.get_cds_legacy(smdobj)[1] mol.stdout.close() assert numpy.linalg.norm(fd_cds - grad_cds) < 1e-8