Skip to content

Commit

Permalink
fix #216: angle parameters are not being read from topology file (#264)
Browse files Browse the repository at this point in the history
* fix #216: angle parameters are not being read from topology file
* refactored heavily to support all angle types
* add tests (with pytest.MonkeyPatch)
* update CHANGES
  • Loading branch information
jandom authored Nov 9, 2023
1 parent fa64735 commit 691393c
Show file tree
Hide file tree
Showing 5 changed files with 216 additions and 100 deletions.
4 changes: 3 additions & 1 deletion CHANGES
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
29 changes: 29 additions & 0 deletions gromacs/fileformats/blocks.py
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@
"""

import logging
from enum import IntEnum


class System(object):
Expand Down Expand Up @@ -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)
Expand Down
190 changes: 93 additions & 97 deletions gromacs/fileformats/top.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------------------------------
Expand All @@ -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"):
"""
Expand Down Expand Up @@ -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"]

Expand Down Expand Up @@ -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"],
Expand Down
4 changes: 2 additions & 2 deletions setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
89 changes: 89 additions & 0 deletions tests/fileformats/top/test_TOP.py
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 691393c

Please sign in to comment.