troubleshooting sifile script

cctk/gaussian_file.py

@@ -530,11 +530,12 @@ def _read_gjf_file(cls, filename, return_lines=False):
         else:
             return f
 
-    def get_molecule(self, num=None):
+    def get_molecule(self, num=None, properties=False):
         """
         Returns the last molecule (from an optimization job) or the only molecule (from other jobs).
 
         If ``num`` is specified, returns ``self.ensemble.molecule_list()[num]``
+        If ``properties`` is True, returns ``(molecule, properties)``.
         """
         # some methods pass num=None, which overrides setting the default above
         if num is None:
@@ -543,7 +544,10 @@ def get_molecule(self, num=None):
         if not isinstance(num, int):
             raise TypeError("num must be int")
 
-        return self.ensemble.molecule_list()[num]
+        if properties:
+            return self.ensemble.molecule_list()[num], self.ensemble.properties_list()[num]
+        else:
+            return self.ensemble.molecule_list()[num]
 
     @classmethod
     def _assign_job_types(cls, header):

cctk/molecule.py

@@ -312,14 +312,14 @@ def formula(self, return_dict=False):
             elements = list(formula_dict.keys())
 
             #### H and C always come first
-            if "H" in elements:
-                elements.remove("H")
-                formula += f"H{formula_dict['H']}"
-
             if "C" in elements:
                 elements.remove("C")
                 formula += f"C{formula_dict['C']}"
 
+            if "H" in elements:
+                elements.remove("H")
+                formula += f"H{formula_dict['H']}"
+
             for element in sorted(elements):
                 formula += f"{element}{formula_dict[element]}"
 

cctk/si_file.py

@@ -30,22 +30,25 @@ def write_file(self, filename, append=False):
             filename (str): path to the new file
             append (Bool): whether or not to append to file
         """
+        first = True
         for title, (molecule, properties) in zip(self.titles, self.ensemble.items()):
             assert isinstance(molecule, cctk.Molecule), "molecule is not a valid Molecule object!"
 
             text = f"{title}\n"
             for key, value in generate_info(molecule, properties).items():
                 text += f"{key}:\t{value}\n"
 
+            text += f"Cartesian Coordinates (Å):\n"
             for index, Z in enumerate(molecule.atomic_numbers, start=1):
                 line = molecule.get_vector(index)
                 text += f"{get_symbol(Z):>2}       {line[0]:>13.8f} {line[1]:>13.8f} {line[2]:>13.8f}\n"
 
-            if append:
+            if first:
+                super().write_file(filename, text)
+                first = False
+            else:
                 text += "\n"
                 super().append_to_file(filename, text)
-            else:
-                super().write_file(filename, text)
 
 
 def generate_info(molecule, properties):
@@ -56,9 +59,18 @@ def generate_info(molecule, properties):
         "Multiplicity": molecule.multiplicity,
     }
 
+    # for now manually handling route card and imaginaries, which typically aren't linked to cctk.Molecule.
+    # long-term would be good to manually pass an extra info_dict from the calling environment
+    # to avoid these ad hoc carveouts. ccw 3.8.21
+
     if "route_card" in properties:
         info["Route Card"] = properties["route_card"]
 
+    if "imaginaries" in properties:
+        info["Imaginary Frequencies (cm-1)"] = properties["imaginaries"]
+    else:
+        info["Imaginary Frequencies (cm-1)"] = "None"
+
     if "energy" in properties:
         info["Energy"] = properties["energy"]
     if "enthalpy" in properties:
@@ -67,6 +79,8 @@ def generate_info(molecule, properties):
         info["Gibbs Free Energy"] = properties["gibbs_free_energy"]
     if "quasiharmonic_gibbs_free_energy" in properties:
         info["Gibbs Free Energy (Quasiharmonic Correction)"] = properties["quasiharmonic_gibbs_free_energy"]
+    if "dipole_moment" in properties:
+        info["Dipole Moment (Debye)"] = properties["dipole_moment"]
 
     return info
 

scripts/generate_si.py

@@ -1,6 +1,23 @@
-import cctk
+import cctk, sys, yaml, tqdm
 
-ensemble = cctk.GaussianFile.read_file("test/static/gaussian_file.out").ensemble
-si_file = cctk.SIFile(ensemble, ["title"] * len(ensemble))
+# usage: python generate_si.py config.txt output.txt
+# config should be tab-separated list of titles and filenames
 
-si_file.write_file("si.txt")
+ensemble = cctk.Ensemble()
+titles = list()
+
+with open(sys.argv[1], "r+") as f:
+    for line in tqdm.tqdm(f.readlines()):
+        try:
+            title, path = line.rsplit(' ', 1)
+            gf = cctk.GaussianFile.read_file(path.strip())
+            m, p = gf.get_molecule(properties=True)
+            p["route_card"] = gf.route_card
+            p["imaginaries"] = gf.imaginaries()
+            ensemble.add_molecule(m, p)
+            titles.append(title.strip())
+        except Exception as e:
+            print(f"skipping {line} - error {e}")
+
+si_file = cctk.SIFile(ensemble, titles)
+si_file.write_file(sys.argv[2])