- Add dynamic extraction of units

- Add units to Gaussian schema
- Remove double Gaussian frequency quantities

electronicparsers/gaussian/metainfo/gaussian.py

@@ -603,45 +603,23 @@ class x_gaussian_section_frequencies(MSection):
         validate=False,
     )
 
-    x_gaussian_frequency_values = Quantity(
-        type=str,
-        shape=['number_of_frequency_rows'],
-        description="""
-        values of frequencies, in cm-1
-        """,
-    )
-
     x_gaussian_frequencies = Quantity(
         type=np.float64,
+        unit='1/m',
         shape=['number_of_frequencies'],
         description="""
-        values of frequencies, in cm-1
+        values of frequencies
         """,
-    )
-
-    x_gaussian_reduced_masses = Quantity(
-        type=np.float64,
-        shape=['number_of_reduced_masses_rows'],
-        description="""
-        values of normal mode reduced masses
-        """,
-    )
+    )  # only store the '--' header, not '---'
 
     x_gaussian_red_masses = Quantity(
         type=np.float64,
+        unit='kg',
         shape=['number_of_frequencies'],
         description="""
         values of normal mode reduced masses
         """,
-    )
-
-    x_gaussian_normal_modes = Quantity(
-        type=str,
-        shape=['number_of_normal_modes_rows'],
-        description="""
-        normal mode vectors
-        """,
-    )
+    )  # only store the '--' header, not '---'
 
     x_gaussian_normal_mode_values = Quantity(
         type=np.float64,

electronicparsers/gaussian/parser.py

@@ -100,6 +100,25 @@ def str_to_force_constants(val_in):
             fc = fc + fc.T - np.diag(fc.diagonal())
             return fc
 
+        def str_to_units(unit: str):
+            """Map native Gaussian units to pint units.
+            Assumes lower case string input."""  # TODO handle compound units recursively
+            unit = unit.lower()
+            if unit == 'cm**-1':
+                return ureg.cm_1
+            elif unit == 'ghz':
+                return ureg.gigahertz
+            elif unit == 'kcal/mol':
+                return ureg.kilocalorie / ureg.mole
+            elif unit == 'kj/mol':
+                return ureg.kilojoule / ureg.mole
+            elif unit == 'j':
+                return ureg.joule
+            elif unit == 'amu':
+                return ureg.amu
+            else:
+                raise ValueError(f'Unknown unit {unit}')
+
         orientation_quantities = [
             Quantity(
                 'standard_orientation',
@@ -390,6 +409,16 @@ def str_to_force_constants(val_in):
                 dtype=float,
                 unit='debye * angstrom**3',
             ),
+            Quantity(
+                'frequency_unit',
+                r'[Hh]armonic frequencies \((\S+)\)',
+                str_operation=str_to_units,
+            ),
+            Quantity(
+                'reduced_mass_unit',
+                r'reduced masses \((\S+)\)',
+                str_operation=str_to_units,
+            ),
             Quantity(
                 'frequencies',
                 r'\n *Frequencies [\-]{2}\s+(.+)',
@@ -1109,18 +1138,19 @@ def parse_energy_corrections(method, iteration=False):
         # vibrational frequencies
         frequencies = section.get('frequencies')
         if frequencies is not None:
-            # frequencies in old parsers are in J, not consistent with metainfo
             sec_frequencies = x_gaussian_section_frequencies()
             sec_run.x_gaussian_section_frequencies.append(sec_frequencies)
-            sec_frequencies.x_gaussian_frequencies = np.hstack(frequencies)
+
+            sec_frequencies.x_gaussian_frequencies = np.hstack(
+                frequencies
+            ) * section.get('frequency_unit', ureg.cm_1)
+
             reduced_masses = section.get('reduced_masses')
             if reduced_masses is not None:
-                reduced_masses = (
-                    np.array(np.hstack(reduced_masses), dtype=np.float64) * ureg.amu
-                )
-                sec_frequencies.x_gaussian_red_masses = reduced_masses.to(
-                    'kg'
-                ).magnitude
+                sec_frequencies.x_gaussian_red_masses = np.hstack(
+                    reduced_masses
+                ) * section.get('reduced_mass_unit', ureg.amu)
+
             normal_modes = section.get('normal_modes')
             if normal_modes is not None:
                 normal_modes = np.hstack(normal_modes)