Continue working on scripts

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

avogadro/qtplugins/forcefield/forcefield.cpp

@@ -51,10 +51,12 @@ 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);
@@ -124,7 +126,8 @@ void Forcefield::optimize()
   m_molecule->undoMolecule()->setInteractive(true);
 
   //cppoptlib::LbfgsSolver<EnergyCalculator> solver;
-  cppoptlib::ConjugatedGradientDescentSolver<EnergyCalculator> solver;
+  //cppoptlib::ConjugatedGradientDescentSolver<EnergyCalculator> solver;
+  cppoptlib::GradientDescentSolver<EnergyCalculator> solver;
 
   int n = m_molecule->atomCount();
   // we have to cast the current 3d positions into a VectorXd
@@ -142,15 +145,15 @@ void Forcefield::optimize()
   // Create a Criteria class so we can update coords every N steps
   cppoptlib::Criteria<Real> crit = cppoptlib::Criteria<Real>::defaults();
 
-  // e.g., every 5 steps, update coordinates
-  crit.iterations = 5;
+  // e.g., every N steps, update coordinates
+  crit.iterations = 10;
   // 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 = "UFF";
+  std::string recommended = "MMFF94";
   qDebug() << "Energy method: " << recommended.c_str();
 
   if (m_method == nullptr) {
@@ -208,7 +211,6 @@ void Forcefield::optimize()
 
     // todo - merge these into one undo step
     if (isFinite) {
-      qDebug() << " finite! ";
       m_molecule->undoMolecule()->setAtomPositions3d(pos,
                                                      tr("Optimize Geometry"));
       m_molecule->setForceVectors(forces);

avogadro/qtplugins/forcefield/obmmenergy.cpp

@@ -22,16 +22,18 @@ class OBMMEnergy::ProcessListener : public QObject
 public:
   ProcessListener(QProcess *proc) : QObject(), m_finished(false), m_process(proc)
   {
-    connect(m_process, SIGNAL(readyRead()), SLOT(readyRead()));
   }
 
   bool waitForOutput(QByteArray output, int msTimeout = 2000)
   {
+    connect(m_process, SIGNAL(readyRead()), SLOT(readyRead()));
     if (!wait(msTimeout))
       return false;
 
     // success!
     output = m_output;
+    disconnect(m_process, nullptr, nullptr, nullptr);
+    m_finished = false;
     return true;
   }
 
@@ -209,6 +211,7 @@ void OBMMEnergy::setMolecule(Core::Molecule* mol)
     // appendError("Error running process.");
     return;
   }
+  qDebug() << "OBMM start: " << result;
 
   // okay, we need to write "load <filename>" to the interpreter
   // and then read the response
@@ -233,9 +236,9 @@ void OBMMEnergy::setMolecule(Core::Molecule* mol)
   input = QByteArray("energy\n");
   result.clear();
   m_process->write(input);
-  if (!listener.waitForOutput(result)) {
-    // appendError("Error running process.");
-    return;
+  result.clear();
+  while (!result.contains("command >")) {
+    result += m_process->readLine();
   }
   qDebug() << "OBMM energy: " << result;
 }
@@ -262,7 +265,6 @@ Real OBMMEnergy::value(const Eigen::VectorXd& x)
   input += "\n";
 
   m_process->write(input);
-  m_process->waitForBytesWritten();
   m_process->waitForReadyRead();
   result.clear();
   while (!result.contains("command >")) {
@@ -275,7 +277,6 @@ Real OBMMEnergy::value(const Eigen::VectorXd& x)
   input = "energy\n\n";
   result.clear();
   m_process->write(input);
-  m_process->waitForBytesWritten();
   m_process->waitForReadyRead();
   while (!result.contains("command >")) {
     result += m_process->readLine();

avogadro/qtplugins/forcefield/scripts/gaff.py

@@ -0,0 +1,88 @@
+#  This source file is part of the Avogadro project.
+#  This source code is released under the 3-Clause BSD License, (see "LICENSE").
+
+import argparse
+import json
+import sys
+
+try:
+    from openbabel import pybel
+    import numpy as np
+
+    imported = True
+except ImportError:
+    imported = False
+
+
+def getMetaData():
+    # before we return metadata, make sure xtb is in the path
+    if not imported:
+        return {}  # Avogadro will ignore us now
+
+    metaData = {
+        "name": "GAFF",
+        "identifier": "GAFF",
+        "description": "Calculate GAFF energies and gradients",
+        "inputFormat": "cml",
+        "elements": "1,6-9,14-17,35,53",
+        "unitCell": False,
+        "gradients": True,
+        "ion": False,
+        "radical": False,
+    }
+    return metaData
+
+
+def run(filename):
+    # we get the molecule from the supplied filename
+    #  in cjson format (it's a temporary file created by Avogadro)
+    mol = next(pybel.readfile("cml", filename))
+
+    ff = pybel._forcefields["gaff"]
+    success = ff.Setup(mol.OBMol)
+    if not success:
+        # should never happen, but just in case
+        sys.exit("GAFF setup failed")
+
+    # we loop forever - Avogadro will kill the process when done
+    num_atoms = len(mol.atoms)
+    while True:
+        # read new coordinates from stdin
+        for i in range(num_atoms):
+            coordinates = np.fromstring(input(), sep=" ")
+            atom = mol.atoms[i]
+            atom.OBAtom.SetVector(coordinates[0], coordinates[1], coordinates[2])
+
+        # update the molecule geometry for the next energy
+        ff.SetCoordinates(mol.OBMol)
+
+        # first print the energy of these coordinates
+        energy = ff.Energy(True)  # in kJ/mol
+        print("AvogadroEnergy:", energy)  # in kJ/mol
+
+        # now print the gradient on each atom
+        print("AvogadroGradient:")
+        for atom in mol.atoms:
+            grad = ff.GetGradient(atom.OBAtom)
+            print(grad.GetX(), grad.GetY(), grad.GetZ())
+
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser("GAFF calculator")
+    parser.add_argument("--display-name", action="store_true")
+    parser.add_argument("--metadata", action="store_true")
+    parser.add_argument("-f", "--file", nargs=1)
+    parser.add_argument("--lang", nargs="?", default="en")
+    args = vars(parser.parse_args())
+
+    if args["metadata"]:
+        print(json.dumps(getMetaData()))
+    elif args["display_name"]:
+        name = getMetaData().get("name")
+        if name:
+            print(name)
+        else:
+            sys.exit("pybel is unavailable")
+    elif args["file"]:
+        run(args["file"][0])

avogadro/qtplugins/forcefield/scripts/gfn1.py

@@ -62,25 +62,27 @@ def run(filename):
 
     # we loop forever - Avogadro will kill the process when done
     while(True):
-        # first print the energy of these coordinates
-        print(res.get_energy()) # in Hartree
-
-        # now print the gradient 
-        # .. we don't want the "[]" in the output
-        output = np.array2string(res.get_gradient())
-        output = output.replace("[", "").replace("]", "")
-        print(output)
-
         # read new coordinates from stdin
         for i in range(len(atoms)):
             coordinates[i] = np.fromstring(input(), sep=" ")
         # .. convert from Angstrom to Bohr
         coordinates /= 0.52917721067
-        
+
         # update the calculator and run a new calculation
         calc.update(coordinates)
         calc.singlepoint(res)
 
+        # first print the energy of these coordinates
+        print("AvogadroEnergy:", res.get_energy())  # in Hartree
+
+        # now print the gradient
+        # .. we don't want the "[]" in the output
+        print("AvogadroGradient:")
+        grad = res.get_gradient() * 4961.475  # convert units
+        output = np.array2string(grad)
+        output = output.replace("[", "").replace("]", "")
+        print(output)
+
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser("GFN1 calculator")

avogadro/qtplugins/forcefield/scripts/gfn2.py

@@ -62,25 +62,27 @@ def run(filename):
 
     # we loop forever - Avogadro will kill the process when done
     while(True):
-        # first print the energy of these coordinates
-        print(res.get_energy()) # in Hartree
-
-        # now print the gradient 
-        # .. we don't want the "[]" in the output
-        output = np.array2string(res.get_gradient())
-        output = output.replace("[", "").replace("]", "")
-        print(output)
-
         # read new coordinates from stdin
         for i in range(len(atoms)):
             coordinates[i] = np.fromstring(input(), sep=" ")
         # .. convert from Angstrom to Bohr
         coordinates /= 0.52917721067
-        
+
         # update the calculator and run a new calculation
         calc.update(coordinates)
         calc.singlepoint(res)
 
+        # first print the energy of these coordinates
+        print("AvogadroEnergy:", res.get_energy())  # in Hartree
+
+        # now print the gradient
+        # .. we don't want the "[]" in the output
+        print("AvogadroGradient:")
+        grad = res.get_gradient() * 4961.475  # convert units
+        output = np.array2string(grad)
+        output = output.replace("[", "").replace("]", "")
+        print(output)
+
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser("GFN2 calculator")

avogadro/qtplugins/forcefield/scripts/gfnff.py

@@ -80,13 +80,10 @@ def run(filename):
         # first print the energy of these coordinates
         print("AvogadroEnergy:", res.get_energy())  # in Hartree
 
-        with open("/Users/ghutchis/gfnff.log", "a") as f:
-            f.write(str(res.get_energy()) + "\n")
-
         # now print the gradient
         # .. we don't want the "[]" in the output
         print("AvogadroGradient:")
-        grad = res.get_gradient() * 4961.475
+        grad = res.get_gradient() * 4961.475  # convert units
         output = np.array2string(grad)
         output = output.replace("[", "").replace("]", "")
         print(output)

avogadro/qtplugins/forcefield/scripts/mmff94.py

@@ -47,14 +47,6 @@ def run(filename):
     # we loop forever - Avogadro will kill the process when done
     num_atoms = len(mol.atoms)
     while True:
-        # first print the energy of these coordinates
-        print(ff.Energy(True))  # in Hartree
-
-        # now print the gradient on each atom
-        for atom in mol.atoms:
-            grad = ff.GetGradient(atom.OBAtom)
-            print(grad.GetX(), grad.GetY(), grad.GetZ())
-
         # read new coordinates from stdin
         for i in range(num_atoms):
             coordinates = np.fromstring(input(), sep=" ")
@@ -64,6 +56,17 @@ def run(filename):
         # update the molecule geometry for the next energy
         ff.SetCoordinates(mol.OBMol)
 
+        # first print the energy of these coordinates
+        energy = ff.Energy(True)  # in kJ/mol
+        print("AvogadroEnergy:", energy)  # in kJ/mol
+
+        # now print the gradient on each atom
+        print("AvogadroGradient:")
+        for atom in mol.atoms:
+            grad = ff.GetGradient(atom.OBAtom)
+            print(grad.GetX(), grad.GetY(), grad.GetZ())
+
+
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser("MMFF94 calculator")

avogadro/qtplugins/forcefield/scripts/uff.py

@@ -53,19 +53,12 @@ def run(filename):
             atom = mol.atoms[i]
             atom.OBAtom.SetVector(coordinates[0], coordinates[1], coordinates[2])
 
-        with open("/Users/ghutchis/uff.log", "a") as f:
-            f.write("Coordinates:\n")
-            for atom in mol.atoms:
-                f.write(str(atom.OBAtom.GetX()) + " " + str(atom.OBAtom.GetY()) + " " + str(atom.OBAtom.GetZ()) + "\n")
-
         # update the molecule geometry for the next energy
         ff.SetCoordinates(mol.OBMol)
 
         # first print the energy of these coordinates
-        energy = ff.Energy(True)  # in Hartree
-        print("AvogadroEnergy:", energy)  # in Hartree
-        with open("/Users/ghutchis/uff.log", "a") as f:
-            f.write(str(energy) + "\n")
+        energy = ff.Energy(True)  # in kJ/mol
+        print("AvogadroEnergy:", energy)  # in kJ/mol
 
         # now print the gradient on each atom
         print("AvogadroGradient:")