fix Charmm NBThole calculation bug and Improper dihedral parameter match issue (#2952)

* nbthole

* Update charmmpsffile.py

* Update charmmpsffile.py

* Update charmmparameterset.py

fix multiple dihedral with the same periodicity and add a warning reminder for improper dihedral

* add test example for NBTHole

add system files of CYT-GUA-CYT DNA segment for testing NBThole

* Update TestCharmmFiles.py

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

@@ -387,7 +387,7 @@ class CharmmParameterSet(object):
                             if dtype != dihedral:
                                 warnings.warn('Replacing dihedral %r with %r' %
                                               (dtype, dihedral))
-                            self.dihedral_types[key]
+                                self.dihedral_types[key][i]=dihedral
                             replaced = True
                             break
                     if not replaced:
@@ -413,10 +413,20 @@ class CharmmParameterSet(object):
                     theteq = tmp
                 except IndexError:
                     pass # Do nothing
-                # Improper types seem not to have the central atom defined in
-                # the first place, so just have the key a fully sorted list. We
-                # still depend on the PSF having properly ordered improper atoms
-                key = tuple(sorted([type1, type2, type3, type4]))
+                if type1 < type4:
+                    key = (type1, type2, type3, type4)
+                elif type1 > type4:
+                    key = (type4, type3, type2, type1)
+                else:
+                    # OK, we need to sort by the middle atoms now
+                    if type2 < type3:
+                        key = (type1, type2, type3, type4)
+                    else:
+                        key = (type4, type3, type2, type1)
+                # check repeat improper dihedral
+                if key in self.improper_types:
+                    warnings.warn('Repeat improper dihedral found %r and New param. k: %f  ,thetaq: %f will use' %
+                                  (key,k,theteq))                        
                 self.improper_types[key] = ImproperType(k, theteq)
                 continue
             if section == 'CMAP':

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

@@ -677,13 +677,12 @@ class CharmmPsfFile(object):
             # Store the atoms
             a1, a2, a3, a4 = imp.atom1, imp.atom2, imp.atom3, imp.atom4
             at1, at2, at3, at4 = a1.attype, a2.attype, a3.attype, a4.attype
-            key = tuple(sorted([at1, at2, at3, at4]))
+            key = min((at1,at2,at3,at4), (at4,at3,at2,at1))
             if not key in parmset.improper_types:
-                # Check for wild-cards
-                for anchor in (at2, at3, at4):
-                    key = tuple(sorted([at1, anchor, 'X', 'X']))
-                    if key in parmset.improper_types:
-                        break # This is the right key
+                key = min((at1,'X', 'X',at4),(at4,'X','X',at1))
+                if not key in parmset.improper_types:
+                    raise MissingParameter('No improper dihedral parameters found for '
+                                           '%r' % imp)
             try:
                 imp.improper_type = parmset.improper_types[key]
             except KeyError:
@@ -1490,6 +1489,31 @@ class CharmmPsfFile(object):
                 ii, jj, q, eps, sig = force.getExceptionParameters(i)
                 nbtforce.addExclusion(ii, jj)
 
+        # Using CustomBondForce to Calculate 1-4 atom pairs' NBThole interaction which have been excluded 
+        if has_drude_particle and has_nbthole_terms:
+            nbt14force=mm.CustomBondForce('-138.935456*charge_prod*(1.0+0.5*screen*r)*exp(-1.0*screen*r)/r')
+            nbt14force.addPerBondParameter("charge_prod")
+            nbt14force.addPerBondParameter("screen")
+            nbt14force.setForceGroup(self.NONBONDED_FORCE_GROUP)
+            num_nbt14=0
+            for dih in self.dihedral_list:
+                a1, a4 = dih.atom1, dih.atom4
+                idx_a1, idx_a4 = a1.idx, a4.idx
+                at1, at4 = self.atom_list[idx_a1].type, self.atom_list[idx_a4].type
+                if at1.nbthole and at4.nbthole:
+                    name_a4 = at4.name
+                    nbt_value = at1.nbthole.get(name_a4,0)
+                    if abs(nbt_value)>TINY:
+                        q1, q4 = a1.charge, a4.charge
+                        alpha1= pow(-1*self.drudeconsts_list[idx_a1][0],-1./6.)
+                        alpha4= pow(-1*self.drudeconsts_list[idx_a4][0],-1./6.)
+                        qij=q1*q4
+                        screen=nbt_value*alpha1*alpha4*10.0
+                        nbt14force.addBond(idx_a1,idx_a4,[qij,screen])
+                        num_nbt14+=1
+            if num_nbt14>0:
+                system.addForce(nbt14force)
+
         # Add GB model if we're doing one
         if implicitSolvent is not None:
             if verbose: print('Adding GB parameters...')

wrappers/python/tests/TestCharmmFiles.py

@@ -147,6 +147,26 @@ class TestCharmmFiles(unittest.TestCase):
         ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
         self.assertAlmostEqual(ene, 15559.71602, delta=0.05)
 
+    def test_NBThole(self):
+        """Tests CHARMM system with NBTHole"""
+        warnings.filterwarnings('ignore', category=CharmmPSFWarning)
+        psf = CharmmPsfFile('systems/cyt-gua-cyt.psf')
+        crd = CharmmCrdFile('systems/cyt-gua-cyt.crd')
+        params = CharmmParameterSet('systems/toppar_drude_master_protein_2013e.str','systems/toppar_drude_nucleic_acid_2017b.str')
+        # Box dimensions (cubic box)
+        psf.setBox(30.0*angstroms, 30.0*angstroms, 30.0*angstroms)
+
+        # Now compute the full energy
+        plat = Platform.getPlatformByName('Reference')
+        system = psf.createSystem(params, nonbondedMethod=PME, ewaldErrorTolerance=0.00005)
+        integrator = DrudeLangevinIntegrator(300*kelvin, 1.0/picosecond, 1*kelvin, 10/picosecond, 0.001*picoseconds)
+        con = Context(system, integrator, plat)
+        con.setPositions(crd.positions)
+
+        state = con.getState(getEnergy=True, enforcePeriodicBox=True)
+        ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+        self.assertAlmostEqual(ene, -292.73015, delta=1.0)
+        
     def test_Drude(self):
         """Test CHARMM systems with Drude force field"""
         warnings.filterwarnings('ignore', category=CharmmPSFWarning)