From 53ea698b32454e5193a0efb8b9198e01991171b5 Mon Sep 17 00:00:00 2001 From: Alexander Date: Thu, 6 Jun 2024 15:32:01 -0400 Subject: [PATCH] Adds option to freeze all torsions and then unfreeze a subset --- optking/addIntcos.py | 87 ++++++++++++++++---------- optking/optimize.py | 2 +- optking/optparams.py | 6 ++ optking/simple.py | 4 +- optking/tests/test_frozen_internals.py | 27 ++++++++ 5 files changed, 89 insertions(+), 37 deletions(-) diff --git a/optking/addIntcos.py b/optking/addIntcos.py index 85638bb..0eea54e 100644 --- a/optking/addIntcos.py +++ b/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) diff --git a/optking/optimize.py b/optking/optimize.py index 3562b6c..c1e2ffd 100644 --- a/optking/optimize.py +++ b/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 diff --git a/optking/optparams.py b/optking/optparams.py index 930da39..cf56540 100644 --- a/optking/optparams.py +++ b/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. diff --git a/optking/simple.py b/optking/simple.py index e66b75e..991212a 100644 --- a/optking/simple.py +++ b/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": diff --git a/optking/tests/test_frozen_internals.py b/optking/tests/test_frozen_internals.py index 9eaca28..c3f9254 100644 --- a/optking/tests/test_frozen_internals.py +++ b/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) +