read vibrational modes

cctk/parse_gaussian.py

@@ -520,6 +520,7 @@ def parse_modes(freq_block, num_atoms, hpmodes=False):
     freqs = list()
     masses = list()
     force_ks = list()
+    intensities = list()
     displacements = list()
 
     if len(freq_block) == 0:
@@ -549,6 +550,7 @@ def parse_modes(freq_block, num_atoms, hpmodes=False):
 
             masses += list(filter(None, re.split(" +", lines[1])))[3:]
             force_ks += list(filter(None, re.split(" +", lines[2])))[3:]
+            intensities += list(filter(None, re.split(" +", lines[3])))[3:]
 
             for line in lines[6:]:
                 fields = re.split(" +", line)
@@ -572,6 +574,7 @@ def parse_modes(freq_block, num_atoms, hpmodes=False):
             freqs += re.split(" +", lines[0])[2:]
             masses += re.split(" +", lines[1])[4:]
             force_ks += re.split(" +", lines[2])[4:]
+            intensities += re.split(" +", lines[3])[4:]
 
             for line in lines[5:]:
                 fields = re.split(" +", line)
@@ -595,15 +598,16 @@ def parse_modes(freq_block, num_atoms, hpmodes=False):
     freqs = [float(x) for x in freqs]
     masses = [float(x) for x in masses]
     force_ks = [float(x) for x in force_ks]
+    intensities = [float(x) for x in intensities]
 
     assert len(freqs) == len(masses)
     assert len(freqs) == len(force_ks)
     assert len(freqs) == len(displacements)
 
     modes = list()
-    for f, m, k, d in zip(freqs, masses, force_ks, displacements):
+    for f, m, k, i, d in zip(freqs, masses, force_ks, intensities, displacements):
         k *= 143.9326 # mdyne Å**-1 to kcal/mol Å**-2
-        modes.append(cctk.VibrationalMode(frequency=f, reduced_mass=m, force_constant=k, displacements=d))
+        modes.append(cctk.VibrationalMode(frequency=f, reduced_mass=m, force_constant=k, intensity=i, displacements=d))
 
     return modes
 

cctk/vibrational_mode.py

@@ -22,10 +22,11 @@ class VibrationalMode:
         frequency (float): frequency, in cm-1
         force_constant (float): force constant, in kcal/mol per Å
         reduced_mass (float): mass, in amus
+        intensity (float): IR intensity
         displacements (cctk.OneIndexedArray): displacements
 
     """
-    def __init__(self, frequency, force_constant, reduced_mass, displacements):
+    def __init__(self, frequency, force_constant, reduced_mass, intensity, displacements):
         assert isinstance(frequency, float)
         self.frequency = frequency
 
@@ -35,6 +36,9 @@ def __init__(self, frequency, force_constant, reduced_mass, displacements):
         assert isinstance(reduced_mass, float)
         self.reduced_mass = reduced_mass
 
+        assert isinstance(intensity, float)
+        self.intensity = intensity
+
         assert isinstance(displacements, cctk.OneIndexedArray)
         self.displacements = displacements
 

test/test_freqs.py

@@ -15,9 +15,14 @@ def test_read(self):
         mol = file.get_molecule()
         self.assertEqual(len(mol.vibrational_modes), 9)
         self.assertEqual(mol.vibrational_modes[-1].frequency, 3119.6807)
+        self.assertEqual(mol.vibrational_modes[-1].intensity, 26.1449)
 
         path = "test/static/methane_hpmodes.out"
         file = cctk.GaussianFile.read_file(path)
+        mol = file.get_molecule()
+        self.assertEqual(len(mol.vibrational_modes), 9)
+        self.assertEqual(mol.vibrational_modes[-1].frequency, 3121.6678)
+        self.assertEqual(mol.vibrational_modes[-1].intensity, 26.1383)
 
     def draw_qho(self):
         from asciichartpy import plot