Add electronic and CD spectra to parsed data

Also confirm total charge, multiplicity, and dipole moment

Signed-off-by: Geoff Hutchison <[email protected]>

avogadro/quantumio/orca.cpp

@@ -78,6 +78,27 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
       molecule.setVibrationRamanIntensities(m_RamanIntensities);
   }
 
+  if (m_electronicTransitions.size() > 0 &&
+      m_electronicTransitions.size() == m_electronicIntensities.size()) {
+    MatrixX electronicData(m_electronicTransitions.size(), 2);
+    for (size_t i = 0; i < m_electronicTransitions.size(); ++i) {
+      electronicData(i, 0) = m_electronicTransitions[i];
+      electronicData(i, 1) = m_electronicIntensities[i];
+    }
+    molecule.setSpectra("Electronic", electronicData);
+    std::cout << "UV/Vis data found." << electronicData.rows() << std::endl;
+    if (m_electronicRotations.size() == m_electronicTransitions.size()) {
+      MatrixX electronicRotations(m_electronicTransitions.size(), 2);
+      for (size_t i = 0; i < m_electronicTransitions.size(); ++i) {
+        electronicRotations(i, 0) = m_electronicTransitions[i];
+        electronicRotations(i, 1) = m_electronicRotations[i];
+      }
+      molecule.setSpectra("CircularDichroism", electronicRotations);
+      std::cout << "CircularDichroism data found." << electronicRotations.rows()
+                << std::endl;
+    }
+  }
+
   // guess bonds and bond orders
   molecule.perceiveBondsSimple();
   molecule.perceiveBondOrders();
@@ -102,6 +123,10 @@ bool ORCAOutput::read(std::istream& in, Core::Molecule& molecule)
   basis->setMolecule(&molecule);
   load(basis);
 
+  molecule.setData("totalCharge", m_charge);
+  molecule.setData("totalSpinMultiplicity", m_spin);
+  molecule.setData("dipoleMoment", m_dipoleMoment);
+
   return true;
 }
 
@@ -157,16 +182,45 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis)
         getline(in, key);
       }
     }
+  } else if (Core::contains(key, "Total Charge")) {
+    list = Core::split(key, ' ');
+    if (list.size() > 4)
+      m_charge = Core::lexicalCast<int>(list[4]);
+  } else if (Core::contains(key, "Multiplicity")) {
+    list = Core::split(key, ' ');
+    if (list.size() > 3)
+      m_spin = Core::lexicalCast<int>(list[3]);
   } 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")) {
     m_currentMode = NotParsing; // no longer reading GTOs
   } else if (Core::contains(key, "Number of Electrons")) {
     list = Core::split(key, ' ');
     m_electrons = Core::lexicalCast<int>(list[5]);
+  } else if (Core::contains(key, "Total Dipole Moment")) {
+    list = Core::split(key, ' ');
+    m_dipoleMoment = Eigen::Vector3d(Core::lexicalCast<double>(list[4]),
+                                     Core::lexicalCast<double>(list[5]),
+                                     Core::lexicalCast<double>(list[6]));
+    // convert from atomic units to Debye
+    // e.g. https://en.wikipedia.org/wiki/Debye
+    m_dipoleMoment *= 2.54174628;
   } else if (Core::contains(key, "Mayer bond orders")) {
     m_currentMode = BondOrders;
     // starts at the next line
+  } else if (Core::contains(
+               key,
+               "ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS")) {
+    m_currentMode = Electronic;
+    for (int i = 0; i < 4; ++i) {
+      getline(in, key); // skip header
+    }
+    // starts at the next line
+  } else if (Core::contains(key, "CD SPECTRUM")) {
+    m_currentMode = ECD;
+    for (int i = 0; i < 4; ++i) {
+      getline(in, key); // skip header
+    }
   } else if (Core::contains(key, "ORBITAL ENERGIES")) {
     m_currentMode = OrbitalEnergies;
     getline(in, key); // skip ------------
@@ -509,6 +563,61 @@ void ORCAOutput::processLine(std::istream& in, GaussianSet* basis)
         m_currentMode = NotParsing;
         break;
       }
+      case Electronic: {
+        if (key.empty())
+          break;
+        list = Core::split(key, ' ');
+        double wavenumbers;
+        while (!key.empty()) {
+          // should have 8 columns
+          if (list.size() != 8) {
+            getline(in, key);
+            key = Core::trimmed(key);
+            list = Core::split(key, ' ');
+            continue; // skip any spin-forbidden transitions
+          }
+
+          wavenumbers = Core::lexicalCast<double>(list[1]);
+          // convert to eV
+          m_electronicTransitions.push_back(wavenumbers / 8065.544);
+          m_electronicIntensities.push_back(Core::lexicalCast<double>(list[3]));
+
+          getline(in, key);
+          key = Core::trimmed(key);
+          list = Core::split(key, ' ');
+          if (list.size() < 2)
+            break; // hit the blank line
+        }
+        m_currentMode = NotParsing;
+      }
+      case ECD: {
+        if (key.empty())
+          break;
+        list = Core::split(key, ' ');
+
+        double wavenumbers;
+        while (!key.empty()) {
+          // should have 7 columns
+          if (list.size() != 7) {
+            getline(in, key);
+            key = Core::trimmed(key);
+            list = Core::split(key, ' ');
+            continue; // skip any spin-forbidden transitions
+          }
+
+          wavenumbers = Core::lexicalCast<double>(list[1]);
+          // convert to eV
+          m_electronicTransitions.push_back(wavenumbers / 8065.544);
+          m_electronicRotations.push_back(Core::lexicalCast<double>(list[3]));
+
+          getline(in, key);
+          key = Core::trimmed(key);
+          list = Core::split(key, ' ');
+          if (list.size() < 2)
+            break; // hit the blank line
+        }
+        m_currentMode = NotParsing;
+      }
       case GTO: {
         //            // should start at the first newGTO
         if (key.empty())