Skip to content
Snippets Groups Projects
Commit d4fabd07 authored by James Graham's avatar James Graham
Browse files

Add test for itp

parent 384cfaee
Branches
No related tags found
No related merge requests found
...@@ -226,26 +226,29 @@ class BondSet: ...@@ -226,26 +226,29 @@ class BondSet:
e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, mol),) e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, mol),)
raise raise
def write_itp(self, filename, mapping): def write_itp(self, filename, mapping, exclude=set()):
""" """
Output a GROMACS .itp file containing atoms/beads and bonded terms. Output a GROMACS .itp file containing atoms/beads and bonded terms.
:param filename: Name of output file :param filename: Name of output file
:param mapping: AA->CG Mapping from which to collect bead properties :param mapping: AA->CG Mapping from which to collect bead properties
:param exclude: Set of molecule names to be excluded from itp
""" """
self._populate_atom_numbers(mapping) self._populate_atom_numbers(mapping)
backup_file(filename, verbose=True) backup_file(filename, verbose=True)
def write_bond_angle_dih(bonds, section_header, itp, fconst=True): def write_bond_angle_dih(bonds, section_header, itp, fconst=True, multiplicity=None):
if bonds: if bonds:
print("\n[ {0:s} ]".format(section_header), file=itp) print("\n[ {0:s} ]".format(section_header), file=itp)
for bond in bonds: for bond in bonds:
# Factor is usually 1, unless doing correction # Factor is usually 1, unless doing correction
eqm = bond.eqm * self._empirical_correction_factor eqm = bond.eqm * self._empirical_correction_factor
line = " ".join(["{0:4d}".format(atnum + 1) for atnum in bond.atom_numbers]) line = " ".join(["{0:4d}".format(atnum + 1) for atnum in bond.atom_numbers])
line += " {0:12.5f}".format(eqm) line += " {0:4d} {1:12.5f}".format(1, eqm)
if fconst: if fconst:
line += " {0:12.5f}".format(bond.fconst) line += " {0:12.5f}".format(bond.fconst)
if multiplicity is not None:
line += " {0:4d}".format(multiplicity)
print(line, file=itp) print(line, file=itp)
with open(filename, "w") as itp: with open(filename, "w") as itp:
...@@ -258,7 +261,7 @@ class BondSet: ...@@ -258,7 +261,7 @@ class BondSet:
print(header, file=itp) print(header, file=itp)
# Print molecule # Print molecule
for mol in self._molecules: for mol in filter(lambda mol: mol not in exclude, self._molecules):
print("\n[ moleculetype ]", file=itp) print("\n[ moleculetype ]", file=itp)
print("{0:4s} {1:4d}".format(mol, 1), file=itp) print("{0:4s} {1:4d}".format(mol, 1), file=itp)
...@@ -271,7 +274,7 @@ class BondSet: ...@@ -271,7 +274,7 @@ class BondSet:
write_bond_angle_dih(self.get_bond_lengths(mol), "bonds", itp) write_bond_angle_dih(self.get_bond_lengths(mol), "bonds", itp)
write_bond_angle_dih(self.get_bond_angles(mol), "angles", itp) write_bond_angle_dih(self.get_bond_angles(mol), "angles", itp)
write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", itp) write_bond_angle_dih(self.get_bond_dihedrals(mol), "dihedrals", itp, multiplicity=1)
write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", itp, fconst=False) write_bond_angle_dih(self.get_bond_length_constraints(mol), "constraints", itp, fconst=False)
def apply(self, frame): def apply(self, frame):
......
;
; Topology prepared automatically using PyCGTOOL
; James Graham <J.A.Graham@soton.ac.uk> 2016
; University of Southampton
; https://github.com/jag1g13/pycgtool
;
[ moleculetype ]
ALLA 1
[ atoms ]
1 P3 1 ALLA C1 1 0.000
2 P3 1 ALLA C2 2 0.000
3 P3 1 ALLA C3 3 0.000
4 P3 1 ALLA C4 4 0.000
5 P2 1 ALLA C5 5 0.000
6 P4 1 ALLA O5 6 0.000
[ bonds ]
2 3 1 0.21998 35162.10022
3 4 1 0.21648 65545.10313
4 5 1 0.26454 19022.17121
5 6 1 0.20452 53457.20659
[ angles ]
1 2 3 1 77.69557 383.16060
2 3 4 1 115.59599 262.66648
3 4 5 1 108.51038 698.56915
4 5 6 1 81.23827 1098.41706
5 6 1 1 146.49898 48.08599
6 1 2 1 99.84045 1571.73513
[ dihedrals ]
1 2 3 4 1 -85.45004 234.10405 1
2 3 4 5 1 65.21378 117.46132 1
3 4 5 6 1 -22.73796 94.20175 1
4 5 6 1 1 55.18087 37.92107 1
5 6 1 2 1 -96.52349 40.73340 1
6 1 2 3 1 73.82477 259.20126 1
[ constraints ]
1 2 1 0.22209
6 1 1 0.17931
import unittest import unittest
import filecmp
from pycgtool.bondset import BondSet from pycgtool.bondset import BondSet
from pycgtool.frame import Frame from pycgtool.frame import Frame
...@@ -103,3 +104,18 @@ class BondSetTest(unittest.TestCase): ...@@ -103,3 +104,18 @@ class BondSetTest(unittest.TestCase):
for bond in bondset.get_bond_lengths("ETH", True): for bond in bondset.get_bond_lengths("ETH", True):
self.assertAlmostEqual(1., bond.eqm) self.assertAlmostEqual(1., bond.eqm)
self.assertEqual(float("inf"), bond.fconst) self.assertEqual(float("inf"), bond.fconst)
def test_full_itp_sugar(self):
measure = BondSet("test/data/sugar.bnd", DummyOptions)
frame = Frame("test/data/sugar.gro", xtc="test/data/sugar.xtc")
mapping = Mapping("test/data/sugar.map", DummyOptions)
cgframe = mapping.apply(frame)
while frame.next_frame():
cgframe = mapping.apply(frame, cgframe=cgframe, exclude={"SOL"})
measure.apply(cgframe)
measure.boltzmann_invert()
measure.write_itp("sugar_out.itp", mapping, exclude={"SOL"})
self.assertTrue(filecmp.cmp("sugar_out.itp", "test/data/sugar_out.itp"))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment