Skip to content

Commit

Permalink
Fix bug in lpb fock routine and update tests (#2963)
Browse files Browse the repository at this point in the history 
* Fix bug in lpb fock routine and update tests

* Slightly more stable tolerances

* fix typo

* And another one
  • Loading branch information
@mfherbst @loriab
mfherbst authored and loriab committed Jul 10, 2023
1 parent 2dc6942 commit da28f08
Show file tree
Hide file tree
Showing 3 changed files with 46 additions and 34 deletions.
4 changes: 2 additions & 2 deletions external/upstream/ddx/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
if(${ENABLE_ddx})
if(NOT (${CMAKE_DISABLE_FIND_PACKAGE_ddx}))
include(FindPythonModule)
find_python_module(pyddx ATLEAST 0.4.0 QUIET)
find_python_module(pyddx ATLEAST 0.4.2 QUIET)
endif()

if(${pyddx_FOUND})
Expand All @@ -21,7 +21,7 @@ if(${ENABLE_ddx})

ExternalProject_Add(ddx_external
BUILD_ALWAYS 1
URL https://github.com/ddsolvation/ddX/archive/v0.3.0.tar.gz
URL https://github.com/ddsolvation/ddX/archive/v0.4.2.tar.gz
CONFIGURE_COMMAND ""
UPDATE_COMMAND ""
BUILD_COMMAND ${Python_EXECUTABLE} setup.py build
Expand Down
12 changes: 8 additions & 4 deletions psi4/driver/procrouting/solvent/ddx.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,7 @@ class DdxInterface:
def __init__(self, molecule, options, basisset):
# verify that the minimal version is used if pyddx is provided
# from outside the Psi4 ecosystem
min_version = "0.4.0"
min_version = "0.4.2"
if parse_version(pyddx.__version__) < parse_version(min_version):
raise ModuleNotFoundError("pyddx version {} is required at least. "
"Version {}"
Expand Down Expand Up @@ -229,13 +229,13 @@ def get_solvation_contributions(self, density_matrix, state=None,
# Compute electronic contributions
psi = self.numints.dd_density_integral(self.scaled_ylms, density_matrix).np.T
dummy_charges = core.Vector.from_array(np.ones(self.model.n_cav))
coords = core.Matrix.from_array(self.model.cavity.T)
phi = self.mints.electrostatic_potential_value(dummy_charges, coords, density_matrix).np
cavcoords = core.Matrix.from_array(self.model.cavity.T)
phi = self.mints.electrostatic_potential_value(dummy_charges, cavcoords, density_matrix).np

elec_field = None
derivative_order = self.model.required_phi_derivative_order(compute_forces=False)
if derivative_order > 0:
elec_field = self.mints.electric_field_value(coords, density_matrix).np.T
elec_field = self.mints.electric_field_value(cavcoords, density_matrix).np.T

# elec_only = True
if not elec_only:
Expand Down Expand Up @@ -279,6 +279,10 @@ def get_solvation_contributions(self, density_matrix, state=None,
for (block, ylm) in zip(self.dftgrid.blocks(), self.scaled_ylms)]
V_ddx = self.numints.potential_integral(eta)

if model.model in ("lpb", ):
V_ind = self.mints.induction_operator(cavcoords, core.Matrix.from_array(state.zeta_dip.T))
V_ddx.add(V_ind)

extern = core.ExternalPotential()
for xi_nj, pos_nj in zip(state.xi, model.cavity.T):
extern.addCharge(xi_nj, pos_nj[0], pos_nj[1], pos_nj[2])
Expand Down
64 changes: 36 additions & 28 deletions tests/pytests/test_ddx.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@
}


def _base_test_fock(fock_term, density_matrix, eps=1e-4, tol=1e-6):
def _base_test_fock(fock_term, density_matrix, eps=1e-4, tol=1e-8):
E, V = fock_term(density_matrix)

perturbation = np.random.random(density_matrix.np.shape)
Expand All @@ -96,21 +96,23 @@ def _base_test_fock(fock_term, density_matrix, eps=1e-4, tol=1e-6):
"dm": core.Matrix.from_array(np.array([[0.]])),
"ddx": {"model": "cosmo", "solvent_epsilon": 1e8, "eta": 0, "radii": [1.0]},
"ref": -0.2645886054599999, # from Gaussian
"tol": 1e-6,
}, id='h'),
pytest.param({
"geom": __geoms["h2"],
"dm": core.Matrix.from_array(0.6682326961201372 * np.ones((2, 2))),
"ddx": {"model": "cosmo", "solvent_epsilon": 1e8, "eta": 0, "radii": [1.5873, 1.5873]},
"ref": -0.0002016948, # from Gaussian
"tol": 1e-6,
}, id='h2'),
}, id='h2cosmo'),
pytest.param({
"geom": __geoms["h2"],
"dm": core.Matrix.from_array(0.6682326961201372 * np.ones((2, 2))),
"ddx": {"model": "pcm", "solvent_epsilon": 80, "radii": [1.5873, 1.5873]},
}, id='h2pcm'),
pytest.param({
"geom": __geoms["h2"],
"dm": core.Matrix.from_array(0.6682326961201372 * np.ones((2, 2))),
"ddx": {"model": "lpb", "solvent_epsilon": 80, "solvent_kappa": 1.5,
"radii": [1.5873, 1.5873]},
"tol": 3e-5,
}, id='h2lpb'),
])
def test_ddx_fock_build(inp):
Expand Down Expand Up @@ -141,7 +143,7 @@ def get_EV(density_matrix):
E, V, _ = ddx_iface.get_solvation_contributions(density_matrix)
return E, V

_base_test_fock(get_EV, inp["dm"], tol=inp["tol"])
_base_test_fock(get_EV, inp["dm"], tol=1e-8)

if "ref" in inp:
E, _ = get_EV(inp["dm"])
Expand Down Expand Up @@ -189,24 +191,25 @@ def ddenergy(model, solvent_kappa=0.0, solvent_epsilon=80.0):
e_cosmo, v_cosmo = ddenergy("cosmo", solvent_epsilon=1e12)

e_lpb_einf, v_lpb_einf = ddenergy("lpb", solvent_kappa=0.1, solvent_epsilon=1e12)
e_lpb0, v_lpb0 = ddenergy("lpb", solvent_kappa=1e-3, solvent_epsilon=80.0)
e_lpbinf, v_lpbinf = ddenergy("lpb", solvent_kappa=8, solvent_epsilon=80.0)

assert abs(e_pcminf - e_cosmo) / abs(e_cosmo) < 1e-5
assert abs(e_pcminf - e_lpb_einf) / abs(e_lpb_einf) < 1e-4
assert abs(e_cosmo - e_lpb_einf) / abs(e_lpb_einf) < 1e-4
assert abs(e_lpbinf - e_lpb_einf) / abs(e_lpb_einf) < 1e-4
assert abs(e_lpbinf - e_cosmo) / abs(e_cosmo) < 1e-4
assert abs(e_lpbinf - e_pcminf) / abs(e_pcminf) < 1e-4
assert abs(e_lpb0 - e_pcm) / abs(e_pcm) < 1e-3

assert np.max(np.abs(v_pcminf.np - v_cosmo.np)) < 1e-5
assert np.max(np.abs(v_pcminf.np - v_lpb_einf.np)) < 1e-4
assert np.max(np.abs(v_cosmo.np - v_lpb_einf.np)) < 1e-4
assert np.max(np.abs(v_lpbinf.np - v_lpb_einf.np)) < 1e-4
assert np.max(np.abs(v_lpbinf.np - v_cosmo.np)) < 1e-4
assert np.max(np.abs(v_lpbinf.np - v_pcminf.np)) < 1e-4
assert np.max(np.abs(v_lpb0.np - v_pcm.np)) < 1e-3
e_lpb0, v_lpb0 = ddenergy("lpb", solvent_kappa=1e-5, solvent_epsilon=80.0)
e_lpbinf, v_lpbinf = ddenergy("lpb", solvent_kappa=10, solvent_epsilon=80.0)

assert abs(e_pcminf - e_cosmo) / abs(e_cosmo) < 1e-10 # noqa: E221
assert abs(e_pcminf - e_lpb_einf) / abs(e_lpb_einf) < 1e-10 # noqa: E221
assert abs(e_cosmo - e_lpb_einf) / abs(e_lpb_einf) < 1e-10 # noqa: E221
assert abs(e_lpb0 - e_pcm) / abs(e_pcm) < 1e-5 # noqa: E221
assert abs(e_lpbinf - e_lpb_einf) / abs(e_lpb_einf) < 1e-4 # noqa: E221
assert abs(e_lpbinf - e_cosmo) / abs(e_cosmo) < 1e-4 # noqa: E221
assert abs(e_lpbinf - e_pcminf) / abs(e_pcminf) < 1e-4 # noqa: E221

assert np.max(np.abs(v_pcminf.np - v_cosmo.np)) < 1e-10 # noqa: E221
assert np.max(np.abs(v_pcminf.np - v_lpb_einf.np)) < 1e-10 # noqa: E221
assert np.max(np.abs(v_cosmo.np - v_lpb_einf.np)) < 1e-10 # noqa: E221
assert np.max(np.abs(v_lpb0.np - v_pcm.np)) < 1e-5 # noqa: E221
assert np.max(np.abs(v_lpbinf.np - v_lpb_einf.np)) < 1e-4 # noqa: E221
assert np.max(np.abs(v_lpbinf.np - v_cosmo.np)) < 1e-4 # noqa: E221
assert np.max(np.abs(v_lpbinf.np - v_pcminf.np)) < 1e-4 # noqa: E221


@pytest.mark.quick
@uusing("ddx")
Expand Down Expand Up @@ -284,13 +287,13 @@ def test_ddx_rhf_reference(inp):
"geom": __geoms["fcm"],
"basis": "cc-pvdz",
"ddx": {"model": "pcm", "solvent": "water", "radii_set": "uff", },
"ref": -597.9718942424215,
"ref": -597.9718943192062,
}, id='fcm-pcm'),
pytest.param({
"geom": __geoms["nh3"],
"basis": "cc-pvdz",
"ddx": {"model": "lpb", "solvent": "water", "radii_set": "uff", "solvent_kappa": 0.11},
"ref": -56.1988043665621,
"ref": -56.1988043810054,
}, id='nh3-lpb'),
])
def test_ddx_rhf_consistency(inp):
Expand Down Expand Up @@ -322,6 +325,11 @@ def test_ddx_eri_algorithms(scf_type):
"ddx_solvent": "water",
"ddx_radii_set": "uff",
})
ref = -56.1715394
if scf_type in ("df", "cd"):
tol = 4
else:
tol = 9

ref = -56.171546236617495
scf_e = psi4.energy('SCF')
assert compare_values(ref, scf_e, 4, "Total SCF energy with DDX versus reference data")
assert compare_values(ref, scf_e, tol, "Total SCF energy with DDX versus reference data")

0 comments on commit da28f08

Please sign in to comment.