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 diff --git a/avogadro/io/xyzformat.cpp b/avogadro/io/xyzformat.cpp index 229d67cac0..59278bf48a 100644 --- a/avogadro/io/xyzformat.cpp +++ b/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(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 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; } diff --git a/avogadro/quantumio/orca.cpp b/avogadro/quantumio/orca.cpp index 387fa589b5..3ce4a32562 100644 --- a/avogadro/quantumio/orca.cpp +++ b/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 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(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 006872a0b9..581937f802 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; @@ -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;