From dfad83196265387f5b1db831b82bd78199be4241 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 14 Oct 2023 13:55:36 -0400 Subject: [PATCH] More cleanups, including a manager class to filter methods Signed-off-by: Geoff Hutchison --- avogadro/calc/CMakeLists.txt | 2 + avogadro/calc/chargemanager.h | 2 +- avogadro/calc/energycalculator.h | 29 ++++ avogadro/calc/energymanager.cpp | 124 +++++++++++++++ avogadro/calc/energymanager.h | 143 ++++++++++++++++++ avogadro/calc/flameminimize.cpp | 113 -------------- avogadro/calc/flameminimize.h | 49 ------ avogadro/calc/lennardjones.cpp | 4 + avogadro/calc/lennardjones.h | 21 ++- avogadro/qtplugins/forcefield/CMakeLists.txt | 12 +- .../qtplugins/forcefield/forcefielddialog.ui | 51 +++++-- 11 files changed, 365 insertions(+), 185 deletions(-) create mode 100644 avogadro/calc/energymanager.cpp create mode 100644 avogadro/calc/energymanager.h delete mode 100644 avogadro/calc/flameminimize.cpp delete mode 100644 avogadro/calc/flameminimize.h diff --git a/avogadro/calc/CMakeLists.txt b/avogadro/calc/CMakeLists.txt index 12b5d06a0e..df6deef248 100644 --- a/avogadro/calc/CMakeLists.txt +++ b/avogadro/calc/CMakeLists.txt @@ -5,6 +5,7 @@ avogadro_headers(Calc chargemanager.h defaultmodel.h energycalculator.h + energymanager.h lennardjones.h ) @@ -13,6 +14,7 @@ target_sources(Calc PRIVATE chargemanager.cpp defaultmodel.cpp energycalculator.cpp + energymanager.cpp lennardjones.cpp ) diff --git a/avogadro/calc/chargemanager.h b/avogadro/calc/chargemanager.h index 7d611c5b8e..13137d8598 100644 --- a/avogadro/calc/chargemanager.h +++ b/avogadro/calc/chargemanager.h @@ -53,7 +53,7 @@ class AVOGADROCALC_EXPORT ChargeManager /** * @brief Register a new charge model with the manager. - * @param format An instance of the format to manage, the manager assumes + * @param model An instance of the model to manage, the manager assumes * ownership of the object passed in. * @return True on success, false on failure. */ diff --git a/avogadro/calc/energycalculator.h b/avogadro/calc/energycalculator.h index 8b65d53345..f59c73546b 100644 --- a/avogadro/calc/energycalculator.h +++ b/avogadro/calc/energycalculator.h @@ -8,6 +8,7 @@ #include "avogadrocalcexport.h" +#include #include #include @@ -46,6 +47,34 @@ class AVOGADROCALC_EXPORT EnergyCalculator : public cppoptlib::Problem */ virtual bool setConfiguration(Core::VariantMap& config) { return true; } + /** + * @brief Indicate if your method only treats a subset of elements + * @return an element mask corresponding to the defined subset + */ + virtual Core::Molecule::ElementMask elements() const = 0; + + /** + * @brief Indicate if your method can handle unit cells + * @return true if unit cells are supported + */ + virtual bool acceptsUnitCell() const { return false; } + + /** + * @brief Indicate if your method can handle ions + * Many methods only treat neutral systems, either + * a neutral molecule or a neutral unit cell. + * + * @return true if ions are supported + */ + virtual bool acceptsIons() const { return false; } + + /** + * @brief Indicate if your method can handle radicals + * Most methods only treat closed-shell molecules. + * @return true if radicals are supported + */ + virtual bool acceptsRadicals() const { return false; } + /** * Calculate the gradients for this method, defaulting to numerical * finite-difference methods diff --git a/avogadro/calc/energymanager.cpp b/avogadro/calc/energymanager.cpp new file mode 100644 index 0000000000..fdd23cd6dd --- /dev/null +++ b/avogadro/calc/energymanager.cpp @@ -0,0 +1,124 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#include "energymanager.h" +#include "energycalculator.h" +#include "lennardjones.h" + +#include +#include + +namespace Avogadro::Calc { + +EnergyManager& EnergyManager::instance() +{ + static EnergyManager instance; + return instance; +} + +void EnergyManager::appendError(const std::string& errorMessage) +{ + m_error += errorMessage + "\n"; +} + +bool EnergyManager::registerModel(EnergyCalculator* model) +{ + return instance().addModel(model); +} + +bool EnergyManager::unregisterModel(const std::string& identifier) +{ + return instance().removeModel(identifier); +} + +bool EnergyManager::addModel(EnergyCalculator* model) +{ + if (model == nullptr) { + appendError("Supplied model was null."); + return false; + } + + if (m_identifiers.find(model->identifier()) != m_identifiers.end()) { + appendError("Model " + model->identifier() + " already loaded."); + return false; + } + + // If we got here then the format is unique enough to be added. + size_t index = m_models.size(); + m_models.push_back(model); + m_identifiers[model->identifier()] = index; + m_identifierToName[model->identifier()] = model->name(); + + return true; +} + +bool EnergyManager::removeModel(const std::string& identifier) +{ + auto ids = m_identifiers[identifier]; + m_identifiers.erase(identifier); + m_identifierToName.erase(identifier); + + auto* model = m_models[ids]; + + if (model != nullptr) { + m_models[ids] = nullptr; + delete model; + } + + return true; +} + +std::string EnergyManager::nameForModel(const std::string& identifier) const +{ + auto it = m_identifierToName.find(identifier); + if (it == m_identifierToName.end()) { + return identifier; + } + return it->second; +} + +EnergyManager::EnergyManager() +{ + // add any default models here (EEM maybe?) +} + +EnergyManager::~EnergyManager() +{ + // Delete the models that were loaded. + for (auto& m_model : m_models) { + delete m_model; + } + m_models.clear(); +} + +std::set EnergyManager::identifiersForMolecule( + const Core::Molecule& molecule) const +{ + std::set identifiers; + + // check our models for compatibility + for (auto m_model : m_models) { + // we can check easy things first + // - is the molecule an ion based on total charge + // - is the molecule a radical based on spin multiplicity + // - does the molecule have a unit cell + if (molecule.totalCharge() != 0 && !m_model->acceptsIons()) + continue; + if (molecule.totalSpinMultiplicity() != 1 && !m_model->acceptsRadicals()) + continue; + if (molecule.unitCell() != nullptr && !m_model->acceptsUnitCell()) + continue; + + // Finally, we check that every element in the molecule + // is handled by the model + auto mask = m_model->elements() & molecule.elements(); + if (mask.count() == molecule.elements().count()) + identifiers.insert(m_model->identifier()); // this one will work + } + + return identifiers; +} + +} // namespace Avogadro::Calc diff --git a/avogadro/calc/energymanager.h b/avogadro/calc/energymanager.h new file mode 100644 index 0000000000..308bd501f4 --- /dev/null +++ b/avogadro/calc/energymanager.h @@ -0,0 +1,143 @@ +/****************************************************************************** + This source file is part of the Avogadro project. + This source code is released under the 3-Clause BSD License, (see "LICENSE"). +******************************************************************************/ + +#ifndef AVOGADRO_CALC_ENERGYMANAGER_H +#define AVOGADRO_CALC_ENERGYMANAGER_H + +#include "avogadrocalcexport.h" + +#include +#include +#include + +#include +#include +#include +#include + +namespace Avogadro { +namespace Core { +class Molecule; +} +namespace Calc { + +class EnergyCalculator; + +/** + * @class EnergyManager chargemanager.h + * + * @brief Class to manage registration, searching and creation of force field + * (energy) calculators. + * @author Geoffrey R. Hutchison + * + * The energy manager is a singleton class that handles the runtime + * registration, search, creation and eventual destruction of calculators + * for geometry optimization and molecular dynamics. + * It can be used to gain a listing of available models, register new + * models, etc. + * + * All energy calculation can take place independent of this code, but for + * automated registration and look up, this is the preferred API. + */ +class AVOGADROCALC_EXPORT EnergyManager +{ +public: + /** + * Get the singleton instance of the energy manager. This instance should + * not be deleted. + */ + static EnergyManager& instance(); + + /** + * @brief Register a new model with the manager. + * @param model An instance of the calculator to manage, the manager assumes + * ownership of the object passed in. + * @return True on success, false on failure. + */ + static bool registerModel(EnergyCalculator* model); + + /** + * @brief Unregister a charge model from the manager. + * @param identifier The identifier for the model to remove. + * @return True on success, false on failure. + */ + static bool unregisterModel(const std::string& identifier); + + /** + * Add the supplied @p model to the manager, registering its ID and other + * relevant data for later lookup. The manager assumes ownership of the + * supplied object. + * @return True on success, false on failure. + */ + bool addModel(EnergyCalculator* model); + + /** + * Remove the model with the identifier @a identifier from the manager. + * @return True on success, false on failure. + */ + bool removeModel(const std::string& identifier); + + /** + * New instance of the model for the specified @p identifier. Ownership + * is passed to the caller. + * @param identifier The unique identifier of the model. + * @return Instance of the model, nullptr if not found. Ownership passes to + * the caller. + */ + EnergyCalculator* newModelFromIdentifier(const std::string& identifier) const; + + /** + * Get a list of all loaded identifiers + */ + std::set identifiers() const; + + /** + * @brief Get a list of models that work for this molecule. + * + * This is probably the method you want to get a list for a user + */ + std::set identifiersForMolecule( + const Core::Molecule& molecule) const; + + /** + * @brief Get the name of the model for the specified identifier. + * + * The name is a user-visible string, and may be translated. + * @param identifier The unique identifier of the model. + * @return The name of the model, or an empty string if not found. + */ + std::string nameForModel(const std::string& identifier) const; + + /** + * Get any errors that have been logged when loading models. + */ + std::string error() const; + +private: + typedef std::map ModelIdMap; + + EnergyManager(); + ~EnergyManager(); + + EnergyManager(const EnergyManager&); // Not implemented. + EnergyManager& operator=(const EnergyManager&); // Not implemented. + + /** + * @brief Append warnings/errors to the error message string. + * @param errorMessage The error message to append. + */ + void appendError(const std::string& errorMessage); + + std::vector m_models; + mutable ModelIdMap m_identifiers; + mutable std::map m_identifierToName; + + std::string m_error; +}; + +} // namespace Calc +} // namespace Avogadro + +#endif // AVOGADRO_CALC_ENERGYMANAGER_H diff --git a/avogadro/calc/flameminimize.cpp b/avogadro/calc/flameminimize.cpp deleted file mode 100644 index 7946855088..0000000000 --- a/avogadro/calc/flameminimize.cpp +++ /dev/null @@ -1,113 +0,0 @@ -/****************************************************************************** - This source file is part of the Avogadro project. - This source code is released under the 3-Clause BSD License, (see "LICENSE"). -******************************************************************************/ - -#include "flameminimize.h" - -#include - -#include - -namespace Avogadro::Calc { - -FlameMinimize::FlameMinimize() - : m_molecule(nullptr), m_calc(nullptr) -{} - -void FlameMinimize::setMolecule(Core::Molecule* mol) -{ - m_molecule = mol; - - m_forces.resize(3 * mol->atomCount()); - m_forces.setZero(); - m_velocities.resize(3 * mol->atomCount()); - m_velocities.setZero(); - m_accel.resize(3 * mol->atomCount()); - m_accel.setZero(); - - m_invMasses.resize(3 * mol->atomCount()); - m_invMasses.setZero(); - for (unsigned int i = 0; i < mol->atomCount(); ++i) { - //@todo should this be set to 1.0 amu? - double scaledMass = log(Core::Elements::mass(mol->atom(i).atomicNumber())); - - m_invMasses[3 * i] = 1.0 / scaledMass; - m_invMasses[3 * i + 1] = 1.0 / scaledMass; - m_invMasses[3 * i + 2] = 1.0 / scaledMass; - } -} - -bool FlameMinimize::minimize(EnergyCalculator& calc, - Eigen::VectorXd& positions) -{ - if (m_molecule == nullptr) - return false; - - m_calc = &calc; - - //@todo - set convergence criteria (e.g., max steps, min gradients, energy, - // etc.) - - double alpha = 0.1; // start - double deltaT = 0.1 * 1.0e-15; // fs - unsigned int positiveSteps = 0; - - m_forces.setZero(); - m_velocities.setZero(); - m_accel.setZero(); - - for (unsigned int i = 0; i < 20; ++i) { - verletIntegrate(positions, deltaT); - //qDebug() << "vvi forces " << m_forces.norm() << " vel " << m_velocities.norm(); - - // Step 1 - double power = m_forces.dot(m_velocities); - - // Step 2 - m_velocities = (1.0 - alpha) * m_velocities + alpha* - m_forces.cwiseProduct(m_velocities.cwiseAbs()); - - if (power > 0.0) { - // Step 3 - positiveSteps++; - if (positiveSteps > 5) { - deltaT = std::min(1.1 * deltaT, 1.0); - alpha = 0.99 * alpha; - } - } else { - // Step 4 - positiveSteps = 0; - deltaT = 0.5 * deltaT; - m_velocities.setZero(); - alpha = 0.1; - } - - double Frms = m_forces.norm() / sqrt(positions.rows()); - if (Frms < 1.0e-5) - break; - } - - return true; -} - -void FlameMinimize::verletIntegrate(Eigen::VectorXd& positions, double deltaT) -{ - // See https://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet - // (as one of many examples) - if (m_molecule == nullptr || m_calc == nullptr) - return; - - positions += deltaT * m_velocities + (deltaT * deltaT / 2.0) * m_accel; - m_calc->gradient(positions, m_forces); - m_forces = -1*m_forces; - // F = m * a ==> a = F/m - // use coefficient-wise product from Eigen - // see http://eigen.tuxfamily.org/dox/group__TutorialArrayClass.html - Eigen::VectorXd newAccel(3 * m_molecule->atomCount()); - newAccel = m_forces.cwiseProduct(m_invMasses); - m_velocities += 0.5 * deltaT * (m_accel + newAccel); - m_accel = newAccel; -} - -} // namespace Avogadro diff --git a/avogadro/calc/flameminimize.h b/avogadro/calc/flameminimize.h deleted file mode 100644 index ddd65daa08..0000000000 --- a/avogadro/calc/flameminimize.h +++ /dev/null @@ -1,49 +0,0 @@ -/****************************************************************************** - This source file is part of the Avogadro project. - This source code is released under the 3-Clause BSD License, (see "LICENSE"). -******************************************************************************/ - -#ifndef AVOGADRO_CALC_FLAMEMINIMIZE_H -#define AVOGADRO_CALC_FLAMEMINIMIZE_H - -#include - -#include "energycalculator.h" - -namespace Avogadro { -namespace Core { -class Molecule; -} - -class FlameMinimize -{ -public: - FlameMinimize(QObject* parent_ = 0); - ~FlameMinimize() {} - - bool minimize(EnergyCalculator &cal, Eigen::VectorXd &positions); - - // @todo probably want a "take N steps" and a way to set convergence - // (e.g., if the forces/gradients are very small) - // ( might also check if energy changes are v. small) - - /** - * Called when the current molecule changes. - */ - void setMolecule(QtGui::Molecule* mol); - -protected: - - void verletIntegrate(Eigen::VectorXd &positions, double deltaT); - - QtGui::Molecule* m_molecule; - EnergyCalculator* m_calc; - Eigen::VectorXd m_invMasses; - Eigen::VectorXd m_forces; - Eigen::VectorXd m_velocities; - Eigen::VectorXd m_accel; -}; - -} // namespace Avogadro - -#endif // AVOGADRO_FLAMEMINIMIZE_H diff --git a/avogadro/calc/lennardjones.cpp b/avogadro/calc/lennardjones.cpp index 63e90e0c58..b2999cba2b 100644 --- a/avogadro/calc/lennardjones.cpp +++ b/avogadro/calc/lennardjones.cpp @@ -14,6 +14,10 @@ namespace Avogadro::Calc { LennardJones::LennardJones() : m_vdw(true), m_depth(100.0), m_exponent(6), m_cell(nullptr) { + // defined for 1-118 + for (unsigned int i = 1; i <= 118; ++i) { + m_elements.set(i); + } } LennardJones::~LennardJones() {} diff --git a/avogadro/calc/lennardjones.h b/avogadro/calc/lennardjones.h index 3364ff3f3e..5fb7466644 100644 --- a/avogadro/calc/lennardjones.h +++ b/avogadro/calc/lennardjones.h @@ -24,23 +24,30 @@ class AVOGADROCALC_EXPORT LennardJones : public EnergyCalculator LennardJones(); ~LennardJones(); - virtual std::string identifier() const override + std::string identifier() const override { return "LJ"; } - virtual std::string name() const override + std::string name() const override { return "Lennard-Jones"; } - virtual std::string description() const override + std::string description() const override { return "Universal Lennard-Jones potential"; } - virtual Real value(const Eigen::VectorXd& x) override; - virtual void gradient(const Eigen::VectorXd& x, + bool acceptsUnitCell() const override { return true; } + + Core::Molecule::ElementMask elements() const override + { + return (m_elements); + } + + Real value(const Eigen::VectorXd& x) override; + void gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) override; /** * Called when the current molecule changes. */ - virtual void setMolecule(Core::Molecule* mol) override; + void setMolecule(Core::Molecule* mol) override; protected: Core::Molecule* m_molecule; @@ -49,6 +56,8 @@ class AVOGADROCALC_EXPORT LennardJones : public EnergyCalculator bool m_vdw; Real m_depth; int m_exponent; + + Core::Molecule::ElementMask m_elements; }; } // namespace Calc diff --git a/avogadro/qtplugins/forcefield/CMakeLists.txt b/avogadro/qtplugins/forcefield/CMakeLists.txt index b3d670d4bf..74fd65392d 100644 --- a/avogadro/qtplugins/forcefield/CMakeLists.txt +++ b/avogadro/qtplugins/forcefield/CMakeLists.txt @@ -12,12 +12,14 @@ avogadro_plugin(Forcefield target_link_libraries(Forcefield PRIVATE Avogadro::Calc) -# Bundled forcefield scripts (open babel) +# Bundled forcefield scripts set(forcefields + scripts/gfn1.py + scripts/gfn2.py + scripts/gfnff.py + scripts/mmff94.py scripts/uff.py - scripts/gaff.py - scripts/mmff.py ) -#install(PROGRAMS ${forcefields} -#DESTINATION "${INSTALL_LIBRARY_DIR}/avogadro2/scripts/forcefields/") +install(PROGRAMS ${forcefields} +DESTINATION "${INSTALL_LIBRARY_DIR}/avogadro2/scripts/forcefields/") diff --git a/avogadro/qtplugins/forcefield/forcefielddialog.ui b/avogadro/qtplugins/forcefield/forcefielddialog.ui index b4622e8153..8c72edc764 100644 --- a/avogadro/qtplugins/forcefield/forcefielddialog.ui +++ b/avogadro/qtplugins/forcefield/forcefielddialog.ui @@ -1,12 +1,12 @@ - Avogadro::QtPlugins::OBForceFieldDialog - + Avogadro::QtPlugins::ForceFieldDialog + 0 0 - 357 + 365 295 @@ -42,11 +42,11 @@ - "Energy" convergence: + Energy convergence: - + Step limit: @@ -72,7 +72,7 @@ - + steps @@ -84,10 +84,36 @@ 100000 + 50 + + 250 + + + + + + Gradient convergence: + + + + + + + units + + + 10^ + + + -10 + + + -1 + - 2500 + -4 @@ -133,14 +159,17 @@ + + 1 + - Steepest Descent + Conjugate Gradient - Conjugate Gradient + BFGS @@ -163,7 +192,7 @@ buttonBox accepted() - Avogadro::QtPlugins::OBForceFieldDialog + Avogadro::QtPlugins::ForceFieldDialog accept() @@ -179,7 +208,7 @@ buttonBox rejected() - Avogadro::QtPlugins::OBForceFieldDialog + Avogadro::QtPlugins::ForceFieldDialog reject()