Fix the AmberCustomGBForce.getStandardParameters and make it work with

mm.GBSAOBCForce as well.

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

@@ -988,23 +988,28 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
             gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         else:
             raise ValueError("Illegal value specified for implicit solvent model")
-        gb_parms = gb.getStandardParameters(gbmodel, elements)
+        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 self._raw_data['SCREEN']]
+            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 self._raw_data['RADII']]
+        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[0]) > 1e-4 or abs(s - gb_parms[1]) > 1e-4:
+            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[0], gb_parms[1] = r, s
+            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

@@ -33,7 +33,7 @@ from __future__ import division, absolute_import
 
 from collections import defaultdict
 import copy
-import simtk.openmm.app.element as E
+from simtk.openmm.app import element as E
 from simtk.openmm import CustomGBForce, Continuous2DFunction
 import simtk.unit as u
 
@@ -226,8 +226,8 @@ def _is_carboxylateO(atom, all_bonds):
 
 def _bondi_radii(topology):
     """ Sets the bondi radii """
-    radii = [0.0 for atom in topology.getAtoms()]
-    for i, atom in enumerate(topology.getAtoms()):
+    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):
@@ -252,9 +252,9 @@ def _bondi_radii(topology):
 
 def _mbondi_radii(topology):
     """ Sets the mbondi radii """
-    radii = [0.0 for atom in topology.getAtoms()]
+    radii = [0.0 for atom in topology.atoms()]
     all_bonds = _get_bonded_atom_list(topology)
-    for i, atom in enumerate(topology.getAtoms()):
+    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]
@@ -288,9 +288,9 @@ def _mbondi_radii(topology):
 
 def _mbondi2_radii(topology):
     """ Sets the mbondi2 radii """
-    radii = [0.0 for atom in topology.getAtoms()]
+    radii = [0.0 for atom in topology.atoms()]
     all_bonds = _get_bonded_atom_list(topology)
-    for i, atom in enumerate(topology.getAtoms()):
+    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]
@@ -324,7 +324,7 @@ 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.getAtoms()):
+    for i, atom in enumerate(topology.atoms()):
         # carboxylate and HH/HE (ARG)
         if _is_carboxylateO(atom, all_bonds):
             radii[i] = 1.4
@@ -448,7 +448,7 @@ class GBSAHCTForce(CustomAmberGBForce):
     @staticmethod
     def getStandardParameters(topology):
         radii = [[x/10] for x in _mbondi_radii(topology)]
-        for i, atom in zip(topology.getAtoms()):
+        for i, atom in enumerate(topology.atoms()):
             radii[i].append(_screen_parameter(atom)[0])
         return radii
 
@@ -477,7 +477,7 @@ class GBSAOBC1Force(CustomAmberGBForce):
     @staticmethod
     def getStandardParameters(topology):
         radii = [[x/10] for x in _mbondi2_radii(topology)]
-        for i, atom in zip(topology.getAtoms()):
+        for i, atom in enumerate(topology.atoms()):
             radii[i].append(_screen_parameter(atom)[0])
         return radii
 
@@ -546,7 +546,7 @@ class GBSAGBnForce(CustomAmberGBForce):
     @staticmethod
     def getStandardParameters(topology):
         radii = [[x/10] for x in _bondi_radii(topology)]
-        for i, atom in zip(topology.getAtoms()):
+        for i, atom in enumerate(topology.atoms()):
             radii[i].append(_screen_parameter(atom)[1])
         return radii
 
@@ -588,7 +588,7 @@ class GBSAGBn2Force(GBSAGBnForce):
     @staticmethod
     def getStandardParameters(topology):
         radii = [[x/10] for x in _mbondi3_radii(topology)]
-        for i, atom in zip(topology.getAtoms()):
+        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])