Add initial obmm support

Signed-off-by: Geoff Hutchison <[email protected]>

avogadro/calc/energycalculator.cpp

@@ -20,13 +20,18 @@ 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])) {
+    if (!std::isfinite(grad[i]) || std::isnan(grad[i])) {
       grad[i] = 0.0;
     }
   }
 
   // freeze any masked atoms or coordinates
-  grad = grad.cwiseProduct(m_mask);
+  /*
+  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

avogadro/qtgui/pythonscript.cpp

@@ -221,37 +221,32 @@ void PythonScript::processFinished(int, QProcess::ExitStatus)
   emit finished();
 }
 
-QByteArray PythonScript::asyncWriteAndResponse(QByteArray input, const int expectedLines)
+void PythonScript::asyncTerminate()
+{
+  if (m_process != nullptr) {
+    m_process->terminate();
+    disconnect(m_process, SIGNAL(finished()), this, SLOT(processFinished()));
+    m_process->deleteLater();
+    m_process = nullptr;
+  }
+}
+
+QByteArray PythonScript::asyncWriteAndResponse(QByteArray input)
 {
   if (m_process == nullptr) {
     return QByteArray(); // wait
   }
 
-  qDebug() << " writing to process: " << input << " " << expectedLines;
-
   m_process->write(input);
   QByteArray buffer;
 
-  if (expectedLines == -1) {
-    m_process->waitForFinished();
-    buffer = m_process->readAll();
-  } else {
-    // only read a certain number of lines
-    // (e.g., energy and gradients)
-    int remainingLines = expectedLines;
-    bool ready = m_process->waitForReadyRead();
-    if (ready) {
-      while (m_process->canReadLine() && remainingLines > 0) {
-        buffer += m_process->readLine();
-        remainingLines--;
-      }
-      // clear the rest of the buffer
-      qDebug() << " remaining " << m_process->readAll();
+  bool ready = m_process->waitForReadyRead();
+  if (ready) {
+    while (m_process->canReadLine()) {
+      buffer += m_process->readLine();
     }
   }
 
-  qDebug() << " read from process: " << buffer;
-
   return buffer;
 }
 

avogadro/qtgui/pythonscript.h

@@ -12,9 +12,9 @@
 #include <avogadro/core/avogadrocore.h>
 
 #include <QtCore/QByteArray>
+#include <QtCore/QProcess>
 #include <QtCore/QString>
 #include <QtCore/QStringList>
-#include <QtCore/QProcess>
 
 namespace Avogadro {
 namespace QtGui {
@@ -92,33 +92,36 @@ class AVOGADROQTGUI_EXPORT PythonScript : public QObject
    * Start a new process to execute asynchronously
    * "<m_pythonInterpreter> <scriptFilePath()> [args ...]",
    * optionally passing scriptStdin to the processes standard input.
-   * 
+   *
    * Will send asyncFinished() signal when finished
    */
   void asyncExecute(const QStringList& args,
-                          const QByteArray& scriptStdin = QByteArray());
+                    const QByteArray& scriptStdin = QByteArray());
 
   /**
    * Write input to the asynchronous process' standard input and return the
    * standard output when ready. Does not wait for the process to terminate
    * before returning (e.g. "server mode").
-   * 
+   *
    * @param input The input to write to the process' standard input
-   * @param expectedLines The number of lines to expect in the output. If -1,
-   *                     the output is read until the process terminates.
    * @return The standard output of the process
+   */
+  QByteArray asyncWriteAndResponse(QByteArray input);
+
+  /**
+   * Terminate the asynchronous process.
   */
-  QByteArray asyncWriteAndResponse(QByteArray input, const int expectedLines = -1);
+  void asyncTerminate();
 
   /**
    * Returns the standard output of the asynchronous process when finished.
    */
   QByteArray asyncResponse();
 
 signals:
-/**
- * The asynchronous execution is finished or timed out
- */
+  /**
+   * The asynchronous execution is finished or timed out
+   */
   void finished();
 
 public slots:

avogadro/qtplugins/forcefield/CMakeLists.txt

@@ -1,5 +1,6 @@
 set(forcefield_srcs
   forcefield.cpp
+  obmmenergy.cpp
   scriptenergy.cpp
 )
 
@@ -18,8 +19,6 @@ set(forcefields
   scripts/gfn1.py
   scripts/gfn2.py
   scripts/gfnff.py
-  scripts/mmff94.py
-  scripts/uff.py
 )
 
 install(PROGRAMS ${forcefields}

avogadro/qtplugins/forcefield/forcefield.cpp

@@ -5,6 +5,7 @@
 ******************************************************************************/
 
 #include "forcefield.h"
+#include "obmmenergy.h"
 #include "scriptenergy.h"
 
 #include <QtCore/QDebug>
@@ -50,7 +51,10 @@ const int constraintAction = 5;
 
 Forcefield::Forcefield(QObject* parent_) : ExtensionPlugin(parent_)
 {
-  refreshScripts();
+  //refreshScripts();
+  Calc::EnergyManager::registerModel(new OBMMEnergy("MMFF94"));
+  Calc::EnergyManager::registerModel(new OBMMEnergy("UFF"));
+  Calc::EnergyManager::registerModel(new OBMMEnergy("GAFF"));
 
   QAction* action = new QAction(this);
   action->setEnabled(true);
@@ -119,14 +123,16 @@ void Forcefield::optimize()
   bool isInteractive = m_molecule->undoMolecule()->isInteractive();
   m_molecule->undoMolecule()->setInteractive(true);
 
-  cppoptlib::LbfgsSolver<EnergyCalculator> solver;
+  //cppoptlib::LbfgsSolver<EnergyCalculator> solver;
+  cppoptlib::ConjugatedGradientDescentSolver<EnergyCalculator> solver;
 
   int n = m_molecule->atomCount();
   // 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 lastPositions = positions;
 
   Eigen::VectorXd gradient = Eigen::VectorXd::Zero(3 * n);
   // just to get the right size / shape
@@ -137,54 +143,93 @@ void Forcefield::optimize()
   cppoptlib::Criteria<Real> crit = cppoptlib::Criteria<Real>::defaults();
 
   // e.g., every 5 steps, update coordinates
-  crit.iterations = m_nSteps;
+  crit.iterations = 5;
   // we don't set function or gradient criteria
   // .. these seem to be broken in the solver code
   // .. so we handle ourselves
   solver.setStopCriteria(crit);
 
   // set the method
-  std::string recommended = recommendedForceField();
+  std::string recommended = "UFF";
   qDebug() << "Energy method: " << recommended.c_str();
 
   if (m_method == nullptr) {
     // we have to create the calculator
     m_method = Calc::EnergyManager::instance().model(recommended);
   }
   m_method->setMolecule(m_molecule);
-  m_method->setMask(m_molecule->frozenAtomMask());
+  auto mask = m_molecule->frozenAtomMask();
+  if (mask.rows() != m_molecule->atomCount()) {
+    mask = Eigen::VectorXd::Zero(3 * m_molecule->atomCount());
+    // set to 1.0
+    for (Index i = 0; i < 3 * m_molecule->atomCount(); ++i) {
+      mask[i] = 1.0;
+    }
+  }
+  m_method->setMask(mask);
 
   Real energy = m_method->value(positions);
+  m_method->gradient(positions, gradient);
+  qDebug() << " initial " << energy 
+  << " gradNorm: " << gradient.norm()
+  << " posNorm: " << positions.norm();
+
+  Real currentEnergy = 0.0;
   for (unsigned int i = 0; i < m_maxSteps / crit.iterations; ++i) {
     solver.minimize(*m_method, positions);
 
-    Real currentEnergy = m_method->value(positions);
+    currentEnergy = m_method->value(positions);
     // get the current gradient for force visualization
     m_method->gradient(positions, gradient);
+    qDebug() << " optimize " << i << currentEnergy
+             << " gradNorm: " << gradient.norm()
+             << " posNorm: " << positions.norm();
 
     // update coordinates
-    const double* d = positions.data();
-    // casting back would be lovely...
-    for (size_t i = 0; i < n; ++i) {
-      pos[i] = Vector3(*(d), *(d + 1), *(d + 2));
-      d += 3;
-
-      forces[i] = -1.0 * Vector3(gradient[3 * i], gradient[3 * i + 1],
-                                 gradient[3 * i + 2]);
+    bool isFinite = std::isfinite(currentEnergy);
+    if (isFinite) {
+      const double* d = positions.data();
+      bool isFinite = true;
+      // casting back would be lovely...
+      for (size_t i = 0; i < n; ++i) {
+        if (!std::isfinite(*d) || !std::isfinite(*(d + 1)) ||
+            !std::isfinite(*(d + 2))) {
+          isFinite = false;
+          break;
+        }
+
+        pos[i] = Vector3(*(d), *(d + 1), *(d + 2));
+        d += 3;
+
+        forces[i] = -1.0 * Vector3(gradient[3 * i], gradient[3 * i + 1],
+                                   gradient[3 * i + 2]);
+      }
     }
+
     // todo - merge these into one undo step
-    m_molecule->undoMolecule()->setAtomPositions3d(pos,
-                                                   tr("Optimize Geometry"));
-    m_molecule->setForceVectors(forces);
-    Molecule::MoleculeChanges changes = Molecule::Atoms | Molecule::Modified;
-    m_molecule->emitChanged(changes);
-
-    // check for convergence
-    if (fabs(gradient.maxCoeff()) < m_gradientTolerance)
-      break;
-    if (fabs(currentEnergy - energy) < m_tolerance)
-      break;
-    energy = currentEnergy;
+    if (isFinite) {
+      qDebug() << " finite! ";
+      m_molecule->undoMolecule()->setAtomPositions3d(pos,
+                                                     tr("Optimize Geometry"));
+      m_molecule->setForceVectors(forces);
+      Molecule::MoleculeChanges changes = Molecule::Atoms | Molecule::Modified;
+      m_molecule->emitChanged(changes);
+      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);
+    }
   }
 
   m_molecule->undoMolecule()->setInteractive(isInteractive);

avogadro/qtplugins/forcefield/forcefield.h

@@ -80,7 +80,7 @@ private slots:
 
   // defaults
   Minimizer m_minimizer = LBFGS;
-  unsigned int m_maxSteps = 250;
+  unsigned int m_maxSteps = 100;
   unsigned int m_nSteps = 5;
   double m_tolerance = 1.0e-6;
   double m_gradientTolerance = 1.0e-4;