Skip to content

GitLab

Commit 2a2fae31 authored by peastman's avatar peastman
Browse files

Merge pull request #549 from swails/charmm_nbfix 

Add NBFIX support to CHARMM parsers
parents 213c77eb 920ddf97
Expand all Show whitespace changes
Inline Side-by-side
Showing with 12187 additions and 33 deletions
wrappers/python/CMakeLists.txt
View file @ 2a2fae31
...@@ -55,13 +55,16 @@ foreach(SUBDIR ${SUBDIRS}) ...@@ -55,13 +55,16 @@ foreach(SUBDIR ${SUBDIRS})
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prm"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.crd"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.parm7" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.parm7"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.rst7" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.rst7"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ncrst" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ncrst"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.dms" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.dms"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.top" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.top"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.par" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.par"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.str"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*psf" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*psf"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/charmm22.*" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/charmm22.*"
) )
......
wrappers/python/simtk/openmm/app/__init__.py
View file @ 2a2fae31
...@@ -28,7 +28,7 @@ from desmonddmsfile import DesmondDMSFile ...@@ -28,7 +28,7 @@ from desmonddmsfile import DesmondDMSFile
from checkpointreporter import CheckpointReporter from checkpointreporter import CheckpointReporter
from charmmcrdfiles import CharmmCrdFile, CharmmRstFile from charmmcrdfiles import CharmmCrdFile, CharmmRstFile
from charmmparameterset import CharmmParameterSet from charmmparameterset import CharmmParameterSet
from charmmpsffile import CharmmPsfFile from charmmpsffile import CharmmPsfFile, CharmmPSFWarning
# Enumerated values # Enumerated values
......
wrappers/python/simtk/openmm/app/charmmparameterset.py
View file @ 2a2fae31
...@@ -13,7 +13,7 @@ Copyright (c) 2014 the Authors ...@@ -13,7 +13,7 @@ Copyright (c) 2014 the Authors
Author: Jason M. Swails Author: Jason M. Swails
Contributors: Contributors:
Date: April 18, 2014 Date: July 17, 2014
Permission is hereby granted, free of charge, to any person obtaining a Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"), copy of this software and associated documentation files (the "Software"),
...@@ -123,6 +123,8 @@ class CharmmParameterSet(object): ...@@ -123,6 +123,8 @@ class CharmmParameterSet(object):
elif arg.endswith('.str'): elif arg.endswith('.str'):
strs.append(arg) strs.append(arg)
elif arg.endswith('.inp'): 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] fname = os.path.split(arg)[1]
if 'par' in fname: if 'par' in fname:
pars.append(arg) pars.append(arg)
...@@ -436,11 +438,29 @@ class CharmmParameterSet(object): ...@@ -436,11 +438,29 @@ class CharmmParameterSet(object):
try: try:
at1 = words[0] at1 = words[0]
at2 = words[1] 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') 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, rmin, emin,
rmin14, emin14)
self.atom_types_str[at2].add_nbfix(at1, rmin, emin,
rmin14, emin14)
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: except IndexError:
raise CharmmFileError('Could not parse NBFIX terms.') raise CharmmFileError('Could not parse NBFIX terms.')
self.nbfix_types[(min(at1, at2), max(at1, at2))] = (emin, rmin) self.nbfix_types[(min(at1, at2), max(at1, at2))] = (rmin, emin)
# Now we're done. Load the nonbonded types into the relevant AtomType # Now we're done. Load the nonbonded types into the relevant AtomType
# instances. In order for this to work, all keys in nonbonded_types # instances. In order for this to work, all keys in nonbonded_types
# must be in the self.atom_types_str dict. Raise a RuntimeError if this # must be in the self.atom_types_str dict. Raise a RuntimeError if this
......
wrappers/python/simtk/openmm/app/charmmpsffile.py
View file @ 2a2fae31
...@@ -995,11 +995,25 @@ class CharmmPsfFile(object): ...@@ -995,11 +995,25 @@ class CharmmPsfFile(object):
u.kilojoule_per_mole) u.kilojoule_per_mole)
ene_conv = dihe_frc_conv 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() system = mm.System()
if verbose: print('Adding particles...') if verbose: print('Adding particles...')
for atom in self.atom_list: for atom in self.atom_list:
typenames.add(atom.type.name)
system.addParticle(atom.mass) 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 # Set up the constraints
if verbose and (constraints is not None and not rigidWater): if verbose and (constraints is not None and not rigidWater):
print('Adding constraints...') print('Adding constraints...')
...@@ -1240,9 +1254,85 @@ class CharmmPsfFile(object): ...@@ -1240,9 +1254,85 @@ class CharmmPsfFile(object):
# Add per-particle nonbonded parameters (LJ params) # Add per-particle nonbonded parameters (LJ params)
sigma_scale = 2**(-1/6) * 2 sigma_scale = 2**(-1/6) * 2
for i, atm in enumerate(self.atom_list): if not has_nbfix_terms:
for atm in self.atom_list:
force.addParticle(atm.charge, sigma_scale*atm.type.rmin*length_conv, force.addParticle(atm.charge, sigma_scale*atm.type.rmin*length_conv,
abs(atm.type.epsilon*ene_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:
rij, wdij, rij14, wdij14 = 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.setUseLongRangeCorrection(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 # Add 1-4 interactions
excluded_atom_pairs = set() # save these pairs so we don't zero them out excluded_atom_pairs = set() # save these pairs so we don't zero them out
...@@ -1283,6 +1373,13 @@ class CharmmPsfFile(object): ...@@ -1283,6 +1373,13 @@ class CharmmPsfFile(object):
continue continue
force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0) force.addException(atom.idx, atom2.idx, 0.0, 0.1, 0.0)
system.addForce(force) 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 # Add GB model if we're doing one
if implicitSolvent is not None: if implicitSolvent is not None:
......
wrappers/python/simtk/openmm/app/internal/charmm/_charmmfile.py
View file @ 2a2fae31
...@@ -47,7 +47,6 @@ class CharmmFile(object): ...@@ -47,7 +47,6 @@ class CharmmFile(object):
""" """
def __init__(self, fname, mode='r'): def __init__(self, fname, mode='r'):
self.closed = False
if mode not in ('r', 'w'): if mode not in ('r', 'w'):
raise ValueError('Cannot open CharmmFile with mode "%s"' % mode) raise ValueError('Cannot open CharmmFile with mode "%s"' % mode)
if mode == 'r': if mode == 'r':
...@@ -58,6 +57,7 @@ class CharmmFile(object): ...@@ -58,6 +57,7 @@ class CharmmFile(object):
self._handle = open(fname, mode) self._handle = open(fname, mode)
except IOError, e: except IOError, e:
raise CharmmFileError(str(e)) raise CharmmFileError(str(e))
self.closed = False
self.line_number = 0 self.line_number = 0
def write(self, *args, **kwargs): def write(self, *args, **kwargs):
......
wrappers/python/simtk/openmm/app/internal/charmm/topologyobjects.py
View file @ 2a2fae31
...@@ -42,6 +42,32 @@ TINY = 1e-8 ...@@ -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): class AtomType(object):
""" """
Atom types can either be compared by indexes or names. Can be assigned with Atom types can either be compared by indexes or names. Can be assigned with
...@@ -61,6 +87,9 @@ class AtomType(object): ...@@ -61,6 +87,9 @@ class AtomType(object):
- _member_number (int, private) : The order in which this atom type - _member_number (int, private) : The order in which this atom type
was 'added' this is used to make sure that atom types added was 'added' this is used to make sure that atom types added
last have priority in assignment in the generated hash tables 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: Example:
>>> at = AtomType('HA', 1, 1.008, 1) >>> at = AtomType('HA', 1, 1.008, 1)
>>> at.name, at.number >>> at.name, at.number
...@@ -93,6 +122,9 @@ class AtomType(object): ...@@ -93,6 +122,9 @@ class AtomType(object):
self.atomic_number = atomic_number self.atomic_number = atomic_number
# We have no LJ parameters as of yet # We have no LJ parameters as of yet
self.epsilon = self.rmin = self.epsilon_14 = self.rmin_14 = None 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): def __eq__(self, other):
""" """
...@@ -123,6 +155,12 @@ class AtomType(object): ...@@ -123,6 +155,12 @@ class AtomType(object):
""" The integer representation of an AtomType is its index """ """ The integer representation of an AtomType is its index """
return self.number 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): def __str__(self):
return self.name return self.name
...@@ -311,7 +349,7 @@ class Atom(object): ...@@ -311,7 +349,7 @@ class Atom(object):
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
class AtomList(list): class AtomList(TrackedList):
""" A list of Atom instances. """ """ A list of Atom instances. """
def unmark(self): def unmark(self):
...@@ -424,7 +462,7 @@ class ResidueList(list): ...@@ -424,7 +462,7 @@ class ResidueList(list):
elif (self._last_residue != (resname, resnum) or elif (self._last_residue != (resname, resnum) or
system != self._last_residue.system): system != self._last_residue.system):
if (self._last_residue.idx == resnum and if (self._last_residue.idx == resnum and
system == self._last_residue.system): self._last_residue.system == system):
lresname = self._last_residue.resname lresname = self._last_residue.resname
warnings.warn('Residue %d split into separate residues %s ' warnings.warn('Residue %d split into separate residues %s '
'and %s' % (resnum, lresname, resname), 'and %s' % (resnum, lresname, resname),
...@@ -1047,28 +1085,6 @@ class _CmapGrid(object): ...@@ -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__': if __name__ == '__main__':
import doctest import doctest
doctest.testmod() doctest.testmod()
wrappers/python/tests/TestCharmmFiles.py
View file @ 2a2fae31
...@@ -100,6 +100,33 @@ class TestCharmmFiles(unittest.TestCase): ...@@ -100,6 +100,33 @@ class TestCharmmFiles(unittest.TestCase):
totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu) totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
self.assertAlmostEqual(totalMass1, totalMass2) self.assertAlmostEqual(totalMass1, totalMass2)
def test_NBFIX(self):
"""Tests CHARMM systems with NBFIX Lennard-Jones modifications"""
import warnings
warnings.filterwarnings('ignore', category=CharmmPSFWarning)
psf = CharmmPsfFile('systems/ala3_solv.psf')
crd = CharmmCrdFile('systems/ala3_solv.crd')
params = CharmmParameterSet('systems/par_all36_prot.prm',
'systems/toppar_water_ions.str')
# Box dimensions (found from bounding box)
psf.setBox(32.7119500*angstroms, 32.9959600*angstroms, 33.0071500*angstroms)
# Turn off charges so we only test the Lennard-Jones energies
for a in psf.atom_list:
a.charge = 0.0
# Now compute the full energy
plat = Platform.getPlatformByName('Reference')
system = psf.createSystem(params, nonbondedMethod=PME,
nonbondedCutoff=8*angstroms)
con = Context(system, VerletIntegrator(2*femtoseconds), plat)
con.setPositions(crd.positions)
state = con.getState(getEnergy=True, enforcePeriodicBox=True)
ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
self.assertAlmostEqual(ene, 15490.0033559, delta=0.05)
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
wrappers/python/tests/systems/ala3_solv.crd 0 → 100644
View file @ 2a2fae31
This source diff could not be displayed because it is too large. You can view the blob instead.
wrappers/python/tests/systems/ala3_solv.psf 0 → 100644
View file @ 2a2fae31
This diff is collapsed.
wrappers/python/tests/systems/par_all36_prot.prm 0 → 100644
View file @ 2a2fae31
This diff is collapsed.
wrappers/python/tests/systems/toppar_water_ions.str 0 → 100644
View file @ 2a2fae31
* Toplogy and parameter information for water and ions.
*
!Testcase
!test_water_ions.inp
! IMPORTANT NOTE: this file contains NBFixes between carboxylates and sodium,
! which will only apply if the main files containing carboxylate atom types
! have been read in first!
!references
!
!TIP3P water model
!
!W.L. Jorgensen; J.Chandrasekhar; J.D. Madura; R.W. Impey;
!M.L. Klein; "Comparison of simple potential functions for
!simulating liquid water", J. Chem. Phys. 79 926-935 (1983).
!
!IONS
!
!Ions from Roux and coworkers
!
!Beglov, D. and Roux, B., Finite Representation of an Infinite
!Bulk System: Solvent Boundary Potential for Computer Simulations,
!Journal of Chemical Physics, 1994, 100: 9050-9063
!
!ZINC
!
!Stote, R.H. and Karplus, M. Zinc Binding in Proteins and
!Solution: A Simple but Accurate Nonbonded Representation, PROTEINS:
!Structure, Function, and Genetics 23:12-31 (1995)
!test "append" to determine if previous toppar files have been read and
!add append to "read rtf card" if true
set nat ?NATC
set app
!We're exploiting what is arguably a bug in the parser. On the left hand side,
!the quotes have priority, so NAT is correctly substituted. On the right hand
!side, the ? has priority and NATC" (sic) is not a valid substitution...
if "@NAT" ne "?NATC" if @nat ne 0 set app append
read rtf card @app
* Topology for water and ions
*
31 1
MASS 1 HT 1.00800 H ! TIPS3P WATER HYDROGEN
MASS 2 HX 1.00800 H ! hydroxide hydrogen
MASS 3 OT 15.99940 O ! TIPS3P WATER OXYGEN
MASS 4 OX 15.99940 O ! hydroxide oxygen
MASS 5 LIT 6.94100 LI ! Lithium ion
MASS 6 SOD 22.98977 NA ! Sodium Ion
MASS 7 MG 24.30500 MG ! Magnesium Ion
MASS 8 POT 39.09830 K ! Potassium Ion
MASS 9 CAL 40.08000 CA ! Calcium Ion
MASS 10 RUB 85.46780 RB ! Rubidium Ion
MASS 11 CES 132.90545 CS ! Cesium Ion
MASS 12 BAR 137.32700 BA ! Barium Ion
MASS 13 ZN 65.37000 ZN ! zinc (II) cation
MASS 14 CAD 112.41100 CD ! cadmium (II) cation
MASS 15 CLA 35.45000 CL ! Chloride Ion
default first none last none
RESI TIP3 0.000 ! tip3p water model, generate using noangle nodihedral
GROUP
ATOM OH2 OT -0.834
ATOM H1 HT 0.417
ATOM H2 HT 0.417
BOND OH2 H1 OH2 H2 H1 H2 ! the last bond is needed for shake
ANGLE H1 OH2 H2 ! required
DONOR H1 OH2
DONOR H2 OH2
ACCEPTOR OH2
PATCHING FIRS NONE LAST NONE
RESI TP3M 0.000 ! "mmff" water model, as an analog of tip3p
GROUP
ATOM OH2 OT -0.834 ! these charges are replaced by the mmff setup
ATOM H1 HT 0.417 ! these charges are replaced by the mmff setup
ATOM H2 HT 0.417 ! these charges are replaced by the mmff setup
BOND OH2 H1 OH2 H2 ! omits the H1-H2 bond, which is needed for shake with tip3p
ANGLE H1 OH2 H2 ! required
DONOR H1 OH2
DONOR H2 OH2
ACCEPTOR OH2
PATCHING FIRS NONE LAST NONE
RESI OH -1.00 ! hydroxide ion by adm.jr.
GROUP
ATOM O1 OX -1.32
ATOM H1 HX 0.32
BOND O1 H1
DONOR H1 O1
ACCEPTOR O1
! Ion parameters from Benoit Roux and Coworkers
! As of 8/10 new NBFIX terms required
!
RESI LIT 1.00 ! Lithium Ion
GROUP
ATOM LIT LIT 1.00
PATCHING FIRST NONE LAST NONE
RESI SOD 1.00 ! Sodium Ion
GROUP
ATOM SOD SOD 1.00
PATCHING FIRST NONE LAST NONE
RESI MG 2.00 ! Magnesium Ion
GROUP
ATOM MG MG 2.00
PATCHING FIRST NONE LAST NONE
RESI POT 1.00 ! Potassium Ion
GROUP
ATOM POT POT 1.00
PATCHING FIRST NONE LAST NONE
RESI CAL 2.00 ! Calcium Ion
GROUP
ATOM CAL CAL 2.00
PATCHING FIRST NONE LAST NONE
RESI RUB 1.00 ! Rubidium Ion
GROUP
ATOM RUB RUB 1.00
PATCHING FIRST NONE LAST NONE
RESI CES 1.00 ! Cesium Ion
GROUP
ATOM CES CES 1.00
PATCHING FIRST NONE LAST NONE
RESI BAR 2.00 ! Barium Ion
GROUP
ATOM BAR BAR 2.00
PATCHING FIRST NONE LAST NONE
RESI ZN2 2.00 ! Zinc (II) cation, Roland Stote
GROUP
ATOM ZN ZN 2.00
PATCHING FIRST NONE LAST NONE
RESI CD2 2.00 ! Cadmium (II) cation
GROUP
ATOM CD CAD 2.00
PATCHING FIRST NONE LAST NONE
RESI CLA -1.00 ! Chloride Ion
GROUP
ATOM CLA CLA -1.00
PATCHING FIRST NONE LAST NONE
END
read para card flex @app
* Parameters for water and ions
*
ATOMS
MASS 1 HT 1.00800 ! TIPS3P WATER HYDROGEN
MASS 2 HX 1.00800 ! hydroxide hydrogen
MASS 3 OT 15.99940 ! TIPS3P WATER OXYGEN
MASS 4 OX 15.99940 ! hydroxide oxygen
MASS 5 LIT 6.94100 ! Lithium ion
MASS 6 SOD 22.98977 ! Sodium Ion
MASS 7 MG 24.30500 ! Magnesium Ion
MASS 8 POT 39.09830 ! Potassium Ion
MASS 9 CAL 40.08000 ! Calcium Ion
MASS 10 RUB 85.46780 ! Rubidium Ion
MASS 11 CES 132.90545 ! Cesium Ion
MASS 12 BAR 137.32700 ! Barium Ion
MASS 13 ZN 65.37000 ! zinc (II) cation
MASS 14 CAD 112.41100 ! cadmium (II) cation
MASS 15 CLA 35.45000 ! Chloride Ion
BONDS
!
!V(bond) = Kb(b - b0)**2
!
!Kb: kcal/mole/A**2
!b0: A
!
!atom type Kb b0
!
HT HT 0.0 1.5139 ! from TIPS3P geometry (for SHAKE w/PARAM)
HT OT 450.0 0.9572 ! from TIPS3P geometry
OX HX 545.0 0.9700 ! hydroxide ion
ANGLES
!
!V(angle) = Ktheta(Theta - Theta0)**2
!
!V(Urey-Bradley) = Kub(S - S0)**2
!
!Ktheta: kcal/mole/rad**2
!Theta0: degrees
!Kub: kcal/mole/A**2 (Urey-Bradley)
!S0: A
!
!atom types Ktheta Theta0 Kub S0
!
HT OT HT 55.0 104.52 ! FROM TIPS3P GEOMETRY
DIHEDRALS
!
!V(dihedral) = Kchi(1 + cos(n(chi) - delta))
!
!Kchi: kcal/mole
!n: multiplicity
!delta: degrees
!
!atom types Kchi n delta
!
!
IMPROPER
!
!V(improper) = Kpsi(psi - psi0)**2
!
!Kpsi: kcal/mole/rad**2
!psi0: degrees
!note that the second column of numbers (0) is ignored
!
!atom types Kpsi psi0
!
NONBONDED nbxmod 5 atom cdiel fshift vatom vdistance vfswitch -
cutnb 14.0 ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac 1.0 wmin 1.5
!TIP3P LJ parameters
HT 0.0 -0.046 0.2245
OT 0.0 -0.1521 1.7682
!for hydroxide
OX 0.000000 -0.120000 1.700000 ! ALLOW POL ION
! JG 8/27/89
HX 0.000000 -0.046000 0.224500 ! ALLOW PEP POL SUL ARO ALC
! same as TIP3P hydrogen, adm jr., 7/20/89
!ions
LIT 0.0 -0.00233 1.2975 ! Lithium
! From S Noskov, target ddG(Li-Na) was 23-26.0 kcal/mol (see JPC B, Lamoureux&Roux,2006)
SOD 0.0 -0.0469 1.41075 ! new CHARMM Sodium
! ddG of -18.6 kcal/mol with K+ from S. Noskov
MG 0.0 -0.0150 1.18500 ! Magnesium
! B. Roux dA = -441.65
POT 0.0 -0.0870 1.76375 ! Potassium
! D. Beglovd and B. Roux, dA=-82.36+2.8 = -79.56 kca/mol
CAL 0.0 -0.120 1.367 ! Calcium
! S. Marchand and B. Roux, dA = -384.8 kcal/mol
RUB 0.0000 -0.15 1.90 ! Rubidium
! delta A with respect to POT is +6.0 kcal/mol in bulk water
CES 0.0 -0.1900 2.100 ! Cesium
! delta A with respect to POT is +12.0 kcal/mol
BAR 0.0 -0.150 1.890 ! Barium
! B. Roux, dA = dA[calcium] + 64.2 kcal/mol
ZN 0.000000 -0.250000 1.090000 ! Zinc
! RHS March 18, 1990
CAD 0.000000 -0.120000 1.357000 ! Cadmium
! S. Marchand and B. Roux, from delta delta G
CLA 0.0 -0.150 2.27 ! Chloride
! D. Beglovd and B. Roux, dA=-83.87+4.46 = -79.40 kcal/mol
NBFIX
! Emin Rmin
! (kcal/mol) (A)
SOD CLA -0.083875 3.731 ! From osmotic pressure calibration, J. Phys.Chem.Lett. 1:183-189
POT CLA -0.114236 4.081 ! From osmotic pressure calibration, J. Phys.Chem.Lett. 1:183-189
END
! The following section contains NBFixes for sodium interacting with
! carboxylate oxygens of various CHARMM force fields. It will generate
! level -1 warnings whenever any of these force fields have not been
! read prior to the current stream file. Since we don't want to force
! the user to always read all the force fields, we're suppressing the
! warnings. The only side effect is that you will have "most severe
! warning was at level 0" at the end of your output. Also note that
! the user is responsible for reading the current file last if they
! want the NBFixes to apply. A more elegant solution would require new
! features to be added to CHARMM.
! parallel fix, to avoid duplicated messages in the log
set para
if ?NUMNODE gt 1 set para node 0
set wrn ?WRNLEV
! Some versions of CHARMM don't seem to initialize wrnlev...
if "@WRN" eq "?WRNLEV" set wrn 5
set bom ?bomlev
WRNLEV -1 @PARA
BOMLEV -1 @PARA
read para card flex append
* NBFix between carboxylate and sodium
*
! These NBFixes will only apply if the main files have been read in first!!!
NBFIX
!new SOD NBFIX values
! Simulations of Anionic Lipid Membranes: Development of Interaction-Specific
! Ion Parameters and Validation using NMR Data.
! Venable, R.M.; Luo, Y,; Gawrisch, K.; Roux, B.; Pastor, R.W.
! J. Phys. Chem. B 2013, 117 (35), pp 10183–10192. DOI: 10.1021/jp401512z
!
SOD OCL -0.07502 3.23 ! osmotic P; carboxylate =O
SOD OBL -0.07502 3.13 ! POPC optim.; ester =O
SOD O2L -0.07502 3.16 ! POPC optim.; phosphate =O
!
SOD OC -0.07502 3.23 ! For prot carboxylate groups
SOD OC2D2 -0.07502 3.23 ! For carb carboxylate groups
SOD OG2D2 -0.07502 3.23 ! For CGenFF carboxylate groups
END
BOMLEV @bom @PARA
WRNLEV @wrn @PARA
return
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