Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Added option for box shape when building solvent boxes #3480

Merged
merged 1 commit into from
Feb 27, 2022
Merged

Added option for box shape when building solvent boxes #3480

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 33 additions & 13 deletions wrappers/python/openmm/app/modeller.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2021 Stanford University and the Authors.
Portions copyright (c) 2012-2022 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
Expand Down Expand Up @@ -48,7 +48,7 @@
import sys
import xml.etree.ElementTree as etree
from copy import deepcopy
from math import ceil, floor
from math import ceil, floor, sqrt
from collections import defaultdict

class Modeller(object):
Expand Down Expand Up @@ -375,8 +375,8 @@ def _addIons(self, forcefield, numWaters, replaceableMols, ionCutoff=0.05*nanome
self.topology = modeller.topology
self.positions = modeller.positions

def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, padding=None, numAdded=None, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar, neutralize=True):
"""Add solvent (both water and ions) to the model to fill a rectangular box.
def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, padding=None, numAdded=None, boxShape='cube', positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar, neutralize=True):
"""Add solvent (both water and ions) to the model to fill a periodic box.
The algorithm works as follows:
Expand All @@ -390,12 +390,18 @@ def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, p
1. You can explicitly give the vectors defining the periodic box to use.
2. Alternatively, for a rectangular box you can simply give the dimensions of the unit cell.
3. You can give a padding distance. The largest dimension of the solute (along the x, y, or z axis) is determined, and a cubic
box of size (largest dimension)+2*padding is used.
4. You can specify the total number of molecules (both waters and ions) to add. A cubic box is then created whose size is
3. You can give a padding distance. A bounding sphere containing the solute is determined, and the box size is
set to (sphere radius)+(padding). This guarantees no atom in the solute will come closer than the padding
distance to any atom of another periodic copy.
4. You can specify the total number of molecules (both waters and ions) to add. A box is then created whose size is
just large enough hold the specified amount of solvent.
5. Finally, if none of the above options is specified, the existing Topology's box vectors are used.
When specifying either a padding distance or a number of molecules, you can specify a shape for the periodic box:
cubic, rhombic dodecahedron, or truncated octahedron. Using a non-rectangular box allows the same distance
between periodic copies to be achieved with a smaller box. The most compact option is a rhombic dodecahedron,
for which the box volume is 70.7% the volume of a cubic box with the same amount of padding.
Parameters
----------
forcefield : ForceField
Expand All @@ -410,6 +416,9 @@ def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, p
the padding distance to use
numAdded : int=None
the total number of molecules (waters and ions) to add
boxShape: str='cube'
the box shape to use. Allowed values are 'cube', 'dodecahedron', and 'octahedron'. If padding and numAdded
are both None, this is ignored.
positiveIon : string='Na+'
the type of positive ion to add. Allowed values are 'Cs+', 'K+', 'Li+', 'Na+', and 'Rb+'
negativeIon : string='Cl-'
Expand Down Expand Up @@ -449,7 +458,7 @@ def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, p
if numAdded is not None:
# Select a padding distance which is guaranteed to give more than the specified number of molecules.

padding = 1.1*(numAdded/((len(pdbResidues)/pdbBoxSize[0]**3)*8))**(1.0/3.0)
padding = 2.2*(numAdded/((len(pdbResidues)/pdbBoxSize[0]**3)*8))**(1.0/3.0)
if padding < 0.5:
padding = 0.5 # Ensure we have enough when adding very small numbers of molecules
if boxVectors is not None:
Expand All @@ -466,12 +475,23 @@ def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, p
if is_quantity(padding):
padding = padding.value_in_unit(nanometer)
if len(self.positions) == 0:
maxSize = 0
radius = 0
else:
positions = self.positions.value_in_unit(nanometer)
minRange = Vec3(*(min((pos[i] for pos in positions)) for i in range(3)))
maxRange = Vec3(*(max((pos[i] for pos in positions)) for i in range(3)))
center = 0.5*(minRange+maxRange)
radius = max(unit.norm(center-pos) for pos in positions)
width = 2*radius+padding
box = width*Vec3(1, 1, 1)
if boxShape == 'cube':
vectors = (Vec3(width, 0, 0), Vec3(0, width, 0), Vec3(0, 0, width))
elif boxShape == 'dodecahedron':
vectors = (Vec3(width, 0, 0), Vec3(0, width, 0), Vec3(0.5, 0.5, 0.5*sqrt(2))*width)
elif boxShape == 'octahedron':
vectors = (Vec3(width, 0, 0), Vec3(1/3, 2*sqrt(2)/3, 0)*width, Vec3(-1/3, sqrt(2)/3, sqrt(6)/3)*width)
else:
maxSize = max(max((pos[i] for pos in self.positions))-min((pos[i] for pos in self.positions)) for i in range(3))
maxSize = maxSize.value_in_unit(nanometer)
box = (maxSize+2*padding)*Vec3(1, 1, 1)
vectors = (Vec3(maxSize+2*padding, 0, 0), Vec3(0, maxSize+2*padding, 0), Vec3(0, 0, maxSize+2*padding))
raise ValueError(f'Illegal box shape: {boxShape}')
else:
box = self.topology.getUnitCellDimensions().value_in_unit(nanometer)
vectors = self.topology.getPeriodicBoxVectors().value_in_unit(nanometer)
Expand Down
32 changes: 23 additions & 9 deletions wrappers/python/tests/TestModeller.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
from collections import defaultdict
import unittest
import math
import sys

from validateModeller import *
from openmm.app import *
Expand Down Expand Up @@ -333,7 +331,6 @@ def test_addSolventPeriodicBox(self):
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)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, boxSize = Vec3(3.6, 4.6, 5.6)*nanometers)
Expand All @@ -345,7 +342,6 @@ def test_addSolventPeriodicBox(self):
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)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, boxVectors = (Vec3(3.4, 0, 0), Vec3(0.5, 4.4, 0), Vec3(-1.0, -1.5, 5.4))*nanometers)
Expand All @@ -357,25 +353,43 @@ def test_addSolventPeriodicBox(self):
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)
modeller.deleteWater()
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))
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(1.924363, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 1.924363, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 1.924363))

# 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)
modeller.deleteWater()
numInitial = len(list(modeller.topology.residues()))
modeller.addSolvent(self.forcefield, numAdded=1000)
self.assertEqual(numInitial+1000, len(list(modeller.topology.residues())))

def test_addSolventBoxShape(self):
"""Test the addSolvent() method; test the different box shapes."""
modeller = Modeller(self.pdb.topology, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, padding=1.0*nanometers, boxShape='cube')
cubeVectors = modeller.getTopology().getPeriodicBoxVectors()
modeller = Modeller(self.pdb.topology, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, padding=1.0*nanometers, boxShape='dodecahedron')
dodecVectors = modeller.getTopology().getPeriodicBoxVectors()
modeller = Modeller(self.pdb.topology, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, padding=1.0*nanometers, boxShape='octahedron')
octVectors = modeller.getTopology().getPeriodicBoxVectors()
cubeVolume = cubeVectors[0][0]*cubeVectors[1][1]*cubeVectors[2][2]/(nanometers**3)
dodecVolume = dodecVectors[0][0]*dodecVectors[1][1]*dodecVectors[2][2]/(nanometers**3)
octVolume = octVectors[0][0]*octVectors[1][1]*octVectors[2][2]/(nanometers**3)
self.assertAlmostEqual(0.707, dodecVolume/cubeVolume, places=3)
self.assertAlmostEqual(0.770, octVolume/cubeVolume, places=3)

def test_addSolventNeutralSolvent(self):
""" Test the addSolvent() method; test adding ions to neutral solvent. """

Expand Down