Skip to content

Commit

Permalink
pricise unit test for solvent models
Browse files Browse the repository at this point in the history
  • Loading branch information
wxj6000 committed Apr 22, 2024
1 parent ac82cfb commit 884c93c
Show file tree
Hide file tree
Showing 6 changed files with 52 additions and 49 deletions.
3 changes: 2 additions & 1 deletion gpu4pyscf/solvent/grad/smd.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
4 changes: 2 additions & 2 deletions gpu4pyscf/solvent/smd.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 = {
Expand Down Expand Up @@ -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
Expand Down
37 changes: 19 additions & 18 deletions gpu4pyscf/solvent/tests/test_pcm.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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')
Expand Down
53 changes: 28 additions & 25 deletions gpu4pyscf/solvent/tests/test_pcm_grad.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down
2 changes: 0 additions & 2 deletions gpu4pyscf/solvent/tests/test_smd.py
Original file line number Diff line number Diff line change
Expand Up @@ -114,15 +114,13 @@ 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):
smdobj = smd.SMD(mol)
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):
Expand Down
2 changes: 1 addition & 1 deletion gpu4pyscf/solvent/tests/test_smd_grad.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down

0 comments on commit 884c93c

Please sign in to comment.