New forcefield framework

avogadro/calc/CMakeLists.txt

@@ -4,15 +4,21 @@ avogadro_headers(Calc
   chargemodel.h
   chargemanager.h
   defaultmodel.h
+  energycalculator.h
+  energymanager.h
+  lennardjones.h
 )
 
 target_sources(Calc PRIVATE
   chargemodel.cpp
   chargemanager.cpp
   defaultmodel.cpp
+  energycalculator.cpp
+  energymanager.cpp
+  lennardjones.cpp
 )
 
 avogadro_add_library(Calc)
 
 target_link_libraries(Calc
-  PUBLIC Avogadro::Core)
+  PUBLIC Avogadro::Core cppoptlib)

avogadro/calc/chargemanager.cpp

@@ -10,7 +10,6 @@
 #include <algorithm>
 #include <memory>
 
-
 namespace Avogadro::Calc {
 
 ChargeManager& ChargeManager::instance()
@@ -88,7 +87,7 @@ ChargeManager::ChargeManager()
 ChargeManager::~ChargeManager()
 {
   // Delete the models that were loaded.
-  for (auto & m_model : m_models) {
+  for (auto& m_model : m_models) {
     delete m_model;
   }
   m_models.clear();
@@ -102,7 +101,6 @@ std::set<std::string> ChargeManager::identifiersForMolecule(
 
   // check our models for compatibility
   for (auto m_model : m_models) {
-
     // We check that every element in the molecule
     // is handled by the model
     auto mask = m_model->elements() & molecule.elements();
@@ -113,8 +111,8 @@ std::set<std::string> ChargeManager::identifiersForMolecule(
   return identifiers;
 }
 
-MatrixX ChargeManager::partialCharges(
-  const std::string& identifier, Core::Molecule& molecule) const
+MatrixX ChargeManager::partialCharges(const std::string& identifier,
+                                      Core::Molecule& molecule) const
 {
   // first check if the type is found in the molecule
   // (i.e., read from a file not computed dynamically)
@@ -181,4 +179,4 @@ Core::Array<double> ChargeManager::potentials(
   return model->potentials(molecule, points);
 }
 
-} // namespace Avogadro
+} // namespace Avogadro::Calc

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.
    */
@@ -80,16 +80,6 @@ class AVOGADROCALC_EXPORT ChargeManager
    */
   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 format.
-   * @return Instance of the format, nullptr if not found. Ownership passes to
-   * the
-   * caller.
-   */
-  ChargeModel* newModelFromIdentifier(const std::string& identifier) const;
-
   /**
    * Get a list of all loaded identifiers
    */
@@ -111,7 +101,7 @@ class AVOGADROCALC_EXPORT ChargeManager
    * @return atomic partial charges for the molecule, or 0.0 if undefined
    */
   MatrixX partialCharges(const std::string& identifier,
-                               Core::Molecule& mol) const;
+                         Core::Molecule& mol) const;
 
   /**
    * @return the potential at the point for the molecule, or 0.0 if the model is

avogadro/calc/chargemodel.cpp

@@ -39,7 +39,7 @@ double ChargeModel::potential(Molecule& mol, const Vector3& point) const
     double distance = (positions[i] - point).norm();
     if (distance > 0.01) {
       // drop small distances to avoid overflow
-      potential += charges(i,0) / distance;
+      potential += charges(i, 0) / distance;
     }
   }
 
@@ -52,12 +52,13 @@ Array<double> ChargeModel::potentials(Core::Molecule& mol,
   // This is naive and slow, but can be re-implemented by methods
   // for batching
   Array<double> potentials(points.size(), 0.0);
-  for(unsigned int i = 0; i < points.size(); ++i)
+  for (unsigned int i = 0; i < points.size(); ++i)
     potentials[i] = potential(mol, points[i]);
   return potentials;
 }
 
-void ChargeModel::appendError(const std::string& errorString, bool newLine) const
+void ChargeModel::appendError(const std::string& errorString,
+                              bool newLine) const
 {
   m_error += errorString;
   if (newLine)

avogadro/calc/energycalculator.cpp

@@ -0,0 +1,36 @@
+/******************************************************************************
+  This source file is part of the Avogadro project.
+  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+******************************************************************************/
+
+#include "energycalculator.h"
+
+#include <iostream>
+
+namespace Avogadro::Calc {
+
+void EnergyCalculator::gradient(const TVector& x, TVector& grad)
+{
+  finiteGradient(x, grad);
+  cleanGradients(grad);
+}
+
+void EnergyCalculator::cleanGradients(TVector& grad)
+{
+  unsigned int size = grad.rows();
+  // check for overflows -- in case of divide by zero, etc.
+  for (unsigned int i = 0; i < size; ++i) {
+    if (!std::isfinite(grad[i]) || std::isnan(grad[i])) {
+      grad[i] = 0.0;
+    }
+  }
+
+  // freeze any masked atoms or coordinates
+  if (m_mask.rows() == size)
+    grad = grad.cwiseProduct(m_mask);
+  else
+    std::cerr << "Error: mask size " << m_mask.rows() << " " << grad.rows()
+              << std::endl;
+}
+
+} // namespace Avogadro::Calc

avogadro/calc/energycalculator.h

@@ -0,0 +1,127 @@
+/******************************************************************************
+  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_ENERGYCALCULATOR_H
+#define AVOGADRO_CALC_ENERGYCALCULATOR_H
+
+#include "avogadrocalcexport.h"
+
+#include <avogadro/core/molecule.h>
+#include <avogadro/core/variantmap.h>
+#include <avogadro/core/vector.h>
+
+#include <cppoptlib/problem.h>
+
+namespace Avogadro {
+namespace Core {
+class Molecule;
+}
+
+namespace Calc {
+
+class AVOGADROCALC_EXPORT EnergyCalculator : public cppoptlib::Problem<Real>
+{
+public:
+  EnergyCalculator() {}
+  ~EnergyCalculator() {}
+
+  /**
+   * Create a new instance of the model. Ownership passes to the
+   * caller.
+   */
+  virtual EnergyCalculator* newInstance() const = 0;
+
+  /**
+   * @return a unique identifier for this calculator.
+   */
+  virtual std::string identifier() const = 0;
+
+  /**
+   * @return A short translatable name for this method (e.g., MMFF94, UFF, etc.)
+   */
+  virtual std::string name() const = 0;
+
+  /**
+   * @return a description of the method
+   */
+  virtual std::string description() const = 0;
+
+  /**
+   * Called to set the configuration (e.g., for a GUI options dialog)
+   */
+  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
+   */
+  virtual void gradient(const TVector& x, TVector& grad) override;
+
+  /**
+   * Called to 'clean' gradients @param grad (e.g., for constraints)
+   */
+  void cleanGradients(TVector& grad);
+
+  /**
+   * Called to update the "frozen" mask (e.g., during editing)
+   */
+  void setMask(TVector mask) { m_mask = mask; }
+
+  /**
+   * @return the frozen atoms mask
+   */
+  TVector mask() const { return m_mask; }
+
+  /**
+   * Called when the current molecule changes.
+   */
+  virtual void setMolecule(Core::Molecule* mol) = 0;
+
+protected:
+  /**
+   * @brief Append an error to the error string for the model.
+   * @param errorString The error to be added.
+   * @param newLine Add a new line after the error string?
+   */
+  void appendError(const std::string& errorString, bool newLine = true) const;
+
+  TVector m_mask; // optimize or frozen atom mask
+
+private:
+  mutable std::string m_error;
+};
+
+} // end namespace Calc
+} // end namespace Avogadro
+
+#endif // AVOGADRO_CALC_ENERGYCALCULATOR_H