"""
This module contains an amber prmtop class that will read in all
parameters and allow users to manipulate that data and write a new
prmtop object.
"""
from __future__ import division
import copy as _copy
from collections import defaultdict
from math import sqrt
from warnings import warn
import numpy as np
from .. import unit as u
from ..constants import (DEG_TO_RAD, IFBOX, IFCAP, IFPERT, MBONA, MBPER, MDPER,
MGPER, MPHIA, MTHETA, NATOM, NATYP, NBONA, NBONH,
NBPER, NCOPY, NDPER, NEXT, NGPER, NHPARM, NMXRS, NNB,
NPARM, NPHB, NPHIA, NPHIH, NPTRA, NRES, NTHETA,
NTHETH, NTYPES, NUMANG, NUMBND, NUMEXTRA, RAD_TO_DEG,
SMALL, TRUNCATED_OCTAHEDRON_ANGLE)
from ..exceptions import AmberError, AmberWarning, MoleculeError
from ..geometry import box_lengths_and_angles_to_vectors
from ..periodic_table import AtomicNum, element_by_mass
from ..residue import ALLION_NAMES, SOLVENT_NAMES
from ..structure import Structure, needs_openmm
from ..topologyobjects import (Angle, AngleType, Atom, AtomList, AtomType,
Bond, BondType, Dihedral, DihedralType,
DihedralTypeList, ExtraPoint, Cmap, CmapType)
from ..utils.six import iteritems, string_types
from ..utils.six.moves import range, zip
from ..vec3 import Vec3
from .amberformat import AmberFormat
from .asciicrd import AmberAsciiRestart
from .netcdffiles import NetCDFRestart
try:
from simtk import openmm as mm
from simtk.openmm import app
except ImportError:
mm = app = None
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[docs]class AmberParm(AmberFormat, Structure):
"""
Amber Topology (parm7 format) class. Gives low, and some high, level access
to topology data. You can interact with the raw data in the topology file
directly or interact with some of the high-level classes comprising the
system topology and parameters.
Parameters
----------
prm_name : str, optional
If provided, this file is parsed and the data structures will be loaded
from the data in this file
xyz : str or array, optional
If provided, the coordinates and unit cell dimensions from the provided
Amber inpcrd/restart file will be loaded into the molecule, or the
coordinates will be loaded from the coordinate array
box : array, optional
If provided, the unit cell information will be set from the provided
unit cell dimensions (a, b, c, alpha, beta, and gamma, respectively)
Attributes
----------
parm_data : dict {str : list}
A dictionary that maps FLAG names to all of the data contained in that
section of the Amber file.
formats : dict {str : FortranFormat}
A dictionary that maps FLAG names to the FortranFormat instance in which
the data is stored in that section
parm_comments : dict {str : list}
A dictionary that maps FLAG names to the list of COMMENT lines that were
stored in the original file
flag_list : list
An ordered list of all FLAG names. This must be kept synchronized with
`parm_data`, `formats`, and `parm_comments` such that every item in
`flag_list` is a key to those 3 dicts and no other keys exist
charge_flag : str='CHARGE'
The name of the name of the FLAG that describes partial atomic charge
data. If this flag is found, then its data are multiplied by the
ELECTROSTATIC_CONSTANT to convert back to fractions of electrons
version : str
The VERSION string from the Amber file
name : str
The file name of the originally parsed file (set to the fname parameter)
LJ_types : dict {str : int}
A mapping whose keys are atom types paired with the nonbonded index of
that type
LJ_radius : list(float)
A list of floats corresponding to the Rmin/2 parameter for every
Lennard-Jones type. The indexes are the nonbonded index (`nb_idx`
attribute of the `Atom` class) minus 1 to account for indexing from 0 in
Python. The values are in Angstroms. To get the radius for a particular
atom type, you can use `LJ_radius[LJ_types["type"]-1]`
LJ_depth : list(float)
A list of Lennard-Jones epsilon parameters laid out the same way as
LJ_radius, described above.
atoms : AtomList(Atom)
List of all atoms in the system
residues : ResidueList(Residue)
List of all residues in the system
bonds : TrackedList(Bond)
List of bonds between two atoms in the system
angles : TrackedList(Angle)
List of angles between three atoms in the system
dihedrals : TrackedList(Angle)
List of all proper and improper torsions between 4 atoms in the system
box : list of 6 floats
Periodic boundary unit cell dimensions and angles
bond_types : TrackedList(BondType)
The bond types containing the parameters for each bond stretching term
angle_types : TrackedList(AngleType)
The angle types containing the parameters for each angle bending term
dihedral_types : TrackedList(DihedralType)
The dihedral types containing the parameters for each torsional term
bonds_inc_h : iterator(Bond)
Read-only generator that loops through all bonds that contain Hydrogen
bonds_without_h : iterator(Bond)
Read-only generator that loops through all bonds that do not contain
Hydrogen
angles_inc_h : iterator(Angle)
Read-only generator that loops through all angles that contain Hydrogen
angles_without_h : iterator(Angle)
Read-only generator that loops through all angles that do not contain
Hydrogen
dihedrals_inc_h : iterator(Dihedral)
Read-only generator that loops through all dihedrals that contain
Hydrogen
dihedrals_without_h : iterator(Dihedral)
Read-only generator that loops through all dihedrals that do not contain
Hydrogen
chamber : bool=False
On AmberParm instances, this is always False to indicate that it is not
a CHAMBER-style topology file
amoeba : bool=False
On AmberParm instances, this is always False to indicate that it is not
an AMOEBA-style topology file
has_cmap : bool=False
On AmberParm instances, this is always False to indicate that it does
not have correction maps (unique to CHARMM force field and chamber
topologies)
"""
_cmap_prefix = ""
#===================================================
[docs] def __init__(self, prm_name=None, xyz=None, box=None, rst7_name=None):
"""
Instantiates an AmberParm object from data in prm_name and establishes
validity based on presence of POINTERS and CHARGE sections. In general,
you should use LoadParm from the readparm module instead. LoadParm will
correctly dispatch the object to the 'correct' flavor of AmberParm
"""
AmberFormat.__init__(self, prm_name)
Structure.__init__(self)
self.hasvels = False
self.hasbox = False
self._box = None
self.crdname = None
if xyz is None and rst7_name is not None:
warn('rst7_name keyword is deprecated. Use xyz instead',
DeprecationWarning)
self.crdname = xyz = rst7_name
elif xyz is not None and rst7_name is not None:
warn('rst7_name keyword is deprecated and ignored in favor of xyz',
DeprecationWarning)
if isinstance(xyz, string_types):
self.crdname = xyz
if prm_name is not None:
self.initialize_topology(xyz, box)
#===================================================
[docs] def initialize_topology(self, xyz=None, box=None):
"""
Initializes topology data structures, like the list of atoms, bonds,
etc., after the topology file has been read.
"""
from parmed import load_file
# We need to handle RESIDUE_ICODE properly since it may have picked up
# some extra values
if 'RESIDUE_ICODE' in self.flag_list:
self._truncate_array('RESIDUE_ICODE',
self.parm_data['POINTERS'][NRES])
# Set up some of the attributes provided
self.pointers = {}
self.LJ_types = {}
self.LJ_radius = []
self.LJ_depth = []
self.load_pointers()
self.fill_LJ()
# The way we check if the combining rules are geometric vs.
# Lorentz-Berthelot is to try and change the combining rules to
# geometric *if we detect NBFIXes* and see if we still have NBFIXes. If
# not, then our combining rules are clearly geometric.
if not self.has_1012() and self.has_NBFIX():
self.combining_rule = 'geometric'
if self.has_NBFIX():
self.combining_rule = 'lorentz'
# Instantiate the Structure data structures
self.load_structure()
if isinstance(xyz, string_types):
f = load_file(xyz, skip_bonds=True)
if not hasattr(f, 'coordinates') or f.coordinates is None:
raise TypeError('%s does not have coordinates' % xyz)
self.coordinates = f.coordinates
if hasattr(f, 'velocities') and f.velocities is not None:
self.velocities = f.velocities
if hasattr(f, 'box') and f.box is not None and box is None:
self.box = f.box
else:
self.coordinates = xyz
if box is not None:
self.box = box
# If all else fails, set the box from the prmtop file
if self.parm_data['POINTERS'][IFBOX] > 0 and self.box is None:
box = self.parm_data['BOX_DIMENSIONS']
self.box = list(box[1:]) + [box[0], box[0], box[0]]
self.hasbox = self.box is not None
#===================================================
[docs] @classmethod
def from_rawdata(cls, rawdata):
"""
Take the raw data from a AmberFormat object and initialize an AmberParm
from that data.
Parameters
----------
rawdata : :class:`AmberFormat`
An AmberFormat instance that has already been instantiated
Returns
-------
parm : :class:`AmberParm`
An instance of this type from the data in rawdata
"""
inst = cls()
inst.name = rawdata.name
inst.version = rawdata.version
inst.formats = rawdata.formats
inst.parm_data = rawdata.parm_data
inst.parm_comments = rawdata.parm_comments
inst.flag_list = rawdata.flag_list
inst.initialize_topology()
# See if the rawdata has any kind of structural attributes, like
# coordinates and an atom list with positions and/or velocities
if hasattr(rawdata, 'coordinates'):
inst.coordinates = _copy.copy(rawdata.coordinates)
if hasattr(rawdata, 'velocities'):
inst.velocities = _copy.copy(rawdata.velocities)
if hasattr(rawdata, 'box'):
inst.box = _copy.copy(rawdata.box)
inst.hasbox = inst.box is not None
inst.hasvels = inst.velocities is not None
n_copy = inst.pointers.get('NCOPY', 1)
if n_copy >= 2:
inst._label_alternates()
return inst
#===================================================
[docs] @classmethod
def from_structure(cls, struct, copy=False):
"""
Take a Structure instance and initialize an AmberParm instance from that
data.
Parameters
----------
struct : :class:`Structure`
The input structure from which to construct an AmberParm instance
copy : bool
If True, the input struct is deep-copied to make sure it does not
share any objects with the original ``struct``. Default is False
Returns
-------
inst : :class:`AmberParm`
The AmberParm instance derived from the input structure
Raises
------
TypeError
If the structure has parameters not supported by the standard Amber
force field (i.e., standard bond, angle, and dihedral types)
Notes
-----
Due to the nature of the prmtop file, struct almost *always* returns a
deep copy. The one exception is when struct is already of type
:class:`AmberParm`, in which case the original object is returned unless
``copy`` is ``True``.
"""
if type(struct) is cls:
if copy:
return _copy.copy(struct)
return struct
if struct.unknown_functional:
raise TypeError('Cannot instantiate an AmberParm from unknown functional')
if (struct.urey_bradleys or struct.impropers or struct.cmaps or
struct.trigonal_angles or struct.pi_torsions or
struct.out_of_plane_bends or struct.stretch_bends or
struct.torsion_torsions or struct.multipole_frames):
if (struct.trigonal_angles or struct.pi_torsions or
struct.out_of_plane_bends or struct.torsion_torsions or
struct.multipole_frames or struct.stretch_bends):
raise TypeError('AmberParm does not support all of the '
'parameters defined in the input Structure')
# Maybe it just has CHARMM parameters?
raise TypeError('AmberParm does not support all of the parameters '
'defined in the input Structure. Try ChamberParm')
# Convert all RB torsions to propers and delete the originals.
for dihedral in struct.rb_torsions:
proper_types = DihedralTypeList.from_rbtorsion(dihedral.type)
for proper_type in proper_types:
proper = _copy.copy(dihedral)
proper.type = proper_type
struct.dihedrals.append(proper)
struct.dihedral_types.append(proper.type)
del struct.rb_torsions[:]
del struct.rb_torsion_types[:]
struct.dihedrals.claim()
struct.dihedral_types.claim()
inst = struct.copy(cls, split_dihedrals=True)
inst.update_dihedral_exclusions()
inst._add_missing_13_14()
del inst.adjusts[:]
inst.pointers = {}
inst.LJ_types = {}
nbfixes = inst._set_nbidx_lj_params()
inst._add_standard_flags()
inst.pointers['NATOM'] = len(inst.atoms)
inst.parm_data['POINTERS'][NATOM] = len(inst.atoms)
# pmemd likes to skip torsions with periodicities of 0, which may be
# present as a way to hack entries into the 1-4 pairlist. See
# https://github.com/ParmEd/ParmEd/pull/145 for discussion. The solution
# here is to simply set that periodicity to 1.
inst._cleanup_dihedrals_with_periodicity_zero()
inst.remake_parm()
inst._set_nonbonded_tables(nbfixes)
n_copy = inst.pointers.get('NCOPY', 1)
if n_copy >= 2:
inst._label_alternates()
# Setting box sets some of the flag data appropriately
inst.box = inst.get_box()
return inst
#===================================================
def _set_nbidx_lj_params(self):
nbfixes = self.atoms.assign_nbidx_from_types()
# Give virtual sites a name that Amber understands
for atom in self.atoms:
atom.type = atom.type if not isinstance(atom, ExtraPoint) else 'EP'
# Fill the Lennard-Jones arrays/dicts
ntyp = 0
for atom in self.atoms:
self.LJ_types[atom.type] = atom.nb_idx
ntyp = max(ntyp, atom.nb_idx)
self.LJ_radius = [0 for i in range(ntyp)]
self.LJ_depth = [0 for i in range(ntyp)]
for atom in self.atoms:
self.LJ_radius[atom.nb_idx-1] = atom.atom_type.rmin
self.LJ_depth[atom.nb_idx-1] = atom.atom_type.epsilon
return nbfixes
#===================================================
def __copy__(self):
""" Needs to copy a few additional data structures """
other = super(AmberParm, self).__copy__()
other.initialize_topology()
# Copy coordinates, box, etc. if applicable
other.coordinates = _copy.copy(self.get_coordinates('all'))
other._box = _copy.copy(self._box)
# Now we should have a full copy
return other
#===================================================
def __getitem__(self, selection):
other = super(AmberParm, self).__getitem__(selection)
if isinstance(other, Atom):
return other
other.pointers = {}
self._copy_lj_data(other)
other._add_standard_flags()
other.remake_parm()
other._set_nonbonded_tables()
# Add back the original L-J tables to restore any NBFIXes that may have
# been there
other.parm_data['LENNARD_JONES_ACOEF'] = \
self.parm_data['LENNARD_JONES_ACOEF'][:]
other.parm_data['LENNARD_JONES_BCOEF'] = \
self.parm_data['LENNARD_JONES_BCOEF'][:]
if 'LENNARD_JONES_CCOEF' in self.parm_data:
other.parm_data['LENNARD_JONES_CCOEF'] = \
self.parm_data['LENNARD_JONES_CCOEF'][:]
return other
def _copy_lj_data(self, other):
""" Copies Lennard-Jones lists and dicts from myself to a copy """
other.LJ_types = self.LJ_types.copy()
other.LJ_radius = _copy.copy(self.LJ_radius)
other.LJ_depth = _copy.copy(self.LJ_depth)
for atom in other.atoms:
other.LJ_radius[atom.nb_idx-1] = atom.atom_type.rmin
other.LJ_depth[atom.nb_idx-1] = atom.atom_type.epsilon
#===================================================
def __imul__(self, ncopies, other=None):
super(AmberParm, self).__imul__(ncopies, other)
self.remake_parm()
return self
def __iadd__(self, other):
if self.has_NBFIX() or self.has_1012():
raise ValueError('Cannot combine Amber systems with NBFIX or 10-12 parameters')
super(AmberParm, self).__iadd__(other)
# Make sure we properly recompute the LJ tables
nbfixes = self._set_nbidx_lj_params()
self.remake_parm()
self._set_nonbonded_tables(nbfixes)
return self
#===================================================
[docs] def load_pointers(self):
"""
Loads the data in POINTERS section into a pointers dictionary with each
key being the pointer name according to http://ambermd.org/formats.html
"""
self.pointers["NATOM"] = self.parm_data["POINTERS"][NATOM]
self.pointers["NTYPES"] = self.parm_data["POINTERS"][NTYPES]
self.pointers["NBONH"] = self.parm_data["POINTERS"][NBONH]
self.pointers["MBONA"] = self.parm_data["POINTERS"][MBONA]
self.pointers["NTHETH"] = self.parm_data["POINTERS"][NTHETH]
self.pointers["MTHETA"] = self.parm_data["POINTERS"][MTHETA]
self.pointers["NPHIH"] = self.parm_data["POINTERS"][NPHIH]
self.pointers["MPHIA"] = self.parm_data["POINTERS"][MPHIA]
self.pointers["NHPARM"] = self.parm_data["POINTERS"][NHPARM]
self.pointers["NPARM"] = self.parm_data["POINTERS"][NPARM]
self.pointers["NEXT"] = self.parm_data["POINTERS"][NEXT]
self.pointers["NNB"] = self.parm_data["POINTERS"][NNB] # alias for above
self.pointers["NRES"] = self.parm_data["POINTERS"][NRES]
self.pointers["NBONA"] = self.parm_data["POINTERS"][NBONA]
self.pointers["NTHETA"] = self.parm_data["POINTERS"][NTHETA]
self.pointers["NPHIA"] = self.parm_data["POINTERS"][NPHIA]
self.pointers["NUMBND"] = self.parm_data["POINTERS"][NUMBND]
self.pointers["NUMANG"] = self.parm_data["POINTERS"][NUMANG]
self.pointers["NPTRA"] = self.parm_data["POINTERS"][NPTRA]
self.pointers["NATYP"] = self.parm_data["POINTERS"][NATYP]
self.pointers["NPHB"] = self.parm_data["POINTERS"][NPHB]
self.pointers["IFPERT"] = self.parm_data["POINTERS"][IFPERT]
self.pointers["NBPER"] = self.parm_data["POINTERS"][NBPER]
self.pointers["NGPER"] = self.parm_data["POINTERS"][NGPER]
self.pointers["NDPER"] = self.parm_data["POINTERS"][NDPER]
self.pointers["MBPER"] = self.parm_data["POINTERS"][MBPER]
self.pointers["MGPER"] = self.parm_data["POINTERS"][MGPER]
self.pointers["MDPER"] = self.parm_data["POINTERS"][MDPER]
self.pointers["IFBOX"] = self.parm_data["POINTERS"][IFBOX]
self.pointers["NMXRS"] = self.parm_data["POINTERS"][NMXRS]
self.pointers["IFCAP"] = self.parm_data["POINTERS"][IFCAP]
self.pointers["NUMEXTRA"] = self.parm_data["POINTERS"][NUMEXTRA]
if self.parm_data['POINTERS'][IFBOX] > 0:
self.pointers['IPTRES'] = self.parm_data['SOLVENT_POINTERS'][0]
self.pointers['NSPM'] = self.parm_data['SOLVENT_POINTERS'][1]
self.pointers['NSPSOL'] = self.parm_data['SOLVENT_POINTERS'][2]
# The next is probably only there for LES-prmtops
try:
self.pointers["NCOPY"] = self.parm_data["POINTERS"][NCOPY]
except:
pass
# If CMAP is not present, don't load the pointers
if self.has_cmap:
self.pointers['CMAP'] = self.parm_data[self._cmap_prefix + 'CMAP_COUNT'][0]
self.pointers['CMAP_TYPES'] = self.parm_data[self._cmap_prefix + 'CMAP_COUNT'][1]
#===================================================
[docs] def load_structure(self):
"""
Loads all of the topology instance variables. This is necessary if we
actually want to modify the topological layout of our system
(like deleting atoms)
"""
self._check_section_lengths()
self._load_atoms_and_residues()
self.load_atom_info()
self._load_bond_info()
self._load_angle_info()
self._load_dihedral_info()
self._load_cmap_info()
self._load_extra_exclusions()
super(AmberParm, self).unchange()
#===================================================
[docs] def load_atom_info(self):
"""
Loads atom properties into the atoms that have been loaded. If any
arrays are too short or too long, an IndexError will be raised
"""
# Collect all of the atom properties present in our topology file
zeros = _zeros(len(self.atoms))
anam = self.parm_data['ATOM_NAME']
chg = self.parm_data['CHARGE']
mass = self.parm_data['MASS']
nbtyp = self.parm_data['ATOM_TYPE_INDEX']
atyp = self.parm_data['AMBER_ATOM_TYPE']
join = self.parm_data['JOIN_ARRAY']
irot = self.parm_data['IROTAT']
tree = self.parm_data['TREE_CHAIN_CLASSIFICATION']
try:
radii = self.parm_data['RADII']
except KeyError:
radii = zeros
try:
screen = self.parm_data['SCREEN']
except KeyError:
screen = zeros
try:
atnum = self.parm_data['ATOMIC_NUMBER']
replace_atnum = True
except KeyError:
atnum = [AtomicNum[element_by_mass(m)] for m in mass]
replace_atnum = False
try:
occu = self.parm_data['ATOM_OCCUPANCY']
except KeyError:
occu = zeros
try:
bfac = self.parm_data['ATOM_BFACTOR']
except KeyError:
bfac = zeros
try:
anum = self.parm_data['ATOM_NUMBER']
except KeyError:
anum = [-1 for atom in self.atoms]
for i, atom in enumerate(self.atoms):
atom.name = anam[i]
atom.charge = chg[i]
atom.mass = mass[i]
atom.nb_idx = nbtyp[i]
atom.type = atyp[i]
atom.join = join[i]
atom.irotat = irot[i]
atom.tree = tree[i]
atom.solvent_radius = radii[i]
atom.screen = screen[i]
if replace_atnum or atom.atomic_number == 0:
if atnum[i] == -1:
atom.atomic_number = AtomicNum[element_by_mass(mass[i])]
else:
atom.atomic_number = atnum[i]
atom.atom_type = AtomType(atyp[i], None, mass[i], atnum[i])
atom.occupancy = occu[i]
atom.bfactor = bfac[i]
atom.number = anum[i]
depth = self.LJ_depth[atom.nb_idx-1]
radius = self.LJ_radius[atom.nb_idx-1]
try:
depth14 = self.LJ_14_depth[atom.nb_idx-1]
radius14 = self.LJ_14_radius[atom.nb_idx-1]
except AttributeError:
depth14 = radius14 = None
atom.atom_type.set_lj_params(depth, radius, depth14, radius14)
#===================================================
[docs] def ptr(self, pointer):
"""
Returns the value of the given pointer, and converts to upper-case so
it's case-insensitive. A non-existent pointer meets with a KeyError
Parameters
----------
pointer : str
The AMBER pointer for which to extract the value
Returns
-------
int
The returned integer is the value of that pointer
"""
return self.pointers[pointer.upper()]
#===================================================
[docs] def write_rst7(self, name, netcdf=None):
"""
Writes a restart file with the current coordinates and velocities and
box info if it's present
Parameters
----------
name : str
Name of the file to write the restart file to
netcdf : bool=False
If True, write a NetCDF-format restart file (requires a NetCDF
backend; scipy, netCDF4, or ScientificPython; to be installed)
Notes
-----
If `netcdf` is not specified and the filename extension given by `name`
is `.ncrst`, the a NetCDF restart file will be written. However, an
explicit value for `netcdf` will override any filename extensions.
"""
# By default, determine file type by extension (.ncrst is NetCDF)
netcdf = netcdf or (netcdf is None and name.endswith('.ncrst'))
# Check that we have a rst7 loaded, then overwrite it with a new one if
# necessary
rst7 = Rst7(natom=len(self.atoms))
# Now fill in the rst7 coordinates
rst7.coordinates = [0.0 for i in range(len(self.atoms)*3)]
if self.velocities is not None:
rst7.vels = [0.0 for i in range(len(self.atoms)*3)]
for i, at in enumerate(self.atoms):
i3 = i * 3
rst7.coordinates[i3 ] = at.xx
rst7.coordinates[i3+1] = at.xy
rst7.coordinates[i3+2] = at.xz
if rst7.hasvels:
rst7.vels[i3 ] = at.vx
rst7.vels[i3+1] = at.vy
rst7.vels[i3+2] = at.vz
rst7.box = _copy.copy(self.box)
# Now write the restart file
rst7.write(name, netcdf)
#===================================================
[docs] def write_parm(self, name):
"""
Writes the current data in parm_data into a new topology file with a given name.
Parameters
----------
name : str or file-like
The name of the file to write the prmtop to or the file object to write to
"""
self.remake_parm()
AmberFormat.write_parm(self, name)
#===================================================
[docs] def remake_parm(self):
"""
Fills :attr:`parm_data` from the data in the parameter and topology
arrays (e.g., :attr:`atoms`, :attr:`bonds`, :attr:`bond_types`, ...)
"""
# Get rid of terms containing deleted atoms and empty residues
self.prune_empty_terms()
self.residues.prune()
self.rediscover_molecules()
# Transfer information from the topology lists
self._xfer_atom_info()
self._xfer_residue_info()
self._xfer_bond_info()
self._xfer_angle_info()
self._xfer_dihedral_info()
self._xfer_cmap_properties()
self._set_ifbox()
# Load the pointers dict
self.load_pointers()
# Mark atom list as unchanged
super(AmberParm, self).unchange()
#===================================================
[docs] def is_changed(self):
"""
Determines if any of the topological arrays have changed since the
last upload
"""
is_changed = super(AmberParm, self).is_changed()
if is_changed and hasattr(self, '_topology'):
del self._topology
return is_changed
#===================================================
[docs] def strip(self, selection):
"""
Deletes a subset of the atoms corresponding to an atom-based selection.
Parameters
----------
selection : AmberMask, str, or iterable of bool
This is the selection of atoms that will be deleted from this
structure. If it is a string, it will be interpreted as an
AmberMask. If it is an AmberMask, it will be converted to a
selection of atoms. If it is an iterable, it must be the same length
as the `atoms` list.
"""
super(AmberParm, self).strip(selection)
self.remake_parm()
#===================================================
[docs] def rediscover_molecules(self, solute_ions=True, fix_broken=True):
"""
This determines the molecularity and sets the ATOMS_PER_MOLECULE and
SOLVENT_POINTERS sections of the prmtops. Returns the new atom sequence
in terms of the 'old' atom indexes if re-ordering was necessary to fix
the tleap bug. Returns None otherwise.
"""
# Bail out of we are not doing a solvated prmtop
if self.parm_data['POINTERS'][IFBOX] == 0 or self.box is None:
return None
owner = set_molecules(self)
all_solvent = SOLVENT_NAMES
if not solute_ions:
all_solvent = all_solvent | ALLION_NAMES
first_solvent = None
for i, res in enumerate(self.residues):
if res.name in all_solvent:
first_solvent = i
break
# Now remake our SOLVENT_POINTERS and ATOMS_PER_MOLECULE section
self.parm_data['SOLVENT_POINTERS'] = [first_solvent, len(owner), 0]
# Find the first solvent molecule if there are any
if first_solvent is None:
self.parm_data['SOLVENT_POINTERS'][0] = len(self.residues)
self.parm_data['SOLVENT_POINTERS'][2] = len(owner) + 1
else:
self.parm_data['SOLVENT_POINTERS'][2] = \
self.residues[first_solvent].atoms[0].marked
# Now set up ATOMS_PER_MOLECULE and catch any errors
self.parm_data['ATOMS_PER_MOLECULE'] = [len(mol) for mol in owner]
# Check that all of our molecules are contiguous, because we have to
# re-order atoms if they're not
try:
for mol in owner:
sortedmol = sorted(list(mol))
for i in range(1, len(sortedmol)):
if sortedmol[i] != sortedmol[i-1] + 1:
raise StopIteration()
except StopIteration:
if not fix_broken:
raise MoleculeError('Molecule atoms are not contiguous!')
# Non-contiguous molecules detected... time to fix (ugh!)
warn('Molecule atoms are not contiguous! I am attempting to '
'reorder the atoms to fix this.', AmberWarning)
# Make sure that no residues are split up by this
for res in self.residues:
molid = res.atoms[0].marked
for atom in res:
if molid != atom.marked:
warn('Residues cannot be part of 2 molecules! Molecule '
'section will not be correctly set. [Offending '
'residue is %d: %r]' % (res.idx, res),
AmberWarning)
return None
new_atoms = AtomList()
for mol in owner:
for idx in sorted(mol):
new_atoms.append(self.atoms[idx])
self.atoms = new_atoms
# Re-sort our residues and residue list for new atom ordering
for residue in self.residues:
residue.sort()
self.residues.sort()
return owner
return None
#===================================================
[docs] def fill_LJ(self):
"""
Calculates the Lennard-Jones parameters (Rmin/2 and epsilon) for each
atom type by computing their values from the A and B coefficients of
each atom interacting with itself.
This fills the :attr:`LJ_radius`, :attr:`LJ_depth`, and :attr:`LJ_types`
data structures.
"""
self.LJ_radius = [] # empty LJ_radii so it can be re-filled
self.LJ_depth = [] # empty LJ_depths so it can be re-filled
self.LJ_types = {} # empty LJ_types so it can be re-filled
one_sixth = 1 / 6 # we need to raise some numbers to the 1/6th power
pd = self.parm_data
acoef = pd['LENNARD_JONES_ACOEF']
bcoef = pd['LENNARD_JONES_BCOEF']
natom = self.pointers['NATOM']
ntypes = self.pointers['NTYPES']
for i in range(natom): # fill the LJ_types array
self.LJ_types[pd["AMBER_ATOM_TYPE"][i]] = pd["ATOM_TYPE_INDEX"][i]
for i in range(ntypes):
lj_index = pd["NONBONDED_PARM_INDEX"][ntypes*i+i] - 1
if lj_index < 0 or acoef[lj_index] < 1.0e-10 or bcoef[lj_index] < 1.0e-10:
self.LJ_radius.append(0)
self.LJ_depth.append(0)
else:
factor = 2 * acoef[lj_index] / bcoef[lj_index]
self.LJ_radius.append(pow(factor, one_sixth) * 0.5)
self.LJ_depth.append(bcoef[lj_index] / 2 / factor)
#===================================================
[docs] def recalculate_LJ(self):
"""
Fills the ``LENNARD_JONES_ACOEF`` and ``LENNARD_JONES_BCOEF`` arrays in
the :attr:`parm_data` raw data dictionary by applying the canonical
Lorentz-Berthelot combining rules to the values in :attr:`LJ_radius` and
:attr:`LJ_depth`.
Notes
-----
This will undo any off-diagonal L-J modifications you may have made, so
call this function with care.
"""
assert self.combining_rule in ('lorentz', 'geometric'), "Unrecognized combining rule"
if self.combining_rule == 'lorentz':
comb_sig = lambda sig1, sig2: 0.5 * (sig1 + sig2)
elif self.combining_rule == 'geometric':
comb_sig = lambda sig1, sig2: sqrt(sig1 * sig2)
pd = self.parm_data
ntypes = self.pointers['NTYPES']
fac = 2**(-1/6) * 2
LJ_sigma = [x*fac for x in self.LJ_radius]
fac = 2**(1/6)
for i in range(ntypes):
for j in range(i, ntypes):
index = pd['NONBONDED_PARM_INDEX'][ntypes*i+j] - 1
if index < 0:
# This indicates *either* a 10-12 potential for this pair
# _or_ it indicates a placeholder for HW atom type
# interactions and serves as a tag for fast water routines
# (or did in the past, anyway). So just skip to the next one
continue
rij = comb_sig(LJ_sigma[i], LJ_sigma[j]) * fac
wdij = sqrt(self.LJ_depth[i] * self.LJ_depth[j])
pd["LENNARD_JONES_ACOEF"][index] = wdij * rij**12
pd["LENNARD_JONES_BCOEF"][index] = 2 * wdij * rij**6
#===================================================
[docs] def has_NBFIX(self):
"""
This routine determines whether there are any off-diagonal Lennard-Jones
modifications (i.e., if any two atoms have a L-J pair interaction that
does not match the combined L-J parameters for that pair).
Returns
-------
nbfix : bool
If True, off-diagonal elements in the combined Lennard-Jones matrix
exist. If False, they do not.
"""
assert self.combining_rule in ('lorentz', 'geometric'), "Unrecognized combining rule"
if self.combining_rule == 'lorentz':
comb_sig = lambda sig1, sig2: 0.5 * (sig1 + sig2)
elif self.combining_rule == 'geometric':
comb_sig = lambda sig1, sig2: sqrt(sig1 * sig2)
fac = 2**(-1/6) * 2
LJ_sigma = [x*fac for x in self.LJ_radius]
pd = self.parm_data
ntypes = self.parm_data['POINTERS'][NTYPES]
fac = 2**(1/6)
for i in range(ntypes):
for j in range(ntypes):
idx = pd['NONBONDED_PARM_INDEX'][ntypes*i+j] - 1
if idx < 0: continue
rij = comb_sig(LJ_sigma[i], LJ_sigma[j]) * fac
wdij = sqrt(self.LJ_depth[i] * self.LJ_depth[j])
a = pd['LENNARD_JONES_ACOEF'][idx]
b = pd['LENNARD_JONES_BCOEF'][idx]
if a == 0 or b == 0:
if a != 0 or b != 0 or (wdij != 0 and rij != 0):
return True
elif (abs((a - (wdij * rij**12)) / a) > 1e-6 or
abs((b - (2 * wdij * rij**6)) / b) > 1e-6):
return True
return False
#===================================================
[docs] def has_1012(self):
"""
This routine determines whether there are any defined 10-12
Lennard-Jones interactions that are non-zero
Returns
-------
has_10_12 : bool
If True, 10-12 interactions *are* defined for this particular system
"""
indices_with_1012 = []
ntypes = self.parm_data['POINTERS'][NTYPES]
for i in range(ntypes):
for j in range(ntypes):
idx = self.parm_data['NONBONDED_PARM_INDEX'][i*ntypes+j] - 1
if idx >= 0: continue
# It was negative, so we should have ADDED 1 to adjust for
# indexing from 0
idx = -idx - 2
a = self.parm_data['HBOND_ACOEF'][idx]
b = self.parm_data['HBOND_BCOEF'][idx]
if a == 0 and b == 0: continue
indices_with_1012.append((i, j))
if not indices_with_1012: return False
# Now make sure that some of the atoms *have* those indices
active_indices = set()
for atom in self.atoms: active_indices.add(atom.nb_idx-1)
for i, j in indices_with_1012:
if i in active_indices and j in active_indices:
return True
return False
#===================================================
[docs] def load_rst7(self, rst7):
""" Loads coordinates into the AmberParm class
Parameters
----------
rst7 : str or :class:`Rst7`
The Amber coordinate file (either ASCII restart or NetCDF restart)
object or filename to assign atomic coordinates from.
"""
if not hasattr(rst7, 'coordinates'):
rst7 = Rst7.open(rst7)
self.coordinates = rst7.coordinates
self.hasvels = rst7.hasvels
self.box = _copy.copy(rst7.box)
self.hasbox = self.box is not None
if self.hasvels:
self.velocities = rst7.vels
#===================================================
[docs] @needs_openmm
def omm_nonbonded_force(self, nonbondedMethod=None,
nonbondedCutoff=8*u.angstroms,
switchDistance=0*u.angstroms,
ewaldErrorTolerance=0.0005,
reactionFieldDielectric=78.5):
"""
Creates the OpenMM NonbondedForce (and CustomNonbondedForce if
necessary) to define the nonbonded interatomic potential for this
system. A CustomNonbondedForce is used for the r^-4 part of the 12-6-4
Lennard-Jones potential as well as any modified off-diagonal (i.e.,
NBFIX) terms
Parameters
----------
nonbondedMethod : cutoff method
This is the cutoff method. It can be either the NoCutoff,
CutoffNonPeriodic, CutoffPeriodic, PME, or Ewald objects from the
simtk.openmm.app namespace
nonbondedCutoff : float or distance Quantity
The nonbonded cutoff must be either a floating point number
(interpreted as nanometers) or a Quantity with attached units. This
is ignored if nonbondedMethod is NoCutoff.
switchDistance : float or distance Quantity
The distance at which the switching function is turned on for van
der Waals interactions. This is ignored when no cutoff is used, and
no switch is used if switchDistance is 0, negative, or greater than
the cutoff
ewaldErrorTolerance : float=0.0005
When using PME or Ewald, the Ewald parameters will be calculated
from this value
reactionFieldDielectric : float=78.5
If the nonbondedMethod is CutoffPeriodic or CutoffNonPeriodic, the
region beyond the cutoff is treated using a reaction field method
with this dielectric constant. It should be set to 1 if another
implicit solvent model is being used (e.g., GB)
Returns
-------
NonbondedForce [, CustomNonbondedForce]
If a CustomNonbondedForce is necessary, the return value is a
2-element tuple of NonbondedForce, CustomNonbondedForce. If only a
NonbondedForce is necessary, that is the return value
"""
if not self.atoms: return None
nonbfrc = super(AmberParm, self).omm_nonbonded_force(
nonbondedMethod, nonbondedCutoff, switchDistance,
ewaldErrorTolerance, reactionFieldDielectric
)
has1012 = self.has_1012()
hasnbfix = self.has_NBFIX()
has1264 = 'LENNARD_JONES_CCOEF' in self.flag_list
if not hasnbfix and not has1264 and not has1012:
if self.chamber:
self._modify_nonb_exceptions(nonbfrc, None)
return nonbfrc
# If we have NBFIX, omm_nonbonded_force returned a tuple
if hasnbfix: nonbfrc = nonbfrc[0]
# If we have a 10-12 potential, toggle hasnbfix since the 12-6 terms
# need to be handled via a lookup table (to avoid counting *both* 10-12
# and 12-6 terms for the same pairs). Do this now, since nonbfrc is only
# a tuple if hasnbfix is True *without* 10-12 detection.
hasnbfix = hasnbfix or has1012
# We need a CustomNonbondedForce... determine what it needs to calculate
if hasnbfix and has1264:
if has1012:
force = mm.CustomNonbondedForce(
'(a/r6)^2-b/r6-c/r4+(ah/r6)^2-bh/(r6*r4); r6=r4*r2;'
'r4=r2^2; r2=r^2;'
'a=acoef(type1, type2);'
'b=bcoef(type1, type2);'
'c=ccoef(type1, type2);'
'ah=ahcoef(type1, type2);'
'bh=bhcoef(type1, type2);'
)
else:
force = mm.CustomNonbondedForce('(a/r6)^2-b/r6-c/r4; r6=r4*r2;'
'r4=r2^2; r2=r^2;'
'a=acoef(type1, type2);'
'b=bcoef(type1, type2);'
'c=ccoef(type1, type2);')
elif hasnbfix:
if has1012:
force = mm.CustomNonbondedForce(
'(a/r6)^2-b/r6+(ah/r6)^2-bh/(r6*r4);'
'r6=r4*r2;r4=r2^2;r2=r^2;'
'a=acoef(type1, type2);'
'b=bcoef(type1, type2);'
'ah=ahcoef(type1, type2);'
'bh=bhcoef(type1, type2);'
)
else:
force = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r2*r2*r2;'
'r2=r^2; a=acoef(type1, type2);'
'b=bcoef(type1, type2);')
elif has1264:
force = mm.CustomNonbondedForce('-c/r^4;c=ccoef(type1, type2);')
# Set up the force with all of the particles
force.addPerParticleParameter('type')
force.setForceGroup(self.NONBONDED_FORCE_GROUP)
for atom in self.atoms: force.addParticle([atom.nb_idx-1])
# Now construct the lookup tables
ene_conv = u.kilocalories.conversion_factor_to(u.kilojoules)
length_conv = u.angstroms.conversion_factor_to(u.nanometers)
ntypes = self.parm_data['POINTERS'][NTYPES]
if hasnbfix:
acoef = [0 for i in range(ntypes*ntypes)]
parm_acoef = self.parm_data['LENNARD_JONES_ACOEF']
bcoef = acoef[:]
parm_bcoef = self.parm_data['LENNARD_JONES_BCOEF']
if has1264:
ccoef = acoef[:]
parm_ccoef = self.parm_data['LENNARD_JONES_CCOEF']
if has1012:
ahcoef = acoef[:]
bhcoef = acoef[:]
parm_ahcoef = self.parm_data['HBOND_ACOEF']
parm_bhcoef = self.parm_data['HBOND_BCOEF']
afac = sqrt(ene_conv) * length_conv**6
bfac = ene_conv * length_conv**6
cfac = ene_conv * length_conv**4
ahfac = sqrt(ene_conv) * length_conv**6
bhfac = ene_conv * length_conv**10
nbidx = self.parm_data['NONBONDED_PARM_INDEX']
for i in range(ntypes):
for j in range(ntypes):
idx = nbidx[ntypes*i+j] - 1
ii = i + ntypes * j
if idx < 0 and has1012:
idx = -idx - 2 # Fix adjust for 0 since it was negative
ahcoef[ii] = sqrt(parm_ahcoef[idx]) * ahfac
bhcoef[ii] = parm_bhcoef[idx] * bhfac
if has1264:
ccoef[ii] = parm_ccoef[idx] * cfac
elif idx >= 0:
acoef[ii] = sqrt(parm_acoef[idx]) * afac
bcoef[ii] = parm_bcoef[idx] * bfac
if has1264:
ccoef[ii] = parm_ccoef[idx] * cfac
force.addTabulatedFunction('acoef', mm.Discrete2DFunction(ntypes, ntypes, acoef))
force.addTabulatedFunction('bcoef', mm.Discrete2DFunction(ntypes, ntypes, bcoef))
if has1264:
force.addTabulatedFunction('ccoef', mm.Discrete2DFunction(ntypes, ntypes, ccoef))
# Our CustomNonbondedForce is taking care of the LJ part of our
# potential, so we need to go through nonbfrc and zero-out the
# Lennard-Jones parameters, but keep the charge parameters in place.
for i in range(nonbfrc.getNumParticles()):
chg, sig, eps = nonbfrc.getParticleParameters(i)
nonbfrc.setParticleParameters(i, chg, 0.5, 0.0)
# We still let the NonbondedForce handle our nonbonded exceptions,
# but we may have to modify the Lennard-Jones part with
# off-diagonal modifications. Offload this to a private method so it
# can be overridden in the ChamberParm class
self._modify_nonb_exceptions(nonbfrc, force)
elif has1264:
# Here we have JUST the r^-4 or r^-10/r^-12 parts, since the
# hasnbfix block above handled the "hasnbfix and has1264" case
if has1264:
ccoef = [0 for i in range(ntypes*ntypes)]
parm_ccoef = self.parm_data['LENNARD_JONES_CCOEF']
cfac = ene_conv * length_conv**4
ahfac = ene_conv * length_conv**12
bhfac = ene_conv * length_conv**10
nbidx = self.parm_data['NONBONDED_PARM_INDEX']
for i in range(ntypes):
for j in range(ntypes):
idx = nbidx[ntypes*i+j] - 1
if idx < 0: continue
ccoef[i+ntypes*j] = parm_ccoef[idx] * cfac
force.addTabulatedFunction('ccoef', mm.Discrete2DFunction(ntypes, ntypes, ccoef))
# Copy the exclusions
for ii in range(nonbfrc.getNumExceptions()):
i, j, qq, ss, ee = nonbfrc.getExceptionParameters(ii)
force.addExclusion(i, j)
if has1012:
force.addTabulatedFunction('ahcoef', mm.Discrete2DFunction(ntypes, ntypes, ahcoef))
force.addTabulatedFunction('bhcoef', mm.Discrete2DFunction(ntypes, ntypes, bhcoef))
# Copy the switching function information to the CustomNonbondedForce
if nonbfrc.getUseSwitchingFunction():
force.setUseSwitchingFunction(True)
force.setSwitchingDistance(nonbfrc.getSwitchingDistance())
# Set the dispersion correction on (by default)
force.setUseLongRangeCorrection(True)
# Determine which nonbonded method we should use and transfer the
# nonbonded cutoff
assert nonbondedMethod in (app.NoCutoff, app.CutoffNonPeriodic, app.PME, app.Ewald,
app.CutoffPeriodic), 'Bad nonbondedMethod'
if nonbondedMethod is app.NoCutoff:
force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff)
elif nonbondedMethod is app.CutoffNonPeriodic:
force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)
elif nonbondedMethod in (app.PME, app.Ewald, app.CutoffPeriodic):
force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffPeriodic)
force.setCutoffDistance(nonbfrc.getCutoffDistance())
return nonbfrc, force
#===================================================
# Iterators for parameters with and without hydrogen
@property
def bonds_inc_h(self):
""" All bonds including hydrogen """
for bond in self.bonds:
if (bond.atom1.atomic_number == 1 or
bond.atom2.atomic_number == 1):
yield bond
@property
def bonds_without_h(self):
""" All bonds without hydrogen """
for bond in self.bonds:
if (bond.atom1.atomic_number == 1 or
bond.atom2.atomic_number == 1):
continue
yield bond
@property
def angles_inc_h(self):
""" All angles including hydrogen """
for angle in self.angles:
if (angle.atom1.atomic_number == 1 or angle.atom2.atomic_number == 1
or angle.atom3.atomic_number == 1):
yield angle
@property
def angles_without_h(self):
""" All angles including hydrogen """
for angle in self.angles:
if (angle.atom1.atomic_number == 1 or angle.atom2.atomic_number == 1
or angle.atom3.atomic_number == 1):
continue
yield angle
@property
def dihedrals_inc_h(self):
""" All dihedrals including hydrogen """
for dihed in self.dihedrals:
if (dihed.atom1.atomic_number == 1
or dihed.atom2.atomic_number == 1
or dihed.atom3.atomic_number == 1
or dihed.atom4.atomic_number == 1):
yield dihed
@property
def dihedrals_without_h(self):
""" All dihedrals including hydrogen """
for dihed in self.dihedrals:
if (dihed.atom1.atomic_number == 1
or dihed.atom2.atomic_number == 1
or dihed.atom3.atomic_number == 1
or dihed.atom4.atomic_number == 1):
continue
yield dihed
#===================================================
@property
def chamber(self):
""" Whether this instance uses the CHARMM force field """
return False
@property
def amoeba(self):
""" Whether this instance uses the Amoeba force field """
return False
@property
def has_cmap(self):
""" Whether this instance has correction map terms or not """
return len(self.cmaps) > 0 or (self._cmap_prefix + 'CMAP_COUNT') in self.parm_data
#=========== PRIVATE INSTANCE METHODS ============
def _truncate_array(self, section, length):
""" Truncates an array to get the given length """
self.parm_data[section] = self.parm_data[section][:length]
#===================================================
def _load_cmap_info(self):
""" Loads the CHARMM CMAP types and array """
if not self.has_cmap: return
del self.cmaps[:]
del self.cmap_types[:]
resolution_key = self._cmap_prefix + 'CMAP_RESOLUTION'
parameter_key = self._cmap_prefix + 'CMAP_PARAMETER_%02d'
for i in range(self.pointers['CMAP_TYPES']):
resolution = self.parm_data[resolution_key][i]
grid = self.parm_data[parameter_key % (i+1)]
cmts = self.parm_comments[parameter_key % (i+1)]
self.cmap_types.append(CmapType(resolution, grid, cmts, list=self.cmap_types))
it = iter(self.parm_data[self._cmap_prefix + 'CMAP_INDEX'])
for i, j, k, l, m, n in zip(it, it, it, it, it, it):
self.cmaps.append(
Cmap(self.atoms[i-1], self.atoms[j-1], self.atoms[k-1],
self.atoms[l-1], self.atoms[m-1], self.cmap_types[n-1])
)
#===================================================
def _check_section_lengths(self):
"""
Checks that all of the raw sections have the appropriate length as
specified by the POINTER section.
Raises
------
AmberError if any of the lengths are incorrect
"""
def check_length(key, length, required=True):
if not required and key not in self.parm_data:
return
if len(self.parm_data[key]) != length:
raise AmberError('FLAG %s has %d elements; expected %d' %
(key, len(self.parm_data[key]), length))
natom = self.ptr('NATOM')
check_length('ATOM_NAME', natom)
check_length('CHARGE', natom)
check_length('MASS', natom)
check_length('ATOM_TYPE_INDEX', natom)
check_length('NUMBER_EXCLUDED_ATOMS', natom)
check_length('JOIN_ARRAY', natom)
check_length('IROTAT', natom)
check_length('RADIUS', natom, False)
check_length('SCREEN', natom, False)
check_length('ATOMIC_NUMBER', natom, False)
ntypes = self.ptr('NTYPES')
check_length('NONBONDED_PARM_INDEX', ntypes*ntypes)
check_length('LENNARD_JONES_ACOEF', ntypes*(ntypes+1)//2)
check_length('LENNARD_JONES_BCOEF', ntypes*(ntypes+1)//2)
check_length('LENNARD_JONES_CCOEF', ntypes*(ntypes+1)//2, False)
nres = self.ptr('NRES')
check_length('RESIDUE_LABEL', nres)
check_length('RESIDUE_POINTER', nres)
check_length('RESIDUE_CHAINID', nres, False)
check_length('RESIDUE_ICODE', nres, False)
check_length('RESIDUE_NUMBER', nres, False)
check_length('BOND_FORCE_CONSTANT', self.ptr('NUMBND'))
check_length('BOND_EQUIL_VALUE', self.ptr('NUMBND'))
check_length('ANGLE_FORCE_CONSTANT', self.ptr('NUMANG'))
check_length('ANGLE_EQUIL_VALUE', self.ptr('NUMANG'))
check_length('DIHEDRAL_FORCE_CONSTANT', self.ptr('NPTRA'))
check_length('DIHEDRAL_PERIODICITY', self.ptr('NPTRA'))
check_length('DIHEDRAL_PHASE', self.ptr('NPTRA'))
check_length('SCEE_SCALE_FACTOR', self.ptr('NPTRA'), False)
check_length('SCNB_SCALE_FACTOR', self.ptr('NPTRA'), False)
check_length('SOLTY', self.ptr('NATYP'))
check_length('BONDS_INC_HYDROGEN', self.ptr('NBONH')*3)
check_length('BONDS_WITHOUT_HYDROGEN', self.ptr('MBONA')*3)
check_length('ANGLES_INC_HYDROGEN', self.ptr('NTHETH')*4)
check_length('ANGLES_WITHOUT_HYDROGEN', self.ptr('NTHETA')*4)
check_length('DIHEDRALS_INC_HYDROGEN', self.ptr('NPHIH')*5)
check_length('DIHEDRALS_WITHOUT_HYDROGEN', self.ptr('NPHIA')*5)
check_length('HBOND_ACOEF', self.ptr('NPHB'))
check_length('HBOND_BCOEF', self.ptr('NPHB'))
check_length('SOLVENT_POINTERS', 3, False)
if 'SOLVENT_POINTERS' in self.parm_data:
check_length('ATOMS_PER_MOLECULE', self.parm_data['SOLVENT_POINTERS'][1], False)
if self.has_cmap:
check_length(self._cmap_prefix + 'CMAP_COUNT', 2)
check_length(self._cmap_prefix + 'CMAP_RESOLUTION', self.pointers['CMAP_TYPES'])
resolution_key = self._cmap_prefix + 'CMAP_RESOLUTION'
parameter_key = self._cmap_prefix + 'CMAP_PARAMETER_%02d'
for i in range(self.pointers['CMAP_TYPES']):
res = self.parm_data[resolution_key][i]
check_length(parameter_key % (i+1), res*res)
#===================================================
def _load_atoms_and_residues(self):
"""
Loads the atoms and residues (which are always done together) into the
data structure
"""
del self.residues[:]
del self.atoms[:]
# Figure out on which atoms the residues start and stop
natom = self.parm_data['POINTERS'][NATOM]
res_ptr = self.parm_data['RESIDUE_POINTER'] + [natom+1]
try:
atnums = self.parm_data['ATOMIC_NUMBER']
except KeyError:
atnums = None
try:
res_icd = self.parm_data['RESIDUE_ICODE']
except KeyError:
res_icd = ['' for i in range(self.parm_data['POINTERS'][NRES])]
try:
res_chn = self.parm_data['RESIDUE_CHAINID']
except KeyError:
res_chn = ['' for i in range(self.parm_data['POINTERS'][NRES])]
for i, resname in enumerate(self.parm_data['RESIDUE_LABEL']):
resstart = res_ptr[i] - 1
resend = res_ptr[i+1] - 1
for j in range(resstart, resend):
if atnums is None:
if self.parm_data['AMBER_ATOM_TYPE'][j] in ('EP', 'LP'):
atom = ExtraPoint()
else:
atom = Atom()
elif atnums[j] == 0:
atom = ExtraPoint()
else:
atom = Atom()
self.add_atom(atom, resname, i, res_chn[i], res_icd[i])
# Residue number may not be unique, since zeros may be assigned to all
# solvent residues by addPDB. So we can't use the RESIDUE_NUMBER section
# to fill in the residue sequence IDs in the add_atom call above. Add it
# in as a post-hoc addition now if that information is present
if 'RESIDUE_NUMBER' in self.parm_data:
for res, num in zip(self.residues, self.parm_data['RESIDUE_NUMBER']):
res.number = num
#===================================================
def _load_extra_exclusions(self):
"""
Look through the exclusion list in the prmtop file and see if any
_additional_ exclusions outside the basic ones defined for bonds,
angles, and dihedrals are specified. If so, add those to the exclusion
list.
This also goes through all atoms and loads the proper bond, angle and
dihedral partners into any extra points. The way extra points are
handled is that they are considered to carry the same topological
connectivity as they actual atom they are bonded to.
"""
num_excluded = self.parm_data['NUMBER_EXCLUDED_ATOMS']
excluded_list = self.parm_data['EXCLUDED_ATOMS_LIST']
first_excl = 0
for i, atom in enumerate(self.atoms):
exclusions = set()
bond_excl = set(atom._bond_partners + atom._angle_partners +
atom._dihedral_partners + atom._tortor_partners +
atom._exclusion_partners)
nexcl = num_excluded[i]
for j in range(first_excl, first_excl+nexcl):
if excluded_list[j]: # Zeros are placeholders
exclusions.add(self.atoms[excluded_list[j]-1])
for eatom in exclusions - bond_excl:
atom.exclude(eatom)
first_excl += nexcl
#===================================================
def _load_bond_info(self):
""" Loads the bond types and bond arrays """
del self.bond_types[:]
del self.bonds[:]
for k, req in zip(self.parm_data['BOND_FORCE_CONSTANT'],
self.parm_data['BOND_EQUIL_VALUE']):
self.bond_types.append(BondType(k, req, self.bond_types))
it = iter(self.parm_data['BONDS_WITHOUT_HYDROGEN'])
for i, j, k in zip(it, it, it):
self.bonds.append(
Bond(self.atoms[i//3], self.atoms[j//3],
self.bond_types[k-1])
)
it = iter(self.parm_data['BONDS_INC_HYDROGEN'])
for i, j, k in zip(it, it, it):
self.bonds.append(
Bond(self.atoms[i//3], self.atoms[j//3],
self.bond_types[k-1])
)
#===================================================
def _load_angle_info(self):
""" Loads the angle types and angle arrays """
del self.angle_types[:]
del self.angles[:]
for k, theteq in zip(self.parm_data['ANGLE_FORCE_CONSTANT'],
self.parm_data['ANGLE_EQUIL_VALUE']):
theteq *= RAD_TO_DEG
self.angle_types.append(AngleType(k, theteq, self.angle_types))
it = iter(self.parm_data['ANGLES_WITHOUT_HYDROGEN'])
for i, j, k, l in zip(it, it, it, it):
self.angles.append(
Angle(self.atoms[i//3], self.atoms[j//3], self.atoms[k//3],
self.angle_types[l-1])
)
it = iter(self.parm_data['ANGLES_INC_HYDROGEN'])
for i, j, k, l in zip(it, it, it, it):
self.angles.append(
Angle(self.atoms[i//3], self.atoms[j//3], self.atoms[k//3],
self.angle_types[l-1])
)
#===================================================
def _load_dihedral_info(self):
""" Loads the dihedral types and dihedral arrays """
del self.dihedral_types[:]
del self.dihedrals[:]
try:
scee = self.parm_data['SCEE_SCALE_FACTOR']
except KeyError:
scee = [1.2 for i in self.parm_data['DIHEDRAL_FORCE_CONSTANT']]
try:
scnb = self.parm_data['SCNB_SCALE_FACTOR']
except KeyError:
scnb = [2.0 for i in self.parm_data['DIHEDRAL_FORCE_CONSTANT']]
for k, per, ph, e, n in zip(self.parm_data['DIHEDRAL_FORCE_CONSTANT'],
self.parm_data['DIHEDRAL_PERIODICITY'],
self.parm_data['DIHEDRAL_PHASE'],
scee, scnb):
ph *= RAD_TO_DEG
self.dihedral_types.append(
DihedralType(k, per, ph, e, n, list=self.dihedral_types)
)
it = iter(self.parm_data['DIHEDRALS_WITHOUT_HYDROGEN'])
for i, j, k, l, m in zip(it, it, it, it, it):
ignore_end = k < 0
improper = l < 0
self.dihedrals.append(
Dihedral(self.atoms[i//3], self.atoms[j//3],
self.atoms[abs(k)//3], self.atoms[abs(l)//3],
improper=improper, ignore_end=ignore_end,
type=self.dihedral_types[m-1])
)
it = iter(self.parm_data['DIHEDRALS_INC_HYDROGEN'])
for i, j, k, l, m in zip(it, it, it, it, it):
ignore_end = k < 0
improper = l < 0
self.dihedrals.append(
Dihedral(self.atoms[i//3], self.atoms[j//3],
self.atoms[abs(k)//3], self.atoms[abs(l)//3],
improper=improper, ignore_end=ignore_end,
type=self.dihedral_types[m-1])
)
#===================================================
def _xfer_atom_info(self):
"""
Sets the various topology file section data from the `atoms` list to the
topology file data in `parm_data`
"""
natom = len(self.atoms)
data = self.parm_data
data['POINTERS'][NATOM] = natom
self.pointers['NATOM'] = natom
data['ATOM_NAME'] = [atom.name[:4] for atom in self.atoms]
data['AMBER_ATOM_TYPE'] = [atom.type[:4] for atom in self.atoms]
data['CHARGE'] = [atom.charge for atom in self.atoms]
data['MASS'] = [atom.mass for atom in self.atoms]
data['ATOM_TYPE_INDEX'] = [atom.nb_idx for atom in self.atoms]
data['JOIN_ARRAY'] = [atom.join for atom in self.atoms]
data['TREE_CHAIN_CLASSIFICATION'] = [atom.tree[:4] for atom in self.atoms]
data['IROTAT'] = [atom.irotat for atom in self.atoms]
data['NUMBER_EXCLUDED_ATOMS'] = [0 for atom in self.atoms]
if 'RADII' in data:
data['RADII'] = [atom.solvent_radius for atom in self.atoms]
if 'SCREEN' in data:
data['SCREEN'] = [atom.screen for atom in self.atoms]
if 'ATOMIC_NUMBER' in data:
data['ATOMIC_NUMBER'] = [atom.atomic_number for atom in self.atoms]
# Do the non-bonded exclusions now
data['EXCLUDED_ATOMS_LIST'] = []
nextra = 0
max_typ = 0
for i, atom in enumerate(self.atoms):
excl = atom.nonbonded_exclusions(index_from=1)
if len(excl) == 0:
excl = [0]
data['EXCLUDED_ATOMS_LIST'] += excl
data['NUMBER_EXCLUDED_ATOMS'][i] = len(excl)
if atom.atomic_number == 0:
nextra += 1
max_typ = max(max_typ, atom.nb_idx)
nnb = len(data['EXCLUDED_ATOMS_LIST'])
data['POINTERS'][NNB] = nnb
self.pointers['NNB'] = self.pointers['NEXT'] = nnb
data['POINTERS'][NUMEXTRA] = nextra
self.pointers['NUMEXTRA'] = nextra
max_typ = max(data['POINTERS'][NTYPES], max_typ)
data['POINTERS'][NTYPES] = max_typ
self.pointers['NTYPES'] = max_typ
#===================================================
def _xfer_residue_info(self):
"""
Sets the various topology file section data from the `residues` list to
the topology file data in `parm_data`
"""
data = self.parm_data
nres = len(self.residues)
data['POINTERS'][NRES] = nres
self.pointers['NRES'] = nres
data['RESIDUE_LABEL'] = [r.name[:4] for r in self.residues]
data['RESIDUE_POINTER'] = [r.atoms[0].idx+1 for r in self.residues]
if 'RESIDUE_NUMBER' in data:
data['RESIDUE_NUMBER'] = [r.number for r in self.residues]
if 'RESIDUE_CHAINID' in data:
data['RESIDUE_CHAINID'] = [res.chain for res in self.residues]
if 'RESIDUE_ICODE' in data:
data['RESIDUE_ICODE'] = [r.insertion_code for r in self.residues]
nmxrs = max([len(res) for res in self.residues]) if self.residues else 0
data['POINTERS'][NMXRS] = nmxrs
self.pointers['NMXRS'] = nmxrs
#===================================================
def _xfer_bond_info(self):
"""
Sets the data for the various bond arrays in the raw data from the
parameter lists
"""
# First do the bond types
data = self.parm_data
for bond_type in self.bond_types:
bond_type.used = False
for bond in self.bonds:
bond.type.used = True
self.bond_types.prune_unused()
data['BOND_FORCE_CONSTANT'] = [type.k for type in self.bond_types]
data['BOND_EQUIL_VALUE'] = [type.req for type in self.bond_types]
data['POINTERS'][NUMBND] = len(self.bond_types)
self.pointers['NUMBND'] = len(self.bond_types)
# Now do the bond arrays
data['BONDS_INC_HYDROGEN'] = bond_array = []
bond_list = list(self.bonds_inc_h)
for bond in bond_list:
bond_array.extend([bond.atom1.idx*3, bond.atom2.idx*3, bond.type.idx+1])
data['POINTERS'][NBONH] = len(bond_list)
self.pointers['NBONH'] = len(bond_list)
data['BONDS_WITHOUT_HYDROGEN'] = bond_array = []
bond_list = list(self.bonds_without_h)
for bond in bond_list:
bond_array.extend([bond.atom1.idx*3, bond.atom2.idx*3, bond.type.idx+1])
data['POINTERS'][MBONA] = data['POINTERS'][NBONA] = len(bond_list)
self.pointers['MBONA'] = self.pointers['NBONA'] = len(bond_list)
#===================================================
def _xfer_angle_info(self):
"""
Sets the data for the various angle arrays in the raw data from the
parameter lists
"""
# First do the angle types
data = self.parm_data
for angle_type in self.angle_types:
angle_type.used = False
for angle in self.angles:
angle.type.used = True
self.angle_types.prune_unused()
data['ANGLE_FORCE_CONSTANT'] = [type.k for type in self.angle_types]
data['ANGLE_EQUIL_VALUE'] = [type.theteq*DEG_TO_RAD for type in self.angle_types]
data['POINTERS'][NUMANG] = len(self.angle_types)
self.pointers['NUMANG'] = len(self.angle_types)
# Now do the angle arrays
data['ANGLES_INC_HYDROGEN'] = angle_array = []
angle_list = list(self.angles_inc_h)
for angle in angle_list:
angle_array.extend([angle.atom1.idx*3, angle.atom2.idx*3,
angle.atom3.idx*3, angle.type.idx+1])
data['POINTERS'][NTHETH] = len(angle_list)
self.pointers['NTHETH'] = len(angle_list)
data['ANGLES_WITHOUT_HYDROGEN'] = angle_array = []
angle_list = list(self.angles_without_h)
for angle in angle_list:
angle_array.extend([angle.atom1.idx*3, angle.atom2.idx*3,
angle.atom3.idx*3, angle.type.idx+1])
data['POINTERS'][NTHETA] = data['POINTERS'][MTHETA] = len(angle_list)
self.pointers['NTHETA'] = self.pointers['MTHETA'] = len(angle_list)
#===================================================
def _xfer_dihedral_info(self):
"""
Sets the data for the various dihedral arrays in the raw data from the
parameter lists
"""
# First do the dihedral types
data = self.parm_data
for dihedral_type in self.dihedral_types:
dihedral_type.used = False
for dihed in self.dihedrals:
dihed.type.used = True
self.dihedral_types.prune_unused()
data['DIHEDRAL_FORCE_CONSTANT'] = [type.phi_k for type in self.dihedral_types]
data['DIHEDRAL_PERIODICITY'] = [type.per for type in self.dihedral_types]
data['DIHEDRAL_PHASE'] = [type.phase*DEG_TO_RAD for type in self.dihedral_types]
if 'SCEE_SCALE_FACTOR' in data:
data['SCEE_SCALE_FACTOR'] = [type.scee for type in self.dihedral_types]
if 'SCNB_SCALE_FACTOR' in data:
data['SCNB_SCALE_FACTOR'] = [type.scnb for type in self.dihedral_types]
data['POINTERS'][NPTRA] = len(self.dihedral_types)
self.pointers['NPTRA'] = len(self.dihedral_types)
# Now do the dihedral arrays
data['DIHEDRALS_INC_HYDROGEN'] = dihed_array = []
dihed_list = list(self.dihedrals_inc_h)
for dihed in dihed_list:
imp_sign = -1 if dihed.improper else 1
end_sign = -1 if dihed.ignore_end else 1
if dihed.atom3.idx == 0 or dihed.atom4.idx == 0:
dihed_array.extend([dihed.atom4.idx*3, dihed.atom3.idx*3,
dihed.atom2.idx*3*end_sign, dihed.atom1.idx*3*imp_sign,
dihed.type.idx+1])
else:
dihed_array.extend([dihed.atom1.idx*3, dihed.atom2.idx*3,
dihed.atom3.idx*3*end_sign, dihed.atom4.idx*3*imp_sign,
dihed.type.idx+1])
data['POINTERS'][NPHIH] = len(dihed_list)
self.pointers['NPHIH'] = len(dihed_list)
data['DIHEDRALS_WITHOUT_HYDROGEN'] = dihed_array = []
dihed_list = list(self.dihedrals_without_h)
for dihed in dihed_list:
imp_sign = -1 if dihed.improper else 1
end_sign = -1 if dihed.ignore_end else 1
if dihed.atom3.idx == 0 or dihed.atom4.idx == 0:
dihed_array.extend([dihed.atom4.idx*3, dihed.atom3.idx*3,
dihed.atom2.idx*3*end_sign, dihed.atom1.idx*3*imp_sign,
dihed.type.idx+1])
else:
dihed_array.extend([dihed.atom1.idx*3, dihed.atom2.idx*3,
dihed.atom3.idx*3*end_sign, dihed.atom4.idx*3*imp_sign,
dihed.type.idx+1])
data['POINTERS'][NPHIA] = data['POINTERS'][MPHIA] = len(dihed_list)
self.pointers['NPHIA'] = self.pointers['MPHIA'] = len(dihed_list)
#===================================================
def _xfer_cmap_properties(self):
""" Sets the topology file section data from the cmap arrays """
# If we have no cmaps, delete all remnants of the CMAP terms in the
# prmtop and bail out
if len(self.cmaps) == 0:
# We have deleted all cmaps. Get rid of them from the parm file and
# bail out. This is probably pretty unlikely, though...
flag_prefix = self._cmap_prefix + 'CMAP'
flags_to_delete = [flag for flag in self.flag_list if flag.startswith(flag_prefix)]
for flag in flags_to_delete:
self.delete_flag(flag)
if 'CMAP' in self.pointers:
del self.pointers['CMAP']
if 'CMAP_TYPES' in self.pointers:
del self.pointers['CMAP_TYPES']
return
# Time to transfer our CMAP types
data = self.parm_data
for ct in self.cmap_types:
ct.used = False
for cmap in self.cmaps:
cmap.type.used = True
self.cmap_types.prune_unused()
# All of our CMAP types are in different topology file sections. We need
# to delete all of the CMAP_PARAMETER_XX sections and then
# recreate them with the correct size and comments. The comments have
# been stored in the CMAP types themselves to prevent them from being
# lost. We will also assume that the Fortran format we're going to use
# is the same for all CMAP types, so just pull it from
# CMAP_PARAMETER_01 (or fall back to 8(F9.5))
parameter_key = self._cmap_prefix + 'CMAP_PARAMETER_%02d'
try:
fmt = str(self.formats[parameter_key % 1])
except KeyError:
fmt = '8(F9.5)'
flags_to_delete = []
for flag in self.flag_list:
if 'CMAP_PARAMETER' in flag:
flags_to_delete.append(flag)
for flag in flags_to_delete:
self.delete_flag(flag)
# Now add them back
after = self._cmap_prefix + 'CMAP_RESOLUTION'
for i, ct in enumerate(self.cmap_types):
newflag = self._cmap_prefix + 'CMAP_PARAMETER_%02d' % (i+1)
self.add_flag(newflag, fmt, data=ct.grid, comments=ct.comments, after=after)
after = newflag
# Now do the CMAP_INDEX section
data[self._cmap_prefix + 'CMAP_INDEX'] = cmap_array = []
for cm in self.cmaps:
cmap_array.extend([cm.atom1.idx+1, cm.atom2.idx+1, cm.atom3.idx+1, cm.atom4.idx+1,
cm.atom5.idx+1, cm.type.idx+1])
data[self._cmap_prefix + 'CMAP_COUNT'] = [len(self.cmaps), len(self.cmap_types)]
data[self._cmap_prefix + 'CMAP_RESOLUTION'] = [ct.resolution for ct in self.cmap_types]
self.pointers['CMAP'] = len(self.cmaps)
self.pointers['CMAP_TYPES'] = len(self.cmap_types)
#===================================================
def _add_standard_flags(self):
""" Adds all of the standard flags to the parm_data array """
self.set_version()
self.add_flag('TITLE', '20a4', num_items=0)
self.add_flag('POINTERS', '10I8', num_items=31)
self.add_flag('ATOM_NAME', '20a4', num_items=0)
self.add_flag('CHARGE', '5E16.8', num_items=0)
self.add_flag('ATOMIC_NUMBER', '10I8', num_items=0)
self.add_flag('MASS', '5E16.8', num_items=0)
self.add_flag('ATOM_TYPE_INDEX', '10I8', num_items=0)
self.add_flag('NUMBER_EXCLUDED_ATOMS', '10I8', num_items=0)
self.add_flag('NONBONDED_PARM_INDEX', '10I8', num_items=0)
self.add_flag('RESIDUE_LABEL', '20a4', num_items=0)
self.add_flag('RESIDUE_POINTER', '10I8', num_items=0)
self.add_flag('BOND_FORCE_CONSTANT', '5E16.8', num_items=0)
self.add_flag('BOND_EQUIL_VALUE', '5E16.8', num_items=0)
self.add_flag('ANGLE_FORCE_CONSTANT', '5E16.8', num_items=0)
self.add_flag('ANGLE_EQUIL_VALUE', '5E16.8', num_items=0)
self.add_flag('DIHEDRAL_FORCE_CONSTANT', '5E16.8', num_items=0)
self.add_flag('DIHEDRAL_PERIODICITY', '5E16.8', num_items=0)
self.add_flag('DIHEDRAL_PHASE', '5E16.8', num_items=0)
self.add_flag('SCEE_SCALE_FACTOR', '5E16.8', num_items=0)
self.add_flag('SCNB_SCALE_FACTOR', '5E16.8', num_items=0)
self.pointers['NATYP'] = self.parm_data['POINTERS'][NATYP] = 1
self.add_flag('SOLTY', '5E16.8', num_items=1)
self.add_flag('LENNARD_JONES_ACOEF', '5E16.8', num_items=0)
self.add_flag('LENNARD_JONES_BCOEF', '5E16.8', num_items=0)
self.add_flag('BONDS_INC_HYDROGEN', '10I8', num_items=0)
self.add_flag('BONDS_WITHOUT_HYDROGEN', '10I8', num_items=0)
self.add_flag('ANGLES_INC_HYDROGEN', '10I8', num_items=0)
self.add_flag('ANGLES_WITHOUT_HYDROGEN', '10I8', num_items=0)
self.add_flag('DIHEDRALS_INC_HYDROGEN', '10I8', num_items=0)
self.add_flag('DIHEDRALS_WITHOUT_HYDROGEN', '10I8', num_items=0)
self.add_flag('EXCLUDED_ATOMS_LIST', '10I8', num_items=0)
self.add_flag('HBOND_ACOEF', '5E16.8', num_items=0)
self.add_flag('HBOND_BCOEF', '5E16.8', num_items=0)
self.add_flag('HBCUT', '5E16.8', num_items=0)
self.add_flag('AMBER_ATOM_TYPE', '20a4', num_items=0)
self.add_flag('TREE_CHAIN_CLASSIFICATION', '20a4', num_items=0)
self.add_flag('JOIN_ARRAY', '10I8', num_items=0)
self.add_flag('IROTAT', '10I8', num_items=0)
if self.has_cmap:
self.add_flag(self._cmap_prefix + 'CMAP_COUNT', '2I8', num_items=2,
comments=['Number of CMAP terms, number of unique CMAP parameters'])
self.add_flag(self._cmap_prefix + 'CMAP_RESOLUTION', '20I4', num_items=0,
comments=['Number of steps along each phi/psi CMAP axis',
'for each CMAP_PARAMETER grid'])
self.add_flag(self._cmap_prefix + 'CMAP_INDEX', '6I8', num_items=0,
comments=['Atom index i,j,k,l,m of the cross term',
'and then pointer to CMAP_PARAMETER_n'])
if self.box is not None:
self.add_flag('SOLVENT_POINTERS', '3I8', num_items=3)
self.add_flag('ATOMS_PER_MOLECULE', '10I8', num_items=0)
self.add_flag('BOX_DIMENSIONS', '5E16.8', num_items=4)
self.add_flag('RADIUS_SET', '1a80', num_items=1)
self.add_flag('RADII', '5E16.8', num_items=0)
self.add_flag('SCREEN', '5E16.8', num_items=0)
self.add_flag('IPOL', '1I8', num_items=1)
#===================================================
def _set_nonbonded_tables(self, nbfixes=None):
"""
Sets the tables of Lennard-Jones nonbonded interaction pairs
"""
ntypes = self.parm_data['POINTERS'][NTYPES]
ntypes2 = ntypes * ntypes
# Set up the index lookup tables (not a unique solution)
self.parm_data['NONBONDED_PARM_INDEX'] = [0 for i in range(ntypes2)]
holder = [0 for i in range(ntypes2)]
idx = 0
for i in range(ntypes):
for j in range(i+1):
idx += 1
holder[ntypes*i+j] = holder[ntypes*j+i] = idx
idx = 0
for i in range(ntypes):
for j in range(ntypes):
self.parm_data['NONBONDED_PARM_INDEX'][idx] = holder[ntypes*i+j]
idx += 1
nttyp = ntypes * (ntypes + 1) // 2
# Now build the Lennard-Jones arrays
self.parm_data['LENNARD_JONES_ACOEF'] = [0 for i in range(nttyp)]
self.parm_data['LENNARD_JONES_BCOEF'] = [0 for i in range(nttyp)]
self.recalculate_LJ()
# Now make any NBFIX modifications we had
if nbfixes is not None:
for i, fix in enumerate(nbfixes):
for terms in fix:
j, rmin, eps, rmin14, eps14 = terms
i, j = min(i, j-1), max(i, j-1)
eps = abs(eps)
eps14 = abs(eps14)
idx = self.parm_data['NONBONDED_PARM_INDEX'][ntypes*i+j] - 1
self.parm_data['LENNARD_JONES_ACOEF'][idx] = eps * rmin**12
self.parm_data['LENNARD_JONES_BCOEF'][idx] = 2*eps * rmin**6
#===================================================
@needs_openmm
def _modify_nonb_exceptions(self, nonbfrc, customforce):
"""
Modifies the nonbonded force exceptions and the custom nonbonded force
exclusions. The exceptions on the nonbonded force might need to be
adjusted if off-diagonal modifications on the L-J matrix are present
"""
# To get into this routine, either NBFIX is present OR this is a chamber
# prmtop and we need to pull the 1-4 L-J parameters from the
# LENNARD_JONES_14_A/BCOEF arrays
length_conv = u.angstroms.conversion_factor_to(u.nanometers)
ene_conv = u.kilocalories.conversion_factor_to(u.kilojoules)
atoms = self.atoms
try:
acoef = self.parm_data['LENNARD_JONES_14_ACOEF']
bcoef = self.parm_data['LENNARD_JONES_14_BCOEF']
except KeyError:
acoef = self.parm_data['LENNARD_JONES_ACOEF']
bcoef = self.parm_data['LENNARD_JONES_BCOEF']
nbidx = self.parm_data['NONBONDED_PARM_INDEX']
ntypes = self.parm_data['POINTERS'][NTYPES]
sigma_scale = 2**(-1/6) * length_conv
for ii in range(nonbfrc.getNumExceptions()):
i, j, qq, ss, ee = nonbfrc.getExceptionParameters(ii)
if qq.value_in_unit(u.elementary_charge**2) == 0 and (
ss.value_in_unit(u.angstroms) == 0 or
ee.value_in_unit(u.kilocalories_per_mole) == 0):
# Copy this exclusion as-is... no need to modify the nonbfrc
# exception parameters
if customforce is not None:
customforce.addExclusion(i, j)
continue
# Figure out what the 1-4 scaling parameters were for this pair...
unscaled_ee = sqrt(self.atoms[i].epsilon_14 * self.atoms[j].epsilon_14) * ene_conv
try:
one_scnb = ee.value_in_unit(u.kilojoules_per_mole) / unscaled_ee
except ZeroDivisionError:
one_scnb = 1
id1 = atoms[i].nb_idx - 1
id2 = atoms[j].nb_idx - 1
idx = nbidx[ntypes*id1+id2] - 1
if idx >= 0:
a = acoef[idx]
b = bcoef[idx]
else:
a = b = 0
if b == 0:
epsilon = 0.0
sigma = 0.5
elif idx >= 0:
# b / a == 2 / r^6 --> (a / b * 2)^(1/6) = rmin
rmin = (a / b * 2)**(1/6)
epsilon = b / (2 * rmin**6) * ene_conv * one_scnb
sigma = rmin * sigma_scale
nonbfrc.setExceptionParameters(ii, i, j, qq, sigma, epsilon)
if customforce is not None:
customforce.addExclusion(i, j)
#===================================================
def _add_missing_13_14(self, ignore_inconsistent_vdw=False):
"""
Uses the bond graph to fill in zero-parameter angles and dihedrals. The
reason this is necessary is that Amber assumes that the list of angles
and dihedrals encompasses *all* 1-3 and 1-4 pairs as determined by the
bond graph, respectively. As a result, Amber programs use the angle and
dihedral lists to set nonbonded exclusions and exceptions.
Parameters
----------
ignore_inconsistent_vdw : bool, optional
If True, do not make inconsistent 1-4 vdW parameters fatal. For
ChamberParm, the 1-4 specific vdW parameters can compensate. For
AmberParm, the 1-4 scaling factor cannot represent arbitrary
exceptions. Default is False (should only be True for ChamberParm)
Returns
-------
n13, n14 : int, int
The number of 1-3 and 1-4 pairs that needed to be added,
respectively. Purely diagnostic
"""
# We need to figure out what 1-4 scaling term to use if we don't have
# explicit exceptions
assert self.combining_rule in ('lorentz', 'geometric'), "Unrecognized combining rule"
if not self.adjusts:
scalings = defaultdict(int)
for dih in self.dihedrals:
if dih.ignore_end or dih.improper: continue
scalings[(dih.type.scee, dih.type.scnb)] += 1
if len(scalings) > 0:
maxkey, maxval = next(iteritems(scalings))
for key, val in iteritems(scalings):
if maxval < val:
maxkey, maxval = key, val
scee, scnb = maxkey
else:
scee = scnb = 1e10
zero_torsion = DihedralType(0, 1, 0, scee, scnb)
else:
# Turn list of exceptions into a dict so we can look it up quickly
adjust_dict = dict()
for pair in self.adjusts:
adjust_dict[tuple(sorted([pair.atom1, pair.atom2]))] = pair
ignored_torsion = None
zero_torsion = None
# Scan through existing dihedrals to make sure the exceptions match
# the dihedral list
if self.combining_rule == 'lorentz':
comb_sig = lambda sig1, sig2: 0.5 * (sig1 + sig2)
elif self.combining_rule == 'geometric':
comb_sig = lambda sig1, sig2: sqrt(sig1 * sig2)
fac = 2**(1/6)
for dihedral in self.dihedrals:
if dihedral.ignore_end: continue
key = tuple(sorted([dihedral.atom1, dihedral.atom4]))
eref = sqrt(dihedral.atom1.epsilon_14 * dihedral.atom4.epsilon_14)
rref = comb_sig(dihedral.atom1.sigma_14, dihedral.atom4.sigma_14) * fac
if key in adjust_dict:
pair = adjust_dict[key]
if pair.type.epsilon == 0:
scnb = 1e10
else:
scnb = eref / pair.type.epsilon
if pair.type.chgscale == 0:
scee = 1e10
else:
scee = 1 / pair.type.chgscale
if ignore_inconsistent_vdw:
scnb = 1.0
elif abs(rref - pair.type.rmin) > SMALL and pair.type.epsilon != 0:
raise TypeError('Cannot translate exceptions')
if (abs(scnb - dihedral.type.scnb) < SMALL and
abs(scee - dihedral.type.scee) < SMALL):
continue
else:
scee = scnb = 1e10
newtype = _copy.copy(dihedral.type)
newtype.scee = scee
newtype.scnb = scnb
dihedral.type = newtype
newtype.list = self.dihedral_types
self.dihedral_types.append(newtype)
zero_angle = AngleType(0, 0)
n13 = n14 = 0
if self.combining_rule == 'lorentz':
comb_sig = lambda sig1, sig2: 0.5 * (sig1 + sig2)
elif self.combining_rule == 'geometric':
comb_sig = lambda sig1, sig2: sqrt(sig1 * sig2)
fac = 2**(1/6)
for atom in self.atoms:
if isinstance(atom, ExtraPoint): continue
for batom in atom.bond_partners:
if isinstance(batom, ExtraPoint): continue
for aatom in batom.bond_partners:
if isinstance(aatom, ExtraPoint) or aatom is atom: continue
for datom in aatom.bond_partners:
if isinstance(datom, ExtraPoint): continue
if (datom in atom.angle_partners + atom.bond_partners +
atom.dihedral_partners or datom is atom):
continue
# Add the missing dihedral
if not self.adjusts:
tortype = zero_torsion
if n14 == 0:
tortype.list = self.dihedral_types
self.dihedral_types.append(tortype)
else:
# Figure out what the scale factors must be
key = tuple(sorted([atom, datom]))
if key not in adjust_dict:
if ignored_torsion is None:
ignored_torsion = DihedralType(0, 1, 0, 1e10, 1e10)
self.dihedral_types.append(ignored_torsion)
ignored_torsion.list = self.dihedral_types
tortype = ignored_torsion
elif 0 in (adjust_dict[key].type.epsilon, adjust_dict[key].type.rmin) \
and adjust_dict[key].type.chgscale == 0:
if ignored_torsion is None:
ignored_torsion = DihedralType(0, 1, 0, 1e10, 1e10,
list=self.dihedral_types)
self.dihedral_types.append(ignored_torsion)
tortype = ignored_torsion
else:
pair = adjust_dict[key]
epsilon = pair.type.epsilon
rmin = pair.type.rmin
# Compare it to the 1-4 parameters that are
# already present
eref = sqrt(pair.atom1.epsilon_14 * pair.atom2.epsilon_14)
if pair.type.epsilon == 0:
scnb = 1e10
else:
scnb = eref / epsilon
if pair.type.chgscale == 0:
scee = 1e10
else:
scee = 1 / pair.type.chgscale
rref = comb_sig(pair.atom1.sigma_14, pair.atom2.sigma_14) * fac
if abs(rmin - rref) > SMALL:
if ignore_inconsistent_vdw:
scnb = 1.0
else:
raise TypeError('Cannot translate exceptions')
tortype = DihedralType(0, 1, 0, scee, scnb,
list=self.dihedral_types)
self.dihedral_types.append(tortype)
dihedral = Dihedral(atom, batom, aatom, datom, ignore_end=False,
improper=False, type=tortype)
self.dihedrals.append(dihedral)
n14 += 1
if aatom in atom.angle_partners + atom.bond_partners:
continue
# Add the missing angle
self.angles.append(Angle(atom, batom, aatom, zero_angle))
n13 += 1
if n13:
self.angle_types.append(zero_angle)
zero_angle.list = self.angle_types
if n14: # See if there is some ambiguity here
if not self.adjusts and len(scalings) > 1:
warn('Multiple 1-4 scaling factors detected. Using the most-used values scee=%f '
'scnb=%f' % (scee, scnb), AmberWarning)
return n13, n14
#===================================================
def _get_atom_collection_for_alternate_labels(self):
atom_collection = [defaultdict(list) for r in self.residues]
for adict, residue in zip(atom_collection, self.residues):
for atom in residue.atoms:
adict[atom.name].append(atom)
return atom_collection
def _label_alternates(self):
atom_collection = self._get_atom_collection_for_alternate_labels()
possible_labels = list('ABCDEFGHIJKLMNOPQRSTUVWXYZ')
for _, adict in enumerate(atom_collection):
for atom_name, atom_list in iteritems(adict):
if len(atom_list) > 1:
for i, atom in enumerate(atom_list):
label = possible_labels[i%len(possible_labels)]
atom.altloc = label
#===================================================
@property
def box(self):
if self._box is not None:
return self._box[0]
return None
@box.setter
def box(self, value):
# Deleting the box is more complicated for AmberParm, so override it
# here
if value is None:
# Delete all of the other box info in the prmtop
self._box = None
for flag in ('IPTRES', 'NSPM', 'NSPSOL'):
if flag in self.pointers:
del self.pointers[flag]
for flag in ('SOLVENT_POINTERS', 'ATOMS_PER_MOLECULE', 'BOX_DIMENSIONS'):
self.delete_flag(flag)
self.hasbox = False
else:
if isinstance(value, np.ndarray):
box = value
else:
box = list(value)
if len(box) != 6:
raise ValueError('Box information must be 6 floats')
if u.is_quantity(box[0]):
box[0] = box[0].value_in_unit(u.angstroms)
if u.is_quantity(box[1]):
box[1] = box[1].value_in_unit(u.angstroms)
if u.is_quantity(box[2]):
box[2] = box[2].value_in_unit(u.angstroms)
if u.is_quantity(box[3]):
box[3] = box[3].value_in_unit(u.degrees)
if u.is_quantity(box[4]):
box[4] = box[4].value_in_unit(u.degrees)
if u.is_quantity(box[5]):
box[5] = box[5].value_in_unit(u.degrees)
box = np.array(box, dtype=np.float64, copy=False, subok=True).reshape((-1, 6))
# We are adding a box for the first time, so make sure we add some flags
if self._box is None:
self._box = box
# We need to add topology information
if 'SOLVENT_POINTERS' not in self.flag_list:
self.add_flag('SOLVENT_POINTERS', '3I8', num_items=3, after='IROTAT')
if 'ATOMS_PER_MOLECULE' not in self.flag_list:
self.add_flag('ATOMS_PER_MOLECULE', '10I8', data=[0], after='SOLVENT_POINTERS')
if 'BOX_DIMENSIONS' not in self.flag_list:
self.add_flag('BOX_DIMENSIONS', '5E16.8', after='ATOMS_PER_MOLECULE',
data=[box[0,3], box[0,0], box[0,1], box[0,2]])
try:
self.rediscover_molecules(fix_broken=False)
except MoleculeError:
# Do not reorder molecules here -- only do that when
# specifically requested. Otherwise we could get out-of-sync
# with coordinates.
pass
self.load_pointers()
else:
self.parm_data['BOX_DIMENSIONS'] = [box[0,3], box[0,0], box[0,1], box[0,2]]
self._box = box
self._set_ifbox()
def _set_ifbox(self):
""" Sets the IFBOX pointers to 1 (ortho), 2 (octahedral) , or 3 (other) """
if self.box is None:
self.parm_data['POINTERS'][IFBOX] = self.pointers['IFBOX'] = 0
elif np.allclose(self.box[3:], 90):
self.parm_data['POINTERS'][IFBOX] = self.pointers['IFBOX'] = 1
elif np.allclose(self.box[3:], TRUNCATED_OCTAHEDRON_ANGLE, atol=0.02):
self.parm_data['POINTERS'][IFBOX] = self.pointers['IFBOX'] = 2
else:
# General triclinic
self.parm_data['POINTERS'][IFBOX] = self.pointers['IFBOX'] = 3
def _cleanup_dihedrals_with_periodicity_zero(self):
"""
For torsions with only a single term and a periodicity set to 0, make sure pmemd still
properly recognizes the necessary exception parameters. update_dihedral_exclusions will
make sure that if a dihedral has a type pn0 *and* ignore_end is set to False (which means
that it is required to specify exclusions), then it is the *only* torsion between those
atoms in the system. This allows us to scan through our dihedrals, look for significant
terms that have pn==0, and simply add another dihedral with pn=1 and k=0 to ensure that
pmemd will always get that exception correct
"""
new_dihedrals = []
for dih in self.dihedrals:
if dih.ignore_end or dih.type.per != 0:
continue
# If we got here, ignore_end must be False and out periodicity must be 0. So add
# another dihedral
dt = DihedralType(0, 1, 0, dih.type.scee, dih.type.scnb, list=self.dihedral_types)
self.dihedral_types.append(dt)
new_dihedrals.append(
Dihedral(dih.atom1, dih.atom2, dih.atom3, dih.atom4, improper=dih.improper,
ignore_end=False, type=dt)
)
# Now that we added the above dihedral, we can start ignoring the end-group interactions
# on this dihedral
dih.ignore_end = True
if new_dihedrals:
self.dihedrals.extend(new_dihedrals)
#===================================================
_AMBERPARM_ATTRS = 'LJ_types LJ_radius LJ_depth parm_data pointers'.split()
def __getstate__(self):
d = Structure.__getstate__(self)
d.update(AmberFormat.__getstate__(self))
for attr in self._AMBERPARM_ATTRS:
if getattr(self, attr, None) is not None:
d[attr] = getattr(self, attr)
return d
def __setstate__(self, d):
AmberFormat.__setstate__(self, d)
Structure.__setstate__(self, d)
for attr in self._AMBERPARM_ATTRS:
if attr in d:
setattr(self, attr, d[attr])
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[docs]class Rst7(object):
"""
Amber input coordinate (or restart coordinate) file. Front-end for the
readers and writers, supports both NetCDF and ASCII restarts.
Parameters
----------
filename : str, optional
If a filename is provided, this file is parsed and the Rst7 data
populated from that file. The format (ASCII or NetCDF) is autodetected
natom : int, optional
If no filename is provided, this value is required. If a filename is
provided, this value is ignored (and instead set to the value of natom
from the coordinate file). This is the number of atoms for which we have
coordinates. If not provided for a new file, it *must* be set later.
title : str, optional
For a file that is to be written, this is the title that will be given
to that file. Default is an empty string
time : float, optional
The time to write to the restart file. This is cosmetic. Default is 0
"""
[docs] def __init__(self, filename=None, natom=None, title='', time=0.0):
"""
Optionally takes a filename to read. This is deprecated, though, as the
alternative constructor "open" should be used instead
"""
self.coordinates = []
self.vels = None
self._box = None
self.natom = natom
self.title = title
self.time = 0
if filename is not None:
self.filename = filename
self._read(filename)
@property
def box(self):
return self._box
@box.setter
def box(self, value):
if value is None:
self._box = None
else:
self._box = np.array(value).reshape((-1, 6))[0]
[docs] @classmethod
def open(cls, filename):
""" Constructor that opens and parses an input coordinate file
Parameters
----------
filename : str
Name of the file to parse
"""
inst = cls()
inst.filename = filename
inst._read(filename)
return inst
def _read(self, filename):
"""
Open and parse an input coordinate file in either ASCII or NetCDF format
"""
try:
f = AmberAsciiRestart(filename, 'r')
self.natom = f.natom
except ValueError:
# Maybe it's a NetCDF file?
try:
f = NetCDFRestart.open_old(filename)
self.natom = f.atom
except (TypeError, RuntimeError):
raise AmberError('Could not parse restart file %s' % filename)
self.coordinates = f.coordinates
if f.hasvels:
self.vels = f.velocities
if f.hasbox:
self.box = f.box
self.title = f.title
self.time = f.time
[docs] @classmethod
def copy_from(cls, thing):
"""
Copies the coordinates, velocities, and box information from another
instance
"""
inst = cls()
inst.natom = thing.natom
inst.title = thing.title
inst.coordinates = thing.coordinates[:]
if hasattr(thing, 'vels'): inst.vels = _copy.deepcopy(thing.vels)
if hasattr(thing, 'box'): inst.box = _copy.deepcopy(thing.box)
inst.time = thing.time
return inst
def __copy__(self):
""" Copy constructor """
return type(self).copy_from(self)
[docs] def write(self, fname, netcdf=False):
""" Writes the coordinates and/or velocities to a restart file """
if netcdf:
if self.natom is None:
raise RuntimeError('Number of atoms must be set for NetCDF '
'Restart files before write time')
f = NetCDFRestart.open_new(fname, self.natom, self.box is not None,
self.vels is not None, self.title)
else:
f = AmberAsciiRestart(fname, 'w', natom=self.natom,
title=self.title)
f.time = self.time
# Now write the coordinates
f.coordinates = self.coordinates
if self.vels is not None:
f.velocities = self.vels
if self.box is not None:
f.box = self.box
f.close()
@property
def positions(self):
""" Atomic coordinates with units """
coordinates = self.coordinates.reshape(self.natom, 3)
return [Vec3(*x) for x in coordinates] * u.angstroms
@property
def velocities(self):
""" Atomic velocities in units of angstroms/picoseconds """
return np.array(self.vels, copy=False).reshape(self.natom, 3)
@property
def box_vectors(self):
""" Unit cell vectors with units """
if self.box is None: return None
return box_lengths_and_angles_to_vectors(*self.box)
@property
def hasbox(self):
""" Whether or not this Rst7 has unit cell information """
return self.box is not None
@property
def hasvels(self):
""" Whether or not this Rst7 has velocities """
return self.vels is not None
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def set_molecules(parm):
"""
Correctly sets the ATOMS_PER_MOLECULE and SOLVENT_POINTERS sections of the
topology file.
"""
from sys import setrecursionlimit, getrecursionlimit
# Since we use a recursive function here, we make sure that the recursion
# limit is large enough to handle the maximum possible recursion depth we'll
# need (NATOM). We don't want to shrink it, though, since we use list
# comprehensions in list constructors in some places that have an implicit
# (shallow) recursion, therefore, reducing the recursion limit too much here
# could raise a recursion depth exceeded exception during a _Type/Atom/XList
# creation. Therefore, set the recursion limit to the greater of the current
# limit or the number of atoms
setrecursionlimit(max(len(parm.atoms), getrecursionlimit()))
# Unmark all atoms so we can track which molecule each goes into
parm.atoms.unmark()
if not parm.ptr('ifbox'):
raise MoleculeError('Only periodic prmtops can have '
'Molecule definitions')
# The molecule "ownership" list
owner = []
# The way I do this is via a recursive algorithm, in which
# the "set_owner" method is called for each bonded partner an atom
# has, which in turn calls set_owner for each of its partners and
# so on until everything has been assigned.
molecule_number = 1 # which molecule number we are on
for i, atom in enumerate(parm.atoms):
# If this atom has not yet been "owned", make it the next molecule
# However, we only increment which molecule number we're on if
# we actually assigned a new molecule (obviously)
if not atom.marked:
tmp = set()
tmp.add(i)
_set_owner(parm, tmp, i, molecule_number)
# Make sure the atom indexes are sorted
owner.append(tmp)
molecule_number += 1
return owner
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _set_owner(parm, owner_array, atm, mol_id):
""" Recursively sets ownership of given atom and all bonded partners """
parm.atoms[atm].marked = mol_id
for partner in parm.atoms[atm].bond_partners:
if not partner.marked:
owner_array.add(partner.idx)
_set_owner(parm, owner_array, partner.idx, mol_id)
elif partner.marked != mol_id:
raise MoleculeError('Atom %d in multiple molecules' %
partner.idx)
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _zeros(length):
""" Returns an array of zeros of the given length """
return [0 for i in range(length)]