From 7bca87bc3dff4013fc7972e778c91facc06fbdfa Mon Sep 17 00:00:00 2001
From: Geoff Hutchison <geoff.hutchison@gmail.com>
Date: Sat, 28 Oct 2023 19:40:55 -0400
Subject: [PATCH] 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 <geoff.hutchison@gmail.com>
---
 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<unsigned char>(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<int>(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<int> m_atomNums;
   std::vector<Eigen::Vector3d> m_atomPos;
 
+  std::vector<std::vector<int>> m_bondOrders;
+
   std::vector<int> shellFunctions;
   std::vector<Core::GaussianSet::orbital> shellTypes;
   std::vector<std::vector<int>> m_orcaNumShells;
@@ -83,6 +85,7 @@ class AVOGADROQUANTUMIO_EXPORT ORCAOutput : public Io::FileFormat
     Raman,
     Electronic,
     NMR,
+    BondOrders,
     NotParsing,
     Unrecognized
   };