New forcefield framework

avogadro/calc/CMakeLists.txt

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

avogadro/calc/energycalculator.cpp

@@ -0,0 +1,48 @@
+/******************************************************************************
+  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"
+
+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])) {
+      grad[i] = 0.0;
+    }
+  }
+
+  // freeze any masked atoms or coordinates
+  grad = grad.cwiseProduct(m_mask);
+}
+
+void EnergyCalculator::freezeAtom(Index atomId)
+{
+  if (atomId * 3 <= m_mask.rows() - 3) {
+    m_mask[atomId*3] = 0.0;
+    m_mask[atomId*3+1] = 0.0;
+    m_mask[atomId*3+2] = 0.0;
+  }
+}
+
+void EnergyCalculator::unfreezeAtom(Index atomId)
+{
+  if (atomId * 3 <= m_mask.rows() - 3) {
+    m_mask[atomId*3] = 1.0;
+    m_mask[atomId*3+1] = 1.0;
+    m_mask[atomId*3+2] = 1.0;
+  }
+}
+
+} // namespace Avogadro

avogadro/calc/energycalculator.h

@@ -0,0 +1,85 @@
+/******************************************************************************
+  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/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() {}
+
+  /**
+   * @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; }
+
+  /**
+   * 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; }
+
+  void freezeAtom(Index atomId);
+  void unfreezeAtom(Index atomId);
+
+  /**
+   * Called when the current molecule changes.
+   */
+  virtual void setMolecule(Core::Molecule* mol) = 0;
+
+protected:
+  TVector m_mask; // optimize or frozen atom mask
+};
+
+} // end namespace Calc
+} // end namespace Avogadro
+
+#endif // AVOGADRO_CALC_ENERGYCALCULATOR_H

avogadro/calc/flameminimize.cpp

@@ -0,0 +1,113 @@
+/******************************************************************************
+  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 <avogadro/core/elements.h>
+
+#include <cmath>
+
+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

avogadro/calc/flameminimize.h

@@ -0,0 +1,49 @@
+/******************************************************************************
+  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 <avogadro/core/vector.h>
+
+#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