mirror of
https://github.com/openmm/openmm
synced 2026-06-03 06:39:48 +09:00
fe0550bb7b
* Create bonds based on chem_comp_bond records * Fixed a test that assumed bonds would be in a particular order
215 lines
9.8 KiB
Python
215 lines
9.8 KiB
Python
import tempfile
|
|
import unittest
|
|
from openmm.app import *
|
|
from openmm import *
|
|
from openmm.unit import *
|
|
import openmm.app.element as elem
|
|
import os
|
|
from io import StringIO
|
|
|
|
class TestPdbxFile(unittest.TestCase):
|
|
"""Test the PDBx/mmCIF file parser"""
|
|
|
|
def test_FormatConversion(self):
|
|
"""Test conversion from PDB to PDBx"""
|
|
|
|
mol = PDBFile('systems/ala_ala_ala.pdb')
|
|
with tempfile.TemporaryDirectory() as tempdir:
|
|
filename = os.path.join(tempdir, 'temp.pdbx')
|
|
PDBxFile.writeFile(mol.topology, mol.positions, filename, keepIds=True)
|
|
try:
|
|
pdbx = PDBxFile(filename)
|
|
except Exception:
|
|
self.fail('Parser failed to read PDBx/mmCIF file')
|
|
|
|
|
|
def test_Triclinic(self):
|
|
"""Test parsing a file that describes a triclinic box."""
|
|
pdb = PDBxFile('systems/triclinic.pdbx')
|
|
self.assertEqual(len(pdb.positions), 8)
|
|
expectedPositions = [
|
|
Vec3(1.744, 2.788, 3.162),
|
|
Vec3(1.048, 0.762, 2.340),
|
|
Vec3(2.489, 1.570, 2.817),
|
|
Vec3(1.027, 1.893, 3.271),
|
|
Vec3(0.937, 0.825, 0.009),
|
|
Vec3(2.290, 1.887, 3.352),
|
|
Vec3(1.266, 1.111, 2.894),
|
|
Vec3(0.933, 1.862, 3.490)]*nanometers
|
|
for (p1, p2) in zip(expectedPositions, pdb.positions):
|
|
self.assertVecAlmostEqual(p1, p2)
|
|
expectedVectors = [
|
|
Vec3(2.5, 0, 0),
|
|
Vec3(0.5, 3.0, 0),
|
|
Vec3(0.7, 0.9, 3.5)]*nanometers
|
|
for (v1, v2) in zip(expectedVectors, pdb.topology.getPeriodicBoxVectors()):
|
|
self.assertVecAlmostEqual(v1, v2, 1e-6)
|
|
self.assertVecAlmostEqual(Vec3(2.5, 3.0, 3.5)*nanometers, pdb.topology.getUnitCellDimensions(), 1e-6)
|
|
for atom in pdb.topology.atoms():
|
|
if atom.index < 4:
|
|
self.assertEqual(elem.chlorine, atom.element)
|
|
self.assertEqual('Cl', atom.name)
|
|
self.assertEqual('Cl', atom.residue.name)
|
|
else:
|
|
self.assertEqual(elem.sodium, atom.element)
|
|
self.assertEqual('Na', atom.name)
|
|
self.assertEqual('Na', atom.residue.name)
|
|
|
|
def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
|
|
unit = p1.unit
|
|
p1 = p1.value_in_unit(unit)
|
|
p2 = p2.value_in_unit(unit)
|
|
scale = max(1.0, norm(p1),)
|
|
for i in range(3):
|
|
diff = abs(p1[i]-p2[i])/scale
|
|
self.assertTrue(diff < tol)
|
|
|
|
def testReporterImplicit(self):
|
|
""" Tests the PDBxReporter without PBC """
|
|
parm = AmberPrmtopFile('systems/alanine-dipeptide-implicit.prmtop')
|
|
system = parm.createSystem()
|
|
sim = Simulation(parm.topology, system, VerletIntegrator(1*femtoseconds),
|
|
Platform.getPlatform('Reference'))
|
|
sim.context.setPositions(PDBFile('systems/alanine-dipeptide-implicit.pdb').getPositions())
|
|
sim.reporters.append(PDBxReporter('test.cif', 1))
|
|
sim.step(10)
|
|
pdb = PDBxFile('test.cif')
|
|
self.assertEqual(len(list(pdb.topology.atoms())), len(list(parm.topology.atoms())))
|
|
self.assertEqual(len(list(pdb.topology.residues())), len(list(parm.topology.residues())))
|
|
for res1, res2 in zip(pdb.topology.residues(), parm.topology.residues()):
|
|
self.assertEqual(res1.name, res2.name)
|
|
for atom1, atom2 in zip(res1.atoms(), res2.atoms()):
|
|
self.assertEqual(atom1.name, atom2.name)
|
|
positions = pdb.getPositions(frame=9)
|
|
self.assertFalse(all(x1 == x2 for x1, x2 in zip(positions, pdb.getPositions(frame=0))))
|
|
# There should only be 10 frames (0 through 9)
|
|
self.assertRaises(IndexError, lambda: pdb.getPositions(frame=10))
|
|
self.assertIs(pdb.topology.getPeriodicBoxVectors(), None)
|
|
del sim
|
|
os.unlink('test.cif')
|
|
|
|
def assertAlmostEqualVec(self, vec1, vec2, *args, **kwargs):
|
|
if is_quantity(vec1):
|
|
vec1 = vec1.value_in_unit_system(md_unit_system)
|
|
if is_quantity(vec2):
|
|
vec2 = vec2.value_in_unit_system(md_unit_system)
|
|
for x, y in zip(vec1, vec2):
|
|
self.assertAlmostEqual(x, y, *args, **kwargs)
|
|
|
|
def testReporterExplicit(self):
|
|
""" Tests the PDBxReporter with PBC """
|
|
parm = AmberPrmtopFile('systems/alanine-dipeptide-explicit.prmtop')
|
|
system = parm.createSystem(nonbondedCutoff=1.0, nonbondedMethod=PME)
|
|
sim = Simulation(parm.topology, system, VerletIntegrator(1*femtoseconds),
|
|
Platform.getPlatform('Reference'))
|
|
orig_pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
|
|
sim.context.setPositions(orig_pdb.getPositions())
|
|
sim.context.setPeriodicBoxVectors(*parm.topology.getPeriodicBoxVectors())
|
|
sim.reporters.append(PDBxReporter('test.cif', 1))
|
|
sim.step(10)
|
|
pdb = PDBxFile('test.cif')
|
|
self.assertEqual(len(list(pdb.topology.atoms())), len(list(parm.topology.atoms())))
|
|
self.assertEqual(len(list(pdb.topology.residues())), len(list(parm.topology.residues())))
|
|
for res1, res2 in zip(pdb.topology.residues(), parm.topology.residues()):
|
|
self.assertEqual(res1.name, res2.name)
|
|
for atom1, atom2 in zip(res1.atoms(), res2.atoms()):
|
|
self.assertEqual(atom1.name, atom2.name)
|
|
positions = pdb.getPositions(frame=9)
|
|
self.assertFalse(all(x1 == x2 for x1, x2 in zip(positions, pdb.getPositions(frame=0))))
|
|
# There should only be 10 frames (0 through 9)
|
|
self.assertRaises(IndexError, lambda: pdb.getPositions(frame=10))
|
|
self.assertAlmostEqualVec(parm.topology.getPeriodicBoxVectors()[0],
|
|
pdb.topology.getPeriodicBoxVectors()[0],
|
|
places=5)
|
|
self.assertAlmostEqualVec(parm.topology.getPeriodicBoxVectors()[1],
|
|
pdb.topology.getPeriodicBoxVectors()[1],
|
|
places=5)
|
|
self.assertAlmostEqualVec(parm.topology.getPeriodicBoxVectors()[2],
|
|
pdb.topology.getPeriodicBoxVectors()[2],
|
|
places=5)
|
|
del sim
|
|
os.unlink('test.cif')
|
|
|
|
def testBonds(self):
|
|
"""Test reading and writing a file that includes bonds."""
|
|
pdb = PDBFile('systems/methanol_ions.pdb')
|
|
output = StringIO()
|
|
PDBxFile.writeFile(pdb.topology, pdb.positions, output)
|
|
input = StringIO(output.getvalue())
|
|
pdbx = PDBxFile(input)
|
|
output.close()
|
|
input.close()
|
|
self.assertEqual(pdb.topology.getNumBonds(), pdbx.topology.getNumBonds())
|
|
for bond1, bond2 in zip(pdb.topology.bonds(), pdbx.topology.bonds()):
|
|
self.assertEqual(bond1[0].name, bond2[0].name)
|
|
self.assertEqual(bond1[1].name, bond2[1].name)
|
|
|
|
def testChemCompBonds(self):
|
|
"""Test creating bonds based on chem_comp_bond records."""
|
|
pdb = PDBxFile('systems/6mvz.cif')
|
|
|
|
def bondCount(res1, atom1, res2, atom2):
|
|
count = 0
|
|
for bond in pdb.topology.bonds():
|
|
if (res1 == bond[0].residue.index and atom1 == bond[0].name) or (res2 == bond[0].residue.index and atom2 == bond[0].name):
|
|
if (res1 == bond[1].residue.index and atom1 == bond[1].name) or (res2 == bond[1].residue.index and atom2 == bond[1].name):
|
|
count += 1
|
|
return count
|
|
|
|
# Check some bonds that ought to be present.
|
|
|
|
self.assertEqual(1, bondCount(0, 'N', 0, 'CA'))
|
|
self.assertEqual(1, bondCount(0, 'C', 1, 'N'))
|
|
self.assertEqual(1, bondCount(1, 'CA', 1, 'CB'))
|
|
self.assertEqual(1, bondCount(1, 'C', 2, 'N'))
|
|
self.assertEqual(1, bondCount(2, 'CB', 2, 'CG'))
|
|
self.assertEqual(1, bondCount(2, 'C', 3, 'N'))
|
|
self.assertEqual(1, bondCount(3, 'C', 3, 'OXT'))
|
|
self.assertEqual(1, bondCount(4, 'O', 4, 'C1'))
|
|
|
|
# Check the types and orders of a few bonds.
|
|
|
|
for bond in pdb.topology.bonds():
|
|
if (bond[0].name == 'C' and bond[1].name == 'O') or (bond[1].name == 'C' and bond[0].name == 'O'):
|
|
self.assertEqual(topology.Double, bond.type)
|
|
self.assertEqual(2, bond.order)
|
|
if (bond[0].name == 'N' or bond[1].name == 'N'):
|
|
if bond[0].residue == bond[1].residue:
|
|
self.assertEqual(topology.Single, bond.type)
|
|
self.assertEqual(1, bond.order)
|
|
else:
|
|
self.assertEqual(None, bond.type)
|
|
self.assertEqual(None, bond.order)
|
|
|
|
def testMultiChain(self):
|
|
"""Test reading and writing a file that includes multiple chains"""
|
|
cif_ori = PDBxFile('systems/multichain.pdbx')
|
|
|
|
output = StringIO()
|
|
PDBxFile.writeFile(cif_ori.topology, cif_ori.positions, output, keepIds=True)
|
|
input = StringIO(output.getvalue())
|
|
cif_new = PDBxFile(input)
|
|
output.close()
|
|
input.close()
|
|
|
|
self.assertEqual(cif_ori.topology.getNumChains(), cif_new.topology.getNumChains())
|
|
|
|
for chain1, chain2 in zip(cif_ori.topology.chains(), cif_new.topology.chains()):
|
|
self.assertEqual(chain1.id, chain2.id)
|
|
|
|
def testInsertionCodes(self):
|
|
"""Test reading a file that uses insertion codes."""
|
|
pdbx = PDBxFile('systems/insertions.pdbx')
|
|
residues = list(pdbx.topology.residues())
|
|
self.assertEqual(7, len(residues))
|
|
names = ['PHE', 'ASP', 'LYS', 'ILE', 'LYS', 'ASN', 'TRP']
|
|
ids = ['59', '60', '60', '60', '60', '60', '61']
|
|
codes = ['', '', 'A', 'B', 'C', 'D', '']
|
|
for res, name, id, code in zip(residues, names, ids, codes):
|
|
self.assertEqual(name, res.name)
|
|
self.assertEqual(id, res.id)
|
|
self.assertEqual(code, res.insertionCode)
|
|
|
|
if __name__ == '__main__':
|
|
unittest.main()
|