Merge pull request #1808 from ghutchis/add-orca-coord-sets

Read ORCA coordinate sets

avogadro/core/variant-inline.h

@@ -55,6 +55,15 @@ inline Variant::Variant(const Vector3f& v) : m_type(Vector)
   m_value.vector = _v;
 }
 
+template <>
+inline Variant::Variant(const std::vector<double>& v) : m_type(Matrix)
+{
+  MatrixX* m = new MatrixX(v.size(), 1);
+  for (size_t i = 0; i < v.size(); ++i)
+    m->coeffRef(i, 0) = v[i];
+  m_value.matrix = m;
+}
+
 inline Variant::Variant(const Variant& variant) : m_type(variant.type())
 {
   if (m_type == String)
@@ -92,6 +101,18 @@ inline bool Variant::setValue(double x, double y, double z)
   return true;
 }
 
+inline bool Variant::setValue(const std::vector<double>& v)
+{
+  clear();
+
+  m_type = Matrix;
+  m_value.matrix = new MatrixX(v.size(), 1);
+  for (size_t i = 0; i < v.size(); ++i)
+    m_value.matrix->coeffRef(i, 0) = v[i];
+
+  return true;
+}
+
 template <typename T>
 inline bool Variant::setValue(T v)
 {
@@ -415,6 +436,19 @@ inline const Vector3& Variant::value() const
   return nullVector;
 }
 
+template <>
+inline std::vector<double> Variant::value() const
+{
+  if (m_type == Matrix && m_value.matrix->cols() == 1) {
+    std::vector<double> list(m_value.matrix->rows());
+    for (int i = 0; i < m_value.matrix->rows(); ++i)
+      list[i] = m_value.matrix->coeff(i, 0);
+    return list;
+  }
+
+  return std::vector<double>();
+}
+
 inline void Variant::clear()
 {
   if (m_type == String) {
@@ -516,6 +550,11 @@ inline Vector3 Variant::toVector3() const
   return value<Vector3>();
 }
 
+inline std::vector<double> Variant::toList() const
+{
+  return value<std::vector<double>>();
+}
+
 // --- Operators ----------------------------------------------------------- //
 inline Variant& Variant::operator=(const Variant& variant)
 {

avogadro/core/variant.h

@@ -71,6 +71,9 @@ class AVOGADROCORE_EXPORT Variant
   /** Sets the value of the variant to a 3D vector */
   bool setValue(double x, double y, double z);
 
+  /** Sets the value of the variant to a vector<double> */
+  bool setValue(const std::vector<double>& v);
+
   /** @return the value of the variant in the type given by \c T. */
   template <typename T>
   T value() const;
@@ -126,6 +129,9 @@ class AVOGADROCORE_EXPORT Variant
   /** @return the value of the variant as a Vector3 */
   inline Vector3 toVector3() const;
 
+  /** @return the value as a vector<double> */
+  inline std::vector<double> toList() const;
+
   /**
    * @return a reference to the value of the variant as a MatrixX.
    * This method will not perform any casting -- if type() is not exactly

avogadro/io/xyzformat.cpp

@@ -37,6 +37,34 @@ using Core::trimmed;
 using std::isalpha;
 #endif
 
+bool findEnergy(const std::string& buffer, double& energyValue)
+{
+  // Check for energy in the comment line
+  // orca uses  E -680.044112849966 (with spaces)
+  // xtb uses energy: -680.044112849966
+  // Open Babel uses Energy: -680.044112849966
+  std::size_t energyStart = buffer.find("energy:");
+  std::size_t offset = 7;
+  if (energyStart == std::string::npos) {
+    energyStart = buffer.find("Energy:");
+  }
+  if (energyStart == std::string::npos) {
+    energyStart = buffer.find(" E ");
+    offset = 3;
+  }
+
+  if (energyStart != std::string::npos) {
+    // find the next whitespace or end of the string
+    std::size_t energyEnd = buffer.find_first_of(" \t", energyStart + offset);
+    if (energyEnd == std::string::npos)
+      energyEnd = buffer.size();
+    std::string energy = buffer.substr(energyStart + offset, energyEnd);
+    energyValue = lexicalCast<double>(energy);
+    return true;
+  }
+  return false;
+}
+
 bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol)
 {
   json opts;
@@ -53,10 +81,17 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol)
 
   string buffer;
   getline(inStream, buffer); // Finish the first line
-  getline(inStream, buffer);
+  getline(inStream, buffer); // comment or name or energy
   if (!buffer.empty())
     mol.setData("name", trimmed(buffer));
 
+  double energy = 0.0;
+  std::vector<double> energies;
+  if (findEnergy(buffer, energy)) {
+    mol.setData("totalEnergy", energy);
+    energies.push_back(energy);
+  }
+
   // check for Lattice= in an extended XYZ from ASE and company
   // e.g. Lattice="H11 H21 H31 H12 H22 H32 H13 H23 H33"
   // https://atomsk.univ-lille.fr/doc/en/format_xyz.html
@@ -171,7 +206,12 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol)
   }
 
   if ((numAtoms2 = lexicalCast<int>(buffer)) && numAtoms == numAtoms2) {
-    getline(inStream, buffer); // Skip the blank
+    getline(inStream, buffer); // comment line
+    // check for properties in the comment line
+    if (findEnergy(buffer, energy)) {
+      energies.push_back(energy);
+    }
+
     mol.setCoordinate3d(mol.atomPositions3d(), 0);
     int coordSet = 1;
     bool done = false;
@@ -240,6 +280,10 @@ bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol)
     mol.setPartialCharges("From File", chargesMatrix);
   }
 
+  if (energies.size() > 1) {
+    mol.setData("energies", energies);
+  }
+
   return true;
 }
 

avogadro/quantumio/orca.cpp

@@ -106,6 +106,19 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
     molecule.setSpectra("NMR", nmrData);
   }
 
+  // this should be the final coordinate set (e.g. the optimized geometry)
+  molecule.setCoordinate3d(molecule.atomPositions3d(), 0);
+  if (m_coordSets.size() > 1) {
+    for (unsigned int i = 0; i < m_coordSets.size(); i++) {
+      Array<Vector3> positions;
+      positions.reserve(molecule.atomCount());
+      for (size_t j = 0; j < molecule.atomCount(); ++j) {
+        positions.push_back(m_coordSets[i][j] * BOHR_TO_ANGSTROM);
+      }
+      molecule.setCoordinate3d(positions, i + 1);
+    }
+  }
+
   // guess bonds and bond orders
   molecule.perceiveBondsSimple();
   molecule.perceiveBondOrders();
@@ -133,6 +146,9 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
   molecule.setData("totalCharge", m_charge);
   molecule.setData("totalSpinMultiplicity", m_spin);
   molecule.setData("dipoleMoment", m_dipoleMoment);
+  molecule.setData("totalEnergy", m_totalEnergy);
+  if (m_energies.size() > 1)
+    molecule.setData("energies", m_energies);
 
   return true;
 }
@@ -152,6 +168,10 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis)
   if (Core::contains(key, "CARTESIAN COORDINATES (A.U.)")) {
     m_coordFactor = 1.; // leave the coords in BOHR ....
     m_currentMode = Atoms;
+    // if there are any current coordinates, push them back
+    if (m_atomPos.size() > 0) {
+      m_coordSets.push_back(m_atomPos);
+    }
     m_atomPos.clear();
     m_atomNums.clear();
     m_atomLabel.clear();
@@ -197,6 +217,11 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis)
     list = Core::split(key, ' ');
     if (list.size() > 3)
       m_spin = Core::lexicalCast<int>(list[3]);
+  } else if (Core::contains(key, "FINAL SINGLE POINT ENERGY")) {
+    list = Core::split(key, ' ');
+    if (list.size() > 4)
+      m_totalEnergy = Core::lexicalCast<double>(list[4]);
+    m_energies.push_back(m_totalEnergy);
   } else if (Core::contains(key, "TOTAL NUMBER OF BASIS SET")) {
     m_currentMode = NotParsing; // no longer reading GTOs
   } else if (Core::contains(key, "NUMBER OF CARTESIAN GAUSSIAN BASIS")) {

avogadro/quantumio/orca.h

@@ -61,6 +61,8 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat
 
   std::vector<int> m_atomNums;
   std::vector<Eigen::Vector3d> m_atomPos;
+  std::vector<std::vector<Eigen::Vector3d>> m_coordSets;
+  std::vector<double> m_energies;
 
   Vector3 m_dipoleMoment;
 
@@ -103,6 +105,7 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat
   int m_homo;
   int m_charge;
   int m_spin;
+  double m_totalEnergy;
 
   int m_currentAtom;
   unsigned int m_numBasisFunctions;