Adds option to freeze all torsions and then unfreeze a subset

optking/addIntcos.py

@@ -652,6 +652,23 @@ def frozen_tors_from_input(frozen_tors_list, o_molsys):
             logger.info("Frozen dihedral not present, so adding it.\n")
             o_molsys.fragments[f].intcos.append(torsAngle)
 
+def freeze_all_torsions(o_molsys, skipped_torsions=[]):
+    """ Just blindly freezes all torsions """
+
+    for frag in o_molsys.fragments:
+        for intco in frag.intcos:
+            if isinstance(intco, tors.Tors):
+                intco.freeze()
+
+        # Seems better to just unfreeze at end rather than check whether each torsion is in this list
+        for indices in skipped_torsions:
+            new_tors = tors.Tors(indices[0] - 1, indices[1] - 1, indices[2] - 1, indices[3] - 1)
+            try:
+                index = frag.intcos.index(new_tors)
+                frag.intcos[index].unfreeze()
+            except ValueError:
+                logger.info("User requested free torsion not present, so adding it.\n")
+                frag.intcos.append(new_tors)
 
 def ranged_tors_from_input(ranged_tors_list, o_molsys):
     """
@@ -881,44 +898,48 @@ def freeze_intrafrag(o_molsys):
         F.freeze()
 
 
-def add_constrained_intcos(o_molsys):
+def add_constrained_intcos(o_molsys, params: op.OptParams):
+
+    if params.freeze_all_dihedrals:
+        freeze_all_torsions(o_molsys, params.skip_freezing_dihedrals)
+
     # Frozen coordinates
-    if op.Params.frozen_distance:
-        frozen_stre_from_input(op.Params.frozen_distance, o_molsys)
-    if op.Params.frozen_bend:
-        frozen_bend_from_input(op.Params.frozen_bend, o_molsys)
-    if op.Params.frozen_dihedral:
-        frozen_tors_from_input(op.Params.frozen_dihedral, o_molsys)
-    if op.Params.frozen_oofp:
-        frozen_oofp_from_input(op.Params.frozen_oofp, o_molsys)
-    if op.Params.frozen_cartesian:
-        frozen_cart_from_input(op.Params.frozen_cartesian, o_molsys)
+    if params.frozen_distance:
+        frozen_stre_from_input(params.frozen_distance, o_molsys)
+    if params.frozen_bend:
+        frozen_bend_from_input(params.frozen_bend, o_molsys)
+    if params.frozen_dihedral:
+        frozen_tors_from_input(params.frozen_dihedral, o_molsys)
+    if params.frozen_oofp:
+        frozen_oofp_from_input(params.frozen_oofp, o_molsys)
+    if params.frozen_cartesian:
+        frozen_cart_from_input(params.frozen_cartesian, o_molsys)
 
     # Ranged coordinates
-    if op.Params.ranged_distance:
-        ranged_stre_from_input(op.Params.ranged_distance, o_molsys)
-    if op.Params.ranged_bend:
-        ranged_bend_from_input(op.Params.ranged_bend, o_molsys)
-    if op.Params.ranged_dihedral:
-        ranged_tors_from_input(op.Params.ranged_dihedral, o_molsys)
-    if op.Params.ranged_oofp:
-        ranged_oofp_from_input(op.Params.ranged_oofp, o_molsys)
-    if op.Params.ranged_cartesian:
-        ranged_cart_from_input(op.Params.ranged_cartesian, o_molsys)
+    if params.ranged_distance:
+        ranged_stre_from_input(params.ranged_distance, o_molsys)
+    if params.ranged_bend:
+        ranged_bend_from_input(params.ranged_bend, o_molsys)
+    if params.ranged_dihedral:
+        ranged_tors_from_input(params.ranged_dihedral, o_molsys)
+    if params.ranged_oofp:
+        ranged_oofp_from_input(params.ranged_oofp, o_molsys)
+    if params.ranged_cartesian:
+        ranged_cart_from_input(params.ranged_cartesian, o_molsys)
 
     # Coordinates with extra forces
-    if op.Params.ext_force_distance:
-        ext_force_stre_from_input(op.Params.ext_force_distance, o_molsys)
-    if op.Params.ext_force_bend:
-        ext_force_bend_from_input(op.Params.ext_force_bend, o_molsys)
-    if op.Params.ext_force_dihedral:
-        ext_force_tors_from_input(op.Params.ext_force_dihedral, o_molsys)
-    if op.Params.ext_force_oofp:
-        ext_force_oofp_from_input(op.Params.ext_force_oofp, o_molsys)
-    if op.Params.ext_force_cartesian:
-        ext_force_cart_from_input(op.Params.ext_force_cartesian, o_molsys)
-
-    if op.Params.freeze_intrafrag:
+    if params.ext_force_distance:
+        ext_force_stre_from_input(params.ext_force_distance, o_molsys)
+    if params.ext_force_bend:
+        ext_force_bend_from_input(params.ext_force_bend, o_molsys)
+    if params.ext_force_dihedral:
+        ext_force_tors_from_input(params.ext_force_dihedral, o_molsys)
+    if params.ext_force_oofp:
+        ext_force_oofp_from_input(params.ext_force_oofp, o_molsys)
+    if params.ext_force_cartesian:
+        ext_force_cart_from_input(params.ext_force_cartesian, o_molsys)
+
+    if params.freeze_intrafrag:
         freeze_intrafrag(o_molsys)
 
 

optking/optimize.py

@@ -643,7 +643,7 @@ def make_internal_coords(o_molsys: Molsys, params: op.OptParams):
         if params.opt_coordinates in ["CARTESIAN", "BOTH"]:
             for F in o_molsys.fragments:
                 F.add_cartesian_intcos()
-    addIntcos.add_constrained_intcos(o_molsys)  # make sure these are in the set
+    addIntcos.add_constrained_intcos(o_molsys, params)  # make sure these are in the set
     return
 
 

optking/optparams.py

@@ -220,6 +220,12 @@ def __init__(self, uod):
         force = uod.get("EXT_FORCE_CARTESIAN", "")
         self.ext_force_cartesian = int_xyz_fx_string(force, 1)
 
+        # Ensures all torsions can be frozen
+        self.freeze_all_dihedrals = uod.get("FREEZE_ALL_DIHEDRALS", False)
+        # skip_dihedrals is only meant for use with all_dihedrals_frozen
+        skip_dihedrals = uod.get("SKIP_FREEZING_DIHEDRALS", "")
+        self.skip_freezing_dihedrals = int_list(tokenize_input_string(skip_dihedrals), 4)
+
         # Should an xyz trajectory file be kept (useful for visualization)?
         # P.print_trajectory_xyz = uod.get('PRINT_TRAJECTORY_XYZ', False)
         # Symmetry tolerance for testing whether a mode is symmetric.

optking/simple.py

@@ -39,10 +39,8 @@ def ranged(self):
         return self.constraint == "ranged"
 
     def freeze(self):
-        if self.constraint == "free":
-            self.constraint = "frozen"
+        self.constraint = "frozen"
         # for now if 'ranged', don't change
-        return
 
     def unfreeze(self):
         if self.constraint == "frozen":

optking/tests/test_frozen_internals.py

@@ -43,3 +43,30 @@ def test_frozen_coords(option, expected, num_steps, check_iter):
 
     assert psi4.compare_values(expected, thisenergy, 6)  # TEST
     utils.compare_iterations(json_output, num_steps, check_iter)
+
+@pytest.mark.parametrize(
+    "option, expected, num_steps", [
+        ({"freeze_all_dihedrals": True}, -154.09270490, 5),
+        ({"freeze_all_dihedrals": True, "skip_freezing_dihedrals": ""}, -154.09270490, 5),
+        ({"freeze_all_dihedrals": True, "skip_freezing_dihedrals": "7 2 6 9 8 2 6 9"}, -154.09270565, 6)]
+)
+def test_freezing_torsions(option, expected, num_steps, check_iter):
+
+    ethanol = psi4.geometry(
+        """
+        0 1
+        C     -2.645200    0.725834    0.013523
+        C     -1.131382    0.457496    0.011963
+        H     -2.885160    1.576571   -0.628480
+        H     -3.184240   -0.148652   -0.357146
+        H     -2.995284    0.943570    1.025249
+        O     -0.442426    1.614674    0.486697
+        H     -0.908114   -0.393441    0.661320
+        H     -0.800148    0.225843   -1.004116
+        H      0.514552    1.429456    0.481888
+        """
+    )
+    psi4.core.clean_options()
+    psi4.set_options({"basis": "cc-PVDZ", "scf_type": "pk", "g_convergence": "gau_tight"})
+    json_output = optking.optimize_psi4("hf", **option)
+