From da5e406dddcb446ce9a3d8e91c2dc470bbad49b0 Mon Sep 17 00:00:00 2001
From: Geoff Hutchison <geoff.hutchison@gmail.com>
Date: Mon, 14 Oct 2024 14:04:09 -0400
Subject: [PATCH] Add initial support for reading molden frequencies

https://discuss.avogadro.cc/t/support-for-molden-frequencies/6331

Signed-off-by: Geoff Hutchison <geoff.hutchison@gmail.com>
---
 avogadro/quantumio/molden.cpp | 94 ++++++++++++++++++++++++++++++++++-
 avogadro/quantumio/molden.h   | 13 ++++-
 2 files changed, 104 insertions(+), 3 deletions(-)

diff --git a/avogadro/quantumio/molden.cpp b/avogadro/quantumio/molden.cpp
index 875bb01fd8..08beeb1b55 100644
--- a/avogadro/quantumio/molden.cpp
+++ b/avogadro/quantumio/molden.cpp
@@ -62,6 +62,24 @@ bool MoldenFile::read(std::istream& in, Core::Molecule& molecule)
   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,6 +104,12 @@ void MoldenFile::processLine(std::istream& in)
     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 {
@@ -93,6 +117,8 @@ void MoldenFile::processLine(std::istream& in)
     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:
@@ -165,7 +191,8 @@ void MoldenFile::processLine(std::istream& in)
         }
 
         // 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;
@@ -176,6 +203,71 @@ void MoldenFile::processLine(std::istream& in)
           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);
+            }
+          }
+
+          // okay, we're either done reading
+          // or we're at the next vibration
+          if (m_vibDisplacements.size() == m_frequencies.size()) {
+            // TODO: might need to go back to read intensities
+            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]));
+
+          getline(in, line);
+          line = Core::trimmed(line);
+        }
         break;
       default:
         break;
diff --git a/avogadro/quantumio/molden.h b/avogadro/quantumio/molden.h
index f06388e0a4..34789700e0 100644
--- a/avogadro/quantumio/molden.h
+++ b/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