From 526d888d9198b450a654a5a9de03429a5c6f2e95 Mon Sep 17 00:00:00 2001 From: ehhartmann <69237857+ehhartmann@users.noreply.github.com> Date: Tue, 27 Feb 2024 19:00:34 +0100 Subject: [PATCH] fix: integration tests (#390) * expect 12 files because of 0_setup * fix: impropers were not correctly removed in break_bond * fuck this bug * remove debug lines --------- Co-authored-by: Eric Hartmann --- src/kimmdy/topology/topology.py | 33 ++++++++++++++++++++------------- tests/test_integration.py | 4 ++-- 2 files changed, 22 insertions(+), 15 deletions(-) diff --git a/src/kimmdy/topology/topology.py b/src/kimmdy/topology/topology.py index c46ddf85..5b288160 100644 --- a/src/kimmdy/topology/topology.py +++ b/src/kimmdy/topology/topology.py @@ -408,8 +408,11 @@ def _get_atom_improper_dihedrals( ai = atom_nr partners = self.atoms[ai].bound_to_nrs if len(partners) >= 3: - for aj, ak, al in combinations(partners, 3): - dihedral_candidate_keys.append((ai, aj, ak, al)) + combs = permutations(partners, 3) + for comb in combs: + aj, ak, al = comb + # center atom is at position 3 + dihedral_candidate_keys.append((aj, ak, ai, al)) # atom in corner of a triangle: for a in self.atoms[atom_nr].bound_to_nrs: @@ -1130,12 +1133,12 @@ def break_bond( reactive_moleculetype.pairs.pop(pairkey, None) # and improper dihedrals - dihedral_k_v = reactive_moleculetype._get_atom_improper_dihedrals( - atompair_nrs[0], self.ff - ) + reactive_moleculetype._get_atom_improper_dihedrals(atompair_nrs[1], self.ff) - for key, _ in dihedral_k_v: + keys_remove = [] + for key in reactive_moleculetype.improper_dihedrals: if all([x in key for x in atompair_nrs]): - reactive_moleculetype.improper_dihedrals.pop(key, None) + keys_remove.append(key) + for key in keys_remove: + reactive_moleculetype.improper_dihedrals.pop(key, None) # warn if atom is part of restraint for atom in atompair: @@ -1221,7 +1224,7 @@ def bind_bond( for a in reactive_moleculetype.atoms.values() if a.resnr == other_atom.resnr ] - type_set = False + name_set = False for key, bond in aa.bonds.items(): if other_atom.atom in key and any( k.upper().startswith("H") for k in key @@ -1237,21 +1240,25 @@ def bind_bond( or h_name == atom.atom ): atom.atom = h_name + name_set = True else: atom.atom = "HX" - type_set = True + name_set = True continue - logging.info(f"HAT hydrogen will be bound to {other_atom}.") + logger.info(f"HAT hydrogen will be bound to {other_atom}.") break else: - if type_set: - logging.warning( + if name_set: + logger.warning( "Found new atomtype but not atomname for HAT hydrogen!" ) else: - logging.warning( + logger.warning( "Found neither new atomtype nor atomname for HAT hydrogen!" ) + if not name_set: + atom.atom = "HX" + logger.warning(f"Named newly bonded hydrogen 'HX'") # update bound_to atompair[0].bound_to_nrs.append(atompair[1].nr) diff --git a/tests/test_integration.py b/tests/test_integration.py index 36724d71..4cf9acec 100644 --- a/tests/test_integration.py +++ b/tests/test_integration.py @@ -104,7 +104,7 @@ def test_integration_homolysis_reaction(arranged_tmp_path): def test_integration_pull(arranged_tmp_path): kimmdy_run() assert "Finished running tasks" in read_last_line(Path("kimmdy.log")) - assert len(list(Path.cwd().glob("kimmdy_001/*"))) == 11 + assert len(list(Path.cwd().glob("kimmdy_001/*"))) == 12 @pytest.mark.require_grappa @@ -117,4 +117,4 @@ def test_integration_pull(arranged_tmp_path): def test_integration_whole_run(arranged_tmp_path): kimmdy_run() assert "Finished running tasks" in read_last_line(Path("kimmdy.log")) - assert len(list(Path.cwd().glob("*/*"))) == 23 + assert len(list(Path.cwd().glob("kimmdy_001/*"))) == 25