From 7bca87bc3dff4013fc7972e778c91facc06fbdfa Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 28 Oct 2023 19:40:55 -0400 Subject: [PATCH 1/3] Parsing orca output would crash when swapping orbitals Probably just masks the underlying problem, but I need more examples https://discuss.avogadro.cc/t/cant-open-a-out-file-in-1-98-0/ Signed-off-by: Geoff Hutchison --- avogadro/quantumio/orca.cpp | 39 ++++++++++++++++++++++++++++++++++++- avogadro/quantumio/orca.h | 3 +++ 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/avogadro/quantumio/orca.cpp b/avogadro/quantumio/orca.cpp index 8349990049..86db7e25b1 100644 --- a/avogadro/quantumio/orca.cpp +++ b/avogadro/quantumio/orca.cpp @@ -78,10 +78,22 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule) molecule.setVibrationRamanIntensities(m_RamanIntensities); } - // Do simple bond perception. + // for now, do simple bond perception molecule.perceiveBondsSimple(); molecule.perceiveBondOrders(); + // add bonds from calculated bond orders + if (m_bondOrders.size() > 0) { + for (unsigned int i = 0; i < m_bondOrders.size(); i++) { + // m_bondOrders[i][0] is the first atom + // m_bondOrders[i][1] is the second atom + // m_bondOrders[i][2] is the bond order + if (m_bondOrders[i].size() > 2) + molecule.addBond(m_bondOrders[i][0], m_bondOrders[i][1], + static_cast(m_bondOrders[i][2])); + } + } + molecule.setBasisSet(basis); basis->setMolecule(&molecule); load(basis); @@ -148,6 +160,9 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) } else if (Core::contains(key, "Number of Electrons")) { list = Core::split(key, ' '); m_electrons = Core::lexicalCast(list[5]); + } else if (Core::contains(key, "Mayer bond orders")) { + m_currentMode = BondOrders; + // starts at the next line } else if (Core::contains(key, "ORBITAL ENERGIES")) { m_currentMode = OrbitalEnergies; getline(in, key); // skip ------------ @@ -266,6 +281,20 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) m_currentMode = NotParsing; break; } + case BondOrders: { + if (key.empty()) + break; + + m_bondOrders.clear(); + while (key[0] == "B") { + // @todo .. parse the bonds based on character position + // e.g. B( 0-Ru, 1-C ) : 0.4881 B( 0-Ru, 4-C ) : 0.6050 + getline(in, key); + key = Core::trimmed(key); + } + + m_currentMode = NotParsing; + } case OrbitalEnergies: { // should start at the first orbital if (!m_readBeta) @@ -579,10 +608,14 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) while (idx < orcaOrbitals.size()) { if (Core::contains(orcaOrbitals.at(idx), "pz")) { for (unsigned int i = 0; i < numColumns; i++) { + if (idx + 1 >= columns[i].size()) + break; std::swap(columns[i].at(idx), columns[i].at(idx + 1)); } idx++; for (unsigned int i = 0; i < numColumns; i++) { + if (idx + 1 >= columns[i].size()) + break; std::swap(columns[i].at(idx), columns[i].at(idx + 1)); } idx++; @@ -654,10 +687,14 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) while (idx < orcaOrbitals.size()) { if (Core::contains(orcaOrbitals.at(idx), "pz")) { for (unsigned int i = 0; i < numColumns; i++) { + if (idx + 1 >= columns[i].size()) + break; std::swap(columns[i].at(idx), columns[i].at(idx + 1)); } idx++; for (unsigned int i = 0; i < numColumns; i++) { + if (idx + 1 >= columns[i].size()) + break; std::swap(columns[i].at(idx), columns[i].at(idx + 1)); } idx++; diff --git a/avogadro/quantumio/orca.h b/avogadro/quantumio/orca.h index 5334d1de33..020b33db89 100644 --- a/avogadro/quantumio/orca.h +++ b/avogadro/quantumio/orca.h @@ -62,6 +62,8 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat std::vector m_atomNums; std::vector m_atomPos; + std::vector> m_bondOrders; + std::vector shellFunctions; std::vector shellTypes; std::vector> m_orcaNumShells; @@ -83,6 +85,7 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat Raman, Electronic, NMR, + BondOrders, NotParsing, Unrecognized }; From de5ee4a6f8ee47a55583a61d2ab7e1824fae7796 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 28 Oct 2023 19:51:29 -0400 Subject: [PATCH 2/3] Fix compile error Signed-off-by: Geoff Hutchison --- avogadro/quantumio/orca.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/avogadro/quantumio/orca.cpp b/avogadro/quantumio/orca.cpp index 86db7e25b1..bc15cecc96 100644 --- a/avogadro/quantumio/orca.cpp +++ b/avogadro/quantumio/orca.cpp @@ -90,7 +90,7 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule) // m_bondOrders[i][2] is the bond order if (m_bondOrders[i].size() > 2) molecule.addBond(m_bondOrders[i][0], m_bondOrders[i][1], - static_cast(m_bondOrders[i][2])); + static_cast(m_bondOrders[i][2])); } } @@ -286,7 +286,7 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) break; m_bondOrders.clear(); - while (key[0] == "B") { + while (key[0] == 'B') { // @todo .. parse the bonds based on character position // e.g. B( 0-Ru, 1-C ) : 0.4881 B( 0-Ru, 4-C ) : 0.6050 getline(in, key); From cc5fa6591593c8421e55d7ec368890a6d629665e Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Sat, 28 Oct 2023 22:11:14 -0400 Subject: [PATCH 3/3] Add support for reading Mayer bond orders >1.6 Should help with adding orders for inorganic / organometallics Signed-off-by: Geoff Hutchison --- avogadro/quantumio/orca.cpp | 59 +++++++++++++++++++++++++++++++++---- 1 file changed, 54 insertions(+), 5 deletions(-) diff --git a/avogadro/quantumio/orca.cpp b/avogadro/quantumio/orca.cpp index bc15cecc96..1748711678 100644 --- a/avogadro/quantumio/orca.cpp +++ b/avogadro/quantumio/orca.cpp @@ -78,19 +78,23 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule) molecule.setVibrationRamanIntensities(m_RamanIntensities); } - // for now, do simple bond perception + // guess bonds and bond orders molecule.perceiveBondsSimple(); molecule.perceiveBondOrders(); - // add bonds from calculated bond orders + // check bonds from calculated bond orders if (m_bondOrders.size() > 0) { for (unsigned int i = 0; i < m_bondOrders.size(); i++) { // m_bondOrders[i][0] is the first atom // m_bondOrders[i][1] is the second atom // m_bondOrders[i][2] is the bond order - if (m_bondOrders[i].size() > 2) - molecule.addBond(m_bondOrders[i][0], m_bondOrders[i][1], - static_cast(m_bondOrders[i][2])); + if (m_bondOrders[i].size() > 2) { + auto bond = molecule.bond(m_bondOrders[i][0], m_bondOrders[i][1]); + if (bond.isValid() && bond.order() != m_bondOrders[i][2]) { + // if the bond order is different, change it + bond.setOrder(static_cast(m_bondOrders[i][2])); + } + } } } @@ -289,6 +293,51 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis) while (key[0] == 'B') { // @todo .. parse the bonds based on character position // e.g. B( 0-Ru, 1-C ) : 0.4881 B( 0-Ru, 4-C ) : 0.6050 + Index firstAtom = Core::lexicalCast(key.substr(2, 3)); + Index secondAtom = Core::lexicalCast(key.substr(9, 3)); + double bondOrder = Core::lexicalCast(key.substr(18, 9)); + + if (bondOrder > 1.6) { + std::vector bond; + bond.push_back(firstAtom); + bond.push_back(secondAtom); + bond.push_back(static_cast(std::round(bondOrder))); + m_bondOrders.push_back(bond); + std::cout << " bond " << firstAtom << " " << secondAtom << " " + << bondOrder << std::endl; + } + + if (key.size() > 54 && key[28] == 'B') { + firstAtom = Core::lexicalCast(key.substr(30, 3)); + secondAtom = Core::lexicalCast(key.substr(37, 3)); + bondOrder = Core::lexicalCast(key.substr(46, 9)); + + if (bondOrder > 1.6) { + std::vector bond; + bond.push_back(firstAtom); + bond.push_back(secondAtom); + bond.push_back(static_cast(std::round(bondOrder))); + std::cout << " bond " << firstAtom << " " << secondAtom << " " + << bondOrder << std::endl; + m_bondOrders.push_back(bond); + } + } + if (key.size() > 82 && key[56] == 'B') { + firstAtom = Core::lexicalCast(key.substr(58, 3)); + secondAtom = Core::lexicalCast(key.substr(65, 3)); + bondOrder = Core::lexicalCast(key.substr(74, 9)); + + if (bondOrder > 1.6) { + std::vector bond; + bond.push_back(firstAtom); + bond.push_back(secondAtom); + bond.push_back(static_cast(std::round(bondOrder))); + std::cout << " bond " << firstAtom << " " << secondAtom << " " + << bondOrder << std::endl; + m_bondOrders.push_back(bond); + } + } + getline(in, key); key = Core::trimmed(key); }