diff --git a/src/python/friendzone/nwx2nwchem/__init__.py b/src/python/friendzone/nwx2nwchem/__init__.py index 609302a..59b33b2 100644 --- a/src/python/friendzone/nwx2nwchem/__init__.py +++ b/src/python/friendzone/nwx2nwchem/__init__.py @@ -15,6 +15,7 @@ from ..friends import is_friend_enabled import pluginplay as pp from simde import TotalEnergy +from simde import EnergyNuclearGradientStdVectorD from ..nwx2qcengine.call_qcengine import call_qcengine @@ -26,18 +27,94 @@ def __init__(self): self.description("Calls NWChem via MolSSI's QCEngine") self.add_input('method') self.add_input("basis set") + self.add_input("keywords").set_default({}) + self.add_input("MPI config").set_default(None) # Kazuumi addition def run_(self, inputs, submods): pt = TotalEnergy() mol, = pt.unwrap_inputs(inputs) method = inputs['method'].value() basis = inputs['basis set'].value() + keywords = inputs['keywords'].value() + MPIconfig = inputs['MPI config'].value() # Kazuumi addition - e = call_qcengine(pt, mol, 'nwchem', method=method, basis=basis) + model = {"method": method, "basis": basis} + e = call_qcengine(pt, + mol, + 'nwchem', + MPIconfig, + model=model, + keywords=keywords) # Kazuumi addition rv = self.results() return pt.wrap_results(rv, e) +class NWChemGradientViaMolSSI(pp.ModuleBase): + + def __init__(self): + pp.ModuleBase.__init__(self) + self.satisfies_property_type(EnergyNuclearGradientStdVectorD()) + self.description("Calls NWChem via MolSSI's QCEngine") + self.add_input('method') + self.add_input("basis set") + self.add_input("keywords").set_default({}) + self.add_input("MPI config").set_default(None) # Kazuumi addition + + def run_(self, inputs, submods): + pt = EnergyNuclearGradientStdVectorD() + # mol, = pt.unwrap_inputs(inputs) # old version + mol, pointset1 = pt.unwrap_inputs(inputs) # new version + method = inputs['method'].value() + basis = inputs['basis set'].value() + keywords = inputs['keywords'].value() + MPIconfig = inputs['MPI config'].value() # Kazuumi addition + + model = {"method": method, "basis": basis} + e, f = call_qcengine(pt, + mol, + 'nwchem', + MPIconfig, + model=model, + keywords=keywords) # Kazuumi addition + f = [c for cs in f for c in cs] # Flatten out the list of lists + rv = self.results() + return pt.wrap_results(rv, f) + + +class NWChemEnergyAndGradientViaMolSSI(pp.ModuleBase): + + def __init__(self): + pp.ModuleBase.__init__(self) + self.satisfies_property_type(EnergyNuclearGradientStdVectorD()) + self.description("Calls NWChem via MolSSI's QCEngine") + self.add_input('method') + self.add_input("basis set") + self.add_input("keywords").set_default({}) + self.add_input("MPI config").set_default(None) # Kazuumi addition + + def run_(self, inputs, submods): + pt = EnergyNuclearGradientStdVectorD() + # mol, = pt.unwrap_inputs(inputs) # old version + mol, pointset1 = pt.unwrap_inputs(inputs) # new version + method = inputs['method'].value() + basis = inputs['basis set'].value() + keywords = inputs['keywords'].value() + MPIconfig = inputs['MPI config'].value() # Kazuumi addition + + model = {"method": method, "basis": basis} + e, f = call_qcengine(pt, + mol, + 'nwchem', + MPIconfig, + model=model, + keywords=keywords) # Kazuumi addition + combined_ef = [c for cs in f + for c in cs] # Flatten out the list of lists + combined_ef.append(e) + rv = self.results() + return pt.wrap_results(rv, combined_ef) + + def load_nwchem_modules(mm): """Loads the collection of modules that wrap NWChem calls. @@ -45,14 +122,28 @@ def load_nwchem_modules(mm): #. NWChem : SCF #. NWChem : MP2 + #. NWChem : B3LYP #. NWChem : CCSD #. NWChem : CCSD(T) - + #. NWChem : SCF Gradient + #. NWChem : MP2 Gradient + #. NWChem : B3LYP Gradient + :param mm: The ModuleManager that the NWChem Modules will be loaded into. :type mm: pluginplay.ModuleManager """ if is_friend_enabled('nwchem'): - for method in ['SCF', 'MP2', 'CCSD', 'CCSD(T)']: + for method in ['SCF', 'MP2', 'B3LYP', 'CCSD', 'CCSD(T)']: mod_key = 'NWChem : ' + method mm.add_module(mod_key, NWChemViaMolSSI()) mm.change_input(mod_key, 'method', method) + + for method in ['SCF', 'MP2', 'B3LYP']: + mod_key = 'NWChem : ' + method + ' Gradient' + mm.add_module(mod_key, NWChemGradientViaMolSSI()) + mm.change_input(mod_key, 'method', method) + + for method in ['SCF', 'MP2', 'B3LYP']: + mod_key = 'NWChem : ' + method + ' EnergyAndGradient' + mm.add_module(mod_key, NWChemEnergyAndGradientViaMolSSI()) + mm.change_input(mod_key, 'method', method) diff --git a/src/python/friendzone/nwx2qcelemental/chemical_system_conversions.py b/src/python/friendzone/nwx2qcelemental/chemical_system_conversions.py index 4e790aa..81bdbe2 100644 --- a/src/python/friendzone/nwx2qcelemental/chemical_system_conversions.py +++ b/src/python/friendzone/nwx2qcelemental/chemical_system_conversions.py @@ -43,7 +43,11 @@ def chemical_system2qc_mol(chem_sys): y = str(atom_i.y * au2ang) z = str(atom_i.z * au2ang) out += symbol + " " + x + " " + y + " " + z + "\n" - return qcel.models.Molecule.from_data(out) + # Fixing CoM and orientation and turning off symmetry to prevent molecular translation + return qcel.models.Molecule.from_data(out, + fix_com=True, + fix_orientation=True, + fix_symmetry="C1") def qc_mol2molecule(qc_mol): diff --git a/src/python/friendzone/nwx2qcengine/call_qcengine.py b/src/python/friendzone/nwx2qcengine/call_qcengine.py index 1947c23..bc150d4 100644 --- a/src/python/friendzone/nwx2qcengine/call_qcengine.py +++ b/src/python/friendzone/nwx2qcengine/call_qcengine.py @@ -18,7 +18,7 @@ from .pt2driver import pt2driver -def call_qcengine(pt, mol, program, **kwargs): +def call_qcengine(pt, mol, program, MPIconfig, **kwargs): """ Wraps calling a program through the QCEngine API. .. note:: @@ -56,7 +56,7 @@ def call_qcengine(pt, mol, program, **kwargs): backend? :type program: str :param kwargs: Key-value pairs which will be forwarded to QCElemental's - ``AtomicInput`` class via the ``model`` key. + ``AtomicInput`` class constructor. :return: The requested property. :rtype: Varies depending on the requested property @@ -64,6 +64,10 @@ def call_qcengine(pt, mol, program, **kwargs): driver = pt2driver(pt) qc_mol = chemical_system2qc_mol(mol) - inp = qcel.models.AtomicInput(molecule=qc_mol, driver=driver, model=kwargs) - results = qcng.compute(inp, program) - return results.return_result + inp = qcel.models.AtomicInput(molecule=qc_mol, driver=driver, **kwargs) + results = qcng.compute(inp, program, task_config=MPIconfig) + if (driver == "gradient" and "qcvars" in results.extras): + return float( + results.extras["qcvars"]["CURRENT ENERGY"]), results.return_result + else: + return results.return_result diff --git a/src/python/friendzone/nwx2qcengine/pt2driver.py b/src/python/friendzone/nwx2qcengine/pt2driver.py index 44215e3..f4b5943 100644 --- a/src/python/friendzone/nwx2qcengine/pt2driver.py +++ b/src/python/friendzone/nwx2qcengine/pt2driver.py @@ -29,7 +29,18 @@ def pt2driver(pt): :raises: Exception if ``pt`` is not a property type which has been registered with this function. """ - if pt.type() == simde.TotalEnergy().type(): - return 'energy' + # Set of Property Types that map to energy driver + energy_pts = { + simde.TotalEnergy().type(), + } + # Set of Property Types that map to gradient driver + gradient_pts = { + simde.EnergyNuclearGradientStdVectorD().type(), + } + + if pt.type() in energy_pts: + return 'energy' + if pt.type() in gradient_pts: + return 'gradient' raise Exception('PropertyType is not registered') diff --git a/tests/python/unit_tests/nwx2nwchem/test_nwchem.py b/tests/python/unit_tests/nwx2nwchem/test_nwchem.py index c8e74bf..ee3b5b8 100644 --- a/tests/python/unit_tests/nwx2nwchem/test_nwchem.py +++ b/tests/python/unit_tests/nwx2nwchem/test_nwchem.py @@ -15,6 +15,8 @@ from pluginplay import ModuleManager from friendzone import friends, load_modules from simde import TotalEnergy +from simde import EnergyNuclearGradientStdVectorD +from chemist import PointSetD from molecules import make_h2 import unittest @@ -49,6 +51,16 @@ def test_ccsd_t(self): egy = self.mm.run_as(TotalEnergy(), key, mol) self.assertAlmostEqual(egy, -1.122251361965036, places=4) + def test_scf_gradient(self): + mol = make_h2() + key = 'NWChem : SCF Gradient' + self.mm.change_input(key, 'basis set', 'sto-3g') + grad = self.mm.run_as(EnergyNuclearGradientStdVectorD(), key, mol, + PointSetD()) + corr = [0.0, 0.0, -0.11827177600466043, 0.0, 0.0, 0.11827177600466043] + for g, c in zip(grad, corr): + self.assertAlmostEqual(g, c, places=4) + def setUp(self): if not friends.is_friend_enabled('nwchem'): self.skipTest("NWChem backend is not enabled!") diff --git a/tests/python/unit_tests/nwx2qcengine/test_pt2driver.py b/tests/python/unit_tests/nwx2qcengine/test_pt2driver.py index 8398004..1bcd073 100644 --- a/tests/python/unit_tests/nwx2qcengine/test_pt2driver.py +++ b/tests/python/unit_tests/nwx2qcengine/test_pt2driver.py @@ -13,7 +13,7 @@ # limitations under the License. from friendzone.nwx2qcengine.pt2driver import pt2driver -from simde import TotalEnergy +import simde import unittest @@ -24,8 +24,18 @@ class NotAPT: class Testpt2driver(unittest.TestCase): def test_pts_that_map_to_energy(self): - for pt in [TotalEnergy()]: + energy_pts = { + simde.TotalEnergy(), + } + for pt in energy_pts: self.assertEqual(pt2driver(pt), 'energy') + def test_pts_that_map_to_gradient(self): + gradient_pts = { + simde.EnergyNuclearGradientStdVectorD(), + } + for pt in gradient_pts: + self.assertEqual(pt2driver(pt), 'gradient') + def test_bad_pt(self): self.assertRaises(Exception, pt2driver, NotAPT())