Add initial support for reading molden frequencies

avogadro/quantumio/molden.cpp

@@ -46,8 +46,9 @@
 {
   // Read the log file line by line, most sections are terminated by an empty
   // line, so they should be retained.
-  while (!in.eof())
+  while (!in.eof() && in.good()) {
     processLine(in);
+  }
 
   auto* basis = new GaussianSet;
 
@@ -62,6 +63,24 @@
   molecule.setBasisSet(basis);
   basis->setMolecule(&molecule);
   load(basis);
+
+  if (m_frequencies.size() > 0 &&
+      m_frequencies.size() == m_vibDisplacements.size()) {
+    molecule.setVibrationFrequencies(m_frequencies);
+    molecule.setVibrationLx(m_vibDisplacements);
+
+    // if we don't have intensities, set them all to zero
+    if (m_IRintensities.size() != m_frequencies.size()) {
+      m_IRintensities.resize(m_frequencies.size());
+      for (unsigned int i = 0; i < m_frequencies.size(); i++)
+        m_IRintensities[i] = 0.0;
+    }
+    molecule.setVibrationIRIntensities(m_IRintensities);
+
+    if (m_RamanIntensities.size())
+      molecule.setVibrationRamanIntensities(m_RamanIntensities);
+  }
+
   return true;
 }
 
@@ -86,100 +105,185 @@
     m_mode = GTO;
   } else if (Core::contains(line, "[MO]")) {
     m_mode = MO;
+  } else if (Core::contains(line, "[FREQ]")) {
+    m_mode = Frequencies;
+  } else if (Core::contains(line, "[FR-NORM-COORD]")) {
+    m_mode = VibrationalModes;
+  } else if (Core::contains(line, "[INT]")) {
+    m_mode = Intensities;
   } else if (Core::contains(line, "[")) { // unknown section
     m_mode = Unrecognized;
   } else {
     // We are in a section, and must parse the lines in that section.
     string shell;
     GaussianSet::orbital shellType;
 
+    std::streampos currentPos = in.tellg();
+
     // Parsing a line of data in a section - what mode are we in?
     switch (m_mode) {
       case Atoms:
         readAtom(list);
         break;
       case GTO: {
         // TODO: detect dead files and make bullet-proof
         int atom = Core::lexicalCast<int>(list[0]);
 
         getline(in, line);
         line = Core::trimmed(line);
         while (!line.empty()) { // Read the shell types in this GTO.
           list = Core::split(line, ' ');
           if (list.size() < 1)
             break;
           shell = list[0];
           shellType = GaussianSet::UU;
           if (shell == "sp")
             shellType = GaussianSet::SP;
           else if (shell == "s")
             shellType = GaussianSet::S;
           else if (shell == "p")
             shellType = GaussianSet::P;
           else if (shell == "d")
             shellType = GaussianSet::D;
           else if (shell == "f")
             shellType = GaussianSet::F;
           else if (shell == "g")
             shellType = GaussianSet::G;
 
           if (shellType != GaussianSet::UU) {
             m_shellTypes.push_back(shellType);
             m_shelltoAtom.push_back(atom);
           } else {
             return;
           }
 
           int numGTOs = Core::lexicalCast<int>(list[1]);
           m_shellNums.push_back(numGTOs);
 
           // Now read all the exponents and contraction coefficients.
           for (int gto = 0; gto < numGTOs; ++gto) {
             getline(in, line);
             line = Core::trimmed(line);
             list = Core::split(line, ' ');
             if (list.size() > 1) {
               m_a.push_back(Core::lexicalCast<double>(list[0]));
               m_c.push_back(Core::lexicalCast<double>(list[1]));
             }
             if (shellType == GaussianSet::SP && list.size() > 2)
               m_csp.push_back(Core::lexicalCast<double>(list[2]));
           }
           // Start reading the next shell.
           getline(in, line);
           line = Core::trimmed(line);
         }
       } break;
 
       case MO:
         // Parse the occupation, spin, energy, etc (Occup, Spin, Ene).
         while (!line.empty() && Core::contains(line, "=")) {
           getline(in, line);
           line = Core::trimmed(line);
           list = Core::split(line, ' ');
           if (Core::contains(line, "Occup"))
             m_electrons += Core::lexicalCast<int>(list[1]);
           else if (Core::contains(line, "Ene"))
             m_orbitalEnergy.push_back(Core::lexicalCast<double>(list[1]));
           // TODO: track alpha beta spin
         }
 
         // Parse the molecular orbital coefficients.
-        while (!line.empty() && !Core::contains(line, "=")) {
+        while (!line.empty() && !Core::contains(line, "=") &&
+               !Core::contains(line, "[")) {
           list = Core::split(line, ' ');
           if (list.size() < 2)
             break;
 
           m_MOcoeffs.push_back(Core::lexicalCast<double>(list[1]));
 
           getline(in, line);
           line = Core::trimmed(line);
           list = Core::split(line, ' ');
         }
+        // go back to previous line
+        in.seekg(currentPos);
+        break;
+
+      case Frequencies:
+        // Parse the frequencies.
+        m_frequencies.clear();
+        while (!line.empty() && !Core::contains(line, "[")) {
+          line = Core::trimmed(line);
+          m_frequencies.push_back(Core::lexicalCast<double>(line));
+          currentPos = in.tellg();
+          getline(in, line);
+        }
+        // go back to previous line
+        in.seekg(currentPos);
+        break;
+
+      case VibrationalModes:
+        // Parse the vibrational modes.
+        // should be "vibration 1" etc.
+        // then the normal mode displacements
+        m_vibDisplacements.clear();
+        // shouldn't be more than the number of frequencies
+        while (!line.empty() && !Core::contains(line, "[")) {
+          if (Core::contains(line, "vibration")) {
+            m_vibDisplacements.push_back(Core::Array<Vector3>());
+            getline(in, line);
+            line = Core::trimmed(line);
+            while (!line.empty() && !Core::contains(line, "[") &&
+                   !Core::contains(line, "vibration")) {
+              list = Core::split(line, ' ');
+              if (list.size() < 3)
+                break;
+
+              m_vibDisplacements.back().push_back(Vector3(
+                Core::lexicalCast<double>(list[0]) * BOHR_TO_ANGSTROM_D,
+                Core::lexicalCast<double>(list[1]) * BOHR_TO_ANGSTROM_D,
+                Core::lexicalCast<double>(list[2]) * BOHR_TO_ANGSTROM_D));
+
+              currentPos = in.tellg();
+              getline(in, line);
+              line = Core::trimmed(line);
+            }
+          } else {
+            // we shouldn't hit this, but better to be safe
+            break;
+          }
+
+          // okay, we're either done reading
+          // or we're at the next vibration
+          if (m_vibDisplacements.size() == m_frequencies.size()) {
+            // reset to make sure we don't miss any other sections
+            // (e.g., intensities)
+            in.seekg(currentPos);
+            break;
+          }
+        }
+        break;
+
+      case Intensities:
+        // could be just IR or two pieces including Raman
+        while (!line.empty() && !Core::contains(line, "[")) {
+          list = Core::split(line, ' ');
+          m_IRintensities.push_back(Core::lexicalCast<double>(list[0]));
+          if (list.size() == 2)
+            m_RamanIntensities.push_back(Core::lexicalCast<double>(list[1]));
+
+          if (m_IRintensities.size() == m_frequencies.size()) {
+            // we're done
+            break;
+          }
+
+          currentPos = in.tellg();
+          getline(in, line);
+          line = Core::trimmed(line);
+        }
         break;
       default:
         break;
     }
   }
 }
 

avogadro/quantumio/molden.h

@@ -7,6 +7,7 @@
 #define AVOGADRO_QUANTUMIO_MOLDEN_H
 
 #include "avogadroquantumioexport.h"
+#include <avogadro/core/array.h>
 #include <avogadro/core/gaussianset.h>
 #include <avogadro/io/fileformat.h>
 
@@ -67,17 +68,25 @@ class AVOGADROQUANTUMIO_EXPORT MoldenFile : public Io::FileFormat
   std::vector<double> m_orbitalEnergy;
   std::vector<double> m_MOcoeffs;
 
+  Core::Array<double> m_frequencies;
+  Core::Array<double> m_IRintensities;
+  Core::Array<double> m_RamanIntensities;
+  Core::Array<Core::Array<Vector3>> m_vibDisplacements;
+
   enum Mode
   {
     Atoms,
     GTO,
     MO,
+    Frequencies,
+    VibrationalModes,
+    Intensities,
     Unrecognized
   };
   Mode m_mode;
 };
 
-} // End namespace
-}
+} // namespace QuantumIO
+} // namespace Avogadro
 
 #endif