Added option for box shape when building solvent boxes (#3480)

wrappers/python/openmm/app/modeller.py

@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 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:
 
@@ -48,7 +48,7 @@ import random
 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):
@@ -375,8 +375,8 @@ class Modeller(object):
             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:
 
@@ -390,12 +390,18 @@ class Modeller(object):
 
         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
@@ -410,6 +416,9 @@ class Modeller(object):
             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-'
@@ -449,7 +458,7 @@ class Modeller(object):
         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:
@@ -466,12 +475,23 @@ class Modeller(object):
             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)

wrappers/python/tests/TestModeller.py

 from collections import defaultdict
 import unittest
-import math
-import sys
 
 from validateModeller import *
 from openmm.app import *
@@ -333,7 +331,6 @@ class TestModeller(unittest.TestCase):
         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)
@@ -345,7 +342,6 @@ class TestModeller(unittest.TestCase):
         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)
@@ -357,25 +353,43 @@ class TestModeller(unittest.TestCase):
         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. """