From 44a76c7666c266b1760324990a75aec5c23308b5 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 14 Oct 2023 16:19:13 -0400 Subject: [PATCH] More work Signed-off-by: Geoff Hutchison --- avogadro/calc/energymanager.cpp | 6 +- avogadro/qtplugins/forcefield/forcefield.cpp | 85 +++++++------- avogadro/qtplugins/forcefield/forcefield.h | 21 ++-- avogadro/qtplugins/forcefield/scripts/gfn1.py | 102 +++++++++++++++++ avogadro/qtplugins/forcefield/scripts/gfn2.py | 102 +++++++++++++++++ .../qtplugins/forcefield/scripts/gfnff.py | 104 ++++++++++++++++++ .../qtplugins/forcefield/scripts/mmff94.py | 85 ++++++++++++++ avogadro/qtplugins/forcefield/scripts/uff.py | 85 ++++++++++++++ 8 files changed, 539 insertions(+), 51 deletions(-) create mode 100644 avogadro/qtplugins/forcefield/scripts/gfn1.py create mode 100644 avogadro/qtplugins/forcefield/scripts/gfn2.py create mode 100644 avogadro/qtplugins/forcefield/scripts/gfnff.py create mode 100644 avogadro/qtplugins/forcefield/scripts/mmff94.py create mode 100644 avogadro/qtplugins/forcefield/scripts/uff.py diff --git a/avogadro/calc/energymanager.cpp b/avogadro/calc/energymanager.cpp index fdd23cd6dd..76e93d78f0 100644 --- a/avogadro/calc/energymanager.cpp +++ b/avogadro/calc/energymanager.cpp @@ -81,7 +81,11 @@ std::string EnergyManager::nameForModel(const std::string& identifier) const EnergyManager::EnergyManager() { - // add any default models here (EEM maybe?) + // add any default models here + + // LJ is the fallback, since it can handle anything + // (maybe not well, but it can handle it) + addModel(new LennardJones); } EnergyManager::~EnergyManager() diff --git a/avogadro/qtplugins/forcefield/forcefield.cpp b/avogadro/qtplugins/forcefield/forcefield.cpp index ca20d539cf..5c9ba4a10e 100644 --- a/avogadro/qtplugins/forcefield/forcefield.cpp +++ b/avogadro/qtplugins/forcefield/forcefield.cpp @@ -10,10 +10,10 @@ #include #include -#include -#include #include #include +#include +#include #include #include @@ -37,22 +37,16 @@ namespace Avogadro { namespace QtPlugins { +using Avogadro::Calc::EnergyCalculator; using Avogadro::QtGui::Molecule; using Avogadro::QtGui::RWMolecule; -using Avogadro::Calc::EnergyCalculator; const int energyAction = 0; const int optimizeAction = 1; const int configureAction = 2; const int freezeAction = 3; -Forcefield::Forcefield(QObject* parent_) - : ExtensionPlugin(parent_) - , m_molecule(nullptr) - , m_minimizer(LBFGS) - , m_method(0) // just Lennard-Jones for now - , m_maxSteps(100) - , m_outputFormat(nullptr) +Forcefield::Forcefield(QObject* parent_) : ExtensionPlugin(parent_) { refreshScripts(); @@ -74,7 +68,7 @@ Forcefield::Forcefield(QObject* parent_) action = new QAction(this); action->setEnabled(true); - action->setText(tr("Freeze Atoms")); // calculate energy + action->setText(tr("Freeze Atoms")); action->setData(freezeAction); connect(action, SIGNAL(triggered()), SLOT(freezeSelected())); m_actions.push_back(action); @@ -106,7 +100,7 @@ void Forcefield::setMolecule(QtGui::Molecule* mol) m_molecule = mol; - // @todo set molecule for a calculator + // @todo set the method } void Forcefield::refreshScripts() @@ -123,9 +117,7 @@ void Forcefield::optimize() bool isInteractive = m_molecule->undoMolecule()->isInteractive(); m_molecule->undoMolecule()->setInteractive(true); - //@todo check m_minimizer for method to use - //cppoptlib::LbfgsSolver solver; - cppoptlib::ConjugatedGradientDescentSolver solver; + cppoptlib::LbfgsSolver solver; int n = m_molecule->atomCount(); Core::Array pos = m_molecule->atomPositions3d(); @@ -138,42 +130,49 @@ void Forcefield::optimize() // Create a Criteria class to adjust stopping criteria cppoptlib::Criteria crit = cppoptlib::Criteria::defaults(); - // @todo allow criteria to be set - crit.iterations = 5; - crit.fDelta = 1.0e-6; - solver.setStopCriteria(crit); // every 5 steps, update coordinates - cppoptlib::Status status = cppoptlib::Status::NotStarted; + + // e.g., every 5 steps, update coordinates + crit.iterations = m_nSteps; + // we don't set function or gradient criteria + // .. these seem to be broken in the solver code + // .. so we handle ourselves + solver.setStopCriteria(crit); // set the method //@todo check m_method for a particular calculator Calc::LennardJones lj; lj.setMolecule(m_molecule); + double energy = lj.value(positions); for (unsigned int i = 0; i < m_maxSteps / crit.iterations; ++i) { solver.minimize(lj, positions); - cppoptlib::Status currentStatus = solver.status(); - // qDebug() << " status: " << (int)currentStatus; - // qDebug() << " energy: " << lj.value(positions); - // check the gradient norm + double currentEnergy = lj.value(positions); lj.gradient(positions, gradient); - // qDebug() << " gradient: " << gradient.norm(); - - // update coordinates - const double* d = positions.data(); - // casting would be lovely... - for (size_t i = 0; i < n; ++i) { - pos[i] = Vector3(*(d), *(d + 1), *(d + 2)); - d += 3; - - forces[i] = Vector3(gradient[3 * i], gradient[3 * i + 1], - gradient[3 * i + 2]); - } - // todo - merge these into one undo step - m_molecule->undoMolecule()->setAtomPositions3d(pos, tr("Optimize Geometry")); - m_molecule->setForceVectors(forces); - Molecule::MoleculeChanges changes = Molecule::Atoms | Molecule::Modified; - m_molecule->emitChanged(changes); + + // update coordinates + const double* d = positions.data(); + // casting would be lovely... + for (size_t i = 0; i < n; ++i) { + pos[i] = Vector3(*(d), *(d + 1), *(d + 2)); + d += 3; + + forces[i] = -1.0 * + Vector3(gradient[3 * i], gradient[3 * i + 1], gradient[3 * i + 2]); + } + // todo - merge these into one undo step + m_molecule->undoMolecule()->setAtomPositions3d(pos, + tr("Optimize Geometry")); + m_molecule->setForceVectors(forces); + Molecule::MoleculeChanges changes = Molecule::Atoms | Molecule::Modified; + m_molecule->emitChanged(changes); + + // check for convergence + if (fabs(gradient.maxCoeff()) < m_gradientTolerance) + break; + if (fabs(currentEnergy - energy) < m_tolerance) + break; + energy = currentEnergy; } m_molecule->undoMolecule()->setInteractive(isInteractive); @@ -212,5 +211,5 @@ void Forcefield::freezeSelected() } } -} // end QtPlugins -} +} // namespace QtPlugins +} // namespace Avogadro diff --git a/avogadro/qtplugins/forcefield/forcefield.h b/avogadro/qtplugins/forcefield/forcefield.h index a96fffaa99..dbe9bcecbb 100644 --- a/avogadro/qtplugins/forcefield/forcefield.h +++ b/avogadro/qtplugins/forcefield/forcefield.h @@ -16,6 +16,11 @@ class QAction; class QDialog; namespace Avogadro { + +namespace Calc { +class EnergyCalculator; +} + namespace QtPlugins { /** @@ -65,18 +70,20 @@ private slots: private: QList m_actions; - QtGui::Molecule* m_molecule; + QtGui::Molecule* m_molecule = nullptr; - Minimizer m_minimizer; - unsigned int m_method; - unsigned int m_maxSteps; + // defaults + Minimizer m_minimizer = LBFGS; + unsigned int m_maxSteps = 250; + unsigned int m_nSteps = 5; + double m_tolerance = 1.0e-6; + double m_gradientTolerance = 1.0e-4; + Calc::EnergyCalculator *m_method = nullptr; // maps program name --> script file path QMap m_forcefieldScripts; - - const Io::FileFormat* m_outputFormat; - QString m_tempFileName; }; + } } diff --git a/avogadro/qtplugins/forcefield/scripts/gfn1.py b/avogadro/qtplugins/forcefield/scripts/gfn1.py new file mode 100644 index 0000000000..3673b1a8ad --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/gfn1.py @@ -0,0 +1,102 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + import numpy as np + from xtb.libxtb import VERBOSITY_MUTED + from xtb.interface import Calculator, Param + imported = True +except ImportError: + imported = False + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "GFN1", + "identifier": "GFN1", + "description": "Calculate GFN1-xtb energies and gradients", + "inputFormat": "cjson", + "elements": "1-86", + "unitCell": False, + "gradients": True, + "ion": True, + "radical": True, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + with open(filename, "r") as f: + mol_cjson = json.load(f) + + # first setup the calculator + atoms = np.array(mol_cjson["atoms"]["elements"]["number"]) + coord_list = mol_cjson["atoms"]["coords"]["3d"] + coordinates = np.array(coord_list, dtype=float).reshape(-1, 3) + # .. we need to convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # check for total charge + # and spin multiplicity + charge = None # neutral + spin = None # singlet + if "properties" in mol_cjson: + if "totalCharge" in mol_cjson["properties"]: + charge = mol_cjson["properties"]["totalCharge"] + if "totalSpinMultiplicity" in mol_cjson["properties"]: + spin = mol_cjson["properties"]["totalSpinMultiplicity"] + + calc = Calculator(Param.GFN1xTB, atoms, coordinates, + charge=charge, uhf=spin) + calc.set_verbosity(VERBOSITY_MUTED) + res = calc.singlepoint() + + # we loop forever - Avogadro will kill the process when done + while(True): + # first print the energy of these coordinates + print(res.get_energy()) # in Hartree + + # now print the gradient + # .. we don't want the "[]" in the output + output = np.array2string(res.get_gradient()) + output = output.replace("[", "").replace("]", "") + print(output) + + # read new coordinates from stdin + for i in range(len(atoms)): + coordinates[i] = np.fromstring(input(), sep=" ") + # .. convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # update the calculator and run a new calculation + calc.update(coordinates) + calc.singlepoint(res) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("GFN1 calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("xtb-python is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/forcefield/scripts/gfn2.py b/avogadro/qtplugins/forcefield/scripts/gfn2.py new file mode 100644 index 0000000000..e7035d147b --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/gfn2.py @@ -0,0 +1,102 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + import numpy as np + from xtb.libxtb import VERBOSITY_MUTED + from xtb.interface import Calculator, Param + imported = True +except ImportError: + imported = False + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "GFN2", + "identifier": "GFN2", + "description": "Calculate GFN2-xtb energies and gradients", + "inputFormat": "cjson", + "elements": "1-86", + "unitCell": False, + "gradients": True, + "ion": True, + "radical": True, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + with open(filename, "r") as f: + mol_cjson = json.load(f) + + # first setup the calculator + atoms = np.array(mol_cjson["atoms"]["elements"]["number"]) + coord_list = mol_cjson["atoms"]["coords"]["3d"] + coordinates = np.array(coord_list, dtype=float).reshape(-1, 3) + # .. we need to convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # check for total charge + # and spin multiplicity + charge = None # neutral + spin = None # singlet + if "properties" in mol_cjson: + if "totalCharge" in mol_cjson["properties"]: + charge = mol_cjson["properties"]["totalCharge"] + if "totalSpinMultiplicity" in mol_cjson["properties"]: + spin = mol_cjson["properties"]["totalSpinMultiplicity"] + + calc = Calculator(Param.GFN2xTB, atoms, coordinates, + charge=charge, uhf=spin) + calc.set_verbosity(VERBOSITY_MUTED) + res = calc.singlepoint() + + # we loop forever - Avogadro will kill the process when done + while(True): + # first print the energy of these coordinates + print(res.get_energy()) # in Hartree + + # now print the gradient + # .. we don't want the "[]" in the output + output = np.array2string(res.get_gradient()) + output = output.replace("[", "").replace("]", "") + print(output) + + # read new coordinates from stdin + for i in range(len(atoms)): + coordinates[i] = np.fromstring(input(), sep=" ") + # .. convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # update the calculator and run a new calculation + calc.update(coordinates) + calc.singlepoint(res) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("GFN2 calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("xtb-python is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/forcefield/scripts/gfnff.py b/avogadro/qtplugins/forcefield/scripts/gfnff.py new file mode 100644 index 0000000000..c6e6223f9c --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/gfnff.py @@ -0,0 +1,104 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + import numpy as np + from xtb.libxtb import VERBOSITY_MUTED + from xtb.interface import Calculator, Param + + imported = True +except ImportError: + imported = False + + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "GFN-FF", + "identifier": "GFN-FF", + "description": "Calculate GFNFF-xtb energies and gradients", + "inputFormat": "cjson", + "elements": "1-86", + "unitCell": False, + "gradients": True, + "ion": True, + "radical": False, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + with open(filename, "r") as f: + mol_cjson = json.load(f) + + # first setup the calculator + atoms = np.array(mol_cjson["atoms"]["elements"]["number"]) + coord_list = mol_cjson["atoms"]["coords"]["3d"] + coordinates = np.array(coord_list, dtype=float).reshape(-1, 3) + # .. we need to convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # check for total charge + # and spin multiplicity + charge = None # neutral + spin = None # singlet + if "properties" in mol_cjson: + if "totalCharge" in mol_cjson["properties"]: + charge = mol_cjson["properties"]["totalCharge"] + if "totalSpinMultiplicity" in mol_cjson["properties"]: + spin = mol_cjson["properties"]["totalSpinMultiplicity"] + + calc = Calculator(Param.GFNFF, atoms, coordinates, + charge=charge, uhf=spin) + calc.set_verbosity(VERBOSITY_MUTED) + res = calc.singlepoint() + + # we loop forever - Avogadro will kill our process when done + while True: + # first print the energy of these coordinates + print(res.get_energy()) # in Hartree + + # now print the gradient + # .. we don't want the "[]" in the output + output = np.array2string(res.get_gradient()) + output = output.replace("[", "").replace("]", "") + print(output) + + # read new coordinates from stdin + for i in range(len(atoms)): + coordinates[i] = np.fromstring(input(), sep=" ") + # .. convert from Angstrom to Bohr + coordinates /= 0.52917721067 + + # update the calculator and run a new calculation + calc.update(coordinates) + calc.singlepoint(res) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("GFN2 calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("xtb-python is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/forcefield/scripts/mmff94.py b/avogadro/qtplugins/forcefield/scripts/mmff94.py new file mode 100644 index 0000000000..e0f826e3da --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/mmff94.py @@ -0,0 +1,85 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + from openbabel import pybel + import numpy as np + + imported = True +except ImportError: + imported = False + + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "MMFF94", + "identifier": "MMFF94", + "description": "Calculate MMFF94 energies and gradients", + "inputFormat": "cml", + "elements": "1,6-9,14-17,35,53", + "unitCell": False, + "gradients": True, + "ion": False, + "radical": False, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + mol = next(pybel.readfile("cml", filename)) + + ff = pybel._forcefields["mmff94"] + success = ff.Setup(mol.OBMol) + if not success: + # should never happen, but just in case + sys.exit("MMFF94 force field setup failed") + + # we loop forever - Avogadro will kill the process when done + num_atoms = len(mol.atoms) + while True: + # first print the energy of these coordinates + print(ff.Energy(True)) # in Hartree + + # now print the gradient on each atom + for atom in mol.atoms: + grad = ff.GetGradient(atom.OBAtom) + print(grad.GetX(), grad.GetY(), grad.GetZ()) + + # read new coordinates from stdin + for i in range(num_atoms): + coordinates = np.fromstring(input(), sep=" ") + atom = mol.atoms[i] + atom.OBAtom.SetVector(coordinates[0], coordinates[1], coordinates[2]) + + # update the molecule geometry for the next energy + ff.SetCoordinates(mol.OBMol) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("MMFF94 calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("pybel is unavailable") + elif args["file"]: + run(args["file"][0]) diff --git a/avogadro/qtplugins/forcefield/scripts/uff.py b/avogadro/qtplugins/forcefield/scripts/uff.py new file mode 100644 index 0000000000..7d35226262 --- /dev/null +++ b/avogadro/qtplugins/forcefield/scripts/uff.py @@ -0,0 +1,85 @@ +# This source file is part of the Avogadro project. +# This source code is released under the 3-Clause BSD License, (see "LICENSE"). + +import argparse +import json +import sys + +try: + from openbabel import pybel + import numpy as np + + imported = True +except ImportError: + imported = False + + +def getMetaData(): + # before we return metadata, make sure xtb is in the path + if not imported: + return {} # Avogadro will ignore us now + + metaData = { + "name": "UFF", + "identifier": "UFF", + "description": "Calculate UFF energies and gradients", + "inputFormat": "cml", + "elements": "1-86", + "unitCell": False, + "gradients": True, + "ion": False, + "radical": False, + } + return metaData + + +def run(filename): + # we get the molecule from the supplied filename + # in cjson format (it's a temporary file created by Avogadro) + mol = next(pybel.readfile("cml", filename)) + + ff = pybel._forcefields["uff"] + success = ff.Setup(mol.OBMol) + if not success: + # should never happen, but just in case + sys.exit("UFF force field setup failed") + + # we loop forever - Avogadro will kill the process when done + num_atoms = len(mol.atoms) + while True: + # first print the energy of these coordinates + print(ff.Energy(True)) # in Hartree + + # now print the gradient on each atom + for atom in mol.atoms: + grad = ff.GetGradient(atom.OBAtom) + print(grad.GetX(), grad.GetY(), grad.GetZ()) + + # read new coordinates from stdin + for i in range(num_atoms): + coordinates = np.fromstring(input(), sep=" ") + atom = mol.atoms[i] + atom.OBAtom.SetVector(coordinates[0], coordinates[1], coordinates[2]) + + # update the molecule geometry for the next energy + ff.SetCoordinates(mol.OBMol) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser("UFF calculator") + parser.add_argument("--display-name", action="store_true") + parser.add_argument("--metadata", action="store_true") + parser.add_argument("-f", "--file", nargs=1) + parser.add_argument("--lang", nargs="?", default="en") + args = vars(parser.parse_args()) + + if args["metadata"]: + print(json.dumps(getMetaData())) + elif args["display_name"]: + name = getMetaData().get("name") + if name: + print(name) + else: + sys.exit("pybel is unavailable") + elif args["file"]: + run(args["file"][0])