Add support for the Amber GB forces. Put in the correct radii determination

(again, code pulled from ParmEd) and update the tests to include this.

wrappers/python/simtk/openmm/app/charmm/psf.py

@@ -797,7 +797,9 @@ class ProteinStructure(object):
 
     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:
@@ -811,6 +813,7 @@ class ProteinStructure(object):
                 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]
@@ -861,9 +864,13 @@ class ProteinStructure(object):
                     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 = [rad * length_conv for rad in self.parm_data['RADII']]
 
         if gb_model is GBn2:
             return zip(radii, screen, alpha, beta, gamma)
@@ -1486,6 +1493,120 @@ def _set_owner(atom_list, owner_array, atm, mol_id):
 
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
+# Routines that set the necessary radii lists based on a list of atoms with
+# their connectivities
+
+def _bondi_radii(atom_list):
+    """ Sets the bondi radii """
+    radii = [0.0 for atom in atom_list]
+    for i, atom in enumerate(atom_list):
+        if atom.type.atomic_number == 6:
+            radii[i] = 1.7
+        elif atom.type.atomic_number == 1:
+            radii[i] = 1.2
+        elif atom.type.atomic_number == 7:
+            radii[i] = 1.55
+        elif atom.type.atomic_number == 8:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 9:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 14:
+            radii[i] = 2.1
+        elif atom.type.atomic_number == 15:
+            radii[i] = 1.85
+        elif atom.type.atomic_number == 16:
+            radii[i] = 1.8
+        elif atom.type.atomic_number == 17:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # converted to nanometers above
+
+def _mbondi_radii(atom_list):
+    """ Sets the mbondi radii """
+    radii = [0.0 for atom in atom_list]
+    for i, atom in enumerate(atom_list):
+        # Radius of H atom depends on element it is bonded to
+        if atom.type.atomic_number == 1:
+            bondeds = list(atom.bond_partners)
+            if bondeds[0].type.atomic_number in (6, 7): # C or N
+                radii[i] = 1.3
+            elif bondeds[0].type.atomic_number in (8, 16): # O or S
+                radii[i] = 0.8
+            else:
+                radii[i] = 1.2
+        # Radius of C atom depends on what type it is
+        elif atom.atomic_number == 6:
+            radii[i] = 1.7
+        # All other elements have fixed radii for all types/partners
+        elif atom.type.atomic_number == 7:
+            radii[i] = 1.55
+        elif atom.type.atomic_number == 8:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 9:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 14:
+            radii[i] = 2.1
+        elif atom.type.atomic_number == 15:
+            radii[i] = 1.85
+        elif atom.type.atomic_number == 16:
+            radii[i] = 1.8
+        elif atom.type.atomic_number == 17:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # converted to nanometers above
+
+def _mbondi2_radii(atom_list):
+    """ Sets the mbondi2 radii """
+    radii = [0.0 for atom in atom_list]
+    for i, atom in enumerate(atom_list):
+        # Radius of H atom depends on element it is bonded to
+        if atom.type.atomic_number == 1:
+            if atom.bond_partners[0].type.atomic_number == 7:
+                radii[i] = 1.3
+            else:
+                radii[i] = 1.2
+        # Radius of C atom depends on what type it is
+        elif atom.type.atomic_number == 6:
+            radii[i] = 1.7
+        # All other elements have fixed radii for all types/partners
+        elif atom.type.atomic_number == 7:
+            radii[i] = 1.55
+        elif atom.type.atomic_number == 8:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 9:
+            radii[i] = 1.5
+        elif atom.type.atomic_number == 14:
+            radii[i] = 2.1
+        elif atom.type.atomic_number == 15:
+            radii[i] = 1.85
+        elif atom.type.atomic_number == 16:
+            radii[i] = 1.8
+        elif atom.type.atomic_number == 17:
+            radii[i] = 1.5
+        else:
+            radii[i] = 1.5
+    return radii  # Converted to nanometers above
+
+def _mbondi3_radii(atom_list):
+    """ Sets the mbondi3 radii """
+    radii = _mbondi2_radii(atom_list)
+    for i, atom in enumerate(atom_list):
+        # Adjust OE (GLU), OD (ASP) and HH/HE (ARG)
+        if atom.residue.resname in ('GLU', 'ASP'):
+            if atom.name.startswith('OE') or atom.name.startswith('OD'):
+                radii[i] = 1.4
+        elif atom.residue.resname == 'ARG':
+            if atom.name.startswith('HH') or atom.name.startswith('HE'):
+                radii[i] = 1.17
+        # Adjust carboxylate O radii on C-termini
+        if atom.name == 'OXT':
+            radii[i] = 1.4
+    return radii  # Converted to nanometers above
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
 if __name__ == '__main__':
     import doctest
     doctest.testmod()

wrappers/python/tests/TestCharmmFiles.py

@@ -62,32 +62,35 @@ class TestCharmmFiles(unittest.TestCase):
             system = self.psf_c.createSystem(removeCMMotion=b)
             self.assertEqual(any(isinstance(f, CMMotionRemover) for f in system.getForces()), b)
 
-#   def test_ImplicitSolvent(self):
-#       """Test implicit solvent using the implicitSolvent parameter.
-
-#       """
-#       system = self.psf_v.createSystem(implicitSolvent=OBC2)
-#       self.assertTrue(any(isinstance(f, CustomGBForce) for f in system.getForces()))
-
-#   def test_ImplicitSolventParameters(self):
-#       """Test that solventDielectric and soluteDielectric are passed correctly.
-
-#       """
-#       system = self.psf_x.createSystem(implicitSolvent=GBn,
-#                                        solventDielectric=50.0, 
-#                                        soluteDielectric = 0.9)
-#       found_matching_solvent_dielectric=False
-#       found_matching_solute_dielectric=False
-#       for force in system.getForces():
-#           if isinstance(force, CustomGBForce):
-#               if force.getSolventDielectric() == 50.0:
-#                   found_matching_solvent_dielectric = True
-#               if force.getSoluteDielectric() == 0.9:
-#                   found_matching_solute_dielectric = True
-#           if isinstance(force, NonbondedForce):
-#               self.assertEqual(force.getReactionFieldDielectric(), 1.0)
-#       self.assertTrue(found_matching_solvent_dielectric and 
-#                       found_matching_solute_dielectric)
+    def test_ImplicitSolvent(self):
+        """Test implicit solvent using the implicitSolvent parameter.
+
+        """
+        system = self.psf_v.createSystem(implicitSolvent=OBC2)
+        self.assertTrue(any(isinstance(f, CustomGBForce) for f in system.getForces()))
+
+    def test_ImplicitSolventParameters(self):
+        """Test that solventDielectric and soluteDielectric are passed correctly.
+
+        """
+        system = self.psf_x.createSystem(implicitSolvent=GBn,
+                                         solventDielectric=50.0, 
+                                         soluteDielectric = 0.9)
+        found_matching_solvent_dielectric=False
+        found_matching_solute_dielectric=False
+        for force in system.getForces():
+            if isinstance(force, CustomGBForce):
+                for i in range(force.getNumGlobalParameters()):
+                    if force.getGlobalParameterName(i) == 'solventDielectric':
+                        if force.getGlobalParameterDefaultValue(i) == 50.0:
+                            found_matching_solvent_dielectric = True
+                    elif force.getGlobalParameterName(i) == 'soluteDielectric':
+                        if force.getGlobalParameterDefaultValue(i) == 0.9:
+                            found_matching_solute_dielectric = True
+            if isinstance(force, NonbondedForce):
+                self.assertEqual(force.getReactionFieldDielectric(), 1.0)
+        self.assertTrue(found_matching_solvent_dielectric and 
+                        found_matching_solute_dielectric)
 
     def test_HydrogenMass(self):
         """Test that altering the mass of hydrogens works correctly."""