diff --git a/avogadro/calc/energymanager.cpp b/avogadro/calc/energymanager.cpp index 53751e45d..925569d04 100644 --- a/avogadro/calc/energymanager.cpp +++ b/avogadro/calc/energymanager.cpp @@ -113,6 +113,9 @@ std::set 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 diff --git a/avogadro/qtplugins/forcefield/CMakeLists.txt b/avogadro/qtplugins/forcefield/CMakeLists.txt index 62079d56a..0ffd8de43 100644 --- a/avogadro/qtplugins/forcefield/CMakeLists.txt +++ b/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/") diff --git a/avogadro/qtplugins/forcefield/forcefield.cpp b/avogadro/qtplugins/forcefield/forcefield.cpp index 9c4092535..9dcf0fcbf 100644 --- a/avogadro/qtplugins/forcefield/forcefield.cpp +++ b/avogadro/qtplugins/forcefield/forcefield.cpp @@ -8,6 +8,10 @@ #include "obmmenergy.h" #include "scriptenergy.h" +#ifdef BUILD_GPL_PLUGINS +#include "obenergy.h" +#endif + #include #include @@ -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); diff --git a/avogadro/qtplugins/forcefield/obenergy.cpp b/avogadro/qtplugins/forcefield/obenergy.cpp new file mode 100644 index 000000000..de28ede83 --- /dev/null +++ b/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 + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +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( + 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(); + 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( + 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" diff --git a/avogadro/qtplugins/forcefield/obenergy.h b/avogadro/qtplugins/forcefield/obenergy.h new file mode 100644 index 000000000..c9a605dce --- /dev/null +++ b/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 + +#include + +#include +#include + +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 diff --git a/cmake/FindOpenBabel3.cmake b/cmake/FindOpenBabel3.cmake new file mode 100644 index 000000000..ba5286ddc --- /dev/null +++ b/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()