Skip to content

GitLab

Commit 5120bb0d authored by peastman's avatar peastman
Browse files

Merge pull request #1349 from jchodera/fix-ionic-strength 

Fix error in ionicStrength definition in Modeller.addSolvent()
parents 6b3a8337 4b2ef27c
Hide whitespace changes
Inline Side-by-side
Showing with 87 additions and 88 deletions
wrappers/python/simtk/openmm/app/modeller.py
View file @ 5120bb0d
......@@ -264,7 +264,7 @@ class Modeller(object):
2. Water molecules are removed if their distance to any solute atom is less than the sum of their van der Waals radii.
3. If the solute is charged and neutralize=True, enough positive or negative ions are added to neutralize it. Each ion is added by
randomly selecting a water molecule and replacing it with the ion.
4. Ion pairs are added to give the requested total ionic strength.
4. Ion pairs are added to give the requested total ionic strength. Note that only monovalent ions are currently supported.
The box size can be specified in any of several ways:
......@@ -298,6 +298,7 @@ class Modeller(object):
ionicStrength : concentration=0*molar
the total concentration of ions (both positive and negative) to add. This
does not include ions that are added to neutralize the system.
Note that only monovalent ions are currently supported.
neutralize : bool=True
whether to add ions to neutralize the system
"""
......@@ -522,7 +523,7 @@ class Modeller(object):
# Add ions based on the desired ionic strength.
numIons = len(addedWaters)*ionicStrength/(55.4*molar) # Pure water is about 55.4 molar (depending on temperature)
numPairs = int(floor(numIons/2+0.5))
numPairs = int(floor(numIons+0.5))
for i in range(numPairs):
addIon(positiveElement)
for i in range(numPairs):
......@@ -1103,7 +1104,7 @@ class Modeller(object):
else:
a2 = newAtoms[matchingAtoms[atom2]]
newTopology.addBond(a1, a2)
for bond in self.topology.bonds():
if bond[0] in newAtoms and bond[1] in newAtoms:
newTopology.addBond(newAtoms[bond[0]], newAtoms[bond[1]])
......@@ -1135,6 +1136,6 @@ class Modeller(object):
delta *= (distance-length)/length
newPositions[atom1] -= weights[0]*delta
newPositions[atom2] += weights[1]*delta
self.topology = newTopology
self.positions = newPositions
wrappers/python/tests/TestModeller.py
View file @ 5120bb0d
......@@ -15,14 +15,14 @@ class TestModeller(unittest.TestCase):
# load the alanine dipeptide pdb file
self.pdb = PDBFile('systems/alanine-dipeptide-explicit.pdb')
self.topology_start = self.pdb.topology
self.positions = self.pdb.positions
self.positions = self.pdb.positions
self.forcefield = ForceField('amber10.xml', 'tip3p.xml')
# load the T4-lysozyme-L99A receptor pdb file
self.pdb2 = PDBFile('systems/lysozyme-implicit.pdb')
self.topology_start2 = self.pdb2.topology
self.positions2 = self.pdb2.positions
# load the metallothionein pdb file
self.pdb3 = PDBFile('systems/1T2Y.pdb')
self.topology_start3 = self.pdb3.topology
......@@ -30,12 +30,12 @@ class TestModeller(unittest.TestCase):
def test_deleteWater(self):
""" Test the deleteWater() method. """
# build the chain dictionary
chain_dict = {0:0}
# 749 water chains are deleted
chain_delta = -749
# Build the residue and atom dictionaries for validate_preserved.
# Also, count the number of deleted residues and atoms.
residues_preserved = 0
......@@ -55,38 +55,38 @@ class TestModeller(unittest.TestCase):
residue_delta -= 1
for atom in residue.atoms():
atom_delta -= 1
modeller = Modeller(self.topology_start, self.positions)
modeller.deleteWater()
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)
validate_deltas(self, self.topology_start, topology_after,
validate_deltas(self, self.topology_start, topology_after,
chain_delta, residue_delta, atom_delta)
def test_delete(self):
""" Test the delete() method. """
modeller = Modeller(self.topology_start, self.positions)
topology_before = modeller.getTopology()
# Create the list of items to be deleted.
# Start with the first 50 water chains
chains = [chain for chain in topology_before.chains()]
toDelete = chains[1:51]
# Next add water residues 103->152 to the list of items to be deleted
residues = [residue for residue in topology_before.residues()]
toDelete.extend(residues[103:153])
# Finally add water atoms 622->771 to the list of items to be deleted
atoms = [atom for atom in topology_before.atoms()]
toDelete.extend(atoms[622:772])
modeller.delete(toDelete)
topology_after = modeller.getTopology()
# build the chain dictionary
chain_dict = {0:0}
for i in range(1,51):
......@@ -95,7 +95,7 @@ class TestModeller(unittest.TestCase):
chain_dict[i+100] = i
for i in range(101, 600):
chain_dict[i+150] = i
# build the residue dictionary
residue_dict = {}
for i in range(3):
......@@ -106,7 +106,7 @@ class TestModeller(unittest.TestCase):
residue_dict[i+100] = i
for i in range(103, 602):
residue_dict[i+150] = i
# build the atom dictionary
atom_dict = {}
for i in range(22):
......@@ -117,37 +117,37 @@ class TestModeller(unittest.TestCase):
atom_dict[i+300] = i
for i in range(322,1819):
atom_dict[i+450] = i
validate_preserved(self, topology_before, topology_after, chain_dict, residue_dict, atom_dict)
chain_delta = -150
residue_delta = -150
atom_delta = -450
validate_deltas(self, topology_before, topology_after, chain_delta, residue_delta, atom_delta)
def test_add(self):
""" Test the add() method. """
# load the methanol-box pdb file
pdb2 = PDBFile('systems/methanol-box.pdb')
topology_toAdd = pdb2.topology
positions_toAdd = pdb2.positions
modeller = Modeller(self.topology_start, self.positions)
modeller.deleteWater()
topology_before = modeller.getTopology()
modeller.add(topology_toAdd, positions_toAdd)
topology_after = modeller.getTopology()
# build the first chain dictionary for the first call of validate_preserved()
chain_counter = 0
chain_dict = {}
for chain in topology_before.chains():
chain_dict[chain.index] = chain_counter
chain_counter += 1
# build the residue and atom dictionaries for the first call of validate_preserved()
residue_counter = 0
residue_dict = {}
......@@ -159,13 +159,13 @@ class TestModeller(unittest.TestCase):
for atom in residue.atoms():
atom_dict[atom.index] = atom_counter
atom_counter += 1
# 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)
# 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.
# build the chain dictionary
chain_delta = 0
chain_dict = {}
......@@ -173,7 +173,7 @@ class TestModeller(unittest.TestCase):
chain_dict[chain.index] = chain_counter
chain_counter += 1
chain_delta += 1
# build the residue and atom dictionaries for the second call of validate_preserved
residue_delta = 0
residue_dict = {}
......@@ -187,26 +187,26 @@ class TestModeller(unittest.TestCase):
atom_dict[atom.index] = atom_counter
atom_counter += 1
atom_delta += 1
# validate that the items in the added topology are preserved
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_deltas(self, topology_before, topology_after, chain_delta, residue_delta, atom_delta)
def test_convertWater(self):
""" Test the convertWater() method. """
for model in ['tip3p', 'spce', 'tip4pew', 'tip5p']:
if model == 'tip5p':
firstmodel = 'tip4pew'
else:
firstmodel = 'tip5p'
modeller = Modeller(self.topology_start, self.positions)
modeller.convertWater(model=firstmodel)
modeller.convertWater(model=model)
topology_after = modeller.getTopology()
for residue in topology_after.residues():
if residue.name == "HOH":
oatom = [atom for atom in residue.atoms() if atom.element == element.oxygen]
......@@ -221,7 +221,7 @@ class TestModeller(unittest.TestCase):
self.assertTrue(len(matoms)==1 and len(m1atoms)==0 and len(m2atoms)==0)
elif model=='tip5p':
self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1)
# build the chain dictionary for validate_preserved
chain_counter = 0
chain_dict = {}
......@@ -229,7 +229,7 @@ class TestModeller(unittest.TestCase):
for chain in self.topology_start.chains():
chain_dict[chain.index] = chain_counter
chain_counter += 1
# build the residue and atom dictionaries for validate_preserved
residue_counter = 0
residue_dict = {}
......@@ -249,16 +249,16 @@ class TestModeller(unittest.TestCase):
if residue.name == 'HOH' and model == 'tip5p':
atom_counter += 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)
validate_deltas(self, self.topology_start, topology_after,
validate_deltas(self, self.topology_start, topology_after,
chain_delta, residue_delta, atom_delta)
def test_addSolventWaterModels(self):
""" Test all addSolvent() method with all possible water models. """
topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers)
for model in ['tip3p', 'spce', 'tip4pew', 'tip5p']:
......@@ -270,16 +270,16 @@ class TestModeller(unittest.TestCase):
# add the solvent to get the "after" topology
modeller.addSolvent(forcefield, model=model)
topology_after = modeller.getTopology()
# First, check that everything that was there before has been preserved.
# build the chain dictionary for validate_preserved
chain_counter = 0
chain_dict = {0:0}
for chain in topology_before.chains():
chain_dict[chain.index] = chain_counter
chain_counter += 1
# build the residue and atom dictionaries for validate_preserved
residue_counter = 0
residue_dict = {}
......@@ -291,10 +291,10 @@ class TestModeller(unittest.TestCase):
for atom in residue.atoms():
atom_dict[atom.index] = atom_counter
atom_counter += 1
# 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)
# Make sure water that was added was the correct model
for residue in topology_after.residues():
if residue.name == 'HOH':
......@@ -310,10 +310,10 @@ class TestModeller(unittest.TestCase):
self.assertTrue(len(matoms)==1 and len(m1atoms)==0 and len(m2atoms)==0)
elif model=='tip5p':
self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1)
def test_addSolventPeriodicBox(self):
""" 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.
topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 4.5, 5.5)*nanometers)
......@@ -322,11 +322,11 @@ class TestModeller(unittest.TestCase):
modeller.addSolvent(self.forcefield)
topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.5, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.5, 0))
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()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
......@@ -334,11 +334,11 @@ class TestModeller(unittest.TestCase):
modeller.addSolvent(self.forcefield, boxSize = Vec3(3.6, 4.6, 5.6)*nanometers)
topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.6, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.6, 0))
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()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
......@@ -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)
topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.4, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0.5, 4.4, 0))
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()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
......@@ -358,11 +358,11 @@ class TestModeller(unittest.TestCase):
modeller.addSolvent(self.forcefield, padding = 1.0*nanometers)
topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(2.8802, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 2.8802, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 2.8802))
# Fifth way: specify a number of molecules to add instead of a box size
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
......@@ -395,7 +395,7 @@ class TestModeller(unittest.TestCase):
total_added = water_count+sodium_count+chlorine_count
self.assertEqual(total_added, 1364)
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(chlorine_count, expected_ions)
......@@ -417,7 +417,7 @@ class TestModeller(unittest.TestCase):
topology_toAdd.addResidue('CL', newChain)
residues = [residue for residue in topology_toAdd.residues()]
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),
Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
modeller.add(topology_toAdd, positions_toAdd)
......@@ -437,7 +437,7 @@ class TestModeller(unittest.TestCase):
total_water_ions = water_count+sodium_count+chlorine_count
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
self.assertEqual(sodium_count, expected_sodium)
self.assertEqual(chlorine_count, expected_chlorine)
......@@ -460,7 +460,7 @@ class TestModeller(unittest.TestCase):
topology_toAdd.addResidue('NA', newChain)
residues = [residue for residue in topology_toAdd.residues()]
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),
Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
......@@ -482,14 +482,14 @@ class TestModeller(unittest.TestCase):
total_water_ions = water_count+sodium_count+chlorine_count
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
self.assertEqual(sodium_count, expected_sodium)
self.assertEqual(chlorine_count, expected_chlorine)
def test_addSolventIons(self):
""" Test the addSolvent() method with all possible choices for positive and negative ions. """
topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers)
......@@ -500,19 +500,19 @@ class TestModeller(unittest.TestCase):
positions_nowater = modeller.getPositions()
expected_ion_fraction = 1.0*molar/(55.4*molar)
for positiveIon in ['Cs+', 'K+', 'Li+', 'Na+', 'Rb+']:
ionName = positiveIon[:-1].upper()
modeller = Modeller(topology_nowater, positions_nowater)
modeller.addSolvent(self.forcefield, positiveIon=positiveIon, ionicStrength=1.0*molar)
topology_after = modeller.getTopology()
water_count = 0
water_count = 0
positive_ion_count = 0
chlorine_count = 0
for residue in topology_after.residues():
if residue.name=='HOH':
water_count += 1
water_count += 1
elif residue.name==ionName:
positive_ion_count += 1
elif residue.name=='CL':
......@@ -520,7 +520,7 @@ class TestModeller(unittest.TestCase):
total_added = water_count+positive_ion_count+chlorine_count
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(chlorine_count, expected_ions)
......@@ -532,7 +532,7 @@ class TestModeller(unittest.TestCase):
topology_after = modeller.getTopology()
water_count = 0
sodium_count = 0
sodium_count = 0
negative_ion_count = 0
for residue in topology_after.residues():
if residue.name=='HOH':
......@@ -541,16 +541,16 @@ class TestModeller(unittest.TestCase):
sodium_count += 1
elif residue.name==ionName:
negative_ion_count += 1
total_added = water_count+sodium_count+negative_ion_count
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(chlorine_count, expected_ions)
def test_addHydrogensPdb2(self):
""" Test the addHydrogens() method on the T4-lysozyme-L99A pdb file. """
# build the Modeller
topology_start = self.topology_start2
positions = self.positions2
......@@ -576,7 +576,7 @@ class TestModeller(unittest.TestCase):
def test_addHydrogensPdb3(self):
""" Test the addHydrogens() method on the metallothionein pdb file. """
# build the Modeller
topology_start = self.topology_start3
positions = self.positions3
......@@ -671,7 +671,7 @@ class TestModeller(unittest.TestCase):
index_list_CYS = [31, 49, 110, 135, 171, 193, 229]
atoms = [atom for atom in topology2.atoms()]
toDelete2 = []
for index in index_list_CYS:
for index in index_list_CYS:
self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index])
modeller2.delete(toDelete2)
......@@ -690,7 +690,7 @@ class TestModeller(unittest.TestCase):
modeller = Modeller(topology_start, positions)
# 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)
# Create a variants list to force the one histidine to be of the right variation.
......@@ -709,13 +709,13 @@ class TestModeller(unittest.TestCase):
modeller.addHydrogens(self.forcefield, variants=variants)
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.
index_list_GLH = [85, 192, 387, 731, 992, 1018, 1718, 2042]
atoms = [atom for atom in topology_GLH.atoms()]
toDelete2 = []
for index in index_list_GLH:
self.assertTrue(atoms[index].element.symbol=='H')
self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index])
modeller.delete(toDelete2)
topology_GLU = modeller.getTopology()
......@@ -811,7 +811,7 @@ class TestModeller(unittest.TestCase):
self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index])
for index in index_list_GLH:
self.assertTrue(atoms[index].element.symbol=='H')
self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index])
modeller.delete(toDelete2)
topology_ASP_GLU = modeller.getTopology()
......@@ -847,13 +847,13 @@ class TestModeller(unittest.TestCase):
index_list_CYS = [31, 49, 110, 135, 171, 193, 229]
atoms = [atom for atom in topology2.atoms()]
toDelete2 = []
for index in index_list_CYS:
for index in index_list_CYS:
self.assertTrue(atoms[index].element.symbol=='H')
toDelete2.append(atoms[index])
modeller2.delete(toDelete2)
topology_after = modeller2.getTopology()
validate_equivalence(self, topology_CYX, topology_after)
validate_equivalence(self, topology_CYX, topology_after)
def test_addHydrogenspH11(self):
""" Test of addHydrogens() with pH=11. """
......@@ -870,7 +870,7 @@ class TestModeller(unittest.TestCase):
modeller.delete(toDelete)
# 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)
# 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
......@@ -982,7 +982,7 @@ class TestModeller(unittest.TestCase):
topology.addAtom('C', element.carbon, residue)
topology.addAtom('N', element.nitrogen, residue)
topology.addAtom('V', element.oxygen, residue)
# Add the virtual sites.
......@@ -1083,5 +1083,3 @@ class TestModeller(unittest.TestCase):
if __name__ == '__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