diff --git a/src/python/structurefinder/pyberny/__init__.py b/src/python/structurefinder/pyberny/__init__.py index f3c005c..8cf0ac2 100644 --- a/src/python/structurefinder/pyberny/__init__.py +++ b/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): diff --git a/tests/python/unit_tests/test_optimizers/test_pybernyop.py b/tests/python/unit_tests/test_optimizers/test_pybernyop.py index f7a948b..22da6c0 100644 --- a/tests/python/unit_tests/test_optimizers/test_pybernyop.py +++ b/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