Adds option to freeze all dihedral angles in molecule

optking/addIntcos.py

@@ -939,6 +939,39 @@ def add_constrained_intcos(o_molsys):
     if op.Params.freeze_intrafrag:
         freeze_intrafrag(o_molsys)
 
+    if op.Params.freeze_all_dihedrals:
+        freeze_all_torsions(o_molsys, op.Params.unfreeze_dihedrals)
+
+def freeze_all_torsions(molsys, skipped_torsions=[]):
+    """ Freeze all intrafragment torsions.
+
+    Parameters
+    ----------
+    molsys: molsys.Molsys
+    skipped_torsions: list[list[int]]
+        each set of integers should be four integers denoting a dihedral angle / torsion
+    """
+
+    for frag in 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 index_set in skipped_torsions:
+        atoms = [index - 1 for index in index_set]
+        f = check_fragment(atoms, molsys)
+        new_tors = tors.Tors(*atoms)
+
+        try:
+            index = molsys.fragments[f].intcos.index(new_tors)
+            frag.intcos[index].unfreeze()
+        except ValueError:
+            logger.info(
+                "dihedral angle %s was unfrozen but was not present - adding it.",
+                new_tors
+            )
+            frag.intcos.append(new_tors)
 
 def add_dimer_frag_intcos(o_molsys):
     # Look for coordinates in the following order:

optking/optparams.py

@@ -190,6 +190,11 @@ def __init__(self, uod):
         frozen = uod.get("FROZEN_CARTESIAN", "")
         self.frozen_cartesian = int_xyz_float_list(tokenize_input_string(frozen), 1, 1, 0)
 
+        self.freeze_all_dihedrals = uod.get("FREEZE_ALL_DIHEDRALS", False)
+        # For use only with `freeze_all_dihedrals` unfreeze a small subset of dihedrals
+        frozen = uod.get("UNFREEZE_DIHEDRALS", "")
+        self.unfreeze_dihedrals = int_list(tokenize_input_string(frozen), 4)
+
         # Specify distance between atoms to be ranged
         ranged = uod.get("RANGED_DISTANCE", "")
         self.ranged_distance = int_float_list(tokenize_input_string(ranged), 2, 2)

optking/tests/test_frozen_internals.py

@@ -15,11 +15,17 @@
 f2 = {"frozen_bend": "1 2 3 2 3 4"}
 f3 = {"frozen_dihedral": "1 2 3 4"}
 
-optking__frozen_coords = [(f1, OH_frozen_stre_rhf, 9), (f2, OOH_frozen_bend_rhf, 6), (f3, HOOH_frozen_dihedral_rhf, 6)]
+optking__frozen_coords = [
+    (f1, OH_frozen_stre_rhf, 9),
+    (f2, OOH_frozen_bend_rhf, 6),
+    (f3, HOOH_frozen_dihedral_rhf, 6),
+]
 
 
 @pytest.mark.parametrize(
-    "option, expected, num_steps", optking__frozen_coords, ids=["frozen stretch", "frozen bend", "frozen dihedral"]
+    "option, expected, num_steps",
+    optking__frozen_coords,
+    ids=["frozen stretch", "frozen bend", "frozen dihedral"],
 )
 def test_frozen_coords(option, expected, num_steps, check_iter):
     # Constrained minimization with frozen bond, bend, and torsion
@@ -34,7 +40,13 @@ def test_frozen_coords(option, expected, num_steps, check_iter):
 
     psi4.core.clean_options()
 
-    psi4_options = {"diis": "false", "basis": "cc-PVDZ", "scf_type": "pk", "print": 4, "g_convergence": "gau_tight"}
+    psi4_options = {
+        "diis": "false",
+        "basis": "cc-PVDZ",
+        "scf_type": "pk",
+        "print": 4,
+        "g_convergence": "gau_tight",
+    }
     psi4.set_options(psi4_options)
     psi4.set_options(option)
 
@@ -43,3 +55,116 @@ 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)
+
+
+def test_butane_frozen(check_iter):
+    _ = psi4.geometry("pubchem:butane")
+
+    psi4.core.clean_options()
+    psi4_options = {
+        "basis": "6-31G",
+        "g_convergence": "gau_tight",
+        "freeze_all_dihedrals": True,
+    }
+    psi4.set_options(psi4_options)
+
+    result = optking.optimize_psi4("scf")
+    E1 = result["energies"][-1]  # TEST
+
+    psi4.core.clean_options()
+    psi4_options = {
+        "basis": "6-31G",
+        "g_convergence": "gau_tight",
+        "frozen_dihedral": """
+            1 2 4 12
+            1 2 4 13
+            1 2 4 14
+            2 1 3 9
+            2 1 3 10
+            2 1 3 11
+            3 1 2 4
+            3 1 2 7
+            3 1 2 8
+            4 2 1 5
+            4 2 1 6
+            5 1 2 7
+            5 1 2 8
+            5 1 3 9
+            5 1 3 10
+            5 1 3 11
+            6 1 2 7
+            6 1 2 8
+            6 1 3 9
+            6 1 3 10
+            6 1 3 11
+            7 2 4 12
+            7 2 4 13
+            7 2 4 14
+            8 2 4 12
+            8 2 4 13
+            8 2 4 14
+        """
+    }
+    psi4.set_options(psi4_options)
+    result = optking.optimize_psi4("scf")
+    E2 = result["energies"][-1]
+
+    assert psi4.compare_values(E1, E2, 8, "RHF energy")  # TEST
+    utils.compare_iterations(result, 5, check_iter)
+
+def test_butane_skip_frozen(check_iter):
+    _ = psi4.geometry("pubchem:butane")
+
+    psi4.core.clean_options()
+    psi4_options = {
+        "basis": "6-31G",
+        "g_convergence": "gau_tight",
+        "freeze_all_dihedrals": True,
+        "unfreeze_dihedrals": """
+            8 2 4 12
+            8 2 4 13
+            8 2 4 14
+            3 1 2 8
+            5 1 2 8
+            6 1 2 8
+        """
+    }
+    psi4.set_options(psi4_options)
+
+    result = optking.optimize_psi4("scf")
+    E1 = result["energies"][-1]  # TEST
+
+    psi4.core.clean_options()
+    psi4_options = {
+        "basis": "6-31G",
+        "g_convergence": "gau_tight",
+        "frozen_dihedral": """
+            1 2 4 12
+            1 2 4 13
+            1 2 4 14
+            2 1 3 9
+            2 1 3 10
+            2 1 3 11
+            3 1 2 4
+            3 1 2 7
+            4 2 1 5
+            4 2 1 6
+            5 1 2 7
+            5 1 3 9
+            5 1 3 10
+            5 1 3 11
+            6 1 2 7
+            6 1 3 9
+            6 1 3 10
+            6 1 3 11
+            7 2 4 12
+            7 2 4 13
+            7 2 4 14
+        """
+    }
+    psi4.set_options(psi4_options)
+    result = optking.optimize_psi4("scf")
+    E2 = result["energies"][-1]
+
+    assert psi4.compare_values(E1, E2, 8, "RHF energy")  # TEST
+    utils.compare_iterations(result, 5, check_iter)