NWChem Gradient Interface

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,33 +27,123 @@ 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.
 
     Currently, the modules in this collection are:
 
     #.  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)

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):

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,14 +56,18 @@ 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
     """
 
     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

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')

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!")

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())