diff --git a/CHANGES b/CHANGES index 2b685dcd..abb880dd 100644 --- a/CHANGES +++ b/CHANGES @@ -3,9 +3,11 @@ ============================== 2023-xx-xx 0.9.0 -orbeckst +orbeckst, jandom * removed support for legacy Python (<= 3.7) (#259) +* fixed GROMACS TOP reader not reading angle parameters from topology + file (#261) 2023-09-16 0.8.5 diff --git a/gromacs/fileformats/blocks.py b/gromacs/fileformats/blocks.py index c52ba1f8..d9a5042b 100644 --- a/gromacs/fileformats/blocks.py +++ b/gromacs/fileformats/blocks.py @@ -56,6 +56,7 @@ """ import logging +from enum import IntEnum class System(object): @@ -441,6 +442,34 @@ def __eq__(self, other): ) +class AngleFunctionType(IntEnum): + HARMONIC = 1 + G96_ANGLE = 2 + CROSS_BOND_BOND = 3 + CROSS_BOND_ANGLE = 4 + UREY_BRADLEY = 5 + QUARTIC_ANGLE = 6 + # There's no function type 7 as per the given code comments + TABULATED_ANGLE = 8 + LINEAR_ANGLE = 9 + RESTRICTED_BENDING = 10 + + @property + def num_params(self): + # Define the number of parameters expected for each function type + return { + self.HARMONIC: 2, + self.G96_ANGLE: 2, + self.CROSS_BOND_BOND: 3, + self.CROSS_BOND_ANGLE: 4, + self.UREY_BRADLEY: 4, + self.QUARTIC_ANGLE: 6, + self.TABULATED_ANGLE: 2, # This might be different depending on your implementation + self.LINEAR_ANGLE: 2, + self.RESTRICTED_BENDING: 2, + }[self] + + class AngleType(Param): def __init__(self, format): super(AngleType, self).__init__(format) diff --git a/gromacs/fileformats/top.py b/gromacs/fileformats/top.py index fce89983..29e087c0 100644 --- a/gromacs/fileformats/top.py +++ b/gromacs/fileformats/top.py @@ -470,7 +470,7 @@ def _add_info(sys_or_mol, section, container): else: raise NotImplementedError - elif curr_sec in ("angletypes", "angles"): + elif curr_sec in {"angletypes", "angles"}: """ section #at fu #param ---------------------------------- @@ -482,110 +482,101 @@ def _add_info(sys_or_mol, section, container): angles 3 6 6 angles 3 8 ?? """ - ai, aj, ak = fields[:3] - fu = int(fields[3]) - assert fu in (1, 2, 3, 4, 5, 6, 8) # no 7 - if fu not in (1, 2, 5): - raise NotImplementedError( - "function {0:d} is not yet supported".format(fu) + fu = blocks.AngleFunctionType(int(fields[3])) + + if len(fields[4:]) != 0 and len(fields[4:]) != fu.num_params: + raise ValueError( + "Expected {num_params} parameters for function type {fu}, got {len(fields[4:])}".format( + num_params=fu.num_params, fu=fu + ) ) ang = blocks.AngleType("gromacs") - if fu == 1: - if curr_sec == "angletypes": - ang.atype1 = ai - ang.atype2 = aj - ang.atype3 = ak - - tetha0, ktetha = list(map(float, fields[4:6])) - ang.gromacs = { - "param": { - "ktetha": ktetha, - "tetha0": tetha0, - "kub": None, - "s0": None, - }, - "func": fu, + ang.gromacs = {"func": fu} + + if curr_sec == "angletypes": + ang.atype1 = ai + ang.atype2 = aj + ang.atype3 = ak + elif curr_sec == "angles": + ai, aj, ak = list(map(int, [ai, aj, ak])) + ang.atom1 = mol.atoms[ai - 1] + ang.atom2 = mol.atoms[aj - 1] + ang.atom3 = mol.atoms[ak - 1] + + # Parse parameters based on the function type + params = list(map(float, fields[4:])) + + # Handle parameters based on function types + if fu.num_params == len(params): + if fu in { + blocks.AngleFunctionType.HARMONIC, + blocks.AngleFunctionType.G96_ANGLE, + }: + ang.gromacs["param"] = { + "tetha0": params[0], + "ktetha": params[1], } - - self.angletypes.append(ang) - _add_info(self, curr_sec, self.angletypes) - - elif curr_sec == "angles": - ai, aj, ak = list(map(int, [ai, aj, ak])) - ang.atom1 = mol.atoms[ai - 1] - ang.atom2 = mol.atoms[aj - 1] - ang.atom3 = mol.atoms[ak - 1] - ang.gromacs["func"] = fu - - mol.angles.append(ang) - _add_info(mol, curr_sec, mol.angles) - - else: - raise ValueError - - elif fu == 2: - if curr_sec == "angletypes": - raise NotImplementedError() - - elif curr_sec == "angles": - ai, aj, ak = list(map(int, [ai, aj, ak])) - ang.atom1 = mol.atoms[ai - 1] - ang.atom2 = mol.atoms[aj - 1] - ang.atom3 = mol.atoms[ak - 1] - ang.gromacs["func"] = fu - - tetha0, ktetha = list(map(float, fields[4:6])) - ang.gromacs = { - "param": { - "ktetha": ktetha, - "tetha0": tetha0, - "kub": None, - "s0": None, - }, - "func": fu, + elif fu == blocks.AngleFunctionType.CROSS_BOND_BOND: + ang.gromacs["param"] = { + "r1e": params[0], + "r2e": params[1], + "krrprime": params[2], } - - mol.angles.append(ang) - _add_info(mol, curr_sec, mol.angles) - - elif fu == 5: - if curr_sec == "angletypes": - ang.atype1 = ai - ang.atype2 = aj - ang.atype3 = ak - tetha0, ktetha, s0, kub = list(map(float, fields[4:8])) - - ang.gromacs = { - "param": { - "ktetha": ktetha, - "tetha0": tetha0, - "kub": kub, - "s0": s0, - }, - "func": fu, + elif fu == blocks.AngleFunctionType.CROSS_BOND_ANGLE: + ang.gromacs["param"] = { + "r1e": params[0], + "r2eprime": params[1], + "r3e": params[2], + "krtheta": params[3], + } + elif fu == blocks.AngleFunctionType.UREY_BRADLEY: + ang.gromacs["param"] = { + "tetha0": params[0], + "ktetha": params[1], + "s0": params[2], + "kub": params[3], + } + elif fu == blocks.AngleFunctionType.QUARTIC_ANGLE: + ang.gromacs["param"] = { + "tetha0": params[0], + "C1": params[1], + "C2": params[2], + "C3": params[3], + "C4": params[4], + "C5": params[5], + } + elif fu == blocks.AngleFunctionType.TABULATED_ANGLE: + ang.gromacs["param"] = { + "table_number": params[0], + "k": params[1], + } # Assuming 'table number' is a parameter here + elif fu == blocks.AngleFunctionType.LINEAR_ANGLE: + ang.gromacs["param"] = { + "a0": params[0], + "klin": params[1], + } + elif fu == blocks.AngleFunctionType.RESTRICTED_BENDING: + ang.gromacs["param"] = { + "tetha0": params[0], + "ktheta": params[1], } - - self.angletypes.append(ang) - _add_info(self, curr_sec, self.angletypes) - - elif curr_sec == "angles": - ai, aj, ak = list(map(int, [ai, aj, ak])) - ang.atom1 = mol.atoms[ai - 1] - ang.atom2 = mol.atoms[aj - 1] - ang.atom3 = mol.atoms[ak - 1] - ang.gromacs["func"] = fu - - mol.angles.append(ang) - _add_info(mol, curr_sec, mol.angles) - else: - raise ValueError + raise NotImplementedError( + "Function type {fu} is not implemented".forma(fu=fu) + ) + # Add the angle to the appropriate list and call _add_info + if curr_sec == "angletypes": + self.angletypes.append(ang) + _add_info(self, curr_sec, self.angletypes) + elif curr_sec == "angles": + mol.angles.append(ang) + _add_info(mol, curr_sec, mol.angles) else: - raise NotImplementedError + raise ValueError("Unknown section while parsing angles") elif curr_sec in ("dihedraltypes", "dihedrals"): """ @@ -1230,8 +1221,8 @@ def _make_angletypes(self, m): ktetha = ang.gromacs["param"]["ktetha"] tetha0 = ang.gromacs["param"]["tetha0"] - kub = ang.gromacs["param"]["kub"] - s0 = ang.gromacs["param"]["s0"] + kub = ang.gromacs["param"].get("kub") + s0 = ang.gromacs["param"].get("s0") fu = ang.gromacs["func"] @@ -1409,7 +1400,12 @@ def _make_angles(self, m): result = [] for ang in m.angles: fu = ang.gromacs["func"] - if ang.gromacs["param"]["ktetha"] and ang.gromacs["param"]["tetha0"]: + has_params = ( + "param" in ang.gromacs + and ang.gromacs["param"]["ktetha"] + and ang.gromacs["param"]["tetha0"] + ) + if has_params: ktetha, tetha0 = ( ang.gromacs["param"]["ktetha"], ang.gromacs["param"]["tetha0"], diff --git a/setup.py b/setup.py index 713de0f1..a597758e 100644 --- a/setup.py +++ b/setup.py @@ -63,9 +63,9 @@ ], }, install_requires=[ - "numpy>=1.0", - "numkit", # numerical helpers "matplotlib", + "numkit", # numerical helpers + "numpy>=1.0", ], tests_require=["pytest", "pandas>=0.17"], zip_safe=True, diff --git a/tests/fileformats/top/test_TOP.py b/tests/fileformats/top/test_TOP.py new file mode 100644 index 00000000..cb28a0d6 --- /dev/null +++ b/tests/fileformats/top/test_TOP.py @@ -0,0 +1,89 @@ +import gromacs as gmx +import io +import random +import pytest +from gromacs.fileformats import blocks + + +# This class simulates a file object +class MockFile(io.StringIO): + def __init__(self, text): + super(MockFile, self).__init__(text) + + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_val, exc_tb): + self.close() + return False + + +def generate_topology_line(func_type): + # Generates random parameters for the different function types + def random_params(num): + return [random.uniform(0.0, 100.0) for _ in range(num)] + + # Define the pattern for each function type + patterns = { + blocks.AngleFunctionType.HARMONIC: "{:5d}{:5d}{:5d} 1{:10.4f}{:10.4f}", + blocks.AngleFunctionType.G96_ANGLE: "{:5d}{:5d}{:5d} 2{:10.4f}{:10.4f}", + blocks.AngleFunctionType.CROSS_BOND_BOND: "{:5d}{:5d}{:5d} 3{:10.4f}{:10.4f}{:10.4f}", + blocks.AngleFunctionType.CROSS_BOND_ANGLE: "{:5d}{:5d}{:5d} 4{:10.4f}{:10.4f}{:10.4f}{:10.4f}", + blocks.AngleFunctionType.UREY_BRADLEY: "{:5d}{:5d}{:5d} 5{:10.4f}{:10.4f}{:10.4f}{:10.4f}", + blocks.AngleFunctionType.QUARTIC_ANGLE: "{:5d}{:5d}{:5d} 6{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}{:10.4f}", + blocks.AngleFunctionType.TABULATED_ANGLE: "{:5d}{:5d}{:5d} 8{:10.4f}{:10.4f}", + blocks.AngleFunctionType.LINEAR_ANGLE: "{:5d}{:5d}{:5d} 9{:10.4f}{:10.4f}", + blocks.AngleFunctionType.RESTRICTED_BENDING: "{:5d}{:5d}{:5d} 10{:10.4f}{:10.4f}", + } + + atoms = [1, 2, 3] + params = random_params(func_type.num_params) + line = patterns[func_type].format(*(atoms + params)) + return line + + +def create_topology_data(func_type): + """ + https://manual.gromacs.org/current/reference-manual/topologies/topology-file-formats.html#tab-topfile2 + """ + line = generate_topology_line(func_type) + return """ +[ moleculetype ] +; Name nrexcl +Example 3 + +[ atoms ] +; nr type resnr residue atom cgnr charge mass typeB cha +; residue 1 ASN rtp ASN q +1.0 + 1 NH3 1 ASN N 1 -0.3 14.007 ; qtot -0.3 + 2 HC 1 ASN H1 2 0.33 1.008 ; qtot 0.03 + 3 HC 1 ASN H2 3 0.33 1.008 ; qtot 0.36 + +[ angles ] +; ai aj ak funct c0 c1 c2 c3 +{line} + +[ molecules ] +; Compound #mols +Example 1 +""".format( + line=line + ) + + +@pytest.mark.parametrize("func_type", blocks.AngleFunctionType) +def test_angles(func_type, monkeypatch): + # Create a custom function that will replace 'open' + def mock_open(*args, **kwargs): + if args[0] == "topol.top": + return MockFile(create_topology_data(func_type)) + else: + return open(*args, **kwargs) + + # Use monkeypatch to replace 'open' with 'mock_open' + monkeypatch.setattr("builtins.open", mock_open) + + topol = gmx.fileformats.top.TOP("topol.top") + [molecule] = topol.molecules + [angle] = molecule.angles + assert len(angle.gromacs["param"].keys()) == func_type.num_params