Skip to content

GitLab

Commit 4b2ef27c authored by jchodera's avatar jchodera
Browse files

Fixed tests

parent 46f62eb2
Hide whitespace changes
Inline Side-by-side
Showing with 82 additions and 84 deletions
wrappers/python/tests/TestModeller.py
View file @ 4b2ef27c
...@@ -15,14 +15,14 @@ class TestModeller(unittest.TestCase): ...@@ -15,14 +15,14 @@ class TestModeller(unittest.TestCase):
# load the alanine dipeptide pdb file # load the alanine dipeptide pdb file
self.pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb') self.pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
self.topology_start = self.pdb.topology self.topology_start = self.pdb.topology
self.positions = self.pdb.positions self.positions = self.pdb.positions
self.forcefield = ForceField('amber10.xml', 'tip3p.xml') self.forcefield = ForceField('amber10.xml', 'tip3p.xml')
# load the T4-lysozyme-L99A receptor pdb file # load the T4-lysozyme-L99A receptor pdb file
self.pdb2 = PDBFile('systems/lysozyme-implicit.pdb') self.pdb2 = PDBFile('systems/lysozyme-implicit.pdb')
self.topology_start2 = self.pdb2.topology self.topology_start2 = self.pdb2.topology
self.positions2 = self.pdb2.positions self.positions2 = self.pdb2.positions
# load the metallothionein pdb file # load the metallothionein pdb file
self.pdb3 = PDBFile('systems/1T2Y.pdb') self.pdb3 = PDBFile('systems/1T2Y.pdb')
self.topology_start3 = self.pdb3.topology self.topology_start3 = self.pdb3.topology
...@@ -30,12 +30,12 @@ class TestModeller(unittest.TestCase): ...@@ -30,12 +30,12 @@ class TestModeller(unittest.TestCase):
def test_deleteWater(self): def test_deleteWater(self):
""" Test the deleteWater() method. """ """ Test the deleteWater() method. """
# build the chain dictionary # build the chain dictionary
chain_dict = {0:0} chain_dict = {0:0}
# 749 water chains are deleted # 749 water chains are deleted
chain_delta = -749 chain_delta = -749
# Build the residue and atom dictionaries for validate_preserved. # Build the residue and atom dictionaries for validate_preserved.
# Also, count the number of deleted residues and atoms. # Also, count the number of deleted residues and atoms.
residues_preserved = 0 residues_preserved = 0
...@@ -55,38 +55,38 @@ class TestModeller(unittest.TestCase): ...@@ -55,38 +55,38 @@ class TestModeller(unittest.TestCase):
residue_delta -= 1 residue_delta -= 1
for atom in residue.atoms(): for atom in residue.atoms():
atom_delta -= 1 atom_delta -= 1
modeller = Modeller(self.topology_start, self.positions) modeller = Modeller(self.topology_start, self.positions)
modeller.deleteWater() modeller.deleteWater()
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
validate_preserved(self, self.topology_start, topology_after, validate_preserved(self, self.topology_start, topology_after,
chain_dict, residue_dict, atom_dict) chain_dict, residue_dict, atom_dict)
validate_deltas(self, self.topology_start, topology_after, validate_deltas(self, self.topology_start, topology_after,
chain_delta, residue_delta, atom_delta) chain_delta, residue_delta, atom_delta)
def test_delete(self): def test_delete(self):
""" Test the delete() method. """ """ Test the delete() method. """
modeller = Modeller(self.topology_start, self.positions) modeller = Modeller(self.topology_start, self.positions)
topology_before = modeller.getTopology() topology_before = modeller.getTopology()
# Create the list of items to be deleted. # Create the list of items to be deleted.
# Start with the first 50 water chains # Start with the first 50 water chains
chains = [chain for chain in topology_before.chains()] chains = [chain for chain in topology_before.chains()]
toDelete = chains[1:51] toDelete = chains[1:51]
# Next add water residues 103->152 to the list of items to be deleted # Next add water residues 103->152 to the list of items to be deleted
residues = [residue for residue in topology_before.residues()] residues = [residue for residue in topology_before.residues()]
toDelete.extend(residues[103:153]) toDelete.extend(residues[103:153])
# Finally add water atoms 622->771 to the list of items to be deleted # Finally add water atoms 622->771 to the list of items to be deleted
atoms = [atom for atom in topology_before.atoms()] atoms = [atom for atom in topology_before.atoms()]
toDelete.extend(atoms[622:772]) toDelete.extend(atoms[622:772])
modeller.delete(toDelete) modeller.delete(toDelete)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
# build the chain dictionary # build the chain dictionary
chain_dict = {0:0} chain_dict = {0:0}
for i in range(1,51): for i in range(1,51):
...@@ -95,7 +95,7 @@ class TestModeller(unittest.TestCase): ...@@ -95,7 +95,7 @@ class TestModeller(unittest.TestCase):
chain_dict[i+100] = i chain_dict[i+100] = i
for i in range(101, 600): for i in range(101, 600):
chain_dict[i+150] = i chain_dict[i+150] = i
# build the residue dictionary # build the residue dictionary
residue_dict = {} residue_dict = {}
for i in range(3): for i in range(3):
...@@ -106,7 +106,7 @@ class TestModeller(unittest.TestCase): ...@@ -106,7 +106,7 @@ class TestModeller(unittest.TestCase):
residue_dict[i+100] = i residue_dict[i+100] = i
for i in range(103, 602): for i in range(103, 602):
residue_dict[i+150] = i residue_dict[i+150] = i
# build the atom dictionary # build the atom dictionary
atom_dict = {} atom_dict = {}
for i in range(22): for i in range(22):
...@@ -117,37 +117,37 @@ class TestModeller(unittest.TestCase): ...@@ -117,37 +117,37 @@ class TestModeller(unittest.TestCase):
atom_dict[i+300] = i atom_dict[i+300] = i
for i in range(322,1819): for i in range(322,1819):
atom_dict[i+450] = i atom_dict[i+450] = i
validate_preserved(self, topology_before, topology_after, chain_dict, residue_dict, atom_dict) validate_preserved(self, topology_before, topology_after, chain_dict, residue_dict, atom_dict)
chain_delta = -150 chain_delta = -150
residue_delta = -150 residue_delta = -150
atom_delta = -450 atom_delta = -450
validate_deltas(self, topology_before, topology_after, chain_delta, residue_delta, atom_delta) validate_deltas(self, topology_before, topology_after, chain_delta, residue_delta, atom_delta)
def test_add(self): def test_add(self):
""" Test the add() method. """ """ Test the add() method. """
# load the methanol-box pdb file # load the methanol-box pdb file
pdb2 = PDBFile('systems/methanol-box.pdb') pdb2 = PDBFile('systems/methanol-box.pdb')
topology_toAdd = pdb2.topology topology_toAdd = pdb2.topology
positions_toAdd = pdb2.positions positions_toAdd = pdb2.positions
modeller = Modeller(self.topology_start, self.positions) modeller = Modeller(self.topology_start, self.positions)
modeller.deleteWater() modeller.deleteWater()
topology_before = modeller.getTopology() topology_before = modeller.getTopology()
modeller.add(topology_toAdd, positions_toAdd) modeller.add(topology_toAdd, positions_toAdd)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
# build the first chain dictionary for the first call of validate_preserved() # build the first chain dictionary for the first call of validate_preserved()
chain_counter = 0 chain_counter = 0
chain_dict = {} chain_dict = {}
for chain in topology_before.chains(): for chain in topology_before.chains():
chain_dict[chain.index] = chain_counter chain_dict[chain.index] = chain_counter
chain_counter += 1 chain_counter += 1
# build the residue and atom dictionaries for the first call of validate_preserved() # build the residue and atom dictionaries for the first call of validate_preserved()
residue_counter = 0 residue_counter = 0
residue_dict = {} residue_dict = {}
...@@ -159,13 +159,13 @@ class TestModeller(unittest.TestCase): ...@@ -159,13 +159,13 @@ class TestModeller(unittest.TestCase):
for atom in residue.atoms(): for atom in residue.atoms():
atom_dict[atom.index] = atom_counter atom_dict[atom.index] = atom_counter
atom_counter += 1 atom_counter += 1
# Validate that the items from the before topology are preserved after addition of items. # Validate that the items from the before topology are preserved after addition of items.
validate_preserved(self, topology_before, topology_after, chain_dict, residue_dict, atom_dict) validate_preserved(self, topology_before, topology_after, chain_dict, residue_dict, atom_dict)
# Next, we build another set of dictionaries to validate that the items added are # Next, we build another set of dictionaries to validate that the items added are
# preserved. Also, we calculate the number of chains, residues, and atoms added. # preserved. Also, we calculate the number of chains, residues, and atoms added.
# build the chain dictionary # build the chain dictionary
chain_delta = 0 chain_delta = 0
chain_dict = {} chain_dict = {}
...@@ -173,7 +173,7 @@ class TestModeller(unittest.TestCase): ...@@ -173,7 +173,7 @@ class TestModeller(unittest.TestCase):
chain_dict[chain.index] = chain_counter chain_dict[chain.index] = chain_counter
chain_counter += 1 chain_counter += 1
chain_delta += 1 chain_delta += 1
# build the residue and atom dictionaries for the second call of validate_preserved # build the residue and atom dictionaries for the second call of validate_preserved
residue_delta = 0 residue_delta = 0
residue_dict = {} residue_dict = {}
...@@ -187,26 +187,26 @@ class TestModeller(unittest.TestCase): ...@@ -187,26 +187,26 @@ class TestModeller(unittest.TestCase):
atom_dict[atom.index] = atom_counter atom_dict[atom.index] = atom_counter
atom_counter += 1 atom_counter += 1
atom_delta += 1 atom_delta += 1
# validate that the items in the added topology are preserved # validate that the items in the added topology are preserved
validate_preserved(self, topology_toAdd, topology_after, chain_dict, residue_dict, atom_dict) validate_preserved(self, topology_toAdd, topology_after, chain_dict, residue_dict, atom_dict)
# validate that the final topology has the correct number of items # validate that the final topology has the correct number of items
validate_deltas(self, topology_before, topology_after, chain_delta, residue_delta, atom_delta) validate_deltas(self, topology_before, topology_after, chain_delta, residue_delta, atom_delta)
def test_convertWater(self): def test_convertWater(self):
""" Test the convertWater() method. """ """ Test the convertWater() method. """
for model in ['tip3p', 'spce', 'tip4pew', 'tip5p']: for model in ['tip3p', 'spce', 'tip4pew', 'tip5p']:
if model == 'tip5p': if model == 'tip5p':
firstmodel = 'tip4pew' firstmodel = 'tip4pew'
else: else:
firstmodel = 'tip5p' firstmodel = 'tip5p'
modeller = Modeller(self.topology_start, self.positions) modeller = Modeller(self.topology_start, self.positions)
modeller.convertWater(model=firstmodel) modeller.convertWater(model=firstmodel)
modeller.convertWater(model=model) modeller.convertWater(model=model)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
for residue in topology_after.residues(): for residue in topology_after.residues():
if residue.name == "HOH": if residue.name == "HOH":
oatom = [atom for atom in residue.atoms() if atom.element == element.oxygen] oatom = [atom for atom in residue.atoms() if atom.element == element.oxygen]
...@@ -221,7 +221,7 @@ class TestModeller(unittest.TestCase): ...@@ -221,7 +221,7 @@ class TestModeller(unittest.TestCase):
self.assertTrue(len(matoms)==1 and len(m1atoms)==0 and len(m2atoms)==0) self.assertTrue(len(matoms)==1 and len(m1atoms)==0 and len(m2atoms)==0)
elif model=='tip5p': elif model=='tip5p':
self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1) self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1)
# build the chain dictionary for validate_preserved # build the chain dictionary for validate_preserved
chain_counter = 0 chain_counter = 0
chain_dict = {} chain_dict = {}
...@@ -229,7 +229,7 @@ class TestModeller(unittest.TestCase): ...@@ -229,7 +229,7 @@ class TestModeller(unittest.TestCase):
for chain in self.topology_start.chains(): for chain in self.topology_start.chains():
chain_dict[chain.index] = chain_counter chain_dict[chain.index] = chain_counter
chain_counter += 1 chain_counter += 1
# build the residue and atom dictionaries for validate_preserved # build the residue and atom dictionaries for validate_preserved
residue_counter = 0 residue_counter = 0
residue_dict = {} residue_dict = {}
...@@ -249,16 +249,16 @@ class TestModeller(unittest.TestCase): ...@@ -249,16 +249,16 @@ class TestModeller(unittest.TestCase):
if residue.name == 'HOH' and model == 'tip5p': if residue.name == 'HOH' and model == 'tip5p':
atom_counter += 2 atom_counter += 2
atom_delta += 2 atom_delta += 2
validate_preserved(self, self.topology_start, topology_after, validate_preserved(self, self.topology_start, topology_after,
chain_dict, residue_dict, atom_dict) chain_dict, residue_dict, atom_dict)
validate_deltas(self, self.topology_start, topology_after, validate_deltas(self, self.topology_start, topology_after,
chain_delta, residue_delta, atom_delta) chain_delta, residue_delta, atom_delta)
def test_addSolventWaterModels(self): def test_addSolventWaterModels(self):
""" Test all addSolvent() method with all possible water models. """ """ Test all addSolvent() method with all possible water models. """
topology_start = self.pdb.topology topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers) topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers)
for model in ['tip3p', 'spce', 'tip4pew', 'tip5p']: for model in ['tip3p', 'spce', 'tip4pew', 'tip5p']:
...@@ -270,16 +270,16 @@ class TestModeller(unittest.TestCase): ...@@ -270,16 +270,16 @@ class TestModeller(unittest.TestCase):
# add the solvent to get the "after" topology # add the solvent to get the "after" topology
modeller.addSolvent(forcefield, model=model) modeller.addSolvent(forcefield, model=model)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
# First, check that everything that was there before has been preserved. # First, check that everything that was there before has been preserved.
# build the chain dictionary for validate_preserved # build the chain dictionary for validate_preserved
chain_counter = 0 chain_counter = 0
chain_dict = {0:0} chain_dict = {0:0}
for chain in topology_before.chains(): for chain in topology_before.chains():
chain_dict[chain.index] = chain_counter chain_dict[chain.index] = chain_counter
chain_counter += 1 chain_counter += 1
# build the residue and atom dictionaries for validate_preserved # build the residue and atom dictionaries for validate_preserved
residue_counter = 0 residue_counter = 0
residue_dict = {} residue_dict = {}
...@@ -291,10 +291,10 @@ class TestModeller(unittest.TestCase): ...@@ -291,10 +291,10 @@ class TestModeller(unittest.TestCase):
for atom in residue.atoms(): for atom in residue.atoms():
atom_dict[atom.index] = atom_counter atom_dict[atom.index] = atom_counter
atom_counter += 1 atom_counter += 1
# validate that the items in the before topology remain after solvent is added # validate that the items in the before topology remain after solvent is added
validate_preserved(self, topology_before, topology_after, chain_dict, residue_dict, atom_dict) validate_preserved(self, topology_before, topology_after, chain_dict, residue_dict, atom_dict)
# Make sure water that was added was the correct model # Make sure water that was added was the correct model
for residue in topology_after.residues(): for residue in topology_after.residues():
if residue.name == 'HOH': if residue.name == 'HOH':
...@@ -310,10 +310,10 @@ class TestModeller(unittest.TestCase): ...@@ -310,10 +310,10 @@ class TestModeller(unittest.TestCase):
self.assertTrue(len(matoms)==1 and len(m1atoms)==0 and len(m2atoms)==0) self.assertTrue(len(matoms)==1 and len(m1atoms)==0 and len(m2atoms)==0)
elif model=='tip5p': elif model=='tip5p':
self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1) self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1)
def test_addSolventPeriodicBox(self): def test_addSolventPeriodicBox(self):
""" Test the addSolvent() method; test that the five ways of passing in the periodic box all work. """ """ Test the addSolvent() method; test that the five ways of passing in the periodic box all work. """
# First way of passing in periodic box vectors: set it in the original topology. # First way of passing in periodic box vectors: set it in the original topology.
topology_start = self.pdb.topology topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 4.5, 5.5)*nanometers) topology_start.setUnitCellDimensions(Vec3(3.5, 4.5, 5.5)*nanometers)
...@@ -322,11 +322,11 @@ class TestModeller(unittest.TestCase): ...@@ -322,11 +322,11 @@ class TestModeller(unittest.TestCase):
modeller.addSolvent(self.forcefield) modeller.addSolvent(self.forcefield)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors() dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.5, 0, 0)) self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.5, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.5, 0)) self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.5, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 5.5)) self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 5.5))
# Second way of passing in the periodic box vectors: with the boxSize parameter to addSolvent() # Second way of passing in the periodic box vectors: with the boxSize parameter to addSolvent()
topology_start = self.pdb.topology topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions) modeller = Modeller(topology_start, self.positions)
...@@ -334,11 +334,11 @@ class TestModeller(unittest.TestCase): ...@@ -334,11 +334,11 @@ class TestModeller(unittest.TestCase):
modeller.addSolvent(self.forcefield, boxSize = Vec3(3.6, 4.6, 5.6)*nanometers) modeller.addSolvent(self.forcefield, boxSize = Vec3(3.6, 4.6, 5.6)*nanometers)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors() dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.6, 0, 0)) self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.6, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.6, 0)) self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.6, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 5.6)) self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 5.6))
# Third way of passing in the periodic box vectors: with the boxVectors parameter to addSolvent() # Third way of passing in the periodic box vectors: with the boxVectors parameter to addSolvent()
topology_start = self.pdb.topology topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions) modeller = Modeller(topology_start, self.positions)
...@@ -346,11 +346,11 @@ class TestModeller(unittest.TestCase): ...@@ -346,11 +346,11 @@ class TestModeller(unittest.TestCase):
modeller.addSolvent(self.forcefield, boxVectors = (Vec3(3.4, 0, 0), Vec3(0.5, 4.4, 0), Vec3(-1.0, -1.5, 5.4))*nanometers) modeller.addSolvent(self.forcefield, boxVectors = (Vec3(3.4, 0, 0), Vec3(0.5, 4.4, 0), Vec3(-1.0, -1.5, 5.4))*nanometers)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors() dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.4, 0, 0)) self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.4, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0.5, 4.4, 0)) self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0.5, 4.4, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(-1.0, -1.5, 5.4)) self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(-1.0, -1.5, 5.4))
# Fourth way of passing in the periodic box vectors: pass a 'padding' value to addSolvent() # Fourth way of passing in the periodic box vectors: pass a 'padding' value to addSolvent()
topology_start = self.pdb.topology topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions) modeller = Modeller(topology_start, self.positions)
...@@ -358,11 +358,11 @@ class TestModeller(unittest.TestCase): ...@@ -358,11 +358,11 @@ class TestModeller(unittest.TestCase):
modeller.addSolvent(self.forcefield, padding = 1.0*nanometers) modeller.addSolvent(self.forcefield, padding = 1.0*nanometers)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors() dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(2.8802, 0, 0)) self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(2.8802, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 2.8802, 0)) self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 2.8802, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 2.8802)) self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 2.8802))
# Fifth way: specify a number of molecules to add instead of a box size # Fifth way: specify a number of molecules to add instead of a box size
topology_start = self.pdb.topology topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions) modeller = Modeller(topology_start, self.positions)
...@@ -395,7 +395,7 @@ class TestModeller(unittest.TestCase): ...@@ -395,7 +395,7 @@ class TestModeller(unittest.TestCase):
total_added = water_count+sodium_count+chlorine_count total_added = water_count+sodium_count+chlorine_count
self.assertEqual(total_added, 1364) self.assertEqual(total_added, 1364)
expected_ion_fraction = 2.0*molar/(55.4*molar) expected_ion_fraction = 2.0*molar/(55.4*molar)
expected_ions = math.floor(total_added*expected_ion_fraction/2+0.5) expected_ions = math.floor(total_added*expected_ion_fraction+0.5)
self.assertEqual(sodium_count, expected_ions) self.assertEqual(sodium_count, expected_ions)
self.assertEqual(chlorine_count, expected_ions) self.assertEqual(chlorine_count, expected_ions)
...@@ -417,7 +417,7 @@ class TestModeller(unittest.TestCase): ...@@ -417,7 +417,7 @@ class TestModeller(unittest.TestCase):
topology_toAdd.addResidue('CL', newChain) topology_toAdd.addResidue('CL', newChain)
residues = [residue for residue in topology_toAdd.residues()] residues = [residue for residue in topology_toAdd.residues()]
for i in range(5): for i in range(5):
topology_toAdd.addAtom('Cl',Element.getBySymbol('Cl'), residues[i]) topology_toAdd.addAtom('Cl',Element.getBySymbol('Cl'), residues[i])
positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0), positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
modeller.add(topology_toAdd, positions_toAdd) modeller.add(topology_toAdd, positions_toAdd)
...@@ -437,7 +437,7 @@ class TestModeller(unittest.TestCase): ...@@ -437,7 +437,7 @@ class TestModeller(unittest.TestCase):
total_water_ions = water_count+sodium_count+chlorine_count total_water_ions = water_count+sodium_count+chlorine_count
expected_ion_fraction = 1.0*molar/(55.4*molar) expected_ion_fraction = 1.0*molar/(55.4*molar)
expected_chlorine = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5 expected_chlorine = math.floor((total_water_ions-10)*expected_ion_fraction+0.5)+5
expected_sodium = expected_chlorine if neutralize else expected_chlorine-5 expected_sodium = expected_chlorine if neutralize else expected_chlorine-5
self.assertEqual(sodium_count, expected_sodium) self.assertEqual(sodium_count, expected_sodium)
self.assertEqual(chlorine_count, expected_chlorine) self.assertEqual(chlorine_count, expected_chlorine)
...@@ -460,7 +460,7 @@ class TestModeller(unittest.TestCase): ...@@ -460,7 +460,7 @@ class TestModeller(unittest.TestCase):
topology_toAdd.addResidue('NA', newChain) topology_toAdd.addResidue('NA', newChain)
residues = [residue for residue in topology_toAdd.residues()] residues = [residue for residue in topology_toAdd.residues()]
for i in range(5): for i in range(5):
topology_toAdd.addAtom('Na',Element.getBySymbol('Na'), residues[i]) topology_toAdd.addAtom('Na',Element.getBySymbol('Na'), residues[i])
positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0), positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
...@@ -482,14 +482,14 @@ class TestModeller(unittest.TestCase): ...@@ -482,14 +482,14 @@ class TestModeller(unittest.TestCase):
total_water_ions = water_count+sodium_count+chlorine_count total_water_ions = water_count+sodium_count+chlorine_count
expected_ion_fraction = 1.0*molar/(55.4*molar) expected_ion_fraction = 1.0*molar/(55.4*molar)
expected_sodium = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5 expected_sodium = math.floor((total_water_ions-10)*expected_ion_fraction+0.5)+5
expected_chlorine = expected_sodium if neutralize else expected_sodium-5 expected_chlorine = expected_sodium if neutralize else expected_sodium-5
self.assertEqual(sodium_count, expected_sodium) self.assertEqual(sodium_count, expected_sodium)
self.assertEqual(chlorine_count, expected_chlorine) self.assertEqual(chlorine_count, expected_chlorine)
def test_addSolventIons(self): def test_addSolventIons(self):
""" Test the addSolvent() method with all possible choices for positive and negative ions. """ """ Test the addSolvent() method with all possible choices for positive and negative ions. """
topology_start = self.pdb.topology topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers) topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers)
...@@ -500,19 +500,19 @@ class TestModeller(unittest.TestCase): ...@@ -500,19 +500,19 @@ class TestModeller(unittest.TestCase):
positions_nowater = modeller.getPositions() positions_nowater = modeller.getPositions()
expected_ion_fraction = 1.0*molar/(55.4*molar) expected_ion_fraction = 1.0*molar/(55.4*molar)
for positiveIon in ['Cs+', 'K+', 'Li+', 'Na+', 'Rb+']: for positiveIon in ['Cs+', 'K+', 'Li+', 'Na+', 'Rb+']:
ionName = positiveIon[:-1].upper() ionName = positiveIon[:-1].upper()
modeller = Modeller(topology_nowater, positions_nowater) modeller = Modeller(topology_nowater, positions_nowater)
modeller.addSolvent(self.forcefield, positiveIon=positiveIon, ionicStrength=1.0*molar) modeller.addSolvent(self.forcefield, positiveIon=positiveIon, ionicStrength=1.0*molar)
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
water_count = 0 water_count = 0
positive_ion_count = 0 positive_ion_count = 0
chlorine_count = 0 chlorine_count = 0
for residue in topology_after.residues(): for residue in topology_after.residues():
if residue.name=='HOH': if residue.name=='HOH':
water_count += 1 water_count += 1
elif residue.name==ionName: elif residue.name==ionName:
positive_ion_count += 1 positive_ion_count += 1
elif residue.name=='CL': elif residue.name=='CL':
...@@ -520,7 +520,7 @@ class TestModeller(unittest.TestCase): ...@@ -520,7 +520,7 @@ class TestModeller(unittest.TestCase):
total_added = water_count+positive_ion_count+chlorine_count total_added = water_count+positive_ion_count+chlorine_count
self.assertEqual(total_added, 1364) self.assertEqual(total_added, 1364)
expected_ions = math.floor(total_added*expected_ion_fraction/2+0.5) expected_ions = math.floor(total_added*expected_ion_fraction+0.5)
self.assertEqual(positive_ion_count, expected_ions) self.assertEqual(positive_ion_count, expected_ions)
self.assertEqual(chlorine_count, expected_ions) self.assertEqual(chlorine_count, expected_ions)
...@@ -532,7 +532,7 @@ class TestModeller(unittest.TestCase): ...@@ -532,7 +532,7 @@ class TestModeller(unittest.TestCase):
topology_after = modeller.getTopology() topology_after = modeller.getTopology()
water_count = 0 water_count = 0
sodium_count = 0 sodium_count = 0
negative_ion_count = 0 negative_ion_count = 0
for residue in topology_after.residues(): for residue in topology_after.residues():
if residue.name=='HOH': if residue.name=='HOH':
...@@ -541,16 +541,16 @@ class TestModeller(unittest.TestCase): ...@@ -541,16 +541,16 @@ class TestModeller(unittest.TestCase):
sodium_count += 1 sodium_count += 1
elif residue.name==ionName: elif residue.name==ionName:
negative_ion_count += 1 negative_ion_count += 1
total_added = water_count+sodium_count+negative_ion_count total_added = water_count+sodium_count+negative_ion_count
self.assertEqual(total_added, 1364) self.assertEqual(total_added, 1364)
expected_ions = math.floor(total_added*expected_ion_fraction/2+0.5) expected_ions = math.floor(total_added*expected_ion_fraction+0.5)
self.assertEqual(positive_ion_count, expected_ions) self.assertEqual(positive_ion_count, expected_ions)
self.assertEqual(chlorine_count, expected_ions) self.assertEqual(chlorine_count, expected_ions)
def test_addHydrogensPdb2(self): def test_addHydrogensPdb2(self):
""" Test the addHydrogens() method on the T4-lysozyme-L99A pdb file. """ """ Test the addHydrogens() method on the T4-lysozyme-L99A pdb file. """
# build the Modeller # build the Modeller
topology_start = self.topology_start2 topology_start = self.topology_start2
positions = self.positions2 positions = self.positions2
...@@ -576,7 +576,7 @@ class TestModeller(unittest.TestCase): ...@@ -576,7 +576,7 @@ class TestModeller(unittest.TestCase):
def test_addHydrogensPdb3(self): def test_addHydrogensPdb3(self):
""" Test the addHydrogens() method on the metallothionein pdb file. """ """ Test the addHydrogens() method on the metallothionein pdb file. """
# build the Modeller # build the Modeller
topology_start = self.topology_start3 topology_start = self.topology_start3
positions = self.positions3 positions = self.positions3
...@@ -671,7 +671,7 @@ class TestModeller(unittest.TestCase): ...@@ -671,7 +671,7 @@ class TestModeller(unittest.TestCase):
index_list_CYS = [31, 49, 110, 135, 171, 193, 229] index_list_CYS = [31, 49, 110, 135, 171, 193, 229]
atoms = [atom for atom in topology2.atoms()] atoms = [atom for atom in topology2.atoms()]
toDelete2 = [] toDelete2 = []
for index in index_list_CYS: for index in index_list_CYS:
self.assertTrue(atoms[index].element.symbol=='H') self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index]) toDelete2.append(atoms[index])
modeller2.delete(toDelete2) modeller2.delete(toDelete2)
...@@ -690,7 +690,7 @@ class TestModeller(unittest.TestCase): ...@@ -690,7 +690,7 @@ class TestModeller(unittest.TestCase):
modeller = Modeller(topology_start, positions) modeller = Modeller(topology_start, positions)
# remove hydrogens from the topology # remove hydrogens from the topology
toDelete = [atom for atom in topology_start.atoms() if atom.element==Element.getBySymbol('H')] toDelete = [atom for atom in topology_start.atoms() if atom.element==Element.getBySymbol('H')]
modeller.delete(toDelete) modeller.delete(toDelete)
# Create a variants list to force the one histidine to be of the right variation. # Create a variants list to force the one histidine to be of the right variation.
...@@ -709,13 +709,13 @@ class TestModeller(unittest.TestCase): ...@@ -709,13 +709,13 @@ class TestModeller(unittest.TestCase):
modeller.addHydrogens(self.forcefield, variants=variants) modeller.addHydrogens(self.forcefield, variants=variants)
topology_GLH = modeller.getTopology() topology_GLH = modeller.getTopology()
# There should be extra hydrogens on the GLU residues. Assert that they exist, # There should be extra hydrogens on the GLU residues. Assert that they exist,
# then we delete them and validate that the topology matches what we started with. # then we delete them and validate that the topology matches what we started with.
index_list_GLH = [85, 192, 387, 731, 992, 1018, 1718, 2042] index_list_GLH = [85, 192, 387, 731, 992, 1018, 1718, 2042]
atoms = [atom for atom in topology_GLH.atoms()] atoms = [atom for atom in topology_GLH.atoms()]
toDelete2 = [] toDelete2 = []
for index in index_list_GLH: for index in index_list_GLH:
self.assertTrue(atoms[index].element.symbol=='H') self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index]) toDelete2.append(atoms[index])
modeller.delete(toDelete2) modeller.delete(toDelete2)
topology_GLU = modeller.getTopology() topology_GLU = modeller.getTopology()
...@@ -811,7 +811,7 @@ class TestModeller(unittest.TestCase): ...@@ -811,7 +811,7 @@ class TestModeller(unittest.TestCase):
self.assertTrue(atoms[index].element.symbol=='H') self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index]) toDelete2.append(atoms[index])
for index in index_list_GLH: for index in index_list_GLH:
self.assertTrue(atoms[index].element.symbol=='H') self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index]) toDelete2.append(atoms[index])
modeller.delete(toDelete2) modeller.delete(toDelete2)
topology_ASP_GLU = modeller.getTopology() topology_ASP_GLU = modeller.getTopology()
...@@ -847,13 +847,13 @@ class TestModeller(unittest.TestCase): ...@@ -847,13 +847,13 @@ class TestModeller(unittest.TestCase):
index_list_CYS = [31, 49, 110, 135, 171, 193, 229] index_list_CYS = [31, 49, 110, 135, 171, 193, 229]
atoms = [atom for atom in topology2.atoms()] atoms = [atom for atom in topology2.atoms()]
toDelete2 = [] toDelete2 = []
for index in index_list_CYS: for index in index_list_CYS:
self.assertTrue(atoms[index].element.symbol=='H') self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index]) toDelete2.append(atoms[index])
modeller2.delete(toDelete2) modeller2.delete(toDelete2)
topology_after = modeller2.getTopology() topology_after = modeller2.getTopology()
validate_equivalence(self, topology_CYX, topology_after) validate_equivalence(self, topology_CYX, topology_after)
def test_addHydrogenspH11(self): def test_addHydrogenspH11(self):
""" Test of addHydrogens() with pH=11. """ """ Test of addHydrogens() with pH=11. """
...@@ -870,7 +870,7 @@ class TestModeller(unittest.TestCase): ...@@ -870,7 +870,7 @@ class TestModeller(unittest.TestCase):
modeller.delete(toDelete) modeller.delete(toDelete)
# Create a variants list to force the one histidine to be of the right variation. # Create a variants list to force the one histidine to be of the right variation.
residues = [residue for residue in topology_start.residues()] residues = [residue for residue in topology_start.residues()]
variants = [None]*len(residues) variants = [None]*len(residues)
# For this protein, when you add hydrogens, the hydrogen is added to the delta nitrogen. # For this protein, when you add hydrogens, the hydrogen is added to the delta nitrogen.
# By setting variants[30] to 'HIE', we force the hydrogen onto the epsilon nitrogen, so # By setting variants[30] to 'HIE', we force the hydrogen onto the epsilon nitrogen, so
...@@ -982,7 +982,7 @@ class TestModeller(unittest.TestCase): ...@@ -982,7 +982,7 @@ class TestModeller(unittest.TestCase):
topology.addAtom('C', element.carbon, residue) topology.addAtom('C', element.carbon, residue)
topology.addAtom('N', element.nitrogen, residue) topology.addAtom('N', element.nitrogen, residue)
topology.addAtom('V', element.oxygen, residue) topology.addAtom('V', element.oxygen, residue)
# Add the virtual sites. # Add the virtual sites.
...@@ -1083,5 +1083,3 @@ class TestModeller(unittest.TestCase): ...@@ -1083,5 +1083,3 @@ class TestModeller(unittest.TestCase):
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment