Finish implementing charmm parser for systems with lonepairs and/or with Drude particles.

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

@@ -12,7 +12,7 @@ the ParmEd program and was ported for use with OpenMM.
 Copyright (c) 2014 the Authors
 
 Author: Jason M. Swails
-Contributors:
+Contributors: Jing Huang
 Date: Sep. 17, 2014
 
 Permission is hereby granted, free of charge, to any person obtaining a
@@ -497,21 +497,22 @@ class CharmmParameterSet(object):
                 except IndexError:
                     raise CharmmFileError('Could not parse NBFIX terms.')
                 self.nbfix_types[(min(at1, at2), max(at1, at2))] = (emin, rmin)
-        # Here parse the possible nbthole section
-        if section == 'NBTHOLE':
-            words = line.split()
-            try:
-                at1 = words[0]
-                at2 = words[1]
-                nbt = abs(conv(words[2], float, 'NBTHOLE a'))
+                continue
+            # Here parse the possible nbthole section
+            if section == 'NBTHOLE':
+                words = line.split()
                 try:
-                    self.atom_types_str[at1].add_nbthole(at2, nbt)
-                    self.atom_types_str[at2].add_nbthole(at1, nbt)
-                except KeyError:
-                    pass
-            except IndexError:
-                raise CharmmFileError('Could not parse NBTHOLE terms.')
-            self.nbthole_types[(min(at1, at2), max(at1, at2))] = (nbt)
+                    at1 = words[0]
+                    at2 = words[1]
+                    nbt = abs(conv(words[2], float, 'NBTHOLE a'))
+                    try:
+                        self.atom_types_str[at1].add_nbthole(at2, nbt)
+                        self.atom_types_str[at2].add_nbthole(at1, nbt)
+                    except KeyError:
+                        pass
+                except IndexError:
+                    raise CharmmFileError('Could not parse NBTHOLE terms.')
+                self.nbthole_types[(min(at1, at2), max(at1, at2))] = (nbt)
         # If there were any CMAP terms stored in the parameter set, the last one
         # defined will not have been added to the set. Add it now.
         if current_cmap is not None:

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

@@ -869,6 +869,13 @@ class CharmmPsfFile(object):
                         raise StopIteration
         except StopIteration:
             pass
+        # test if the system containing the Drude particles
+        has_drude_particle = False
+        try:
+            if self.drudeconsts_list:
+                has_drude_particle = True
+        except AttributeError:
+            pass
 
         # Set up the constraints
         if verbose and (constraints is not None and not rigidWater):
@@ -1219,7 +1226,7 @@ class CharmmPsfFile(object):
                 cforce.addParticle((i - 1,)) # adjust for indexing from 0
 
         # Add NBTHOLE terms
-        if IsDrudePSF and has_nbthole_terms:
+        if has_drude_particle and has_nbthole_terms:
             nbt_idx_list = [0 for atom in self.atom_list]
             nbt_alpha_list = [0 for atom in self.atom_list] # only save alpha for NBThole pairs
             num_nbt_types = 0 
@@ -1314,7 +1321,7 @@ class CharmmPsfFile(object):
             idxa = lpsite[0]
             parent_exclude_list[idx].append(idxa)
             force.addException(idx, idxa, 0.0, 0.1, 0.0)
-        if IsDrudePSF:
+        if has_drude_particle:
             for pair in self.drudepair_list:
                 idx = pair[0]
                 idxa = pair[1]
@@ -1346,7 +1353,7 @@ class CharmmPsfFile(object):
         system.addForce(force)
 
         # Add Drude particles (Drude force)
-        if IsDrudePSF:
+        if has_drude_particle:
             if verbose: print('Adding Drude force and Thole screening...')
             drudeforce = mm.DrudeForce()
             drudeforce.setForceGroup(7)

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

@@ -179,8 +179,8 @@ class AtomType(object):
         self.nbfix[typename] = (rmin, epsilon, rmin14, epsilon14)
 
     def add_nbthole(self, typename, nbt):
-            """ Adds a new NBTHOLE screening factor for this atom """
-            self.nbthole[typename] = (nbt)
+        """ Adds a new NBTHOLE screening factor for this atom """
+        self.nbthole[typename] = (nbt)
 
     def __str__(self):
         return self.name