Implement a script interface - currently just energies

Also implement frozen atoms / coordinates in the molecule class

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

avogadro/calc/energycalculator.cpp

@@ -27,22 +27,4 @@ void EnergyCalculator::cleanGradients(TVector& grad)
   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

@@ -27,6 +27,12 @@ class AVOGADROCALC_EXPORT EnergyCalculator : public cppoptlib::Problem<Real>
   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.
    */
@@ -96,16 +102,23 @@ class AVOGADROCALC_EXPORT EnergyCalculator : public cppoptlib::Problem<Real>
    */
   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:
+  /**
+   * @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

avogadro/calc/energymanager.cpp

@@ -81,7 +81,11 @@ std::string EnergyManager::nameForModel(const std::string& identifier) const
 
 EnergyManager::EnergyManager()
 {
-  // add any default models here (EEM maybe?)
+  // add any default models here
+
+  // LJ is the fallback, since it can handle anything
+  // (maybe not well, but it can handle it)
+  addModel(new LennardJones);
 }
 
 EnergyManager::~EnergyManager()

avogadro/calc/lennardjones.h

@@ -14,7 +14,7 @@ namespace Avogadro {
 namespace Core {
 class Molecule;
 class UnitCell;
-}
+} // namespace Core
 
 namespace Calc {
 
@@ -24,25 +24,23 @@ class AVOGADROCALC_EXPORT LennardJones : public EnergyCalculator
   LennardJones();
   ~LennardJones();
 
-  std::string identifier() const override
-  { return "LJ"; }
+  LennardJones* newInstance() const override { return new LennardJones; }
 
-  std::string name() const override
-  { return "Lennard-Jones"; }
+  std::string identifier() const override { return "LJ"; }
+
+  std::string name() const override { return "Lennard-Jones"; }
 
   std::string description() const override
-  { return "Universal Lennard-Jones potential"; }
+  {
+    return "Universal Lennard-Jones potential";
+  }
 
   bool acceptsUnitCell() const override { return true; }
 
-  Core::Molecule::ElementMask elements() const override
-  {
-    return (m_elements);
-  }
+  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;
+  void gradient(const Eigen::VectorXd& x, Eigen::VectorXd& grad) override;
 
   /**
    * Called when the current molecule changes.

avogadro/core/molecule.cpp

@@ -49,6 +49,7 @@ Molecule::Molecule(const Molecule& other)
     m_residues(other.m_residues), m_hallNumber(other.m_hallNumber),
     m_graph(other.m_graph), m_bondOrders(other.m_bondOrders),
     m_atomicNumbers(other.m_atomicNumbers),
+    m_frozenAtomMask(other.m_frozenAtomMask),
     m_layers(LayerManager::getMoleculeLayer(this))
 {
   // Copy over any meshes
@@ -90,6 +91,7 @@ Molecule::Molecule(Molecule&& other) noexcept
     m_residues(other.m_residues), m_hallNumber(other.m_hallNumber),
     m_graph(other.m_graph), m_bondOrders(other.m_bondOrders),
     m_atomicNumbers(other.m_atomicNumbers),
+    m_frozenAtomMask(other.m_frozenAtomMask),
     m_layers(LayerManager::getMoleculeLayer(this))
 {
   m_basisSet = other.m_basisSet;
@@ -133,6 +135,7 @@ Molecule& Molecule::operator=(const Molecule& other)
     m_bondOrders = other.m_bondOrders;
     m_atomicNumbers = other.m_atomicNumbers;
     m_hallNumber = other.m_hallNumber;
+    m_frozenAtomMask = other.m_frozenAtomMask;
 
     clearMeshes();
 
@@ -193,6 +196,7 @@ Molecule& Molecule::operator=(Molecule&& other) noexcept
     m_bondOrders = other.m_bondOrders;
     m_atomicNumbers = other.m_atomicNumbers;
     m_hallNumber = other.m_hallNumber;
+    m_frozenAtomMask = other.m_frozenAtomMask;
 
     clearMeshes();
     m_meshes = std::move(other.m_meshes);
@@ -267,6 +271,59 @@ std::set<std::string> Molecule::partialChargeTypes() const
   return types;
 }
 
+void Molecule::setFrozenAtom(Index atomId, bool frozen)
+{
+  if (atomId >= m_atomicNumbers.size())
+    return;
+
+  // check if we need to resize
+  unsigned int size = m_frozenAtomMask.rows();  
+  if (m_frozenAtomMask.rows() != 3*m_atomicNumbers.size())
+    m_frozenAtomMask.conservativeResize(3*m_atomicNumbers.size());
+
+  // do we need to initialize new values?
+  if (m_frozenAtomMask.rows() > size)
+    for (unsigned int i = size; i < m_frozenAtomMask.rows(); ++i)
+      m_frozenAtomMask[i] = 1.0;
+
+  float value = frozen ? 0.0 : 1.0;
+  if (atomId * 3 <= m_frozenAtomMask.rows() - 3) {
+    m_frozenAtomMask[atomId*3] = value;
+    m_frozenAtomMask[atomId*3+1] = value;
+    m_frozenAtomMask[atomId*3+2] = value;
+  }
+}
+
+bool Molecule::frozenAtom(Index atomId) const
+{
+  bool frozen = false;
+  if (atomId * 3 <= m_frozenAtomMask.rows() - 3) {
+    frozen = (m_frozenAtomMask[atomId*3] == 0.0 &&
+              m_frozenAtomMask[atomId*3+1] == 0.0 &&
+              m_frozenAtomMask[atomId*3+2] == 0.0);
+  }
+  return frozen;
+}
+
+void Molecule::setFrozenAtomAxis(Index atomId, int axis, bool frozen)
+{
+  // check if we need to resize
+  unsigned int size = m_frozenAtomMask.rows();  
+  if (m_frozenAtomMask.rows() != 3*m_atomicNumbers.size())
+    m_frozenAtomMask.conservativeResize(3*m_atomicNumbers.size());
+
+  // do we need to initialize new values?
+  if (m_frozenAtomMask.rows() > size)
+    for (unsigned int i = size; i < m_frozenAtomMask.rows(); ++i)
+      m_frozenAtomMask[i] = 1.0;
+
+  float value = frozen ? 0.0 : 1.0;
+  if (atomId * 3 <= m_frozenAtomMask.rows() - 3) {
+    m_frozenAtomMask[atomId*3 + axis] = value;
+  }
+}
+
+
 void Molecule::setData(const std::string& name, const Variant& value)
 {
   m_data.setValue(name, value);

avogadro/core/molecule.h

@@ -693,6 +693,26 @@ class AVOGADROCORE_EXPORT Molecule
    */
   bool setAtomicNumber(Index atomId, unsigned char atomicNumber);
 
+  /**
+   * Freeze or unfreeze an atom for optimization
+  */
+  void setFrozenAtom(Index atomId, bool frozen);
+
+  /**
+   * Get the frozen status of an atom
+  */
+  bool frozenAtom(Index atomId) const;
+
+  /**
+   * Freeze or unfreeze X, Y, or Z coordinate of an atom for optimization
+   * @param atomId The index of the atom to modify.
+   * @param axis The axis to freeze (0, 1, or 2 for X, Y, or Z)
+   * @param frozen True to freeze, false to unfreeze
+  */
+  void setFrozenAtomAxis(Index atomId, int axis, bool frozen);
+
+  Eigen::VectorXd frozenAtomMask() const { return m_frozenAtomMask; }
+
   /**
    * @return a map of components and count.
    */
@@ -763,6 +783,8 @@ class AVOGADROCORE_EXPORT Molecule
   // This will be stored from the last space group operation
   unsigned short m_hallNumber = 0;
 
+  Eigen::VectorXd m_frozenAtomMask;
+
 private:
   mutable Graph m_graph; // A transformation of the molecule to a graph.
   // edge information

avogadro/qtgui/pythonscript.cpp

@@ -221,6 +221,34 @@ void PythonScript::processFinished(int, QProcess::ExitStatus)
   emit finished();
 }
 
+QByteArray PythonScript::asyncWriteAndResponse(QByteArray input, const int expectedLines)
+{
+  if (m_process == nullptr) {
+    return QByteArray(); // wait
+  }
+
+  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--;
+      }
+    }
+  }
+
+  return buffer;
+}
+
 QByteArray PythonScript::asyncResponse()
 {
   if (m_process == nullptr || m_process->state() == QProcess::Running) {

avogadro/qtgui/pythonscript.h

@@ -98,6 +98,18 @@ class AVOGADROQTGUI_EXPORT PythonScript : public QObject
   void asyncExecute(const QStringList& args,
                           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, const int expectedLines = -1);
+
   /**
    * Returns the standard output of the asynchronous process when finished.
    */

avogadro/qtplugins/forcefield/CMakeLists.txt

@@ -1,5 +1,6 @@
 set(forcefield_srcs
   forcefield.cpp
+  scriptenergy.cpp
 )
 
 avogadro_plugin(Forcefield
@@ -22,4 +23,4 @@ set(forcefields
 )
 
 install(PROGRAMS ${forcefields}
-DESTINATION "${INSTALL_LIBRARY_DIR}/avogadro2/scripts/forcefields/")
+DESTINATION "${INSTALL_LIBRARY_DIR}/avogadro2/scripts/energy/")