Add an optional calculator that links Open Babel (#1591)

* Optional: Use OB code if GPL plugins are used. Otherwise scripts & obmm

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

avogadro/calc/energymanager.cpp

@@ -113,6 +113,9 @@ std::set<std::string> EnergyManager::identifiersForMolecule(
 
   // check our models for compatibility
   for (auto m_model : m_models) {
+    if (m_model == nullptr)
+      continue;
+
     // we can check easy things first
     // - is the molecule an ion based on total charge
     // - is the molecule a radical based on spin multiplicity

avogadro/qtplugins/forcefield/CMakeLists.txt

@@ -5,6 +5,18 @@ set(forcefield_srcs
   scriptenergy.cpp
 )
 
+if (BUILD_GPL_PLUGINS)
+  find_package(OpenBabel3)
+  if (OpenBabel3_LIBRARY)
+    list(APPEND forcefield_srcs
+      obenergy.cpp
+    )
+    add_definitions(-DBUILD_GPL_PLUGINS)
+    include_directories(${OpenBabel3_INCLUDE_DIRS} ${OpenBabel3_INCLUDE_DIR}
+    ${AvogadroLibs_BINARY_DIR}/../prefix/include/openbabel3)
+  endif()
+endif()
+
 avogadro_plugin(Forcefield
   "Force field optimization and dynamics"
   ExtensionPlugin
@@ -16,16 +28,26 @@ avogadro_plugin(Forcefield
 
 target_link_libraries(Forcefield PRIVATE Avogadro::Calc)
 
+if (BUILD_GPL_PLUGINS AND OpenBabel3_LIBRARY)
+  target_link_libraries(Forcefield PRIVATE OpenBabel3)
+endif()
+
 # Bundled forcefield scripts
 set(forcefields
   scripts/ani2x.py
-  scripts/gaff.py
   scripts/gfn1.py
   scripts/gfn2.py
   scripts/gfnff.py
-  scripts/mmff94.py
-  scripts/uff.py
 )
 
+if (NOT BUILD_GPL_PLUGINS)
+  # install the OB / Pybel forcefield scripts
+  list(APPEND forcefields
+    scripts/gaff.py
+    scripts/mmff94.py
+    scripts/uff.py
+  )
+endif()
+
 install(PROGRAMS ${forcefields}
 DESTINATION "${INSTALL_LIBRARY_DIR}/avogadro2/scripts/energy/")

avogadro/qtplugins/forcefield/forcefield.cpp

@@ -8,6 +8,10 @@
 #include "obmmenergy.h"
 #include "scriptenergy.h"
 
+#ifdef BUILD_GPL_PLUGINS
+#include "obenergy.h"
+#endif
+
 #include <QtCore/QDebug>
 #include <QtCore/QSettings>
 
@@ -61,12 +65,22 @@ Forcefield::Forcefield(QObject* parent_)
   m_gradientTolerance = settings.value("gradientTolerance", 1.0e-4).toDouble();
   settings.endGroup();
 
+  // add the openbabel calculators in case they don't exist
+#ifdef BUILD_GPL_PLUGINS
+  // These directly use Open Babel and are fast
+  Calc::EnergyManager::registerModel(new OBEnergy("MMFF94"));
+  Calc::EnergyManager::registerModel(new OBEnergy("UFF"));
+  Calc::EnergyManager::registerModel(new OBEnergy("GAFF"));
+#endif
+
   refreshScripts();
 
-  // add the openbabel calculators in case they don't exist
+#ifndef BUILD_GPL_PLUGINS
+  // These call obmm and can be slow
   Calc::EnergyManager::registerModel(new OBMMEnergy("MMFF94"));
   Calc::EnergyManager::registerModel(new OBMMEnergy("UFF"));
   Calc::EnergyManager::registerModel(new OBMMEnergy("GAFF"));
+#endif
 
   QAction* action = new QAction(this);
   action->setEnabled(true);

avogadro/qtplugins/forcefield/obenergy.cpp

@@ -0,0 +1,185 @@
+/******************************************************************************
+  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/babelconfig.h>
+
+#include <openbabel/atom.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

cmake/FindOpenBabel3.cmake

@@ -0,0 +1,31 @@
+# Find the OpenBabel3 library
+#
+# Defines:
+#
+#  OpenBabel3_FOUND        - system has OpenBabel
+#  OpenBabel3_INCLUDE_DIRS - the OpenBabel include directories
+#  OpenBabel3_LIBRARY      - The OpenBabel library
+#
+find_path(OpenBabel3_INCLUDE_DIR openbabel3/openbabel/babelconfig.h)
+if(OPENBABEL3_INCLUDE_DIR)
+  set(OPENBABEL3_INCLUDE_DIR ${OPENBABEL3_INCLUDE_DIR}/openbabel3)
+endif()
+find_library(OpenBabel3_LIBRARY NAMES openbabel openbabel3 openbabel-3)
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(OpenBabel3 DEFAULT_MSG OpenBabel3_INCLUDE_DIR
+                                 OpenBabel3_LIBRARY)
+
+mark_as_advanced(OpenBabel3_INCLUDE_DIR OpenBabel3_LIBRARY)
+
+if(OpenBabel3_FOUND)
+  set(OpenBabel3_INCLUDE_DIRS "${OpenBabel3_INCLUDE_DIR}")
+
+  if(NOT TARGET OpenBabel3)
+    add_library(OpenBabel3 SHARED IMPORTED GLOBAL)
+    set_target_properties(OpenBabel3 PROPERTIES
+      IMPORTED_LOCATION "${OpenBabel3_LIBRARY}"
+      IMPORTED_IMPLIB "${OpenBabel3_LIBRARY}"
+      INTERFACE_INCLUDE_DIRECTORIES "${OpenBabel3_INCLUDE_DIR}")
+  endif()
+endif()