From d4fabd07ebf81641f5afb3cb8c6acc377a5c8192 Mon Sep 17 00:00:00 2001 From: James Graham <J.A.Graham@soton.ac.uk> Date: Thu, 11 Aug 2016 12:24:43 +0100 Subject: [PATCH] Add test for itp --- pycgtool/bondset.py | 13 ++++++++----- test/data/sugar_out.itp | 43 +++++++++++++++++++++++++++++++++++++++++ test/test_bondset.py | 16 +++++++++++++++ 3 files changed, 67 insertions(+), 5 deletions(-) create mode 100644 test/data/sugar_out.itp diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 129ab5a..11f0dc7 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -226,26 +226,29 @@ class BondSet: e.args = ("Bead(s) {0} do(es) not exist in residue {1}".format(missing, mol),) 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. :param filename: Name of output file :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) 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: print("\n[ {0:s} ]".format(section_header), file=itp) for bond in bonds: # Factor is usually 1, unless doing correction eqm = bond.eqm * self._empirical_correction_factor 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: line += " {0:12.5f}".format(bond.fconst) + if multiplicity is not None: + line += " {0:4d}".format(multiplicity) print(line, file=itp) with open(filename, "w") as itp: @@ -258,7 +261,7 @@ class BondSet: print(header, file=itp) # 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("{0:4s} {1:4d}".format(mol, 1), file=itp) @@ -271,7 +274,7 @@ class BondSet: 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_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) def apply(self, frame): diff --git a/test/data/sugar_out.itp b/test/data/sugar_out.itp new file mode 100644 index 0000000..1fc051b --- /dev/null +++ b/test/data/sugar_out.itp @@ -0,0 +1,43 @@ +; +; 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 diff --git a/test/test_bondset.py b/test/test_bondset.py index faf529e..5028dc8 100644 --- a/test/test_bondset.py +++ b/test/test_bondset.py @@ -1,4 +1,5 @@ import unittest +import filecmp from pycgtool.bondset import BondSet from pycgtool.frame import Frame @@ -103,3 +104,18 @@ class BondSetTest(unittest.TestCase): for bond in bondset.get_bond_lengths("ETH", True): self.assertAlmostEqual(1., bond.eqm) 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")) -- GitLab