From 3352365e4a4b40f06d5a98dce45c6ef62905b370 Mon Sep 17 00:00:00 2001 From: Alexander Heide Date: Wed, 24 Jul 2024 13:32:31 -0400 Subject: [PATCH] Adds option to freeze all dihedral angles in molecule --- optking/addIntcos.py | 33 +++++++ optking/optparams.py | 5 + optking/tests/test_frozen_internals.py | 131 ++++++++++++++++++++++++- 3 files changed, 166 insertions(+), 3 deletions(-) diff --git a/optking/addIntcos.py b/optking/addIntcos.py index d73a74e..dfa0f0b 100644 --- a/optking/addIntcos.py +++ b/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: diff --git a/optking/optparams.py b/optking/optparams.py index 74cdbcc..965ab23 100644 --- a/optking/optparams.py +++ b/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) diff --git a/optking/tests/test_frozen_internals.py b/optking/tests/test_frozen_internals.py index 9eaca28..69c43b7 100644 --- a/optking/tests/test_frozen_internals.py +++ b/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)