Merge remote-tracking branch 'upstream/master' into forcefield

wrappers/python/simtk/openmm/app/charmmpsffile.py

@@ -32,9 +32,7 @@ DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
 OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
 USE OR OTHER DEALINGS IN THE SOFTWARE.
 """
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
+from __future__ import division, absolute_import, print_function
 
 from functools import wraps
 from math import pi, cos, sin, sqrt
@@ -666,88 +664,6 @@ class CharmmPsfFile(object):
 
         return topology
 
-    def _get_gb_params(self, gb_model=HCT):
-        """ Gets the GB parameters. Need this method to special-case GB neck """
-        screen = [0 for atom in self.atom_list]
-        if gb_model is GBn:
-            radii = _bondi_radii(self.atom_list)
-            screen = [0.5 for atom in self.atom_list]
-            for i, atom in enumerate(self.atom_list):
-                if atom.type.atomic_number == 6:
-                    screen[i] = 0.48435382330
-                elif atom.type.atomic_number == 1:
-                    screen[i] = 1.09085413633
-                elif atom.type.atomic_number == 7:
-                    screen[i] = 0.700147318409
-                elif atom.type.atomic_number == 8:
-                    screen[i] = 1.06557401132
-                elif atom.type.atomic_number == 16:
-                    screen[i] = 0.602256336067
-        elif gb_model is GBn2:
-            radii = _mbondi3_radii(self.atom_list)
-            # Add non-optimized values as defaults
-            alpha = [1.0 for i in self.atom_list]
-            beta = [0.8 for i in self.atom_list]
-            gamma = [4.85 for i in self.atom_list]
-            screen = [0.5 for i in self.atom_list]
-            for i, atom in enumerate(self.atom_list):
-                if atom.type.atomic_number == 6:
-                    screen[i] = 1.058554
-                    alpha[i] = 0.733756
-                    beta[i] = 0.506378
-                    gamma[i] = 0.205844
-                elif atom.type.atomic_number == 1:
-                    screen[i] = 1.425952
-                    alpha[i] = 0.788440
-                    beta[i] = 0.798699
-                    gamma[i] = 0.437334
-                elif atom.type.atomic_number == 7:
-                    screen[i] = 0.733599
-                    alpha[i] = 0.503364
-                    beta[i] = 0.316828
-                    gamma[i] = 0.192915
-                elif atom.type.atomic_number == 8:
-                    screen[i] = 1.061039
-                    alpha[i] = 0.867814
-                    beta[i] = 0.876635
-                    gamma[i] = 0.387882
-                elif atom.type.atomic_number == 16:
-                    screen[i] = -0.703469
-                    alpha[i] = 0.867814
-                    beta[i] = 0.876635
-                    gamma[i] = 0.387882
-        else:
-            # Set the default screening parameters
-            for i, atom in enumerate(self.atom_list):
-                if atom.type.atomic_number == 1:
-                    screen[i] = 0.85
-                elif atom.type.atomic_number == 6:
-                    screen[i] = 0.72
-                elif atom.type.atomic_number == 7:
-                    screen[i] = 0.79
-                elif atom.type.atomic_number == 8:
-                    screen[i] = 0.85
-                elif atom.type.atomic_number == 9:
-                    screen[i] = 0.88
-                elif atom.type.atomic_number == 15:
-                    screen[i] = 0.86
-                elif atom.type.atomic_number == 16:
-                    screen[i] = 0.96
-                else:
-                    screen[i] = 0.8
-            # Determine which radii set we need
-            if gb_model is OBC1 or gb_model is OBC2:
-                radii = _mbondi2_radii(self.atom_list)
-            elif gb_model is HCT:
-                radii = _mbondi_radii(self.atom_list)
-
-        length_conv = u.angstrom.conversion_factor_to(u.nanometer)
-        radii = [x * length_conv for x in radii]
-
-        if gb_model is GBn2:
-            return zip(radii, screen, alpha, beta, gamma)
-        return zip(radii, screen)
-
     def createSystem(self, params, nonbondedMethod=ff.NoCutoff,
                      nonbondedCutoff=1.0*u.nanometer,
                      switchDistance=0.0*u.nanometer,
@@ -1249,8 +1165,6 @@ class CharmmPsfFile(object):
         # Add GB model if we're doing one
         if implicitSolvent is not None:
             if verbose: print('Adding GB parameters...')
-            gb_parms = self._get_gb_params(implicitSolvent)
-
             # If implicitSolventKappa is None, compute it from salt
             # concentration
             if implicitSolventKappa is None:
@@ -1289,6 +1203,7 @@ class CharmmPsfFile(object):
             elif implicitSolvent is GBn2:
                 gb = GBSAGBn2Force(solventDielectric, soluteDielectric, None,
                                    cutoff, kappa=implicitSolventKappa)
+            gb_parms = gb.getStandardParameters(self.topology)
             for atom, gb_parm in zip(self.atom_list, gb_parms):
                 gb.addParticle([atom.charge] + list(gb_parm))
             # Set cutoff method

wrappers/python/simtk/openmm/app/gromacstopfile.py

@@ -259,7 +259,7 @@ class GromacsTopFile(object):
         """Process the [ defaults ] line."""
         fields = line.split()
         if len(fields) < 4:
-            raise ValueError('Too few fields in [ defaults ] line: '+line);
+            raise ValueError('Too few fields in [ defaults ] line: '+line)
         if fields[0] != '1':
             raise ValueError('Unsupported nonbonded type: '+fields[0])
         if fields[1] != '2':
@@ -272,7 +272,7 @@ class GromacsTopFile(object):
         """Process a line in the [ moleculetypes ] category."""
         fields = line.split()
         if len(fields) < 1:
-            raise ValueError('Too few fields in [ moleculetypes ] line: '+line);
+            raise ValueError('Too few fields in [ moleculetypes ] line: '+line)
         type = GromacsTopFile._MoleculeType()
         self._moleculeTypes[fields[0]] = type
         self._currentMoleculeType = type
@@ -281,7 +281,7 @@ class GromacsTopFile(object):
         """Process a line in the [ molecules ] category."""
         fields = line.split()
         if len(fields) < 2:
-            raise ValueError('Too few fields in [ molecules ] line: '+line);
+            raise ValueError('Too few fields in [ molecules ] line: '+line)
         self._molecules.append((fields[0], int(fields[1])))
 
     def _processAtom(self, line):
@@ -290,7 +290,7 @@ class GromacsTopFile(object):
             raise ValueError('Found [ atoms ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 5:
-            raise ValueError('Too few fields in [ atoms ] line: '+line);
+            raise ValueError('Too few fields in [ atoms ] line: '+line)
         self._currentMoleculeType.atoms.append(fields)
 
     def _processBond(self, line):
@@ -299,9 +299,9 @@ class GromacsTopFile(object):
             raise ValueError('Found [ bonds ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 3:
-            raise ValueError('Too few fields in [ bonds ] line: '+line);
+            raise ValueError('Too few fields in [ bonds ] line: '+line)
         if fields[2] != '1':
-            raise ValueError('Unsupported function type in [ bonds ] line: '+line);
+            raise ValueError('Unsupported function type in [ bonds ] line: '+line)
         self._currentMoleculeType.bonds.append(fields)
 
     def _processAngle(self, line):
@@ -310,9 +310,9 @@ class GromacsTopFile(object):
             raise ValueError('Found [ angles ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 4:
-            raise ValueError('Too few fields in [ angles ] line: '+line);
+            raise ValueError('Too few fields in [ angles ] line: '+line)
         if fields[3] not in ('1', '5'):
-            raise ValueError('Unsupported function type in [ angles ] line: '+line);
+            raise ValueError('Unsupported function type in [ angles ] line: '+line)
         self._currentMoleculeType.angles.append(fields)
 
     def _processDihedral(self, line):
@@ -321,9 +321,9 @@ class GromacsTopFile(object):
             raise ValueError('Found [ dihedrals ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 5:
-            raise ValueError('Too few fields in [ dihedrals ] line: '+line);
+            raise ValueError('Too few fields in [ dihedrals ] line: '+line)
         if fields[4] not in ('1', '2', '3', '4', '9'):
-            raise ValueError('Unsupported function type in [ dihedrals ] line: '+line);
+            raise ValueError('Unsupported function type in [ dihedrals ] line: '+line)
         self._currentMoleculeType.dihedrals.append(fields)
 
     def _processExclusion(self, line):
@@ -332,7 +332,7 @@ class GromacsTopFile(object):
             raise ValueError('Found [ exclusions ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 2:
-            raise ValueError('Too few fields in [ exclusions ] line: '+line);
+            raise ValueError('Too few fields in [ exclusions ] line: '+line)
         self._currentMoleculeType.exclusions.append(fields)
 
     def _processPair(self, line):
@@ -341,9 +341,9 @@ class GromacsTopFile(object):
             raise ValueError('Found [ pairs ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 3:
-            raise ValueError('Too few fields in [ pairs ] line: '+line);
+            raise ValueError('Too few fields in [ pairs ] line: '+line)
         if fields[2] != '1':
-            raise ValueError('Unsupported function type in [ pairs ] line: '+line);
+            raise ValueError('Unsupported function type in [ pairs ] line: '+line)
         self._currentMoleculeType.pairs.append(fields)
 
     def _processCmap(self, line):
@@ -352,14 +352,14 @@ class GromacsTopFile(object):
             raise ValueError('Found [ cmap ] section before [ moleculetype ]')
         fields = line.split()
         if len(fields) < 6:
-            raise ValueError('Too few fields in [ cmap ] line: '+line);
+            raise ValueError('Too few fields in [ cmap ] line: '+line)
         self._currentMoleculeType.cmaps.append(fields)
 
     def _processAtomType(self, line):
         """Process a line in the [ atomtypes ] category."""
         fields = line.split()
         if len(fields) < 6:
-            raise ValueError('Too few fields in [ atomtypes ] line: '+line);
+            raise ValueError('Too few fields in [ atomtypes ] line: '+line)
         if len(fields[3]) == 1:
             # Bonded type and atomic number are both missing.
             fields.insert(1, None)
@@ -377,27 +377,27 @@ class GromacsTopFile(object):
         """Process a line in the [ bondtypes ] category."""
         fields = line.split()
         if len(fields) < 5:
-            raise ValueError('Too few fields in [ bondtypes ] line: '+line);
+            raise ValueError('Too few fields in [ bondtypes ] line: '+line)
         if fields[2] != '1':
-            raise ValueError('Unsupported function type in [ bondtypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ bondtypes ] line: '+line)
         self._bondTypes[tuple(fields[:2])] = fields
 
     def _processAngleType(self, line):
         """Process a line in the [ angletypes ] category."""
         fields = line.split()
         if len(fields) < 6:
-            raise ValueError('Too few fields in [ angletypes ] line: '+line);
+            raise ValueError('Too few fields in [ angletypes ] line: '+line)
         if fields[3] not in ('1', '5'):
-            raise ValueError('Unsupported function type in [ angletypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ angletypes ] line: '+line)
         self._angleTypes[tuple(fields[:3])] = fields
 
     def _processDihedralType(self, line):
         """Process a line in the [ dihedraltypes ] category."""
         fields = line.split()
         if len(fields) < 7:
-            raise ValueError('Too few fields in [ dihedraltypes ] line: '+line);
+            raise ValueError('Too few fields in [ dihedraltypes ] line: '+line)
         if fields[4] not in ('1', '2', '3', '4', '9'):
-            raise ValueError('Unsupported function type in [ dihedraltypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ dihedraltypes ] line: '+line)
         key = tuple(fields[:5])
         if fields[4] == '9' and key in self._dihedralTypes:
             # There are multiple dihedrals defined for these atom types.
@@ -409,25 +409,25 @@ class GromacsTopFile(object):
         """Process a line in the [ implicit_genborn_params ] category."""
         fields = line.split()
         if len(fields) < 6:
-            raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line);
+            raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line)
         self._implicitTypes[fields[0]] = fields
 
     def _processPairType(self, line):
         """Process a line in the [ pairtypes ] category."""
         fields = line.split()
         if len(fields) < 5:
-            raise ValueError('Too few fields in [ pairtypes] line: '+line);
+            raise ValueError('Too few fields in [ pairtypes] line: '+line)
         if fields[2] != '1':
-            raise ValueError('Unsupported function type in [ pairtypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ pairtypes ] line: '+line)
         self._pairTypes[tuple(fields[:2])] = fields
 
     def _processCmapType(self, line):
         """Process a line in the [ cmaptypes ] category."""
         fields = line.split()
         if len(fields) < 8 or len(fields) < 8+int(fields[6])*int(fields[7]):
-            raise ValueError('Too few fields in [ cmaptypes ] line: '+line);
+            raise ValueError('Too few fields in [ cmaptypes ] line: '+line)
         if fields[5] != '1':
-            raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line);
+            raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line)
         self._cmapTypes[tuple(fields[:5])] = fields
 
     def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None, includeDir=None, defines=None):

wrappers/python/simtk/openmm/app/internal/amber_file_parser.py

@@ -34,8 +34,7 @@ DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
 OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
 USE OR OTHER DEALINGS IN THE SOFTWARE.
 """
-from __future__ import absolute_import
-from __future__ import print_function
+from __future__ import absolute_import, print_function
 
 #=============================================================================================
 # GLOBAL IMPORTS
@@ -553,87 +552,6 @@ class PrmtopLoader(object):
             self._excludedAtoms.append(atomList)
         return self._excludedAtoms
 
-    def getGBParms(self, gbmodel, elements):
-        """Return list giving GB params, Radius and screening factor"""
-        gb_List=[]
-        radii=[float(r) for r in self._raw_data["RADII"]]
-        screen=[float(s) for s in self._raw_data["SCREEN"]]
-        if gbmodel == 'GBn2':
-            alpha = [0 for i in self._raw_data['RADII']]
-            beta = [0 for i in self._raw_data['RADII']]
-            gamma = [0 for i in self._raw_data['RADII']]
-        # Update screening parameters for GBn if specified
-        if gbmodel == 'GBn':
-            if elements is None:
-                raise Exception('GBn model requires element information')
-            for i, element in enumerate(elements):
-                if element is elem.carbon:
-                    screen[i] = 0.48435382330
-                elif element is elem.hydrogen:
-                    screen[i] = 1.09085413633
-                elif element is elem.nitrogen:
-                    screen[i] = 0.700147318409
-                elif element is elem.oxygen:
-                    screen[i] = 1.06557401132
-                elif element is elem.sulfur:
-                    screen[i] = 0.602256336067
-                else:
-                    screen[i] = 0.5
-                # radii is currently in Angstroms right now. GBn lookup tables
-                # only support radii between 1.0 and 2.0
-                if radii[i] < 1.0 or radii[i] > 2.0:
-                    raise ValueError('GBn requires intrinsic radii between 1 and '
-                                     '2 Angstroms (%.3f found for atom %d)' %
-                                     (radii[i], i))
-        if gbmodel == 'GBn2':
-            if elements is None:
-                raise Exception('GBn2 model requires element information')
-            for i, element in enumerate(elements):
-                if element is elem.carbon:
-                    screen[i] = 1.058554
-                    alpha[i] = 0.733756
-                    beta[i] = 0.506378
-                    gamma[i] = 0.205844
-                elif element is elem.hydrogen:
-                    screen[i] = 1.425952
-                    alpha[i] = 0.788440
-                    beta[i] = 0.798699
-                    gamma[i] = 0.437334
-                elif element is elem.nitrogen:
-                    screen[i] = 0.733599
-                    alpha[i] = 0.503364
-                    beta[i] = 0.316828
-                    gamma[i] = 0.192915
-                elif element is elem.oxygen:
-                    screen[i] = 1.061039
-                    alpha[i] = 0.867814
-                    beta[i] = 0.876635
-                    gamma[i] = 0.387882
-                elif element is elem.sulfur:
-                    screen[i] = -0.703469
-                    alpha[i] = 0.867814
-                    beta[i] = 0.876635
-                    gamma[i] = 0.387882
-                else: # not optimized
-                    screen[i] = 0.5
-                    alpha[i] = 1.0
-                    beta[i] = 0.8
-                    gamma[i] = 4.85
-                # radii is currently in Angstroms right now. GBn lookup tables
-                # only support radii between 1.0 and 2.0
-                if radii[i] < 1.0 or radii[i] > 2.0:
-                    raise ValueError('GBn2 requires intrinsic radii between 1 and '
-                                     '2 Angstroms (%.3f found for atom %d)' %
-                                     (radii[i], i))
-        lengthConversionFactor = units.angstrom.conversion_factor_to(units.nanometer)
-        if gbmodel == 'GBn2':
-            for rad, scr, alp, bet, gam in zip(radii, screen, alpha, beta, gamma):
-                gb_List.append((rad*lengthConversionFactor, scr, alp, bet, gam))
-        else:
-            for rad, scr in zip(radii, screen):
-                gb_List.append((rad*lengthConversionFactor, scr))
-        return gb_List
-
     def getBoxBetaAndDimensions(self):
         """Return periodic boundary box beta angle and dimensions"""
         beta=float(self._raw_data["BOX_DIMENSIONS"][0])
@@ -1053,7 +971,6 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
             cutoff = nonbondedCutoff
             if units.is_quantity(cutoff):
                 cutoff = cutoff.value_in_unit(units.nanometers)
-        gb_parms = prmtop.getGBParms(gbmodel, elements)
         if gbmodel == 'HCT':
             gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         elif gbmodel == 'OBC1':
@@ -1070,7 +987,29 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
         elif gbmodel == 'GBn2':
             gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         else:
-            raise Exception("Illegal value specified for implicit solvent model")
+            raise ValueError("Illegal value specified for implicit solvent model")
+        if isinstance(gb, mm.GBSAOBCForce):
+            # Built-in GBSAOBCForce does not have getStandardParameters, so use
+            # the one from the equivalent CustomGBForce
+            gb_parms = customgb.GBSAOBC2Force.getStandardParameters(topology)
+        else:
+            gb_parms = type(gb).getStandardParameters(topology)
+        # Replace radii and screen, but screen *only* gets replaced by the
+        # prmtop contents for HCT, OBC1, and OBC2. GBn and GBn2 both override
+        # the prmtop screen factors from LEaP in sander and pmemd
+        if gbmodel in ('HCT', 'OBC1', 'OBC2'):
+            screen = [float(s) for s in prmtop._raw_data['SCREEN']]
+        else:
+            screen = [gb_parm[1] for gb_parm in gb_parms]
+        radii = [float(r)/10 for r in prmtop._raw_data['RADII']]
+        warned = False
+        for i, (r, s) in enumerate(zip(radii, screen)):
+            if abs(r - gb_parms[i][0]) > 1e-4 or abs(s - gb_parms[i][1]) > 1e-4:
+                if not warned:
+                    warnings.warn('Non-optimal GB parameters detected for GB '
+                                  'model %s' % gbmodel)
+                    warned = True
+            gb_parms[i][0], gb_parms[i][1] = r, s
 
         for charge, gb_parm in zip(charges, gb_parms):
             if gbmodel == 'OBC2' and implicitSolventKappa == 0:

wrappers/python/simtk/openmm/app/internal/customgbforces.py

@@ -29,10 +29,11 @@ OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
 USE OR OTHER DEALINGS IN THE SOFTWARE.
 """
 
-from __future__ import division
-from __future__ import absolute_import
+from __future__ import division, absolute_import
 
+from collections import defaultdict
 import copy
+from simtk.openmm.app import element as E
 from simtk.openmm import CustomGBForce, Continuous2DFunction
 import simtk.unit as u
 
@@ -196,6 +197,141 @@ for i in range (len(d0)):
     d0[i]=d0[i]/10
     m0[i]=m0[i]*10
 
+def _get_bonded_atom_list(topology):
+    """ Returns a list of atoms bonded to each other atom in a dict """
+    bondeds = defaultdict(list)
+    for a1, a2 in topology.bonds():
+        bondeds[a1].append(a2)
+        bondeds[a2].append(a1)
+    return bondeds
+
+def _is_carboxylateO(atom, all_bonds):
+    if atom is not E.oxygen: return False
+    bondeds = all_bonds[atom]
+    if len(bondeds) != 1:
+        return False
+    if bondeds[0].element is not E.carbon:
+        return False
+    bondedsC = all_bonds[bondeds[0]]
+    if len(bondedsC) != 3:
+        return False
+    for a3 in bondedsC:
+        if a3 is atom: continue
+        if a3.element is E.oxygen:
+            break
+    else:
+        return False
+    # If we got here, must be a carboxylate
+    return True
+
+def _bondi_radii(topology):
+    """ Sets the bondi radii """
+    radii = [0.0 for atom in topology.atoms()]
+    for i, atom in enumerate(topology.atoms()):
+        if atom.element is E.carbon:
+            radii[i] = 1.7
+        elif atom.element in (E.hydrogen, E.deuterium):
+            radii[i] = 1.2
+        elif atom.element is E.nitrogen:
+            radii[i] = 1.55
+        elif atom.element is E.oxygen:
+            radii[i] = 1.5
+        elif atom.element is E.fluorine:
+            radii[i] = 1.5
+        elif atom.element is E.silicon:
+            radii[i] = 2.1
+        elif atom.element is E.phosphorus:
+            radii[i] = 1.85
+        elif atom.element is E.sulfur:
+            radii[i] = 1.8
+        elif atom.element is E.chlorine:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # converted to nanometers above
+
+def _mbondi_radii(topology):
+    """ Sets the mbondi radii """
+    radii = [0.0 for atom in topology.atoms()]
+    all_bonds = _get_bonded_atom_list(topology)
+    for i, atom in enumerate(topology.atoms()):
+        # Radius of H atom depends on element it is bonded to
+        if atom.element in (E.hydrogen, E.deuterium):
+            bondeds = all_bonds[atom]
+            if bondeds[0].element in (E.carbon, E.nitrogen):
+                radii[i] = 1.3
+            elif bondeds[0].type.atomic_number in (E.oxygen, E.sulfur):
+                radii[i] = 0.8
+            else:
+                radii[i] = 1.2
+        # Radius of C atom depends on what type it is
+        elif atom.element is E.carbon:
+            radii[i] = 1.7
+        # All other elements have fixed radii for all types/partners
+        elif atom.element is E.nitrogen:
+            radii[i] = 1.55
+        elif atom.element is E.oxygen:
+            radii[i] = 1.5
+        elif atom.element is E.fluorine:
+            radii[i] = 1.5
+        elif atom.element is E.silicon:
+            radii[i] = 2.1
+        elif atom.element is E.phosphorus:
+            radii[i] = 1.85
+        elif atom.element is E.sulfur:
+            radii[i] = 1.8
+        elif atom.element is E.chlorine:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # converted to nanometers above
+
+def _mbondi2_radii(topology):
+    """ Sets the mbondi2 radii """
+    radii = [0.0 for atom in topology.atoms()]
+    all_bonds = _get_bonded_atom_list(topology)
+    for i, atom in enumerate(topology.atoms()):
+        # Radius of H atom depends on element it is bonded to
+        if atom.element in (E.hydrogen, E.deuterium):
+            bondeds = all_bonds[atom]
+            if bondeds[0].element is E.nitrogen:
+                radii[i] = 1.3
+            else:
+                radii[i] = 1.2
+        # Radius of C atom depends on what type it is
+        elif atom.element is E.carbon:
+            radii[i] = 1.7
+        # All other elements have fixed radii for all types/partners
+        elif atom.element is E.nitrogen:
+            radii[i] = 1.55
+        elif atom.element is E.oxygen:
+            radii[i] = 1.5
+        elif atom.element is E.fluorine:
+            radii[i] = 1.5
+        elif atom.element is E.silicon:
+            radii[i] = 2.1
+        elif atom.element == E.phosphorus:
+            radii[i] = 1.85
+        elif atom.element == E.sulfur:
+            radii[i] = 1.8
+        elif atom.element == E.chlorine:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # Converted to nanometers above
+
+def _mbondi3_radii(topology):
+    """ Sets the mbondi3 radii """
+    radii = _mbondi2_radii(topology)
+    all_bonds = _get_bonded_atom_list(topology)
+    for i, atom in enumerate(topology.atoms()):
+        # carboxylate and HH/HE (ARG)
+        if _is_carboxylateO(atom, all_bonds):
+            radii[i] = 1.4
+        elif atom.residue.name == 'ARG':
+            if atom.name.startswith('HH') or atom.name.startswith('HE'):
+                radii[i] = 1.17
+    return radii  # Converted to nanometers above
 
 def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, kappa, offset):
     # Add the energy terms to the CustomGBForce.  These are identical for all the GB models.
@@ -231,6 +367,22 @@ def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, k
             force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*(1/f-"+str(1/cutoff)+");"
                                 "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
 
+_SCREEN_PARAMETERS = { # normal, GBn, GBn2
+        E.hydrogen : (0.85, 1.09085413633, 1.425952),
+        E.carbon : (0.72, 0.48435382330, 1.058554),
+        E.nitrogen : (0.79, 0.700147318409, 0.733599),
+        E.oxygen : (0.85, 1.06557401132, 1.061039),
+        E.fluorine : (0.88, 0.5, 0.5),
+        E.phosphorus : (0.86, 0.5, 0.5),
+        E.sulfur : (0.96, 0.602256336067, -0.703469),
+        None : (0.8, 0.5, 0.5) # default
+}
+_SCREEN_PARAMETERS[E.deuterium] = _SCREEN_PARAMETERS[E.hydrogen]
+def _screen_parameter(atom):
+    if atom.element in _SCREEN_PARAMETERS:
+        return _SCREEN_PARAMETERS[atom.element]
+    return _SCREEN_PARAMETERS[None]
+
 class CustomAmberGBForce(CustomGBForce):
 
     OFFSET = 0.009
@@ -254,6 +406,24 @@ class CustomAmberGBForce(CustomGBForce):
         CustomGBForce.addParticle(self, params)
         return params
 
+    @staticmethod
+    def getStandardParameters(topology):
+        """ Gets list of standard parameters for this GB model based on an input Topology
+
+        Parameters
+        ----------
+        topology : simtk.openmm.app.Topology
+            Topology of the system to get parameters for
+
+        Returns
+        -------
+        list of float
+            List of all parameters needed for this GB model. These can be passed
+            to addParticle or setParticleParameters after the charge is inserted
+            at the beginning of the list
+        """
+        raise NotImplementedError('Must override getStandardParameters in subclasses')
+
 """
 Amber Equivalent: igb = 1
 """
@@ -275,6 +445,13 @@ class GBSAHCTForce(CustomAmberGBForce):
         self.addComputedValue("B", "1/(1/or-I)", CustomGBForce.SingleParticle)
         _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
 
+    @staticmethod
+    def getStandardParameters(topology):
+        radii = [[x/10] for x in _mbondi_radii(topology)]
+        for i, atom in enumerate(topology.atoms()):
+            radii[i].append(_screen_parameter(atom)[0])
+        return radii
+
 """
 Amber Equivalents: igb = 2
 """
@@ -297,10 +474,17 @@ class GBSAOBC1Force(CustomAmberGBForce):
                                    "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
         _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
 
+    @staticmethod
+    def getStandardParameters(topology):
+        radii = [[x/10] for x in _mbondi2_radii(topology)]
+        for i, atom in enumerate(topology.atoms()):
+            radii[i].append(_screen_parameter(atom)[0])
+        return radii
+
 """
 Amber Equivalents: igb = 5
 """
-class GBSAOBC2Force(CustomAmberGBForce):
+class GBSAOBC2Force(GBSAOBC1Force):
 
     def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
                  cutoff=None, kappa=0.0):
@@ -359,6 +543,13 @@ class GBSAGBnForce(CustomAmberGBForce):
         if parameters[1] < 0.1 or parameters[1] > 0.2:
             raise ValueError('Radii must be between 1 and 2 Angstroms for neck lookup')
 
+    @staticmethod
+    def getStandardParameters(topology):
+        radii = [[x/10] for x in _bondi_radii(topology)]
+        for i, atom in enumerate(topology.atoms()):
+            radii[i].append(_screen_parameter(atom)[1])
+        return radii
+
 """
 Amber Equivalents: igb = 8
 """
@@ -393,3 +584,22 @@ class GBSAGBn2Force(GBSAGBnForce):
         self.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
                                      "psi=I*or; radius=or+offset; offset=0.0195141", CustomGBForce.SingleParticle)
         _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.0195141)
+
+    @staticmethod
+    def getStandardParameters(topology):
+        radii = [[x/10] for x in _mbondi3_radii(topology)]
+        for i, atom in enumerate(topology.atoms()):
+            radii[i].append(_screen_parameter(atom)[2])
+            if atom.element in (E.hydrogen, E.deuterium):
+                radii[i].extend([0.788440, 0.798699, 0.437334])
+            elif atom.element is E.carbon:
+                radii[i].extend([0.733756, 0.506378, 0.205844])
+            elif atom.element is E.nitrogen:
+                radii[i].extend([0.503364, 0.316828, 0.192915])
+            elif atom.element is E.oxygen:
+                radii[i].extend([0.867814, 0.876635, 0.387882])
+            elif atom.element is E.sulfur:
+                radii[i].extend([0.867814, 0.876635, 0.387882])
+            else:
+                radii[i].extend([0.8, 4.85, 0.5])
+        return radii

wrappers/python/simtk/openmm/app/modeller.py

@@ -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/simtk/openmm/app/topology.py

@@ -71,20 +71,17 @@ class Topology(object):
     def getNumAtoms(self):
         """Return the number of atoms in the Topology.
         """
-        natom = self._numAtoms
-        return natom
+        return self._numAtoms
 
     def getNumResidues(self):
         """Return the number of residues in the Topology.
         """
-        nres = self._numResidues
-        return nres
+        return self._numResidues
 
     def getNumChains(self):
         """Return the number of chains in the Topology.
         """
-        nchain = len(self._chains)
-        return nchain
+        return len(self._chains)
 
     def addChain(self, id=None):
         """Create a new Chain and add it to the Topology.

wrappers/python/tests/TestModeller.py

@@ -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()
-
-