Pyberny prop update

src/python/structurefinder/pyberny/__init__.py

@@ -13,7 +13,12 @@
 # limitations under the License.
 
 import pluginplay as pp
-from simde import EnergyNuclearGradientStdVectorD, TotalEnergy, MoleculeFromString
+from simde import (
+    EnergyNuclearGradientStdVectorD,
+    TotalEnergyNuclearOptimization,
+    MoleculeFromString,
+    TotalEnergy,
+)
 from berny import Berny, geomlib
 import chemist
 import numpy as np
@@ -23,65 +28,52 @@ class GeomoptViaPyberny(pp.ModuleBase):
 
     def __init__(self):
         pp.ModuleBase.__init__(self)
-        self.satisfies_property_type(TotalEnergy())
+        self.satisfies_property_type(TotalEnergyNuclearOptimization())
         self.description("Performs PyBerny optimization")
         self.add_submodule(TotalEnergy(), "Energy")
         self.add_submodule(EnergyNuclearGradientStdVectorD(), "Gradient")
         self.add_submodule(MoleculeFromString(), "StringConv")
 
     def run_(self, inputs, submods):
-        pt = TotalEnergy()
-        mol, = pt.unwrap_inputs(inputs)
-        molecule = mol.molecule
+        pt = TotalEnergyNuclearOptimization()
+        sys, points = pt.unwrap_inputs(inputs)
+        molecule = sys.molecule
 
         # Convert Chemist Chemical System to XYZ
         xyz = ""
-        xyz += (str(molecule.size()) + "\n\n")
-        for i in range(molecule.size()):
-            xyz += (molecule.at(i).name + " " + str(molecule.at(i).x) + " " +
-                    str(molecule.at(i).y) + " " + str(molecule.at(i).z) + "\n")
+        xyz += str(molecule.size()) + "\n\n"
+
+        #TODO ensure points == molecule.nuclei.charges.point_set
+        for i in range(points.size()):
+            xyz += (molecule.at(i).name + " " + str(points.at(i).x) + " " +
+                    str(points.at(i).y) + " " + str(points.at(i).z) + "\n")
 
         # Loads the geometry string into the Berny optimizer
         # object.
-        optimizer = Berny(geomlib.loads(xyz, fmt='xyz'))
+        optimizer = Berny(geomlib.loads(xyz, fmt="xyz"))
 
         for geom in optimizer:
-
             # Converts the "Berny" geometry object to Chemical System
-            geom2xyz = geom.dumps('xyz')
-            print('Berny Geom to XYZ value: \n' + geom2xyz + '\n')
-            lines = geom2xyz.split('\n')
-            print('Lines of geom2xyz: \n' + str(lines) + '\n')
-            mol_string = '\n'.join(lines[2:])
-            print('Lines to string: \n' + mol_string + '\n')
+            geom2xyz = geom.dumps("xyz")
+            lines = geom2xyz.split("\n")
+            mol_string = "\n".join(lines[2:])
             xyz2chem_mol = submods["StringConv"].run_as(
                 MoleculeFromString(), mol_string)
-            print('String conversion from xyz to chem sys: \n' +
-                  str(xyz2chem_mol.nuclei) + '\n')
             geom = chemist.ChemicalSystem(xyz2chem_mol)
-            print('Chemical system of xyz2chem_mol: \n' +
-                  str(geom.molecule.nuclei) + '\n')
             geom_nuclei = geom.molecule.nuclei.as_nuclei()
-            geom_points = geom_nuclei.charges.point_set
+            geom_points = geom_nuclei.charges.point_set.as_point_set()
 
             # Main optimizer operation
             energy = submods["Energy"].run_as(TotalEnergy(), geom)
-            print('Interim energy: \n' + str(energy) + '\n')
             gradients = submods["Gradient"].run_as(
-                EnergyNuclearGradientStdVectorD(), geom,
-                geom_points.as_point_set())
-            print('Interim gradient: \n' + str(gradients) + '\n')
+                EnergyNuclearGradientStdVectorD(), geom, geom_points)
             optimizer.send((np.array(energy).item(), gradients))
 
-        opt_geom = geom.molecule.nuclei
-        print(
-            'Resulting relaxed geometry (assigned to variable opt_geom): \n' +
-            str(opt_geom))
-        # Optimized energy is of type "float"
-        e = energy
-        print(e)
+        opt_geom_nuclei = geom.molecule.nuclei.as_nuclei()
+        opt_geom_points = opt_geom_nuclei.charges.point_set.as_point_set()
+
         rv = self.results()
-        return pt.wrap_results(rv, e)
+        return pt.wrap_results(rv, energy, opt_geom_points)
 
 
 def load_pyberny_modules(mm):

tests/python/unit_tests/test_optimizers/test_pybernyop.py

@@ -18,7 +18,15 @@
 import pluginplay as pp
 import chemist
 import unittest
-from simde import TotalEnergy
+from simde import TotalEnergyNuclearOptimization
+
+
+def diatomic_bond_distance(coords):
+    val = 0
+    for i in range(int(len(coords) / 2)):
+        val += (coords[i] - coords[i + 3])**2
+    distance = np.sqrt(val)
+    return distance
 
 
 class Test_optimize_pyberny(unittest.TestCase):
@@ -27,19 +35,38 @@ def test_optimize_pyberny(self):
         mm = pp.ModuleManager()
         nwchemex.load_modules(mm)
         structurefinder.load_modules(mm)
-        mm.change_input("NWChem : SCF", "basis set", "sto-3g")
-        mm.change_input("NWChem : SCF Gradient", "basis set", "sto-3g")
-        mm.change_submod("PyBerny", "Gradient", "NWChem : SCF Gradient")
-        mm.change_submod("PyBerny", "Energy", "NWChem : SCF")
-        mm.change_submod("Pyberny", "StringConv",
-                         "ChemicalSystem via QCElemental")
-        egy = mm.run_as(TotalEnergy(), "PyBerny",
-                        chemist.ChemicalSystem(self.mol))
-        print("Energy = " + str(egy))
-        print(np.array(egy).item())
-        self.assertAlmostEqual(np.array(egy).item(), -1.117505879316, 10)
+
+        pyberny_mod = mm.at("PyBerny")
+        nwchem_scf_mod = mm.at("NWChem : SCF")
+        nwchem_grad_mod = mm.at("NWChem : SCF Gradient")
+        string_conv_mod = mm.at("ChemicalSystem via QCElemental")
+
+        nwchem_scf_mod.change_input("basis set", "sto-3g")
+        nwchem_grad_mod.change_input("basis set", "sto-3g")
+        pyberny_mod.change_submod("Energy", nwchem_scf_mod)
+        pyberny_mod.change_submod("Gradient", nwchem_grad_mod)
+        pyberny_mod.change_submod("StringConv", string_conv_mod)
+
+        energy, points = pyberny_mod.run_as(TotalEnergyNuclearOptimization(),
+                                            self.sys, self.point_set_i)
+
+        coords = []
+        for atom in range(points.size()):
+            for coord in range(3):
+                coords.append(points.at(atom).coord(coord))
+
+        distance = diatomic_bond_distance(coords)
+
+        energy = np.array(energy).item()
+        self.assertAlmostEqual(energy, self.energy, 10)
+        self.assertAlmostEqual(distance, self.distance, 8)
 
     def setUp(self):
         self.mol = chemist.Molecule()
         self.mol.push_back(chemist.Atom("H", 1, 1.0079, 0.0, 0.0, 0.0))
         self.mol.push_back(chemist.Atom("H", 1, 1.0079, 0.0, 0.0, 1.0))
+        self.sys = chemist.ChemicalSystem(self.mol)
+        self.nuclei = self.mol.nuclei.as_nuclei()
+        self.point_set_i = self.nuclei.charges.point_set.as_point_set()
+        self.distance = 1.34606231
+        self.energy = -1.117505879316