Reading Gaussian fchk vibrations when present

avogadro/quantumio/gaussianfchk.cpp

@@ -57,6 +57,17 @@
                             m_aPos[i + 1] * BOHR_TO_ANGSTROM,
                             m_aPos[i + 2] * BOHR_TO_ANGSTROM));
   }
+
+  if (m_frequencies.size() > 0 &&
+      m_frequencies.size() == m_vibDisplacements.size() &&
+      m_frequencies.size() == m_IRintensities.size()) {
+    molecule.setVibrationFrequencies(m_frequencies);
+    molecule.setVibrationIRIntensities(m_IRintensities);
+    molecule.setVibrationLx(m_vibDisplacements);
+    if (m_RamanIntensities.size())
+      molecule.setVibrationRamanIntensities(m_RamanIntensities);
+  }
+
   // Do simple bond perception.
   molecule.perceiveBondsSimple();
   molecule.perceiveBondOrders();
@@ -84,14 +95,15 @@
 
   string tmp = line.substr(43);
   vector<string> list = Core::split(tmp, ' ');
+  std::vector<double> tmpVec;
 
   // Big switch statement checking for various things we are interested in
   if (Core::contains(key, "RHF")) {
     m_scftype = Rhf;
   } else if (Core::contains(key, "UHF")) {
     m_scftype = Uhf;
   } else if (key == "Number of atoms" && list.size() > 1) {
-    cout << "Number of atoms = " << Core::lexicalCast<int>(list[1]) << endl;
+    m_numAtoms = Core::lexicalCast<int>(list[1]);
   } else if (key == "Charge" && list.size() > 1) {
     m_charge = Core::lexicalCast<signed char>(list[1]);
   } else if (key == "Multiplicity" && list.size() > 1) {
@@ -104,13 +116,11 @@
     m_electronsBeta = Core::lexicalCast<int>(list[1]);
   } else if (key == "Number of basis functions" && list.size() > 1) {
     m_numBasisFunctions = Core::lexicalCast<int>(list[1]);
-    cout << "Number of basis functions = " << m_numBasisFunctions << endl;
+    // cout << "Number of basis functions = " << m_numBasisFunctions << endl;
   } else if (key == "Atomic numbers" && list.size() > 2) {
     m_aNums = readArrayI(in, Core::lexicalCast<int>(list[2]));
     if (static_cast<int>(m_aNums.size()) != Core::lexicalCast<int>(list[2]))
       cout << "Reading atomic numbers failed.\n";
-    else
-      cout << "Reading atomic numbers succeeded.\n";
   }
   // Now we get to the meat of it - coordinates of the atoms
   else if (key == "Current cartesian coordinates" && list.size() > 2) {
@@ -134,15 +144,16 @@
   } else if (key == "Alpha Orbital Energies") {
     if (m_scftype == Rhf) {
       m_orbitalEnergy = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
-      cout << "MO energies, n = " << m_orbitalEnergy.size() << endl;
+      // cout << "MO energies, n = " << m_orbitalEnergy.size() << endl;
     } else if (m_scftype == Uhf) {
       m_alphaOrbitalEnergy =
         readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
-      cout << "Alpha MO energies, n = " << m_alphaOrbitalEnergy.size() << endl;
+      // cout << "Alpha MO energies, n = " << m_alphaOrbitalEnergy.size() <<
+      // endl;
     }
   } else if (key == "Beta Orbital Energies") {
     if (m_scftype != Uhf) {
-      cout << "UHF detected. Reassigning Alpha properties." << endl;
+      // cout << "UHF detected. Reassigning Alpha properties." << endl;
       m_scftype = Uhf;
       m_alphaOrbitalEnergy = m_orbitalEnergy;
       m_orbitalEnergy = vector<double>();
@@ -152,37 +163,63 @@
     }
 
     m_betaOrbitalEnergy = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
-    cout << "Beta MO energies, n = " << m_betaOrbitalEnergy.size() << endl;
+    // cout << "Beta MO energies, n = " << m_betaOrbitalEnergy.size() << endl;
   } else if (key == "Alpha MO coefficients" && list.size() > 2) {
     if (m_scftype == Rhf) {
       m_MOcoeffs = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
-      if (static_cast<int>(m_MOcoeffs.size()) ==
-          Core::lexicalCast<int>(list[2]))
-        cout << "MO coefficients, n = " << m_MOcoeffs.size() << endl;
     } else if (m_scftype == Uhf) {
       m_alphaMOcoeffs = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
-      if (static_cast<int>(m_alphaMOcoeffs.size()) ==
-          Core::lexicalCast<int>(list[2]))
-        cout << "Alpha MO coefficients, n = " << m_alphaMOcoeffs.size() << endl;
     } else {
       cout << "Error, alpha MO coefficients, n = " << m_MOcoeffs.size() << endl;
     }
   } else if (key == "Beta MO coefficients" && list.size() > 2) {
     m_betaMOcoeffs = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
-    if (static_cast<int>(m_betaMOcoeffs.size()) ==
-        Core::lexicalCast<int>(list[2]))
-      cout << "Beta MO coefficients, n = " << m_betaMOcoeffs.size() << endl;
   } else if (key == "Total SCF Density" && list.size() > 2) {
-    if (readDensityMatrix(in, Core::lexicalCast<int>(list[2]), 16))
-      cout << "SCF density matrix read in " << m_density.rows() << endl;
-    else
+    if (!readDensityMatrix(in, Core::lexicalCast<int>(list[2]), 16))
       cout << "Error reading in the SCF density matrix.\n";
   } else if (key == "Spin SCF Density" && list.size() > 2) {
-    if (readSpinDensityMatrix(in, Core::lexicalCast<int>(list[2]), 16))
-      cout << "SCF spin density matrix read in " << m_spinDensity.rows()
-           << endl;
-    else
+    if (!readSpinDensityMatrix(in, Core::lexicalCast<int>(list[2]), 16))
       cout << "Error reading in the SCF spin density matrix.\n";
+  } else if (key == "Number of Normal Modes" && list.size() > 1) {
+    m_normalModes = Core::lexicalCast<int>(list[1]);
+  } else if (key == "Vib-E2" && list.size() > 2) {
+    m_frequencies.clear();
+    m_IRintensities.clear();
+    m_RamanIntensities.clear();
+
+    unsigned threeN = m_numAtoms * 3; // degrees of freedom
+    tmpVec = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
+
+    // read in the first 3N-6 elements as frequencies
+    for (unsigned int i = 0; i < m_normalModes; ++i) {
+      m_frequencies.push_back(tmpVec[i]);
+    }
+    // skip to after threeN elements then read IR intensities
+    for (unsigned int i = threeN; i < threeN + m_normalModes; ++i) {
+      m_IRintensities.push_back(tmpVec[i]);
+    }
+    // now check if we have Raman intensities
+    if (tmpVec[threeN + m_normalModes] != 0.0) {
+      for (unsigned int i = threeN + m_normalModes;
+           i < threeN + 2 * m_normalModes; ++i) {
+        m_RamanIntensities.push_back(tmpVec[i]);
+      }
+    }
+  } else if (key == "Vib-Modes" && list.size() > 2) {
+    tmpVec = readArrayD(in, Core::lexicalCast<int>(list[2]), 16);
+    m_vibDisplacements.clear();
+    if (tmpVec.size() == m_numAtoms * 3 * m_normalModes) {
+      for (unsigned int i = 0; i < m_normalModes; ++i) {
+        Core::Array<Vector3> mode;
+        for (unsigned int j = 0; j < m_numAtoms; ++j) {
+          Vector3 v(tmpVec[i * m_numAtoms * 3 + j * 3],
+                    tmpVec[i * m_numAtoms * 3 + j * 3 + 1],
+                    tmpVec[i * m_numAtoms * 3 + j * 3 + 2]);
+          mode.push_back(v);
+        }
+        m_vibDisplacements.push_back(mode);
+      }
+    }
   }
 }
 

avogadro/quantumio/gaussianfchk.h

@@ -7,6 +7,8 @@
 #define AVOGADRO_QUANTUMIO_GAUSSIANFCHK_H
 
 #include "avogadroquantumioexport.h"
+
+#include <avogadro/core/array.h>
 #include <avogadro/core/gaussianset.h>
 #include <avogadro/io/fileformat.h>
 
@@ -67,6 +69,8 @@
   int m_electrons;
   int m_electronsAlpha;
   int m_electronsBeta;
+  int m_normalModes;
+  int m_numAtoms;
   unsigned char m_spin;
   signed char m_charge;
   unsigned int m_numBasisFunctions;
@@ -87,8 +91,14 @@
   MatrixX m_density;     /// Total density matrix
   MatrixX m_spinDensity; /// Spin density matrix
   Core::ScfType m_scftype;
+
+  Core::Array<double> m_frequencies;
+  Core::Array<double> m_IRintensities;
+  Core::Array<double> m_RamanIntensities;
+  Core::Array<Core::Array<Vector3>> m_vibDisplacements;
 };
-}
-}
+
+} // namespace QuantumIO
+} // namespace Avogadro
 
 #endif