22 include constraints and ghost atoms from gaussian output files in the new com files

aqme/argument_parser.py

@@ -70,6 +70,7 @@
     "gen_atoms": [],
     "bs_nogen": "",
     "bs_gen": "",
+    "freeze": [],
     "lowest_only": False,
     "lowest_n": None,
     "e_threshold_qprep": None,

aqme/qcorr.py

@@ -687,6 +687,10 @@ def qcorr_fixing(self, cclib_data, file, atom_types, cartesians):
             program = "gaussian"
         elif cclib_data["metadata"]["QM program"].lower().find("orca") > -1:
             program = "orca"
+
+        # frozen flag treatment-- Heidi Klem
+        if "frozenatoms" in cclib_data["metadata"].keys():
+            self.args.freeze = cclib_data["metadata"]["frozenatoms"]
 
         if program in ["gaussian", "orca"]:
             qprep(
@@ -706,6 +710,7 @@ def qcorr_fixing(self, cclib_data, file, atom_types, cartesians):
                 bs_gen=self.args.bs_gen,
                 bs_nogen=self.args.bs_nogen,
                 gen_atoms=self.args.gen_atoms,
+                freeze=self.args.freeze,
                 create_dat=False,
             )
         else:

aqme/qcorr_utils.py

@@ -436,6 +436,23 @@ def get_json_data(self, file, cclib_data):
                 cclib_data["metadata"]["solvation"] = qm_solv
                 cclib_data["metadata"]["dispersion model"] = qm_disp
                 cclib_data["metadata"]["ground or transition state"] = calc_type
+
+            # Parse for frozen atom flags-- see Heidi Klem
+            elif "Symbolic Z-matrix" in line:
+                # A pull request for cclib has been submitted to parse frozen flags
+                #  so if that cclib version exists AQME will use that instead of the except loop.
+                try:  
+                    cclib_data["metadata"]["frozenatoms"] = cclib_data["optimization"]["frozen atom indices"]
+                except:
+                    if len(outlines[i+2].split()) == 5:
+                        atom_idx = -1 # to make atoms 0 indexed
+                        cclib_data["metadata"]["frozenatoms"] = []
+                        for j in range(i + 2, i + 2 + 10000): # to make sure it goes over all atoms
+                            atom_idx += 1 
+                            if len(outlines[j].split()) == 0:
+                                break
+                            elif outlines[j].split()[1] == '-1':
+                                cclib_data["metadata"]["frozenatoms"].append(atom_idx)
 
             # Basis set name
             elif line[1:15] == "Standard basis":

aqme/qprep.py

@@ -51,6 +51,8 @@
    e_threshold_qprep : float, default=None
       Only create inputs for conformers below the energy threshold (to the lowest conformer)
       of the SDF file
+   freeze : list of int, default=[]
+      Atom indices (zero indexed) to receive -1 flag to be frozen during gaussian optimization 
 """
 ######################################################.
 #        This file stores the QPREP class            #
@@ -108,6 +110,11 @@ def __init__(self, create_dat=True, **kwargs):
             self.args.log.finalize()
             sys.exit()
 
+        if self.args.freeze != [] and self.args.program != "gaussian":
+            self.args.log.write('\nx  Freeze option is only supported for Gaussian input creation currently.')
+            self.args.log.finalize()
+            sys.exit()     
+
         if self.args.destination is None:
             destination = self.args.initial_dir.joinpath("QCALC")
         elif self.args.initial_dir.joinpath(self.args.destination).exists():
@@ -349,6 +356,7 @@ def get_tail(self, qprep_data):
 
                 txt += f"{qm_end_local}\n\n"
 
+
             if self.args.gen_atoms != [] and len(self.args.gen_atoms) > 0:
                 # writes part for Gen/GenECP
                 ecp_used, ecp_not_used, gen_type = [], [], "gen"
@@ -388,6 +396,7 @@ def write(self, qprep_data):
 
         if self.args.program.lower() == "gaussian":
             extension = "com"
+
         elif self.args.program.lower() == "orca":
             extension = "inp"
 
@@ -404,9 +413,18 @@ def write(self, qprep_data):
         fileout.write(header)
 
         for atom_idx in range(0, len(qprep_data["atom_types"])):
+            if self.args.freeze != []:
+                # writes atom flags
+                if atom_idx in self.args.freeze:
+                    atom_flag='-1'
+                else:
+                    atom_flag='0'
+            else:
+                atom_flag=''
             fileout.write(
-                "{0:>2} {1:12.8f} {2:12.8f} {3:12.8f}".format(
+                "{0:>2} {1:>3} {2:12.8f} {3:12.8f} {4:12.8f}".format(
                     qprep_data["atom_types"][atom_idx],
+                    atom_flag,
                     qprep_data["cartesians"][atom_idx][0],
                     qprep_data["cartesians"][atom_idx][1],
                     qprep_data["cartesians"][atom_idx][2],

aqme/utils.py

@@ -340,6 +340,7 @@ def command_line_args():
         "constraints_dist",
         "constraints_angle",
         "constraints_dihedral",
+        "freeze",
         "atom_types",
         "cartesians",
         "nmr_atoms",
@@ -877,4 +878,4 @@ def get_files(value):
                 new_value.append(val)
         else:
             new_value.append(val.as_posix())
-    return new_value
+    return new_value