Skip to content

GitLab

Unverified Commit e936acac authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2497 from z-gong/psf 

Fixed a potential bug in charmmpsffile parser when building 1-3 and 1-4 exclusions
parents 2ca749a1 71624456
Expand all Show whitespace changes
Inline Side-by-side
Showing with 562 additions and 59 deletions
wrappers/python/simtk/openmm/app/charmmpsffile.py
View file @ e936acac
......@@ -164,6 +164,7 @@ class CharmmPsfFile(object):
CMAP_FORCE_GROUP = 5
NONBONDED_FORCE_GROUP = 6
GB_FORCE_GROUP = 6
DRUDE_FORCE_GROUP = 7
@_catchindexerror
def __init__(self, psf_name, periodicBoxVectors=None, unitCellDimensions=None):
......@@ -466,6 +467,37 @@ class CharmmPsfFile(object):
else:
self.box_vectors = periodicBoxVectors
def _build_exclusion_list(self):
pair_12_set = set()
pair_13_set = set()
pair_14_set = set()
for bond in self.bond_list:
a1, a2 = bond.atom1, bond.atom2
pair = (min(a1.idx, a2.idx), max(a1.idx, a2.idx),)
pair_12_set.add(pair)
for bond in self.bond_list:
a2, a3 = bond.atom1, bond.atom2
for a1 in a2.bond_partners:
pair = (min(a1.idx, a3.idx), max(a1.idx, a3.idx),)
if a1 != a3:
pair_13_set.add(pair)
for a4 in a3.bond_partners:
pair = (min(a2.idx, a4.idx), max(a2.idx, a4.idx),)
if a2 != a4:
pair_13_set.add(pair)
for bond in self.bond_list:
a2, a3 = bond.atom1, bond.atom2
for a1 in a2.bond_partners:
for a4 in a3.bond_partners:
pair = (min(a1.idx, a4.idx), max(a1.idx, a4.idx),)
if a1 != a4:
pair_14_set.add(pair)
# in case there are 3,4,5-member rings
self.pair_12_list = list(sorted(pair_12_set))
self.pair_13_list = list(sorted(pair_13_set - pair_12_set))
self.pair_14_list = list(sorted(pair_14_set - pair_13_set.union(pair_12_set)))
@staticmethod
def _convert(string, type, message):
"""Converts a string to a specific type, making sure to raise
......@@ -1320,31 +1352,29 @@ class CharmmPsfFile(object):
# now, add the actual force to the system
system.addForce(nbtforce)
# build 1-2, 1-3 and 1-4 pairs from connectivity
if verbose:
print('Build exclusion list...')
self._build_exclusion_list()
if verbose:
print('\tNumber of 1-2 pairs: %i' % len(self.pair_12_list))
print('\tNumber of 1-3 pairs: %i' % len(self.pair_13_list))
print('\tNumber of 1-4 pairs: %i' % len(self.pair_14_list))
# Add 1-4 interactions
excluded_atom_pairs = set() # save these pairs so we don't zero them out
sigma_scale = 2**(-1/6)
nbxmod = abs(params.nbxmod)
if nbxmod == 4:
for ia1, ia4 in self.pair_14_list:
force.addException(ia1, ia4, 0.0, 0.1, 0.0)
if nbxmod == 5:
for tor in self.dihedral_parameter_list:
# First check to see if atoms 1 and 4 are already excluded because
# they are 1-2 or 1-3 pairs (would happen in 6-member rings or
# fewer). Then check that they're not already added as exclusions
if tor.atom1 in tor.atom4.bond_partners: continue
if tor.atom1 in tor.atom4.angle_partners: continue
key = min((tor.atom1.idx, tor.atom4.idx),
(tor.atom4.idx, tor.atom1.idx))
if key in excluded_atom_pairs: continue # multiterm...
charge_prod = (tor.atom1.charge * tor.atom4.charge)
epsilon = (sqrt(abs(tor.atom1.type.epsilon_14) * ene_conv *
abs(tor.atom4.type.epsilon_14) * ene_conv))
sigma = (tor.atom1.type.rmin_14 + tor.atom4.type.rmin_14) * (
length_conv * sigma_scale)
force.addException(tor.atom1.idx, tor.atom4.idx,
charge_prod, sigma, epsilon)
excluded_atom_pairs.add(
min((tor.atom1.idx, tor.atom4.idx),
(tor.atom4.idx, tor.atom1.idx))
)
for ia1, ia4 in self.pair_14_list:
atom1 = self.atom_list[ia1]
atom4 = self.atom_list[ia4]
charge_prod = (atom1.charge * atom4.charge)
epsilon = sqrt(abs(atom1.type.epsilon_14 * atom4.type.epsilon_14)) * ene_conv
sigma = (atom1.type.rmin_14 + atom4.type.rmin_14) * (length_conv * sigma_scale)
force.addException(ia1, ia4, charge_prod, sigma, epsilon)
# Add excluded atoms
# Drude and lonepairs will be excluded based on their parent atoms
......@@ -1368,33 +1398,24 @@ class CharmmPsfFile(object):
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 1-2 and 1-3 pairs as well as the lonepair/Drude attached onto them
if nbxmod > 1:
for atom2 in atom.bond_partners:
if atom2.idx > atom.idx:
for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
for ia1, ia2 in self.pair_12_list:
for excludeatom in [ia1]+parent_exclude_list[ia1]:
for excludeatom2 in [ia2]+parent_exclude_list[ia2]:
force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
if nbxmod > 2:
for atom2 in atom.angle_partners:
if atom2.idx > atom.idx:
for excludeatom in [atom.idx]+parent_exclude_list[atom.idx]:
for excludeatom2 in [atom2.idx]+parent_exclude_list[atom2.idx]:
for ia1, ia3 in self.pair_13_list:
for excludeatom in [ia1]+parent_exclude_list[ia1]:
for excludeatom2 in [ia3]+parent_exclude_list[ia3]:
force.addException(excludeatom, excludeatom2, 0.0, 0.1, 0.0)
if nbxmod > 3:
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 has_drude_particle:
if verbose: print('Adding Drude force and Thole screening...')
drudeforce = mm.DrudeForce()
drudeforce.setForceGroup(7)
drudeforce.setForceGroup(self.DRUDE_FORCE_GROUP)
for pair in self.drudepair_list:
parentatom=pair[0]
drudeatom=pair[1]
......@@ -1424,23 +1445,15 @@ class CharmmPsfFile(object):
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]
# Apply thole screening for 1-2 and 1-3 pairs
for ia1, ia2 in self.pair_12_list + self.pair_13_list:
alpha1 = self.drudeconsts_list[ia1][0]
alpha2 = self.drudeconsts_list[ia2][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
thole1 = self.drudeconsts_list[ia1][1]
thole2 = self.drudeconsts_list[ia2][1]
drude1 = ia1 + 1 # CHARMM psf has hard-coded rule that the Drude is next to its parent
drude2 = ia2 + 1
drudeforce.addScreenedPair(particleMap[drude1], particleMap[drude2], thole1+thole2)
# If we needed a CustomNonbondedForce, map all of the exceptions from
......
wrappers/python/tests/TestCharmmFiles.py
View file @ e936acac
......@@ -349,6 +349,18 @@ class TestCharmmFiles(unittest.TestCase):
energy = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilocalories_per_mole)
self.assertAlmostEqual(energy, modeEnergy[abs(nbxmod)], delta=1e-3*abs(energy))
def test_Nonbonded_Exclusion(self):
"""Test that the 1-2, 1-3 and 1-4 pairs are correctly excluded or scaled."""
psf = CharmmPsfFile('systems/MoS2.psf')
pdb = PDBFile('systems/MoS2.pdb')
params = CharmmParameterSet('systems/MoS2.prm')
system = psf.createSystem(params, nonbondedMethod=NoCutoff)
context = Context(system, VerletIntegrator(1*femtoseconds), Platform.getPlatformByName('Reference'))
context.setPositions(pdb.positions)
energy = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilocalories_per_mole)
# Compare with value computed with NAMD.
self.assertAlmostEqual(energy, -2154.5539, delta=1e-3*abs(energy))
if __name__ == '__main__':
unittest.main()
......
wrappers/python/tests/systems/MoS2-namd.inp 0 → 100644
View file @ e936acac
timestep 1.0 # fs
# initial config
coordinates MoS2.pdb
temperature 300K
outputname out
# force field params
structure MoS2.psf
parameters MoS2.prm
paraTypeCharmm on
vdwGeometricSigma no
LJcorrection no
exclude scaled1-4
1-4scaling 1.0
switching off
cutoff 100 # A. actually no cutoff
pairlistdist 100
run 20
wrappers/python/tests/systems/MoS2.pdb 0 → 100644
View file @ e936acac
TITLE created by fftool
REMARK SIMBOX
CRYST1 50.000 50.000 50.000 90.00 90.00 90.00 P 1 1
ATOM 1 Mo01 MoS2 1 36.632 11.332 15.080 1.00 0.00 Mo
ATOM 2 Mo02 MoS2 1 48.105 10.837 4.221 1.00 0.00 Mo
ATOM 3 Mo03 MoS2 1 48.501 1.004 16.589 1.00 0.00 Mo
ATOM 4 Mo04 MoS2 1 39.006 9.266 15.382 1.00 0.00 Mo
ATOM 5 Mo05 MoS2 1 45.811 10.936 6.393 1.00 0.00 Mo
ATOM 6 Mo06 MoS2 1 48.422 2.971 14.116 1.00 0.00 Mo
ATOM 7 Mo07 MoS2 1 48.184 8.870 6.695 1.00 0.00 Mo
ATOM 8 Mo08 MoS2 1 38.927 11.233 12.909 1.00 0.00 Mo
ATOM 9 Mo09 MoS2 1 46.127 3.070 16.287 1.00 0.00 Mo
ATOM 10 Mo10 MoS2 1 41.380 7.201 15.684 1.00 0.00 Mo
ATOM 11 Mo11 MoS2 1 43.516 11.035 8.565 1.00 0.00 Mo
ATOM 12 Mo12 MoS2 1 48.343 4.937 11.642 1.00 0.00 Mo
ATOM 13 Mo13 MoS2 1 48.264 6.904 9.168 1.00 0.00 Mo
ATOM 14 Mo14 MoS2 1 41.221 11.134 10.737 1.00 0.00 Mo
ATOM 15 Mo15 MoS2 1 43.754 5.135 15.986 1.00 0.00 Mo
ATOM 16 Mo16 MoS2 1 36.553 13.298 12.607 1.00 0.00 Mo
ATOM 17 Mo17 MoS2 1 36.237 21.164 2.712 1.00 0.00 Mo
ATOM 18 Mo18 MoS2 1 45.732 12.902 3.919 1.00 0.00 Mo
ATOM 19 Mo19 MoS2 1 41.301 9.167 13.210 1.00 0.00 Mo
ATOM 20 Mo20 MoS2 1 45.890 8.969 8.867 1.00 0.00 Mo
ATOM 21 Mo21 MoS2 1 46.048 5.036 13.814 1.00 0.00 Mo
ATOM 22 Mo22 MoS2 1 43.674 7.102 13.512 1.00 0.00 Mo
ATOM 23 Mo23 MoS2 1 43.595 9.068 11.038 1.00 0.00 Mo
ATOM 24 Mo24 MoS2 1 45.969 7.003 11.340 1.00 0.00 Mo
ATOM 25 Mo25 MoS2 1 36.474 15.265 10.133 1.00 0.00 Mo
ATOM 26 Mo26 MoS2 1 38.610 19.099 3.014 1.00 0.00 Mo
ATOM 27 Mo27 MoS2 1 43.437 13.001 6.091 1.00 0.00 Mo
ATOM 28 Mo28 MoS2 1 43.358 14.968 3.618 1.00 0.00 Mo
ATOM 29 Mo29 MoS2 1 36.316 19.198 5.186 1.00 0.00 Mo
ATOM 30 Mo30 MoS2 1 38.848 13.199 10.435 1.00 0.00 Mo
ATOM 31 Mo31 MoS2 1 36.395 17.231 7.660 1.00 0.00 Mo
ATOM 32 Mo32 MoS2 1 40.984 17.033 3.316 1.00 0.00 Mo
ATOM 33 Mo33 MoS2 1 41.142 13.100 8.263 1.00 0.00 Mo
ATOM 34 Mo34 MoS2 1 38.769 15.166 7.961 1.00 0.00 Mo
ATOM 35 Mo35 MoS2 1 38.689 17.132 5.488 1.00 0.00 Mo
ATOM 36 Mo36 MoS2 1 41.063 15.067 5.790 1.00 0.00 Mo
ATOM 37 S01 MoS2 1 48.074 13.199 4.080 1.00 0.00 S
ATOM 38 S02 MoS2 1 36.601 13.694 14.939 1.00 0.00 S
ATOM 39 S03 MoS2 1 36.205 23.526 2.571 1.00 0.00 S
ATOM 40 S04 MoS2 1 46.502 11.163 2.512 1.00 0.00 S
ATOM 41 S05 MoS2 1 35.029 11.658 13.371 1.00 0.00 S
ATOM 42 S06 MoS2 1 34.633 21.490 1.003 1.00 0.00 S
ATOM 43 S07 MoS2 1 45.700 15.264 3.778 1.00 0.00 S
ATOM 44 S08 MoS2 1 38.895 13.595 12.767 1.00 0.00 S
ATOM 45 S09 MoS2 1 36.284 21.560 5.045 1.00 0.00 S
ATOM 46 S10 MoS2 1 44.128 13.228 2.210 1.00 0.00 S
ATOM 47 S11 MoS2 1 37.323 11.559 11.199 1.00 0.00 S
ATOM 48 S12 MoS2 1 34.712 19.524 3.477 1.00 0.00 S
ATOM 49 S13 MoS2 1 36.522 15.660 12.466 1.00 0.00 S
ATOM 50 S14 MoS2 1 45.779 13.298 6.252 1.00 0.00 S
ATOM 51 S15 MoS2 1 38.579 21.461 2.873 1.00 0.00 S
ATOM 52 S16 MoS2 1 34.950 13.624 10.898 1.00 0.00 S
ATOM 53 S17 MoS2 1 44.207 11.262 4.684 1.00 0.00 S
ATOM 54 S18 MoS2 1 37.007 19.425 1.305 1.00 0.00 S
ATOM 55 S19 MoS2 1 43.326 17.330 3.477 1.00 0.00 S
ATOM 56 S20 MoS2 1 41.190 13.496 10.596 1.00 0.00 S
ATOM 57 S21 MoS2 1 36.363 19.593 7.519 1.00 0.00 S
ATOM 58 S22 MoS2 1 41.754 15.294 1.909 1.00 0.00 S
ATOM 59 S23 MoS2 1 39.618 11.460 9.028 1.00 0.00 S
ATOM 60 S24 MoS2 1 34.791 17.557 5.951 1.00 0.00 S
ATOM 61 S25 MoS2 1 36.442 17.627 9.992 1.00 0.00 S
ATOM 62 S26 MoS2 1 43.484 13.397 8.424 1.00 0.00 S
ATOM 63 S27 MoS2 1 40.952 19.395 3.175 1.00 0.00 S
ATOM 64 S28 MoS2 1 34.870 15.591 8.424 1.00 0.00 S
ATOM 65 S29 MoS2 1 41.912 11.361 6.856 1.00 0.00 S
ATOM 66 S30 MoS2 1 39.380 17.359 1.607 1.00 0.00 S
ATOM 67 S31 MoS2 1 48.153 11.232 6.554 1.00 0.00 S
ATOM 68 S32 MoS2 1 48.469 3.366 16.448 1.00 0.00 S
ATOM 69 S33 MoS2 1 38.974 11.628 15.241 1.00 0.00 S
ATOM 70 S34 MoS2 1 46.581 9.196 4.986 1.00 0.00 S
ATOM 71 S35 MoS2 1 46.897 1.331 14.880 1.00 0.00 S
ATOM 72 S36 MoS2 1 37.402 9.593 13.673 1.00 0.00 S
ATOM 73 S37 MoS2 1 43.405 15.363 5.950 1.00 0.00 S
ATOM 74 S38 MoS2 1 38.816 15.561 10.294 1.00 0.00 S
ATOM 75 S39 MoS2 1 38.658 19.494 5.347 1.00 0.00 S
ATOM 76 S40 MoS2 1 41.833 13.327 4.382 1.00 0.00 S
ATOM 77 S41 MoS2 1 37.244 13.525 8.726 1.00 0.00 S
ATOM 78 S42 MoS2 1 37.086 17.458 3.779 1.00 0.00 S
ATOM 79 S43 MoS2 1 41.032 17.429 5.648 1.00 0.00 S
ATOM 80 S44 MoS2 1 41.111 15.462 8.122 1.00 0.00 S
ATOM 81 S45 MoS2 1 38.737 17.528 7.820 1.00 0.00 S
ATOM 82 S46 MoS2 1 39.460 15.393 4.080 1.00 0.00 S
ATOM 83 S47 MoS2 1 39.539 13.426 6.554 1.00 0.00 S
ATOM 84 S48 MoS2 1 37.165 15.492 6.252 1.00 0.00 S
ATOM 85 S49 MoS2 1 48.232 9.266 9.027 1.00 0.00 S
ATOM 86 S50 MoS2 1 46.096 5.432 16.146 1.00 0.00 S
ATOM 87 S51 MoS2 1 41.269 11.529 13.069 1.00 0.00 S
ATOM 88 S52 MoS2 1 46.660 7.230 7.459 1.00 0.00 S
ATOM 89 S53 MoS2 1 44.524 3.396 14.578 1.00 0.00 S
ATOM 90 S54 MoS2 1 39.697 9.494 11.501 1.00 0.00 S
ATOM 91 S55 MoS2 1 41.348 9.563 15.543 1.00 0.00 S
ATOM 92 S56 MoS2 1 48.390 5.333 13.974 1.00 0.00 S
ATOM 93 S57 MoS2 1 45.858 11.331 8.726 1.00 0.00 S
ATOM 94 S58 MoS2 1 39.776 7.527 13.975 1.00 0.00 S
ATOM 95 S59 MoS2 1 46.818 3.297 12.406 1.00 0.00 S
ATOM 96 S60 MoS2 1 44.286 9.296 7.158 1.00 0.00 S
ATOM 97 S61 MoS2 1 48.311 7.299 11.501 1.00 0.00 S
ATOM 98 S62 MoS2 1 43.722 7.497 15.845 1.00 0.00 S
ATOM 99 S63 MoS2 1 43.564 11.430 10.897 1.00 0.00 S
ATOM 100 S64 MoS2 1 46.739 5.264 9.933 1.00 0.00 S
ATOM 101 S65 MoS2 1 42.150 5.462 14.277 1.00 0.00 S
ATOM 102 S66 MoS2 1 41.992 9.395 9.329 1.00 0.00 S
ATOM 103 S67 MoS2 1 45.937 9.365 11.199 1.00 0.00 S
ATOM 104 S68 MoS2 1 46.017 7.398 13.673 1.00 0.00 S
ATOM 105 S69 MoS2 1 43.643 9.464 13.371 1.00 0.00 S
ATOM 106 S70 MoS2 1 44.365 7.329 9.631 1.00 0.00 S
ATOM 107 S71 MoS2 1 44.445 5.363 12.105 1.00 0.00 S
ATOM 108 S72 MoS2 1 42.071 7.428 11.803 1.00 0.00 S
END
wrappers/python/tests/systems/MoS2.prm 0 → 100644
View file @ e936acac
* CHARMM FORCE FIELD GENERATED BY fftool
*
ATOMS
MASS 1 MoS 95.9370
MASS 2 SMo 32.0640
BONDS
MoS SMo 51.422084 2.410000
ANGLES
SMo MoS SMo 141.933556 83.800000
MoS SMo MoS 244.980880 83.800000
DIHEDRALS
NONBONDED
!VLJ = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
MoS 0.000000 -0.115918 2.486253 0.000000 -0.057959 2.486253
SMo 0.000000 -0.498327 1.874512 0.000000 -0.249163 1.874512
END
wrappers/python/tests/systems/MoS2.psf 0 → 100644
View file @ e936acac
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment