Update the python api for parsing CHARMM par/top/psf files.

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

@@ -114,6 +114,7 @@ class CharmmParameterSet(object):
         self.improper_types = dict()
         self.cmap_types = dict()
         self.nbfix_types = dict()
+        self.nbthole_types = dict()
         self.parametersets = []
 
         # Load all of the files
@@ -261,6 +262,9 @@ class CharmmParameterSet(object):
             if line.startswith('NBFIX'):
                 section = 'NBFIX'
                 continue
+            if line.startswith('THOLE'):
+                section = 'NBTHOLE'
+                continue
             if line.startswith('HBOND'):
                 section = None
                 continue
@@ -493,6 +497,21 @@ 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'))
+                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

@@ -57,7 +57,7 @@ from simtk.openmm.app.internal.charmm.exceptions import (
 import warnings
 
 TINY = 1e-8
-WATNAMES = ('WAT', 'HOH', 'TIP3', 'TIP4', 'TIP5', 'SPCE', 'SPC')
+WATNAMES = ('WAT', 'HOH', 'TIP3', 'TIP4', 'TIP5', 'SPCE', 'SPC', 'SWM4', 'SWM6')
 if sys.version_info >= (3, 0):
     xrange = range
 
@@ -236,19 +236,31 @@ class CharmmPsfFile(object):
             atom = residue_list.add_atom(system, resid, resname, name,
                             attype, charge, mass, inscode, props)
             atom_list.append(atom)
+            if IsDrudePSF:
+                alpha = conv(words[9], float, 'alpha')
+                thole = conv(words[10], float, 'thole')
+		drudeconsts_list.append([alpha, thole])
         atom_list.assign_indexes()
         atom_list.changed = False
         # Now get the number of bonds
         nbond = conv(psfsections['NBOND'][0], int, 'number of bonds')
         holder = psfsections['NBOND'][1]
         bond_list = TrackedList()
+        if IsDrudePSF:
+            drudepair_list = TrackedList()
         if len(holder) != nbond * 2:
             raise CharmmPSFError('Got %d indexes for %d bonds' %
                                  (len(holder), nbond))
         for i in range(nbond):
             id1 = holder[2*i  ] - 1
             id2 = holder[2*i+1] - 1
-            bond_list.append(Bond(atom_list[id1], atom_list[id2]))
+            # ignore any bond pair involving Drude or lonepairs: possible using atom's prop
+            if (atom_list[id1].name[0]=='D' or atom_list[id2].name[0]=='D'):
+                drudepair_list.append([min(id1,id2), max(id1,id2)])
+            elif (atom_list[id1].name[0:2]=='LP' or atom_list[id2].name[0:2]=='LP' or atom_list[id1].name=='OM' or atom_list[id2].name=='OM'):
+                pass
+            else:
+                bond_list.append(Bond(atom_list[id1], atom_list[id2]))
         bond_list.changed = False
         # Now get the number of angles and the angle list
         ntheta = conv(psfsections['NTHETA'][0], int, 'number of angles')
@@ -349,12 +361,48 @@ class CharmmPsfFile(object):
             if molecule_list != holder:
                 warnings.warn('Detected PSF molecule section that is WRONG. '
                               'Resetting molecularity.', CharmmPSFWarning)
-            # We have a CHARMM PSF file; now do NUMLP/NUMLPH sections
-            numlp, numlph = psfsections['NUMLP NUMLPH'][0]
-            if numlp != 0 or numlph != 0:
-                raise NotImplementedError('Cannot currently handle PSFs with lone '
-                                          'pairs defined in the NUMLP/NUMLPH '
-                                          'section.')
+        # We have a CHARMM PSF file; now do NUMLP/NUMLPH sections
+        numlp, numlph = psfsections['NUMLP NUMLPH'][0]
+        holder = psfsections['NUMLP NUMLPH'][1]
+        lonepair_list = TrackedList()
+        if numlp != 0 or numlph != 0:
+            lp_distance_list=[]
+            lp_angle_list=[]
+            lp_dihe_list=[]
+            for i in range(numlp):
+                lpline = holder[i].split()
+                if len(lpline)!=6 or lpline[0] != '3' or lpline[2] != 'F' or int(lpline[1]) != 4*i+1 :
+                    raise CharmmPSFError('Lonepair format error')
+                else:
+                    lp_distance_list.append(float(lpline[3]))
+                    lp_angle_list.append(float(lpline[4]))
+                    lp_dihe_list.append(float(lpline[5]))
+            for i in range(numlp):
+                lpline = holder[int(i/2)+numlp].split()
+                icolumn = (i%2) * 4
+                id1=int(lpline[icolumn])  -1
+                id2=int(lpline[icolumn+1])-1
+                id3=int(lpline[icolumn+2])-1
+                id4=int(lpline[icolumn+3])-1
+                lonepair_list.append([id1, id2, id3, id4, lp_distance_list[i], lp_angle_list[i], lp_dihe_list[i]])
+        # In Drude psf, here comes anisotropic section
+        if IsDrudePSF:
+            numaniso = psfsections['NUMANISO'][0]
+            holder = psfsections['NUMANISO'][1]
+            aniso_list = TrackedList()
+            if numaniso != 0:
+               k_list=[]
+               for i in range(numaniso):
+                   lpline = holder[i].split()
+                   k_list.append([float(lpline[0]),float(lpline[1]),float(lpline[2])])
+               for i in range(numaniso):
+                    lpline = holder[int(i/2)+numaniso].split()
+                    icolumn = (i%2) * 4
+                    id1=int(lpline[icolumn])  -1
+                    id2=int(lpline[icolumn+1])-1
+                    id3=int(lpline[icolumn+2])-1
+                    id4=int(lpline[icolumn+3])-1
+                    aniso_list.append([id1, id2, id3, id4, k_list[i][0], k_list[i][1], k_list[i][2]])
         # Now do the CMAPs
         ncrterm = conv(psfsections['NCRTERM'][0], int, 'Number of cross-terms')
         holder = psfsections['NCRTERM'][1]
@@ -385,6 +433,11 @@ class CharmmPsfFile(object):
         self.dihedral_list = dihedral_list
         self.dihedral_parameter_list = TrackedList()
         self.improper_list = improper_list
+        self.lonepair_list = lonepair_list
+        if IsDrudePSF:
+            self.drudeconsts_list = drudeconsts_list
+            self.drudepair_list = drudepair_list
+            self.aniso_list = aniso_list
         self.cmap_list = cmap_list
         self.donor_list = donor_list
         self.acceptor_list = acceptor_list
@@ -459,8 +512,8 @@ class CharmmPsfFile(object):
             # blank line) as well as any sections that have 0 members.
             line = psf.readline().strip()
         data = []
-        if title == 'NATOM' or title == 'NTITLE':
-            # Store these two sections as strings (ATOM section we will parse
+        if title == 'NATOM' or title == 'NTITLE' or title == 'NUMLP NUMLPH' or title == 'NUMANISO':
+            # Store these four sections as strings (ATOM section we will parse
             # later). The rest of the sections are integer pointers
             while line:
                 data.append(line)
@@ -787,7 +840,7 @@ class CharmmPsfFile(object):
                             u.kilojoule_per_mole)
         ene_conv = dihe_frc_conv
 
-        # Create the system and determine if any of our atoms have NBFIX (and
+        # Create the system and determine if any of our atoms have NBFIX or NBTHOLE (and
         # therefore requires a CustomNonbondedForce instead)
         typenames = set()
         system = mm.System()
@@ -806,6 +859,17 @@ class CharmmPsfFile(object):
                         raise StopIteration
         except StopIteration:
             pass
+        has_nbthole_terms = False
+        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.nbthole:
+                        has_nbthole_terms = True
+                        raise StopIteration
+        except StopIteration:
+            pass
+
         # Set up the constraints
         if verbose and (constraints is not None and not rigidWater):
             print('Adding constraints...')
@@ -832,6 +896,25 @@ class CharmmPsfFile(object):
                     bond.atom2.residue.resname in WATNAMES):
                     system.addConstraint(bond.atom1.idx, bond.atom2.idx,
                                          bond.bond_type.req*length_conv)
+
+        # Add virtual sites
+        if verbose: print('Adding lonepairs...')
+        for lpsite in self.lonepair_list:
+            index=lpsite[0]
+            atom1=lpsite[1]
+            atom2=lpsite[2]
+            atom3=lpsite[3]
+            if lpsite[4] > 0 : #relative
+                r = lpsite[4] /10.0 # in nanometer
+                xweights = [-1.0, 0.0, 1.0]
+            elif lpsite[4] < 0: # bisector
+                r = lpsite[4] / (-10.0)
+                xweights = [-1.0, 0.5, 0.5]
+            theta = lpsite[5] * pi / 180.0
+            phi = (180.0 - lpsite[6]) * pi / 180.0
+            p = [r*cos(theta), r*sin(theta)*cos(phi), r*sin(theta)*sin(phi)]
+            p = [x if abs(x) > 1e-10 else 0 for x in p] # Avoid tiny numbers caused by roundoff error
+            system.setVirtualSite(index, mm.LocalCoordinatesSite(atom1, atom3, atom2, mm.Vec3(1.0, 0.0, 0.0), mm.Vec3(xweights[0],xweights[1],xweights[2]), mm.Vec3(0.0, -1.0, 1.0), mm.Vec3(p[0],p[1],p[2])))
         # Add Bond forces
         if verbose: print('Adding bonds...')
         force = mm.HarmonicBondForce()
@@ -1079,8 +1162,8 @@ class CharmmPsfFile(object):
                     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
+                    elif not atom.nbfix and not atom.nbthole:
+                        # Only non-NBFIXed and non-NBTholed atom types can be compressed
                         ljtype2 = (atom2.rmin, atom2.epsilon)
                         if ljtype == ljtype2:
                             lj_idx_list[j] = num_lj_types
@@ -1135,6 +1218,68 @@ class CharmmPsfFile(object):
             for i in lj_idx_list:
                 cforce.addParticle((i - 1,)) # adjust for indexing from 0
 
+        # Add NBTHOLE terms
+        if IsDrudePSF 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 
+            nbt_type_list = []
+            nbt_set_list = []
+            for i, atom in enumerate(self.atom_list):
+                atom = atom.type
+                if not atom.nbthole: continue # get them as zero
+                if nbt_idx_list[i]: continue # already assigned
+                num_nbt_types += 1
+                nbt_idx_list[i] = num_nbt_types
+                nbt_idx_list[i+1] = num_nbt_types
+                nbt_alpha_list[i] = pow(-1*self.drudeconsts_list[i][0],-1./6.)
+                nbt_alpha_list[i+1] = pow(-1*self.drudeconsts_list[i][0],-1./6.)
+                nbt_type_list.append(atom)
+                nbt_set_list.append([i,i+1])
+                for j in range(i+1, len(self.atom_list)):
+                    atom2 = self.atom_list[j].type
+                    if nbt_idx_list[j] > 0: continue # already assigned
+                    if atom2 is atom:
+               	       nbt_idx_list[j] = num_nbt_types
+               	       nbt_idx_list[j+1] = num_nbt_types
+                       nbt_alpha_list[j] = pow(-1*self.drudeconsts_list[j][0],-1./6.)
+                       nbt_alpha_list[j+1] = pow(-1*self.drudeconsts_list[j][0],-1./6.)
+                       nbt_set_list[num_nbt_types-1].append(j)
+                       nbt_set_list[num_nbt_types-1].append(j+1)
+            num_total_nbt=num_nbt_types+1 # use zero index for all the atoms with no nbthole
+            nbt_interset_list=[]
+            # need to get all other particles as an independent group, so in total num_nbt_types+1
+            coef = [0 for i in range(num_total_nbt*num_total_nbt)]
+            for i in range(num_nbt_types):
+                for j in range(num_nbt_types):
+                    namej = nbt_type_list[j].name
+                    nbt_value = nbt_type_list[i].nbthole.get(namej,0)
+                    if abs(nbt_value)>TINY and i<j :  nbt_interset_list.append([i+1,j+1])
+                    coef[i+1+num_total_nbt*(j+1)]=nbt_value
+            nbtforce = mm.CustomNonbondedForce('-138.935456*charge1*charge2*(1.0+0.5*screen*r)*exp(-1.0*screen*r)/r; screen=coef(type1, type2) * alpha1*alpha2*10.0')
+            nbtforce.addTabulatedFunction('coef', mm.Discrete2DFunction(num_total_nbt, num_total_nbt, coef))
+            nbtforce.addPerParticleParameter('charge')
+            nbtforce.addPerParticleParameter('alpha')
+            nbtforce.addPerParticleParameter('type')
+            nbtforce.setForceGroup(self.NONBONDED_FORCE_GROUP)
+            # go through all the particles to set up per-particle parameters
+            for i in range(system.getNumParticles()):
+                c=force.getParticleParameters(i)
+                cc=c[0]/u.elementary_charge
+                aa=nbt_alpha_list[i]
+                ti=nbt_idx_list[i]
+                nbtforce.addParticle([cc,aa,ti])
+            # set interaction group
+            for a in nbt_interset_list:
+                ai=a[0]
+                aj=a[1]
+                nbtforce.addInteractionGroup(nbt_set_list[ai-1],nbt_set_list[aj-1])
+            nbtforce.setNonbondedMethod(nbtforce.CutoffPeriodic)
+            nbtforce.setCutoffDistance(0.5*u.nanometer)
+            nbtforce.setUseSwitchingFunction(False)
+            # now, add the actual force to the system
+            system.addForce(nbtforce)
+
         # Add 1-4 interactions
         excluded_atom_pairs = set() # save these pairs so we don't zero them out
         sigma_scale = 2**(-1/6)
@@ -1160,20 +1305,99 @@ class CharmmPsfFile(object):
             )
 
         # Add excluded atoms
+        # Drude and lonepairs will be excluded based on their parent atoms
+        parent_exclude_list=[]
+        for atom in self.atom_list:
+            parent_exclude_list.append([])
+        for lpsite in self.lonepair_list:
+            idx = lpsite[1]
+            idxa = lpsite[0]
+            parent_exclude_list[idx].append(idxa)
+            force.addException(idx, idxa, 0.0, 0.1, 0.0)
+        if IsDrudePSF:
+            for pair in self.drudepair_list:
+                idx = pair[0]
+                idxa = pair[1]
+                parent_exclude_list[idx].append(idxa)
+                force.addException(idx, idxa, 0.0, 0.1, 0.0)
+            # If lonepairs and Drude particles are bonded to the same parent atom, add exception
+            for excludeterm in parent_exclude_list:
+                if(len(excludeterm) >= 2):
+                    for i in range(len(excludeterm)):
+                        for j in range(i):
+                            force.addException(excludeterm[j], excludeterm[i], 0.0, 0.1, 0.0)
+        # Exclude all bonds and angles, as well as the lonepair/Drude attached onto them
         for atom in self.atom_list:
-            # Exclude all bonds and angles
             for atom2 in atom.bond_partners:
                 if atom2.idx > atom.idx:
-                    force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
+                    for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
+                        for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
+                            force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
             for atom2 in atom.angle_partners:
                 if atom2.idx > atom.idx:
-                    force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
+                    for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
+                        for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
+                            force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
             for atom2 in atom.dihedral_partners:
                 if atom2.idx <= atom.idx: continue
                 if ((atom.idx, atom2.idx) in excluded_atom_pairs):
                     continue
                 force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
         system.addForce(force)
+
+        # Add Drude particles (Drude force)
+        if IsDrudePSF:
+            if verbose: print('Adding Drude force and Thole screening...')
+            drudeforce = mm.DrudeForce()
+            drudeforce.setForceGroup(7)
+            for pair in self.drudepair_list:
+                parentatom=pair[0]
+                drudeatom=pair[1]
+                p=[-1, -1, -1]
+                a11 = 0 
+                a22 = 0
+                charge = self.atom_list[drudeatom].charge
+                polarizability = self.drudeconsts_list[parentatom][0]/(-1000.0)
+                for aniso in self.aniso_list: # not smart to do another loop, but the number of aniso is small anyway
+                    if (parentatom==aniso[0]):
+                        p[0]=aniso[1]
+                        p[1]=aniso[2]
+                        p[2]=aniso[3]
+                        k11=aniso[4]
+                        k22=aniso[5]
+                        k33=aniso[6]
+                        # solve out DrudeK, which should equal 500.0
+                        a = k11+k22+3*k33
+                        b = 2*k11*k22+4*k11*k33+4*k22*k33+6*k33*k33
+                        c = 3*k33*(k11+k33)*(k22+k33)
+                        DrudeK = (sqrt(b*b-4*a*c)-b)/2/a
+                        a11=round(DrudeK/(k11+k33+DrudeK),5)
+                        a22=round(DrudeK/(k22+k33+DrudeK),5)
+                drudeforce.addParticle(drudeatom, parentatom, p[0], p[1], p[2], charge, polarizability, a11, a22 )
+            system.addForce(drudeforce)
+            particleMap = {}
+            for i in range(drudeforce.getNumParticles()):
+                particleMap[drudeforce.getParticleParameters(i)[0]] = i
+            
+            for bond in self.bond_list:
+                alpha1 = self.drudeconsts_list[bond.atom1.idx][0]
+                alpha2 = self.drudeconsts_list[bond.atom2.idx][0] 
+                if abs(alpha1) > TINY and abs(alpha2) > TINY: # both are Drude parent atoms
+                    thole1 = self.drudeconsts_list[bond.atom1.idx][1]
+                    thole2 = self.drudeconsts_list[bond.atom2.idx][1]
+                    drude1 = bond.atom1.idx + 1 # CHARMM psf has hard-coded rule that the Drude is next to its parent
+                    drude2 = bond.atom2.idx + 1
+                    drudeforce.addScreenedPair(particleMap[drude1], particleMap[drude2], thole1+thole2)
+            for ang in self.angle_list:
+                alpha1 = self.drudeconsts_list[ang.atom1.idx][0]
+                alpha2 = self.drudeconsts_list[ang.atom3.idx][0] 
+                if abs(alpha1) > TINY and abs(alpha2) > TINY: # both are Drude parent atoms
+                    thole1 = self.drudeconsts_list[ang.atom1.idx][1]
+                    thole2 = self.drudeconsts_list[ang.atom3.idx][1]
+                    drude1 = ang.atom1.idx + 1 # CHARMM psf has hard-coded rule that the Drude is next to its parent
+                    drude2 = ang.atom3.idx + 1
+                    drudeforce.addScreenedPair(particleMap[drude1], particleMap[drude2], thole1+thole2)
+
         # If we needed a CustomNonbondedForce, map all of the exceptions from
         # the NonbondedForce to the CustomNonbondedForce
         if has_nbfix_terms:

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

@@ -100,6 +100,9 @@ class AtomType(object):
     nbfix : dict
         Dictionary that maps nbfix terms with other atom types. Dict entries
         are (rmin, epsilon) -- precombined values for that particular atom pair
+    nbthole : dict
+            Dictionary that maps nbthole terms with other atom types. Dict entries
+            are the value of a -- screening factor for that pair
 
     Examples
     --------
@@ -137,6 +140,8 @@ class AtomType(object):
         # 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()
+        # Likewise, store each NBTHOLE term as a dict 
+        self.nbthole = dict()
 
     def __eq__(self, other):
         """
@@ -173,6 +178,10 @@ class AtomType(object):
         if epsilon14 is None: epsilon14 = epsilon
         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)
+
     def __str__(self):
         return self.name