Fix openbabel forcefield integration

avogadro/qtplugins/forcefield/forcefield.cpp

@@ -31,9 +31,6 @@
 
 #include <cppoptlib/meta.h>
 #include <cppoptlib/problem.h>
-#include <cppoptlib/solver/bfgssolver.h>
-#include <cppoptlib/solver/conjugatedgradientdescentsolver.h>
-#include <cppoptlib/solver/gradientdescentsolver.h>
 #include <cppoptlib/solver/lbfgssolver.h>
 
 namespace Avogadro {
@@ -49,6 +46,7 @@ const int configureAction = 2;
 const int freezeAction = 3;
 const int unfreezeAction = 4;
 const int constraintAction = 5;
+const int forcesAction = 6;
 
 Forcefield::Forcefield(QObject* parent_)
   : ExtensionPlugin(parent_), m_method(nullptr)
@@ -63,12 +61,12 @@ Forcefield::Forcefield(QObject* parent_)
   m_gradientTolerance = settings.value("gradientTolerance", 1.0e-4).toDouble();
   settings.endGroup();
 
-  refreshScripts();
-  /* @todo - finish OBMM interface
+  // add the openbabel calculators in case they don't exist
   Calc::EnergyManager::registerModel(new OBMMEnergy("MMFF94"));
   Calc::EnergyManager::registerModel(new OBMMEnergy("UFF"));
   Calc::EnergyManager::registerModel(new OBMMEnergy("GAFF"));
-  */
+
+  refreshScripts();
 
   QAction* action = new QAction(this);
   action->setEnabled(true);
@@ -86,6 +84,14 @@ Forcefield::Forcefield(QObject* parent_)
   connect(action, SIGNAL(triggered()), SLOT(energy()));
   m_actions.push_back(action);
 
+  action = new QAction(this);
+  action->setEnabled(true);
+  action->setText(tr("Forces")); // calculate gradients
+  action->setData(forcesAction);
+  action->setProperty("menu priority", 910);
+  connect(action, SIGNAL(triggered()), SLOT(forces()));
+  m_actions.push_back(action);
+
   action = new QAction(this);
   action->setEnabled(true);
   action->setText(tr("Configure…"));
@@ -123,7 +129,11 @@ QList<QAction*> Forcefield::actions() const
 QStringList Forcefield::menuPath(QAction* action) const
 {
   QStringList path;
-  path << tr("&Extensions") << tr("&Calculate");
+  if (action->data().toInt() == optimizeAction)
+    path << tr("&Extensions");
+  else
+    path << tr("&Extensions") << tr("&Calculate");
+
   return path;
 }
 
@@ -220,7 +230,6 @@ void Forcefield::optimize()
   m_molecule->undoMolecule()->setInteractive(true);
 
   cppoptlib::LbfgsSolver<EnergyCalculator> solver;
-  // cppoptlib::ConjugatedGradientDescentSolver<EnergyCalculator> solver;
 
   int n = m_molecule->atomCount();
 
@@ -252,7 +261,7 @@ void Forcefield::optimize()
   cppoptlib::Criteria<Real> crit = cppoptlib::Criteria<Real>::defaults();
 
   // e.g., every N steps, update coordinates
-  crit.iterations = 5;
+  crit.iterations = 2;
   // we don't set function or gradient criteria
   // .. these seem to be broken in the solver code
   // .. so we handle ourselves
@@ -295,6 +304,11 @@ void Forcefield::optimize()
         forces[i] = -0.1 * Vector3(gradient[3 * i], gradient[3 * i + 1],
                                    gradient[3 * i + 2]);
       }
+    } else {
+      // reset to last positions
+      positions = lastPositions;
+      gradient = Eigen::VectorXd::Zero(3 * n);
+      break;
     }
 
     // todo - merge these into one undo step
@@ -307,18 +321,12 @@ void Forcefield::optimize()
       lastPositions = positions;
 
       // check for convergence
-      /*
       if (fabs(gradient.maxCoeff()) < m_gradientTolerance)
         break;
       if (fabs(currentEnergy - energy) < m_tolerance)
         break;
-      */
 
       energy = currentEnergy;
-    } else {
-      // reset to last positions
-      positions = lastPositions;
-      gradient = Eigen::VectorXd::Zero(3 * n);
     }
   }
 
@@ -345,6 +353,52 @@ void Forcefield::energy()
   QMessageBox::information(nullptr, tr("Avogadro"), msg);
 }
 
+void Forcefield::forces()
+{
+  if (m_molecule == nullptr || m_method == nullptr)
+    return;
+
+  int n = m_molecule->atomCount();
+
+  // double-check the mask
+  auto mask = m_molecule->frozenAtomMask();
+  if (mask.rows() != 3 * n) {
+    mask = Eigen::VectorXd::Zero(3 * n);
+    // set to 1.0
+    for (Index i = 0; i < 3 * n; ++i) {
+      mask[i] = 1.0;
+    }
+  }
+  m_method->setMolecule(m_molecule);
+  m_method->setMask(mask);
+
+  // we have to cast the current 3d positions into a VectorXd
+  Core::Array<Vector3> pos = m_molecule->atomPositions3d();
+  double* p = pos[0].data();
+  Eigen::Map<Eigen::VectorXd> map(p, 3 * n);
+  Eigen::VectorXd positions = map;
+
+  Eigen::VectorXd gradient = Eigen::VectorXd::Zero(3 * n);
+  // just to get the right size / shape
+  // we'll use this to draw the force arrows
+  Core::Array<Vector3> forces = m_molecule->atomPositions3d();
+
+  m_method->gradient(positions, gradient);
+
+  for (size_t i = 0; i < n; ++i) {
+    forces[i] =
+      -0.1 * Vector3(gradient[3 * i], gradient[3 * i + 1], gradient[3 * i + 2]);
+  }
+
+  m_molecule->setForceVectors(forces);
+  Molecule::MoleculeChanges changes = Molecule::Atoms | Molecule::Modified;
+  m_molecule->emitChanged(changes);
+
+  QString msg(
+    tr("%1 Force Norm = %L2").arg(m_methodName.c_str()).arg(gradient.norm()));
+  QMessageBox::information(nullptr, tr("Avogadro"), msg);
+}
+
 std::string Forcefield::recommendedForceField() const
 {
   // if we have a unit cell, we need to use the LJ calculator

avogadro/qtplugins/forcefield/forcefield.h

@@ -71,6 +71,7 @@ public slots:
 
 private slots:
   void energy();
+  void forces();
   void optimize();
   void freezeSelected();
   void unfreezeSelected();