From e0a8804dfd65f0c1db5bab8beafb770c24fdcb6e Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Wed, 20 Nov 2024 17:41:19 -0500 Subject: [PATCH 1/5] Fix some minor bugs with vector support in the Variant class Signed-off-by: Geoff Hutchison Signed-off-by: Geoff Hutchison --- avogadro/core/variant-inline.h | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/avogadro/core/variant-inline.h b/avogadro/core/variant-inline.h index 77379df89c..cc884d6aba 100644 --- a/avogadro/core/variant-inline.h +++ b/avogadro/core/variant-inline.h @@ -411,7 +411,8 @@ inline const Vector3& Variant::value() const if (m_type == Vector) return *m_value.vector; - return Vector3::Zero(); + static Vector3 nullVector(0, 0, 0); + return nullVector; } inline void Variant::clear() @@ -422,6 +423,9 @@ inline void Variant::clear() } else if (m_type == Matrix) { delete m_value.matrix; m_value.matrix = 0; + } else if (m_type == Vector) { + delete m_value.vector; + m_value.vector = 0; } m_type = Null; @@ -527,6 +531,8 @@ inline Variant& Variant::operator=(const Variant& variant) m_value.string = new std::string(variant.toString()); else if (m_type == Matrix) m_value.matrix = new MatrixX(*variant.m_value.matrix); + else if (m_type == Vector) + m_value.vector = new Vector3(*variant.m_value.vector); else if (m_type != Null) m_value = variant.m_value; } From 125e7a8308267abb087a562e682be2f8c77eacea Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Thu, 21 Nov 2024 11:03:24 -0500 Subject: [PATCH 2/5] Add support for reading energies in XYZ trajectories Signed-off-by: Geoff Hutchison --- avogadro/io/xyzformat.cpp | 48 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 46 insertions(+), 2 deletions(-) diff --git a/avogadro/io/xyzformat.cpp b/avogadro/io/xyzformat.cpp index 229d67cac0..8cb4f5129d 100644 --- a/avogadro/io/xyzformat.cpp +++ b/avogadro/io/xyzformat.cpp @@ -37,6 +37,34 @@ using Core::trimmed; using std::isalpha; #endif +bool findEnergy(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; + std::vector energies; + 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(energy); + return true; + } + return false; // didn't find an energy +} + 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 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(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; } From 4496561abe995040e0517babf5211a355756d345 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Thu, 21 Nov 2024 11:03:49 -0500 Subject: [PATCH 3/5] Initial support for reading coordinate sets from ORCA Needs more testing Signed-off-by: Geoff Hutchison --- avogadro/quantumio/orca.cpp | 16 ++++++++++++++++ avogadro/quantumio/orca.h | 2 ++ 2 files changed, 18 insertions(+) diff --git a/avogadro/quantumio/orca.cpp b/avogadro/quantumio/orca.cpp index 387fa589b5..605a952c94 100644 --- a/avogadro/quantumio/orca.cpp +++ b/avogadro/quantumio/orca.cpp @@ -106,6 +106,18 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule) molecule.setSpectra("NMR", nmrData); } + molecule.setCoordinate3d(molecule.atomPositions3d(), 0); + if (m_coordSets.size() > 1) { + for (unsigned int i = 1; i < m_coordSets.size(); i++) { + Array positions; + positions.reserve(molecule.atomCount()); + for (size_t j = 0; j < molecule.atomCount(); ++j) { + positions.push_back(m_atomPos[j] * BOHR_TO_ANGSTROM); + } + molecule.setCoordinate3d(positions, i); + } + } + // guess bonds and bond orders molecule.perceiveBondsSimple(); molecule.perceiveBondOrders(); @@ -152,6 +164,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(); diff --git a/avogadro/quantumio/orca.h b/avogadro/quantumio/orca.h index 006872a0b9..6d9eed5282 100644 --- a/avogadro/quantumio/orca.h +++ b/avogadro/quantumio/orca.h @@ -61,6 +61,8 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat std::vector m_atomNums; std::vector m_atomPos; + std::vector> m_coordSets; + std::vector m_energies; Vector3 m_dipoleMoment; From 28f510b654942784f0c07220c0b012cbe1028414 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Fri, 22 Nov 2024 13:15:31 -0500 Subject: [PATCH 4/5] Add support for a variant from vector Signed-off-by: Geoff Hutchison --- avogadro/core/variant-inline.h | 39 ++++++++++++++++++++++++++++++++++ avogadro/core/variant.h | 6 ++++++ 2 files changed, 45 insertions(+) diff --git a/avogadro/core/variant-inline.h b/avogadro/core/variant-inline.h index cc884d6aba..84d6986822 100644 --- a/avogadro/core/variant-inline.h +++ b/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& 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& 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 inline bool Variant::setValue(T v) { @@ -415,6 +436,19 @@ inline const Vector3& Variant::value() const return nullVector; } +template <> +inline std::vector Variant::value() const +{ + if (m_type == Matrix && m_value.matrix->cols() == 1) { + std::vector 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(); +} + inline void Variant::clear() { if (m_type == String) { @@ -516,6 +550,11 @@ inline Vector3 Variant::toVector3() const return value(); } +inline std::vector Variant::toList() const +{ + return value>(); +} + // --- Operators ----------------------------------------------------------- // inline Variant& Variant::operator=(const Variant& variant) { diff --git a/avogadro/core/variant.h b/avogadro/core/variant.h index c7b1d09f53..56b33c7a8b 100644 --- a/avogadro/core/variant.h +++ b/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 */ + bool setValue(const std::vector& v); + /** @return the value of the variant in the type given by \c T. */ template 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 */ + inline std::vector 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 From 6aa00f30d5fd323d3f03056471cfbb6c934b3c95 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Fri, 22 Nov 2024 13:22:51 -0500 Subject: [PATCH 5/5] Ensure energies for each frame is reserved Signed-off-by: Geoff Hutchison --- avogadro/io/xyzformat.cpp | 6 +++--- avogadro/quantumio/orca.cpp | 15 ++++++++++++--- avogadro/quantumio/orca.h | 1 + 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/avogadro/io/xyzformat.cpp b/avogadro/io/xyzformat.cpp index 8cb4f5129d..59278bf48a 100644 --- a/avogadro/io/xyzformat.cpp +++ b/avogadro/io/xyzformat.cpp @@ -37,7 +37,7 @@ using Core::trimmed; using std::isalpha; #endif -bool findEnergy(std::string buffer, double energyValue) +bool findEnergy(const std::string& buffer, double& energyValue) { // Check for energy in the comment line // orca uses E -680.044112849966 (with spaces) @@ -45,7 +45,6 @@ bool findEnergy(std::string buffer, double energyValue) // Open Babel uses Energy: -680.044112849966 std::size_t energyStart = buffer.find("energy:"); std::size_t offset = 7; - std::vector energies; if (energyStart == std::string::npos) { energyStart = buffer.find("Energy:"); } @@ -53,6 +52,7 @@ bool findEnergy(std::string buffer, double energyValue) 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); @@ -62,7 +62,7 @@ bool findEnergy(std::string buffer, double energyValue) energyValue = lexicalCast(energy); return true; } - return false; // didn't find an energy + return false; } bool XyzFormat::read(std::istream& inStream, Core::Molecule& mol) diff --git a/avogadro/quantumio/orca.cpp b/avogadro/quantumio/orca.cpp index 605a952c94..3ce4a32562 100644 --- a/avogadro/quantumio/orca.cpp +++ b/avogadro/quantumio/orca.cpp @@ -106,15 +106,16 @@ 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 = 1; i < m_coordSets.size(); i++) { + for (unsigned int i = 0; i < m_coordSets.size(); i++) { Array positions; positions.reserve(molecule.atomCount()); for (size_t j = 0; j < molecule.atomCount(); ++j) { - positions.push_back(m_atomPos[j] * BOHR_TO_ANGSTROM); + positions.push_back(m_coordSets[i][j] * BOHR_TO_ANGSTROM); } - molecule.setCoordinate3d(positions, i); + molecule.setCoordinate3d(positions, i + 1); } } @@ -145,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; } @@ -213,6 +217,11 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) list = Core::split(key, ' '); if (list.size() > 3) m_spin = Core::lexicalCast(list[3]); + } else if (Core::contains(key, "FINAL SINGLE POINT ENERGY")) { + list = Core::split(key, ' '); + if (list.size() > 4) + m_totalEnergy = Core::lexicalCast(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")) { diff --git a/avogadro/quantumio/orca.h b/avogadro/quantumio/orca.h index 6d9eed5282..581937f802 100644 --- a/avogadro/quantumio/orca.h +++ b/avogadro/quantumio/orca.h @@ -105,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;