Implement NBFIX in the CHARMM parsers.

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

@@ -13,7 +13,7 @@ Copyright (c) 2014 the Authors
 
 Author: Jason M. Swails
 Contributors:
-Date: April 18, 2014
+Date: July 17, 2014
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -123,6 +123,8 @@ class CharmmParameterSet(object):
             elif arg.endswith('.str'):
                 strs.append(arg)
             elif arg.endswith('.inp'):
+                # Only consider the file name (since the directory is likely
+                # "toppar" and will screw up file type detection)
                 fname = os.path.split(arg)[1]
                 if 'par' in fname:
                     pars.append(arg)
@@ -436,8 +438,26 @@ class CharmmParameterSet(object):
                 try:
                     at1 = words[0]
                     at2 = words[1]
-                    emin = conv(words[2], float, 'NBFIX Emin')
+                    emin = abs(conv(words[2], float, 'NBFIX Emin'))
                     rmin = conv(words[3], float, 'NBFIX Rmin')
+                    try:
+                        emin14 = abs(conv(words[4], float, 'NBFIX Emin 1-4'))
+                        rmin14 = conv(words[5], float, 'NBFIX Rmin 1-4')
+                    except IndexError:
+                        emin14 = rmin14 = None
+                    try:
+                        self.atom_types_str[at1].add_nbfix(at2, emin, rmin,
+                                                           emin14, rmin14)
+                        self.atom_types_str[at2].add_nbfix(at1, emin, rmin,
+                                                           emin14, rmin14)
+                    except KeyError:
+                        # Some stream files define NBFIX terms with an atom that
+                        # is defined in another toppar file that does not
+                        # necessarily have to be loaded. As a result, not every
+                        # NBFIX found here will necessarily need to be applied.
+                        # If we can't find a particular atom type, don't bother
+                        # adding that nbfix and press on
+                        pass
                 except IndexError:
                     raise CharmmFileError('Could not parse NBFIX terms.')
                 self.nbfix_types[(min(at1, at2), max(at1, at2))] = (emin, rmin)

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

@@ -995,11 +995,25 @@ class CharmmPsfFile(object):
                             u.kilojoule_per_mole)
         ene_conv = dihe_frc_conv
       
-        # Create the system
+        # Create the system and determine if any of our atoms have NBFIX (and
+        # therefore requires a CustomNonbondedForce instead)
+        typenames = set()
         system = mm.System()
         if verbose: print('Adding particles...')
         for atom in self.atom_list:
+            typenames.add(atom.type.name)
             system.addParticle(atom.mass)
+        has_nbfix_terms = False
+        typenames = list(typenames)
+        try:
+            for i, typename in enumerate(typenames):
+                typ = params.atom_types_str[typename]
+                for j in range(i, len(typenames)):
+                    if typenames[j] in typ.nbfix:
+                        has_nbfix_terms = True
+                        raise StopIteration
+        except StopIteration:
+            pass
         # Set up the constraints
         if verbose and (constraints is not None and not rigidWater):
             print('Adding constraints...')
@@ -1240,9 +1254,85 @@ class CharmmPsfFile(object):
 
         # Add per-particle nonbonded parameters (LJ params)
         sigma_scale = 2**(-1/6) * 2
-        for i, atm in enumerate(self.atom_list):
-            force.addParticle(atm.charge, sigma_scale*atm.type.rmin*length_conv,
-                              abs(atm.type.epsilon*ene_conv))
+        if not has_nbfix_terms:
+            for atm in self.atom_list:
+                force.addParticle(atm.charge, sigma_scale*atm.type.rmin*length_conv,
+                                  abs(atm.type.epsilon*ene_conv))
+        else:
+            for atm in self.atom_list:
+                force.addParticle(atm.charge, 1.0, 0.0)
+            # Now add the custom nonbonded force that implements NBFIX. First
+            # thing we need to do is condense our number of types
+            lj_idx_list = [0 for atom in self.atom_list]
+            lj_radii, lj_depths = [], []
+            num_lj_types = 0
+            lj_type_list = []
+            for i, atom in enumerate(self.atom_list):
+                atom = atom.type
+                if lj_idx_list[i]: continue # already assigned
+                num_lj_types += 1
+                lj_idx_list[i] = num_lj_types
+                ljtype = (atom.rmin, atom.epsilon)
+                lj_type_list.append(atom)
+                lj_radii.append(atom.rmin)
+                lj_depths.append(atom.epsilon)
+                for j in range(i+1, len(self.atom_list)):
+                    atom2 = self.atom_list[j].type
+                    if lj_idx_list[j] > 0: continue # already assigned
+                    if atom2 is atom:
+                        lj_idx_list[j] = num_lj_types
+                    elif not atom.nbfix:
+                        # Only non-NBFIXed atom types can be compressed
+                        ljtype2 = (atom2.rmin, atom2.epsilon)
+                        if ljtype == ljtype2:
+                            lj_idx_list[j] = num_lj_types
+            # Now everything is assigned. Create the A-coefficient and
+            # B-coefficient arrays
+            acoef = [0 for i in range(num_lj_types*num_lj_types)]
+            bcoef = acoef[:]
+            for i in range(num_lj_types):
+                for j in range(num_lj_types):
+                    namej = lj_type_list[j].name
+                    try:
+                        wdij, rij, wdij14, rij14 = lj_type_list[i].nbfix[namej]
+                    except KeyError:
+                        rij = (lj_radii[i] + lj_radii[j]) * length_conv
+                        wdij = sqrt(lj_depths[i] * lj_depths[j]) * ene_conv
+                    else:
+                        rij *= length_conv
+                        wdij *= ene_conv
+                    acoef[i+num_lj_types*j] = sqrt(wdij) * rij**6
+                    bcoef[i+num_lj_types*j] = 2 * wdij * rij**6
+            cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6;'
+                                             'a=acoef(type1, type2);'
+                                             'b=bcoef(type1, type2)')
+            cforce.addTabulatedFunction('acoef',
+                    mm.Discrete2DFunction(num_lj_types, num_lj_types, acoef))
+            cforce.addTabulatedFunction('bcoef',
+                    mm.Discrete2DFunction(num_lj_types, num_lj_types, bcoef))
+            cforce.addPerParticleParameter('type')
+            cforce.setForceGroup(self.NONBONDED_FORCE_GROUP)
+            if (nonbondedMethod is ff.PME or nonbondedMethod is ff.Ewald or
+                        nonbondedMethod is ff.CutoffPeriodic):
+                cforce.setNonbondedMethod(cforce.CutoffPeriodic)
+                cforce.setCutoffDistance(nonbondedCutoff)
+                cforce.setLongRangeCorrection(True)
+            elif nonbondedMethod is ff.NoCutoff:
+                cforce.setNonbondedMethod(cforce.NoCutoff)
+            elif nonbondedMethod is ff.CutoffNonPeriodic:
+                cforce.setNonbondedMethod(cforce.CutoffNonPeriodic)
+                cforce.setCutoffDistance(nonbondedCutoff)
+            else:
+                raise ValueError('Unrecognized nonbonded method')
+            if switchDistance and nonbondedMethod is not ff.NoCutoff:
+                # make sure it's legal
+                if switchDistance >= nonbondedCutoff:
+                    raise ValueError('switchDistance is too large compared '
+                                     'to the cutoff!')
+                    cforce.setUseSwitchingFunction(True)
+                    cforce.setSwitchingDistance(switchDistance)
+            for i in lj_idx_list:
+                cforce.addParticle(i - 1) # adjust for indexing from 0
 
         # Add 1-4 interactions
         excluded_atom_pairs = set() # save these pairs so we don't zero them out
@@ -1283,6 +1373,13 @@ class CharmmPsfFile(object):
                     continue
                 force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
         system.addForce(force)
+        # If we needed a CustomNonbondedForce, map all of the exceptions from
+        # the NonbondedForce to the CustomNonbondedForce
+        if has_nbfix_terms:
+            for i in range(force.getNumExceptions()):
+                ii, jj, q, eps, sig = force.getExceptionParameters(i)
+                cforce.addExclusion(ii, jj)
+            system.addForce(cforce)
 
         # Add GB model if we're doing one
         if implicitSolvent is not None:

wrappers/python/simtk/openmm/app/internal/charmm/_charmmfile.py

@@ -47,7 +47,6 @@ class CharmmFile(object):
     """
 
     def __init__(self, fname, mode='r'):
-        self.closed = False
         if mode not in ('r', 'w'):
             raise ValueError('Cannot open CharmmFile with mode "%s"' % mode)
         if mode == 'r':
@@ -58,6 +57,7 @@ class CharmmFile(object):
             self._handle = open(fname, mode)
         except IOError, e:
             raise CharmmFileError(str(e))
+        self.closed = False
         self.line_number = 0
 
     def write(self, *args, **kwargs):

wrappers/python/simtk/openmm/app/internal/charmm/topologyobjects.py

@@ -42,6 +42,32 @@ TINY = 1e-8
 
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
+def _tracking(fcn):
+    """ Decorator to indicate the list has changed """
+    def new_fcn(self, *args):
+        self.changed = True
+        return fcn(self, *args)
+    return new_fcn
+
+class TrackedList(list):
+    """
+    This creates a list type that allows you to see if anything has changed
+    """
+    def __init__(self, arg=[]):
+        self.changed = False
+        list.__init__(self, arg)
+
+    __delitem__ = _tracking(list.__delitem__)
+    append = _tracking(list.append)
+    extend = _tracking(list.extend)
+    __setitem__ = _tracking(list.__setitem__)
+
+# Python 3 does not have __delslice__, but make sure we override it for Python 2
+if hasattr(TrackedList, '__delslice__'):
+    TrackedList.__delslice__ = _tracking(TrackedList.__delslice__)
+
+# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
 class AtomType(object):
     """
     Atom types can either be compared by indexes or names. Can be assigned with
@@ -61,6 +87,9 @@ class AtomType(object):
         - _member_number (int, private) : The order in which this atom type
                 was 'added' this is used to make sure that atom types added
                 last have priority in assignment in the generated hash tables
+        - nbfix (dict) : Dictionary that maps nbfix terms with other atom types.
+                         Dict entries are (rmin, epsilon) -- precombined values
+                         for that particular atom pair
     Example:
     >>> at = AtomType('HA', 1, 1.008, 1)
     >>> at.name, at.number
@@ -93,6 +122,9 @@ class AtomType(object):
         self.atomic_number = atomic_number
         # We have no LJ parameters as of yet
         self.epsilon = self.rmin = self.epsilon_14 = self.rmin_14 = None
+        # Store each NBFIX term as a dict with the atom type string matching to
+        # a 2-element tuple that is rmin, epsilon
+        self.nbfix = dict()
 
     def __eq__(self, other):
         """
@@ -123,6 +155,12 @@ class AtomType(object):
         """ The integer representation of an AtomType is its index """
         return self.number
 
+    def add_nbfix(self, typename, rmin, epsilon, rmin14, epsilon14):
+        """ Adds a new NBFIX exclusion for this atom """
+        if rmin14 is None: rmin14 = rmin
+        if epsilon14 is None: epsilon14 = epsilon
+        self.nbfix[typename] = (rmin, epsilon, rmin14, epsilon14)
+
     def __str__(self):
         return self.name
 
@@ -311,7 +349,7 @@ class Atom(object):
 
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-class AtomList(list):
+class AtomList(TrackedList):
     """ A list of Atom instances.  """
 
     def unmark(self):
@@ -424,7 +462,7 @@ class ResidueList(list):
         elif (self._last_residue != (resname, resnum) or
               system != self._last_residue.system):
             if (self._last_residue.idx == resnum and
-                system == self._last_residue.system):
+                self._last_residue.system == system):
                 lresname = self._last_residue.resname
                 warnings.warn('Residue %d split into separate residues %s '
                               'and %s' % (resnum, lresname, resname),
@@ -1047,28 +1085,6 @@ class _CmapGrid(object):
 
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-def _tracking(fcn):
-    """ Decorator to indicate the list has changed """
-    def new_fcn(self, *args):
-        self.changed = True
-        return fcn(self, *args)
-    return new_fcn
-
-class TrackedList(list):
-    """
-    This creates a list type that allows you to see if anything has changed
-    """
-    def __init__(self, arg=[]):
-        self.changed = False
-        list.__init__(self, arg)
-
-    __delitem__ = _tracking(list.__delitem__)
-    append = _tracking(list.append)
-    extend = _tracking(list.extend)
-    __setitem__ = _tracking(list.__setitem__)
-
-# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
 if __name__ == '__main__':
     import doctest
     doctest.testmod()