"""
This module contains classes for parsing and processing CHARMM parameter,
topology, and stream files.
Author: Jason M. Swails
Contributors:
Date: Mar. 26, 2017
"""
from __future__ import absolute_import, division, print_function
import os
import re
import warnings
from copy import copy as _copy
from collections import OrderedDict
from itertools import combinations
from ..constants import TINY
from ..exceptions import CharmmError, ParameterWarning
from ..modeller import PatchTemplate, ResidueTemplate
from ..parameters import ParameterSet
from ..periodic_table import AtomicNum, element_by_mass
from ..topologyobjects import (AngleType, Atom, AtomType, BondType, CmapType,
DihedralType, DihedralTypeList, ImproperType,
NoUreyBradley, DrudeAtom, DrudeAnisotropy)
from ..utils.io import genopen
from ..utils.six import integer_types, iteritems, string_types
from ..utils.six.moves import zip
from ._charmmfile import CharmmFile, CharmmStreamFile
_penaltyre = re.compile(r'penalty\s*=\s*([\d\.]+)')
class _EmptyStringIterator(object):
""" Always yields an empty string """
def __iter__(self):
while True:
yield ''
def __getitem__(self, idx):
return ''
def _typeconv(name):
if isinstance(name, integer_types):
return name
if name.upper() == name:
return name.replace('*', 'STR').replace('+', 'P').replace('-', 'M')[:6]
# Lowercase letters present -- decorate the type name with LTU --
# Lower To Upper
return ('%sLTU' % name.upper()).replace('*', 'STR').replace(
'+', 'P').replace('-', 'M')[:6]
[docs]class CharmmImproperMatchingMixin(object):
""" Implements CHARMM-style improper matching """
[docs] def match_improper_type(self, a1, a2, a3, a4):
""" Matches an improper type based on atom type names """
typ = self._match_improper_with_typemap(self.improper_types, a1, a2, a3, a4)
if typ is None:
typ = self._match_improper_with_typemap(self.improper_periodic_types, a1, a2, a3, a4)
return typ
def _match_improper_with_typemap(self, typemap, a1, a2, a3, a4):
if (a1, a2, a3, a4) in typemap: return typemap[(a1, a2, a3, a4)]
if (a4, a3, a2, a1) in typemap: return typemap[(a4, a3, a2, a1)]
# Now try any of the sortings. The documented CHARMM ordering does not seem to work for
# all systems CHARMM supports :(
key = tuple(sorted([a1, a2, a3, a4]))
if self._improper_key_map.get(key, None) in typemap:
return typemap[self._improper_key_map[key]]
for exact1, exact2, exact3 in combinations((a1, a2, a3, a4), 3):
key = tuple(sorted([exact1, exact2, exact3, 'X']))
if self._improper_key_map.get(key, None) in typemap:
return typemap[self._improper_key_map[key]]
for exact1, exact2 in combinations((a1, a2, a3, a4), 2):
key = tuple(sorted([exact1, exact2, 'X', 'X']))
if self._improper_key_map.get(key, None) in typemap:
return typemap[self._improper_key_map[key]]
return None
[docs]class CharmmParameterSet(ParameterSet, CharmmImproperMatchingMixin):
"""
Stores a parameter set defined by CHARMM files. It stores the equivalent of
the information found in the MASS section of the CHARMM topology file
(TOP/RTF) and all of the information in the parameter files (PAR)
Parameters
----------
*filenames : variable length arguments of str
The list of topology, parameter, and stream files to load into the
parameter set. The following file type suffixes are recognized:
.rtf, .top -- Residue topology file
.par, .prm -- Parameter file
.str -- Stream file
.inp -- If "par" is in the file name, it is a parameter file. If
"top" is in the file name, it is a topology file.
Otherwise, ValueError is raised.
See Also
--------
:class:`parmed.parameters.ParameterSet`
"""
def __copy__(self):
other = super(CharmmParameterSet, self).__copy__()
other._declared_nbrules = self._declared_nbrules
return other
@staticmethod
def _convert(data, type, msg='', line_index=None, line=None):
"""
Converts a data type to a desired type, raising CharmmError if it
fails
"""
try:
return type(data)
except ValueError:
msg = 'Could not convert %s to %s\n' % (msg, type)
if line_index is not None:
msg += 'input line %d\n' % line_index
if line is not None:
msg += 'input line: %s\n' % line
raise CharmmError(msg)
[docs] def __init__(self, *args):
# Instantiate the list types
super(CharmmParameterSet, self).__init__()
self.parametersets = []
self._declared_nbrules = False
# Load all of the files
tops, pars, strs = [], [], []
for arg in args:
if arg.endswith('.rtf') or arg.endswith('.top'):
tops.append(arg)
elif arg.endswith('.par') or arg.endswith('.prm'):
pars.append(arg)
elif arg.endswith('.str'):
strs.append(arg)
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]
if 'par' in fname:
pars.append(arg)
elif 'top' in fname:
tops.append(arg)
else:
raise ValueError('Unrecognized file type: %s' % arg)
else:
raise ValueError('Unrecognized file type: %s' % arg)
for top in tops:
self.read_topology_file(top)
for par in pars:
self.read_parameter_file(par)
for strf in strs:
self.read_stream_file(strf)
[docs] @classmethod
def from_parameterset(cls, params, copy=False):
"""
Instantiates a CharmmParameterSet from another ParameterSet (or
subclass). The main thing this feature is responsible for is converting
lower-case atom type names into all upper-case and decorating the name
to ensure each atom type name is unique.
Parameters
----------
params : :class:`parmed.parameters.ParameterSet`
ParameterSet containing the list of parameters to be converted to a
CHARMM-compatible set
copy : bool, optional
If True, the returned parameter set is a deep copy of ``params``. If
False, the returned parameter set is a shallow copy, and the
original set may be modified if any lower-case atom type names are
present. Default is False.
Returns
-------
new_params : CharmmParameterSet
The parameter set whose atom type names are converted to all
upper-case
"""
new_params = cls()
if copy:
do_copy = lambda x: _copy(x)
else:
do_copy = lambda x: x
# Convert all parameters
id_typemap = dict()
def copy_paramtype(key, typ, dict):
if isinstance(key, string_types):
key = _typeconv(key)
elif isinstance(key, tuple):
key = tuple(_typeconv(k) for k in key)
# NoUreyBradley should never be copied
if typ is NoUreyBradley:
dict[key] = NoUreyBradley
elif id(typ) in id_typemap:
dict[key] = id_typemap[id(typ)]
else:
newtype = do_copy(typ)
id_typemap[id(typ)] = newtype
dict[key] = newtype
for key, atom_type in iteritems(params.atom_types_tuple):
atom_type.name = _typeconv(atom_type.name)
copy_paramtype(key, atom_type, new_params.atom_types_tuple)
for typename, atom_type in iteritems(params.atom_types):
atom_type.name = _typeconv(atom_type.name)
copy_paramtype(typename, atom_type, new_params.atom_types)
for idx, atom_type in iteritems(params.atom_types_int):
atom_type.name = _typeconv(atom_type.name)
copy_paramtype(idx, atom_type, new_params.atom_types_int)
for key, typ in iteritems(params.bond_types):
copy_paramtype(key, typ, new_params.bond_types)
for key, typ in iteritems(params.angle_types):
copy_paramtype(key, typ, new_params.angle_types)
for key, typ in iteritems(params.urey_bradley_types):
copy_paramtype(key, typ, new_params.urey_bradley_types)
for key, typ in iteritems(params.dihedral_types):
copy_paramtype(key, typ, new_params.dihedral_types)
for key, typ in iteritems(params.improper_periodic_types):
copy_paramtype(key, typ, new_params.improper_periodic_types)
for key, typ in iteritems(params.improper_types):
copy_paramtype(key, typ, new_params.improper_types)
for key, typ in iteritems(params.cmap_types):
assert len(key) == 8, '%d-key cmap type detected!' % len(key)
copy_paramtype(key, typ, new_params.cmap_types)
for key, typ in iteritems(params.nbfix_types):
copy_paramtype(key, typ, new_params.nbfix_types)
return new_params
[docs] @classmethod
def from_structure(cls, struct):
""" Extracts known parameters from a Structure instance
Parameters
----------
struct : :class:`parmed.structure.Structure`
The parametrized ``Structure`` instance from which to extract
parameters into a ParameterSet
Returns
-------
params : :class:`ParameterSet`
The parameter set with all parameters defined in the Structure
Notes
-----
The parameters here are copies of the ones in the Structure, so
modifying the generated ParameterSet will have no effect on ``struct``.
Furthermore, the *first* occurrence of each parameter will be used. If
future ones differ, they will be silently ignored, since this is
expected behavior in some instances (like with Gromacs topologies in the
ff99sb-ildn force field).
Dihedrals are a little trickier. They can be multi-term, which can be
represented either as a *single* entry in dihedrals with a type of
DihedralTypeList or multiple entries in dihedrals with a DihedralType
parameter type. In this case, the parameter is constructed from either
the first DihedralTypeList found or the first DihedralType of each
periodicity found if no matching DihedralTypeList is found.
"""
return cls.from_parameterset(
ParameterSet.from_structure(struct, allow_unequal_duplicates=False)
)
[docs] @classmethod
def load_set(cls, tfile=None, pfile=None, sfiles=None):
"""
Instantiates a CharmmParameterSet from a Topology file and a Parameter
file (or just a Parameter file if it has all information)
Parameters
----------
tfile : str
The name of the Topology (RTF/TOP) file to parse
pfile : str
The name of the Parameter (PAR/PRM) file to parse
sfiles : list(str)
Iterable of stream (STR) file names
Returns
-------
New CharmmParameterSet populated with parameters found in the provided
files
Notes
-----
The RTF file is read first (if provided), followed by the PAR file,
followed by the list of stream files (in the order they are provided).
Parameters in each stream file will overwrite those that came before (or
simply append to the existing set if they are different)
"""
inst = cls()
if tfile is not None:
inst.read_topology_file(tfile)
if pfile is not None:
inst.read_parameter_file(pfile)
if isinstance(sfiles, string_types):
# The API docstring requests a list, but allow for users to pass a
# string with a single filename instead
inst.read_stream_file(sfiles)
elif sfiles is not None:
for sfile in sfiles:
inst.read_stream_file(sfile)
return inst
[docs] def read_parameter_file(self, pfile, comments=None):
"""
Reads all of the parameters from a parameter file. Versions 36 and
later of the CHARMM force field files have an ATOMS section defining
all of the atom types. Older versions need to load this information
from the RTF/TOP files.
Parameters
----------
pfile : str or list of lines
Name of the CHARMM parameter file to read or list of lines to parse
as a file
comments : list of str, optional
List of comments on each of the pfile lines (if pfile is a list of
lines)
Notes
-----
The atom types must all be loaded by the end of this routine. Either
supply a PAR file with atom definitions in them or read in a RTF/TOP
file first. Failure to do so will result in a raised RuntimeError.
"""
conv = CharmmParameterSet._convert
if isinstance(pfile, str):
own_handle = True
f = CharmmFile(pfile)
else:
own_handle = False
f = pfile
if not isinstance(f, CharmmFile) and comments is None:
comments = _EmptyStringIterator()
# What section are we parsing?
section = None
# The current cmap we are building (these span multiple lines)
current_cmap = None
current_cmap2 = None
current_cmap_data = []
current_cmap_res = 0
nonbonded_types = dict() # Holder
parameterset = None
declared_geometric = False
for i, line in enumerate(f):
line = line.strip()
try:
comment = f.comment
except AttributeError:
comment = comments[i]
if not line:
# This is a blank line
continue
if parameterset is None and line.strip().startswith('*>>'):
parameterset = line.strip()[1:78]
continue
# Set section if this is a section header
if line.upper().startswith('ATOM'):
section = 'ATOMS'
continue
if line.upper().startswith('BOND'):
section = 'BONDS'
continue
if line.upper().startswith('ANGLE') or line.upper().startswith('THETA'):
section = 'ANGLES'
continue
if line.upper().startswith('DIHE') or line.upper().startswith('PHI'):
section = 'DIHEDRALS'
continue
if line.upper().startswith('IMPROPER') or line.upper().startswith('IMPHI'):
section = 'IMPROPER'
continue
if line.upper().startswith('CMAP'):
section = 'CMAP'
continue
if line.upper().startswith('NONBONDED'):
read_first_nonbonded = declared_geometric = False
section = 'NONBONDED'
# Get nonbonded keywords
words = line.split()[1:]
scee = None
for i, word in enumerate(words):
if word.upper() == 'E14FAC':
try:
scee = 1 / float(words[i+1])
except (ValueError, IndexError):
raise CharmmError('Could not parse 1-4 electrostatic scaling factor '
'from NONBONDED card')
if self._declared_nbrules:
if len(self.dihedral_types) > 0:
# We already specified it -- make sure it's the same
# as the one we specified before
_, dt0 = next(iteritems(self.dihedral_types))
diff = abs(dt0[0].scee - scee)
if diff > TINY:
raise CharmmError('Inconsistent 1-4 scalings')
else:
for key, dtl in iteritems(self.dihedral_types):
for dt in dtl:
dt.scee = scee
elif word.upper().startswith('GEOM'):
if self._declared_nbrules and self.combining_rule != 'geometric':
raise CharmmError('Cannot combine parameter files with different '
'combining rules')
self.combining_rule = 'geometric'
declared_geometric = True
continue
if line.upper().startswith('NBFIX'):
section = 'NBFIX'
continue
if line.upper().startswith('HBOND'):
section = None
continue
if line.upper().startswith('THOLE'):
section = None
continue
# It seems like files? sections? can be terminated with 'END'
if line[:3].upper() == 'END':
section = None
continue
# If we have no section, skip
if section is None: continue
# See if our comments define a penalty for this line
pens = _penaltyre.findall(comment)
if len(pens) == 1:
penalty = float(pens[0])
else:
penalty = None
# Now handle each section specifically
if section.upper() == 'ATOMS':
if not line.upper().startswith('MASS'): continue # Should this happen?
words = line.split()
if words[0].upper() == 'END':
continue
try:
idx = conv(words[1], int, 'atom type', line_index=i, line=line)
name = words[2].upper()
mass = conv(words[3], float, 'atom mass', line_index=i, line=line)
except IndexError:
raise CharmmError('Could not parse MASS section.')
# The parameter file might or might not have an element name
try:
elem = words[4].upper()
if len(elem) == 2:
elem = elem[0] + elem[1].lower()
atomic_number = AtomicNum[elem]
except (IndexError, KeyError):
# Figure it out from the mass
atomic_number = AtomicNum[element_by_mass(mass)]
atype = AtomType(name=name, number=idx, mass=mass, atomic_number=atomic_number)
self.atom_types_str[atype.name] = atype
self.atom_types_int[atype.number] = atype
self.atom_types_tuple[(atype.name, atype.number)] = atype
continue
if section.upper() == 'BONDS':
words = line.split()
if words[0].upper() == 'END':
continue
try:
type1 = words[0].upper()
type2 = words[1].upper()
k = conv(words[2], float, 'bond force constant', line_index=i, line=line)
req = conv(words[3], float, 'bond equilibrium dist', line_index=i, line=line)
except IndexError:
raise CharmmError('Could not parse bonds.')
key = (min(type1, type2), max(type1, type2))
bond_type = BondType(k, req)
if key in self.bond_types:
# See if existing bond type has a different value and replaces it with a warning
if self.bond_types[key] != bond_type:
# Replace. Warn if they are different
warnings.warn('Replacing bond %r, %r with %r' %
(key, self.bond_types[key], bond_type), ParameterWarning)
self.bond_types[(type1, type2)] = bond_type
self.bond_types[(type2, type1)] = bond_type
else: # key not present
self.bond_types[(type1, type2)] = bond_type
self.bond_types[(type2, type1)] = bond_type
bond_type.penalty = penalty
continue
if section.upper() == 'ANGLES':
words = line.split()
if words[0].upper() == 'END':
continue
try:
type1 = words[0].upper()
type2 = words[1].upper()
type3 = words[2].upper()
k = conv(words[3], float, 'angle force constant', line_index=i, line=line)
theteq = conv(words[4], float, 'angle equilibrium value', line_index=i, line=line)
except IndexError:
raise CharmmError('Could not parse angles.')
angle_type = AngleType(k, theteq)
key = (type1, type2, type3)
if key in self.angle_types:
# See if the existing angle type list has a different value
# and replaces it with a warning
if self.angle_types[key] != angle_type:
# Replace. Warn if they are different
warnings.warn('Replacing angle %r, %r with %r' %
(key, self.angle_types[key], angle_type), ParameterWarning)
self.angle_types[(type1, type2, type3)] = angle_type
self.angle_types[(type3, type2, type1)] = angle_type
else: # key not present
self.angle_types[(type1, type2, type3)] = angle_type
self.angle_types[(type3, type2, type1)] = angle_type
# See if we have a urey-bradley
try:
ubk = conv(words[5], float, 'Urey-Bradley force constant', line_index=i, line=line)
ubeq = conv(words[6], float, 'Urey-Bradley equil. value', line_index=i, line=line)
ubtype = BondType(ubk, ubeq)
ubtype.penalty = penalty
except IndexError:
ubtype = NoUreyBradley
self.urey_bradley_types[(type1, type2, type3)] = ubtype
self.urey_bradley_types[(type3, type2, type1)] = ubtype
angle_type.penalty = penalty
continue
if section.upper() == 'DIHEDRALS':
words = line.split()
if words[0].upper == 'END':
continue
try:
type1 = words[0].upper()
type2 = words[1].upper()
type3 = words[2].upper()
type4 = words[3].upper()
k = conv(words[4], float, 'dihedral force constant', line_index=i, line=line)
n = conv(words[5], float, 'dihedral periodicity', line_index=i, line=line)
phase = conv(words[6], float, 'dihedral phase', line_index=i, line=line)
except IndexError:
raise CharmmError('Could not parse dihedrals.')
key = (type1, type2, type3, type4)
# See if this is a second (or more) term of the dihedral group
# that's already present.
dihedral = DihedralType(k, n, phase)
dihedral.penalty = penalty
if key in self.dihedral_types:
# See if the existing dihedral type list has a term with
# the same periodicity -- If so, replace it
replaced = False
for i, dtype in enumerate(self.dihedral_types[key]):
if dtype.per == dihedral.per:
# Replace. Warn if they are different
if dtype != dihedral:
warnings.warn('Replacing dihedral %r with %r' % (dtype, dihedral),
ParameterWarning)
self.dihedral_types[key][i] = dihedral
replaced = True
break
if not replaced:
self.dihedral_types[key].append(dihedral)
else: # key not present
dtl = DihedralTypeList()
dtl.append(dihedral)
self.dihedral_types[(type1, type2, type3, type4)] = dtl
self.dihedral_types[(type4, type3, type2, type1)] = dtl
continue
if section.upper() == 'IMPROPER':
words = line.split()
if words[0].upper() == 'END':
continue
try:
type1 = words[0].upper()
type2 = words[1].upper()
type3 = words[2].upper()
type4 = words[3].upper()
k = conv(words[4], float, 'improper force constant', line_index=i, line=line)
theteq = conv(words[5], float, 'improper equil. value', line_index=i, line=line)
except IndexError:
raise CharmmError('Could not parse dihedrals.')
# If we have a 7th column, that is the real psi0 (and the 6th
# is the multiplicity, which will indicate this is a periodic
# improper torsion (so it needs to be added to the
# improper_periodic_types list)
try:
tmp = conv(words[6], float, 'improper equil. value', line_index=i, line=line)
except IndexError:
per = 0
else:
per = int(theteq)
theteq = tmp
# Improper types seem not to always 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 = (type1, type2, type3, type4)
self._improper_key_map[tuple(sorted(key))] = key
if per == 0:
improp = ImproperType(k, theteq)
self.improper_types[key] = improp
else:
improp = DihedralType(k, per, theteq)
self.improper_periodic_types[key] = improp
improp.improper = True
improp.penalty = penalty
continue
if section.upper() == 'CMAP':
# This is the most complicated part, since cmap parameters span
# many lines. We won't do much error catching here.
words = line.split()
if words[0].upper() == 'END':
continue
try:
holder = [float(w) for w in words]
current_cmap_data.extend(holder)
except ValueError:
# We assume this is a definition of a new CMAP, so
# terminate the last CMAP if applicable
if current_cmap is not None:
# We have a map to terminate
ty = CmapType(current_cmap_res, current_cmap_data)
self.cmap_types[current_cmap] = ty
self.cmap_types[current_cmap2] = ty
try:
type1 = words[0].upper()
type2 = words[1].upper()
type3 = words[2].upper()
type4 = words[3].upper()
type5 = words[4].upper()
type6 = words[5].upper()
type7 = words[6].upper()
type8 = words[7].upper()
res = conv(words[8], int, 'CMAP resolution', line_index=i, line=line)
except IndexError:
raise CharmmError('Could not parse CMAP data.')
# order the torsions independently
k1 = [type1, type2, type3, type4, type5, type6, type7, type8]
k2 = [type8, type7, type6, type5, type4, type3, type2, type1]
current_cmap = tuple(min(k1, k2))
current_cmap2 = tuple(max(k1, k2))
current_cmap_res = res
current_cmap_data = []
continue
if section.upper() == 'NONBONDED':
# Now get the nonbonded values
words = line.split()
if words[0].upper == 'END':
continue
try:
atype = words[0].upper()
# 1st column is ignored
epsilon = conv(words[2], float, 'vdW epsilon term', line_index=i, line=line)
rmin = conv(words[3], float, 'vdW Rmin/2 term', line_index=i, line=line)
except (IndexError, CharmmError):
# If we haven't read our first nonbonded term yet, we may
# just be parsing the settings that should be used. So
# soldier on
if read_first_nonbonded: raise
for i, word in enumerate(words):
if word.upper() == 'E14FAC':
try:
scee = 1 / float(words[i+1])
except (ValueError, IndexError):
raise CharmmError('Could not parse electrostatic scaling constant')
if self._declared_nbrules:
# We already specified it -- make sure it's the
# same as the one we specified before
_, dt0 = next(iteritems(self.dihedral_types))
diff = abs(dt0[0].scee - scee)
if diff > TINY:
raise CharmmError('Inconsistent 1-4 scalings')
else:
for key, dtl in iteritems(self.dihedral_types):
for dt in dtl:
dt.scee = scee
elif word.upper().startswith('GEOM'):
if self._declared_nbrules and self.combining_rule != 'geometric':
raise CharmmError('Cannot combine parameter files with different '
'combining rules')
self.combining_rule = 'geometric'
declared_geometric = True
continue
else:
# OK, we've read our first nonbonded section for sure now.
# Make sure we did not try to read in a str file that did
# not define GEOM if a previous file did, since
# Lorentz-Berthelot and geometric combining rules are
# incompatible
if (self._declared_nbrules and self.combining_rule == 'geometric' and
not declared_geometric):
raise CharmmError('Cannot combine parameter files with '
'different combining rules')
read_first_nonbonded = True
self._declared_nbrules = True
# See if we have 1-4 parameters
try:
# 4th column is ignored
eps14 = conv(words[5], float, '1-4 vdW epsilon term', line_index=i, line=line)
rmin14 = conv(words[6], float, '1-4 vdW Rmin/2 term', line_index=i, line=line)
except IndexError:
eps14 = rmin14 = None
nonbonded_types[atype] = [epsilon, rmin, eps14, rmin14]
continue
if section.upper() == 'NBFIX':
words = line.split()
if words[0].upper() == 'END':
continue
try:
at1 = words[0].upper()
at2 = words[1].upper()
emin = abs(conv(words[2], float, 'NBFIX Emin', line_index=i, line=line))
rmin = conv(words[3], float, 'NBFIX Rmin', line_index=i, line=line)
try:
emin14 = abs(conv(words[4], float, 'NBFIX Emin 1-4', line_index=i, line=line))
rmin14 = conv(words[5], float, 'NBFIX Rmin 1-4', line_index=i, line=line)
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:
raise CharmmError('Could not parse NBFIX terms.')
self.nbfix_types[(min(at1, at2), max(at1, at2))] = (emin, rmin)
# If we had any CMAP terms, then the last one will not have been added
# yet. Add it here
if current_cmap is not None:
typ = CmapType(current_cmap_res, current_cmap_data)
self.cmap_types[current_cmap] = typ
self.cmap_types[current_cmap2] = typ
# Now we're done. Load the nonbonded types into the relevant AtomType
# 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
# is not satisfied
try:
for key in nonbonded_types:
self.atom_types_str[key].set_lj_params(*nonbonded_types[key])
except KeyError:
warnings.warn('Atom type %s not present in AtomType list' % key, ParameterWarning)
if parameterset is not None:
self.parametersets.append(parameterset)
if own_handle:
f.close()
[docs] def read_topology_file(self, tfile):
"""
Reads _only_ the atom type definitions from a topology file. This is
unnecessary for versions 36 and later of the CHARMM force field.
Parameters
----------
tfile : str
Name of the CHARMM topology file to read
"""
conv = CharmmParameterSet._convert
if isinstance(tfile, str):
own_handle = True
f = iter(CharmmFile(tfile))
else:
own_handle = False
f = tfile
hpatch = tpatch = None # default Head and Tail patches
residues = OrderedDict()
patches = OrderedDict()
hpatches = OrderedDict()
tpatches = OrderedDict()
line = next(f)
line_index = 0
skip_adding_residue = False
try:
while line:
line = line.strip()
if line[:4].upper() == 'MASS':
words = line.split()
try:
idx = conv(words[1], int, 'atom type', line_index=line_index, line=line)
name = words[2].upper()
mass = conv(words[3], float, 'atom mass', line_index=line_index, line=line)
except IndexError:
raise CharmmError('Could not parse MASS section of %s' % tfile)
# The parameter file might or might not have an element name
try:
elem = words[4].upper()
if len(elem) == 2:
elem = elem[0] + elem[1].lower()
atomic_number = AtomicNum[elem]
except (IndexError, KeyError):
# Figure it out from the mass
atomic_number = AtomicNum[element_by_mass(mass)]
atype = AtomType(name=name, number=idx, mass=mass, atomic_number=atomic_number)
self.atom_types_str[atype.name] = atype
self.atom_types_int[atype.number] = atype
self.atom_types_tuple[(atype.name, atype.number)] = atype
elif line[:4].upper() == 'DECL':
pass # Not really sure what this means
elif line[:4].upper() == 'DEFA':
words = line.split()
if len(words) < 5:
warnings.warn('DEFA line has %d tokens; expected 5' % len(words))
else:
it = iter(words[1:5])
for tok, val in zip(it, it):
if val.upper() == 'NONE':
val = None
if tok.upper().startswith('FIRS'):
hpatch = val
elif tok.upper() == 'LAST':
tpatch = val
else:
warnings.warn('DEFA patch %s unknown' % val)
elif line[:4].upper() in ('RESI', 'PRES'):
restype = line[:4].upper()
# Get the residue definition
words = line.split()
resname = words[1].upper()
if resname in self.residues:
warnings.warn('Replacing residue {}'.format(resname), ParameterWarning)
# Assign default patches
hpatches[resname] = hpatch
tpatches[resname] = tpatch
try:
charge = float(words[2])
except (IndexError, ValueError):
warnings.warn('No charge for %s' % resname)
if restype == 'RESI':
res = ResidueTemplate(resname)
elif restype == 'PRES':
res = PatchTemplate(resname)
else:
assert False, 'restype != RESI or PRES'
skip_adding_residue = False
line = next(f)
group = []
ictable = []
while line:
line = line.lstrip()
if line[:5].upper() == 'GROUP':
if group:
res.groups.append(group)
group = []
elif line[:4].upper() == 'ATOM':
words = line.split()
name = words[1].upper()
type = words[2].upper()
charge = float(words[3])
if 'ALPHA' in words:
# This is a polarizable atom.
alpha = float(words[words.index('ALPHA')+1])
thole = 1.3
drude_type = 'DRUD'
if 'THOLE' in words:
thole = float(words[words.index('THOLE')+1])
if 'TYPE' in words:
drude_type = words[words.index('TYPE')+1]
atom = DrudeAtom(name=name, type=type, charge=charge, alpha=alpha, thole=thole, drude_type=drude_type)
else:
atom = Atom(name=name, type=type, charge=charge)
group.append(atom)
res.add_atom(atom)
elif line[:6].upper() == 'DELETE':
words = line.split()
name = words[2].upper()
entity_type = words[1].upper()
if entity_type == 'ATOM':
res.delete_atoms.append(name)
elif entity_type == 'IMPR':
res.delete_impropers.append(words[2:5])
else:
warnings.warn('WARNING: Ignoring "%s" because entity type %s not '
'used.' % (line.strip(), entity_type))
elif line.strip().upper() and line.split()[0].upper() in ('BOND', 'DOUBLE'):
it = iter([w.upper() for w in line.split()[1:]])
for a1, a2 in zip(it, it):
if restype == 'PRES':
# Patches can have bonds that refer to atoms not in the patch, so store these in a list of tuples
order = 1
if line.split()[0].upper() == 'DOUBLE':
order = 2
res.add_bonds.append( (a1, a2, order) )
continue
if a1.startswith('-'):
res.head = res[a2]
continue
if a2.startswith('-'):
res.head = res[a1]
continue
if a1.startswith('+'):
res.tail = res[a2]
continue
if a2.startswith('+'):
res.tail = res[a1]
continue
res.add_bond(a1, a2)
elif line[:4].upper() == 'CMAP':
pass
elif line[:5].upper() == 'DONOR':
pass
elif line[:6].upper() == 'ACCEPT':
pass
elif line[:8].upper() == 'LONEPAIR':
# See: https://www.charmm.org/charmm/documentation/by-version/c40b1/params/doc/lonepair/
# TODO: This currently doesn't handle some formats, like Note 3 in the above URL
words = line.split()
lptype_keyword = words[1][0:4].upper()
if not skip_adding_residue and lptype_keyword not in ['BISE', 'RELA']:
warnings.warn('LONEPAIR type %s not supported; only BISEctor and '
'RELAtive supported' % words[1])
skip_adding_residue = True
break
a1, a2, a3, a4 = words[2:6]
keywords = {words[index][0:4].upper() : float(words[index+1])
for index in range(6,len(words),2) }
r = keywords['DIST'] # angstrom
theta = keywords['ANGL'] # degrees
phi = keywords['DIHE'] # degrees
lptypes = { 'BISE' : 'bisector', 'RELA' : 'relative' }
lonepair = (lptypes[lptype_keyword], a1, a2, a3, a4, r, theta, phi) # TODO: Define a LonePair object?
res.lonepairs.append(lonepair)
elif line[:2].upper() == 'IC':
words = line.split()[1:]
ictable.append(([w.upper() for w in words[:4]],
[float(w) for w in words[4:]]))
elif line[:3].upper() == 'END':
break
elif line[:5].upper() == 'PATCH':
it = iter(line.split()[1:])
for tok, val in zip(it, it):
if val.upper() == 'NONE': val = None
if tok.upper().startswith('FIRS'):
hpatches[resname] = val
elif tok.upper().startswith('LAST'):
tpatches[resname] = val
elif line[:4].upper() in ('IMPR', 'IMPH'):
it = iter(w.upper() for w in line.split()[1:])
for a1, a2, a3, a4 in zip(it, it, it, it):
res._impr.append((a1, a2, a3, a4))
if a2[0] == '-' or a3[0] == '-' or a4 == '-':
res.head = res[a1]
elif line[:10].upper() == 'ANISOTROPY':
words = line.split()
atoms = [res[name] for name in words[1:5]]
keywords = {words[index].upper() : float(words[index+1])
for index in range(5,len(words),2)}
a11 = float(keywords['A11'])
a22 = float(keywords['A22'])
atoms[0].anisotropy = DrudeAnisotropy(*atoms, a11=a11, a22=a22)
elif line[:4].upper() in ('RESI', 'PRES', 'MASS'):
# Back up a line and bail
break
line = next(f)
if group: res.groups.append(group)
_fit_IC_table(res, ictable)
if skip_adding_residue:
# Do not add this residue to the lookup library
continue
elif restype == 'RESI':
residues[resname] = res
elif restype == 'PRES':
patches[resname] = res
else:
assert False, 'restype != RESI or PRES'
# We parsed a line we need to look at. So don't update the
# iterator
continue
# Get the next line and cycle through
line = next(f)
line_index += 1
except StopIteration:
pass
# Go through the patches and add the appropriate one
self.patches.update(patches)
for resname, res in iteritems(residues):
patch_name = hpatches[resname]
if patch_name is not None:
try:
res.first_patch = self.patches[patch_name]
except KeyError:
warnings.warn('Patch %s not found' % patch_name)
patch_name = tpatches[resname]
if patch_name is not None:
try:
res.last_patch = self.patches[patch_name]
except KeyError:
warnings.warn('Patch %s not found' % patch_name)
# Now update the residues and patches with the ones we parsed here
self.residues.update(residues)
if own_handle: f.close()
[docs] def read_stream_file(self, sfile):
"""
Reads RTF and PAR sections from a stream file and dispatches the
sections to read_topology_file or read_parameter_file
Parameters
----------
sfile : str or CharmmStreamFile
Stream file to parse
"""
if isinstance(sfile, CharmmStreamFile):
f = sfile
else:
f = CharmmStreamFile(sfile)
title, section, comments = f.next_section()
while title is not None and section is not None:
words = title.lower().split()
if words[1] == 'rtf':
# This is a Residue Topology File section.
self.read_topology_file(iter(section))
elif words[1].startswith('para'):
# This is a Parameter file section
self.read_parameter_file(section, comments)
title, section, comments = f.next_section()
[docs] def write(self, top=None, par=None, str=None):
""" Write a CHARMM parameter set to a file
Parameters
----------
top : str or file-like object, optional
If provided, the atom types will be written to this file in RTF
format.
par : str or file-like object, optional
If provided, the parameters will be written to this file in PAR
format. Either this or the ``str`` argument *must* be provided
str : str or file-like object, optional
If provided, the atom types and parameters will be written to this
file as separate RTF and PAR cards that can be read as a CHARMM
stream file. Either this or the ``par`` argument *must* be provided
Raises
------
ValueError if both par and str are None
"""
if par is None and str is None:
raise ValueError('Must specify either par *or* str')
if top is not None:
if isinstance(top, string_types):
f = genopen(top, 'w')
ownhandle = True
else:
f = top
ownhandle = False
f.write('*>>>> CHARMM Topology file generated by ParmEd <<<<\n')
f.write('*\n')
self._write_top_to(f, True)
if ownhandle: f.close()
if par is not None:
if isinstance(par, string_types):
f = genopen(par, 'w')
ownhandle = True
else:
f = par
ownhandle = False
f.write('*>>>> CHARMM Parameter file generated by ParmEd <<<<\n')
f.write('*\n\n')
self._write_par_to(f)
if ownhandle: f.close()
if str is not None:
if isinstance(str, string_types):
f = genopen(str, 'w')
ownhandle = True
else:
f = str
ownhandle = False
self._write_str_to(f)
if ownhandle: f.close()
def _write_str_to(self, f):
""" Private method to write stream items to open file object """
f.write('read rtf card\n* Topology generated by ParmEd\n*\n')
self._write_top_to(f, True)
f.write('\nread para card\n* Parameters generated by ParmEd\n*\n')
self._write_par_to(f)
def _write_top_to(self, f, write_version):
""" Private method to write topology items to open file object """
if write_version:
# This version is known to work
f.write('36 1\n')
f.write('\n')
for i, (_, atom) in enumerate(iteritems(self.atom_types)):
f.write('MASS %5d %-6s %9.5f\n' % (i+1, atom.name, atom.mass))
if write_version:
f.write('\nEND\n')
def _write_par_to(self, f):
""" Private method to write parameter items to open file object """
# Find out what the 1-4 electrostatic scaling factors and the 1-4
# van der Waals scaling factors are
scee, scnb = set(), set()
for _, typ in iteritems(self.dihedral_types):
for t in typ:
if t.scee: scee.add(t.scee)
if t.scnb: scnb.add(t.scnb)
if len(scee) > 1 or len(scnb) > 1:
raise ValueError('Mixed 1-4 scaling not supported')
scee = 1.0 if not scee else scee.pop()
scnb = 1.0 if not scnb else scnb.pop()
f.write('ATOMS\n')
self._write_top_to(f, False)
f.write('\nBONDS\n')
written = set()
for key, typ in iteritems(self.bond_types):
if key in written: continue
written.add(key); written.add(tuple(reversed(key)))
f.write('%-6s %-6s %7.2f %10.4f\n' %
(key[0], key[1], typ.k, typ.req))
f.write('\nANGLES\n')
written = set()
for key, typ in iteritems(self.angle_types):
if key in written: continue
written.add(key)
written.add(tuple(reversed(key)))
f.write('%-6s %-6s %-6s %7.2f %8.2f\n' % (key[0], key[1], key[2], typ.k, typ.theteq))
f.write('\nDIHEDRALS\n')
written = set()
for key, typ in iteritems(self.dihedral_types):
if key in written: continue
written.add(key)
written.add(tuple(reversed(key)))
for tor in typ:
f.write('%-6s %-6s %-6s %-6s %11.4f %2d %8.2f\n' %
(key[0], key[1], key[2], key[3], tor.phi_k, int(tor.per), tor.phase))
f.write('\nIMPROPERS\n')
written = set()
for key, typ in sorted(iteritems(self.improper_periodic_types), key=lambda x: x[0]):
f.write('%-6s %-6s %-6s %-6s %11.4f %2d %8.2f\n' %
(key[0], key[1], key[2], key[3], typ.phi_k, int(typ.per), typ.phase))
for key, typ in iteritems(self.improper_types):
f.write('%-6s %-6s %-6s %-6s %11.4f %2d %8.2f\n' %
(key[0], key[1], key[2], key[3], typ.psi_k, 0, typ.psi_eq))
if self.cmap_types:
f.write('\nCMAPS\n')
written = set()
for key, typ in iteritems(self.cmap_types):
if key in written: continue
written.add(key); written.add(tuple(reversed(key)))
f.write('%-6s %-6s %-6s %-6s %-6s %-6s %-6s %-6s %5d\n\n' %
(key[0], key[1], key[2], key[3], key[4], key[5], key[6],
key[7], typ.resolution))
i = 0
for val in typ.grid:
if i:
if i % 5 == 0:
f.write('\n')
if i % typ.resolution == 0:
f.write('\n')
i = 0
elif i % typ.resolution == 0:
f.write('\n\n')
i = 0
i += 1
f.write(' %13.6f' % val)
f.write('\n\n\n')
f.write('\nNONBONDED nbxmod 5 atom cdiel fshift vatom vdistance vfswitch -\ncutnb 14.0 '
'ctofnb 12.0 ctonnb 10.0 eps 1.0 e14fac %s wmin 1.5%s\n\n' %
(1/scee, ' GEOM' if self.combining_rule == 'geometric' else ''))
for key, typ in iteritems(self.atom_types):
f.write('%-6s %14.6f %10.6f %14.6f' % (key, 0.0, -abs(typ.epsilon), typ.rmin))
if typ.epsilon == typ.epsilon_14 and typ.rmin == typ.rmin_14:
f.write('%10.6f %10.6f %14.6f\n' % (0, -abs(typ.epsilon)/scnb, typ.rmin))
else:
f.write('%10.6f %10.6f %14.6f\n' % (0, -abs(typ.epsilon_14), typ.rmin_14))
f.write('\nEND\n')
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _fit_IC_table(res, ictable):
"""
Determines cartesian coordinates from an internal coordinate table stored in
CHARMM residue topology files
Parameters
----------
res : ResidueTemplate
The residue template for which coordinates are being determined
ictable : list[tuple(atoms, measurements)]
The internal coordinate table
Notes
-----
This method assigns an xx, xy, and xz attribute to ``res``. For the time
being, this is just a placeholder, as its functionality has not yet been
implemented (CHARMM does not use a 'traditional' Z-matrix, and I don't know
of any existing code that will compute the proper cartesian coordinates from
the form used by CHARMM)
"""
for atom in res:
atom.xx = atom.xy = atom.xz = 0.0