Add obenergy code

Signed-off-by: Geoff Hutchison <[email protected]>

avogadro/qtplugins/forcefield/obenergy.cpp

@@ -0,0 +1,184 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+  It links to the Open Babel library, which is released under the GNU GPL v2.
+******************************************************************************/
+
+#include "obenergy.h"
+
+#include <avogadro/core/molecule.h>
+
+#include <openbabel/atom.h>
+#include <openbabel/babelconfig.h>
+#include <openbabel/base.h>
+#include <openbabel/forcefield.h>
+#include <openbabel/math/vector3.h>
+#include <openbabel/mol.h>
+#include <openbabel/obconversion.h>
+#include <openbabel/obiter.h>
+
+#include <QDebug>
+#include <QDir>
+
+using namespace OpenBabel;
+
+namespace Avogadro::QtPlugins {
+
+class OBEnergy::Private
+{
+public:
+  // OBMol and OBForceField are owned by this class
+  OBMol* m_obmol = nullptr;
+  OBForceField* m_forceField = nullptr;
+
+  ~Private()
+  {
+    delete m_obmol;
+    delete m_forceField;
+  }
+};
+
+OBEnergy::OBEnergy(const std::string& method)
+  : m_identifier(method), m_name(method), m_molecule(nullptr)
+{
+  d = new Private;
+  // Ensure the plugins are loaded
+  OBPlugin::LoadAllPlugins();
+
+  d->m_forceField = static_cast<OBForceField*>(
+    OBPlugin::GetPlugin("forcefields", method.c_str()));
+
+  if (method == "UFF") {
+    m_description = tr("Universal Force Field");
+    m_elements.reset();
+    for (unsigned int i = 1; i < 102; ++i)
+      m_elements.set(i);
+  } else if (method == "GAFF") {
+    m_description = tr("Generalized Amber Force Field");
+
+    // H, C, N, O, F, P, S, Cl, Br, and I
+    m_elements.set(1);
+    m_elements.set(6);
+    m_elements.set(7);
+    m_elements.set(8);
+    m_elements.set(9);
+    m_elements.set(15);
+    m_elements.set(16);
+    m_elements.set(17);
+    m_elements.set(35);
+    m_elements.set(53);
+  } else if (method == "MMFF94") {
+    m_description = tr("Merck Molecular Force Field 94");
+    m_elements.reset();
+
+    // H, C, N, O, F, Si, P, S, Cl, Br, and I
+    m_elements.set(1);
+    m_elements.set(6);
+    m_elements.set(7);
+    m_elements.set(8);
+    m_elements.set(9);
+    m_elements.set(14);
+    m_elements.set(15);
+    m_elements.set(16);
+    m_elements.set(17);
+    m_elements.set(35);
+    m_elements.set(53);
+  }
+}
+
+OBEnergy::~OBEnergy() {}
+
+Calc::EnergyCalculator* OBEnergy::newInstance() const
+{
+  return new OBEnergy(m_name);
+}
+
+void OBEnergy::setMolecule(Core::Molecule* mol)
+{
+  m_molecule = mol;
+
+  if (mol == nullptr || mol->atomCount() == 0) {
+    return; // nothing to do
+  }
+
+  // set up our internal OBMol
+  d->m_obmol = new OBMol;
+  // copy the atoms, bonds, and coordinates
+  d->m_obmol->BeginModify();
+  for (size_t i = 0; i < mol->atomCount(); ++i) {
+    const Core::Atom& atom = mol->atom(i);
+    OBAtom* obAtom = d->m_obmol->NewAtom();
+    obAtom->SetAtomicNum(atom.atomicNumber());
+    auto pos = atom.position3d().cast<double>();
+    obAtom->SetVector(pos.x(), pos.y(), pos.z());
+  }
+  for (size_t i = 0; i < mol->bondCount(); ++i) {
+    const Core::Bond& bond = mol->bond(i);
+    d->m_obmol->AddBond(bond.atom1().index() + 1, bond.atom2().index() + 1,
+                        bond.order());
+  }
+  d->m_obmol->EndModify();
+
+  // make sure we can set up the force field
+  if (d->m_forceField != nullptr) {
+    d->m_forceField->Setup(*d->m_obmol);
+  } else {
+    d->m_forceField = static_cast<OBForceField*>(
+      OBPlugin::GetPlugin("forcefields", m_identifier.c_str()));
+    if (d->m_forceField != nullptr) {
+      d->m_forceField->Setup(*d->m_obmol);
+    }
+  }
+}
+
+Real OBEnergy::value(const Eigen::VectorXd& x)
+{
+  if (m_molecule == nullptr || m_molecule->atomCount() == 0)
+    return 0.0; // nothing to do
+
+  // update coordinates in our private OBMol
+  for (size_t i = 0; i < m_molecule->atomCount(); ++i) {
+    Eigen::Vector3d pos(x[i * 3], x[i * 3 + 1], x[i * 3 + 2]);
+    d->m_obmol->GetAtom(i + 1)->SetVector(pos.x(), pos.y(), pos.z());
+  }
+
+  double energy = 0.0;
+  if (d->m_forceField != nullptr) {
+    d->m_forceField->SetCoordinates(*d->m_obmol);
+    energy = d->m_forceField->Energy(false);
+  }
+  return energy;
+}
+
+void OBEnergy::gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad)
+{
+  if (m_molecule == nullptr || m_molecule->atomCount() == 0)
+    return;
+
+  // update coordinates in our private OBMol
+  for (size_t i = 0; i < m_molecule->atomCount(); ++i) {
+    Eigen::Vector3d pos(x[i * 3], x[i * 3 + 1], x[i * 3 + 2]);
+    d->m_obmol->GetAtom(i + 1)->SetVector(pos.x(), pos.y(), pos.z());
+  }
+
+  if (d->m_forceField != nullptr) {
+    d->m_forceField->SetCoordinates(*d->m_obmol);
+
+    // make sure gradients are calculated
+    double energy = d->m_forceField->Energy(true);
+    for (size_t i = 0; i < m_molecule->atomCount(); ++i) {
+      OBAtom* atom = d->m_obmol->GetAtom(i + 1);
+      OpenBabel::vector3 obGrad = d->m_forceField->GetGradient(atom);
+      grad[3 * i] = obGrad.x();
+      grad[3 * i + 1] = obGrad.y();
+      grad[3 * i + 2] = obGrad.z();
+    }
+
+    grad *= -1; // OpenBabel outputs forces, not grads
+    cleanGradients(grad);
+  }
+}
+
+} // namespace Avogadro::QtPlugins
+
+#include "obenergy.moc"

avogadro/qtplugins/forcefield/obenergy.h

@@ -0,0 +1,66 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+  It links to the Open Babel library, which is released under the GNU GPL v2.
+******************************************************************************/
+
+#ifndef AVOGADRO_QTPLUGINS_OBENERGY_H
+#define AVOGADRO_QTPLUGINS_OBENERGY_H
+
+#include <avogadro/calc/energycalculator.h>
+
+#include <avogadro/core/avogadrocore.h>
+
+#include <QtCore/QCoreApplication>
+#include <QtCore/QString>
+
+namespace Avogadro {
+
+namespace QtPlugins {
+
+class OBEnergy : public Avogadro::Calc::EnergyCalculator
+{
+  Q_DECLARE_TR_FUNCTIONS(OBEnergy)
+
+public:
+  OBEnergy(const std::string& method = "");
+  ~OBEnergy() override;
+
+  std::string method() const { return m_identifier; }
+  void setupProcess();
+
+  Calc::EnergyCalculator* newInstance() const override;
+
+  std::string identifier() const override { return m_identifier; }
+  std::string name() const override { return m_name; }
+  std::string description() const override
+  {
+    return m_description.toStdString();
+  }
+
+  Core::Molecule::ElementMask elements() const override { return (m_elements); }
+
+  // This will check if the molecule is valid for this script
+  // and then start the external process
+  void setMolecule(Core::Molecule* mol) override;
+  // energy
+  Real value(const Eigen::VectorXd& x) override;
+  // gradient (which may be unsupported and fall back to numeric)
+  void gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) override;
+
+private:
+  class Private;
+
+  Core::Molecule* m_molecule;
+  Private* d;
+
+  Core::Molecule::ElementMask m_elements;
+  std::string m_identifier;
+  std::string m_name;
+  QString m_description;
+};
+
+} // namespace QtPlugins
+} // namespace Avogadro
+
+#endif // AVOGADRO_QTPLUGINS_OBENERGY_H