"""
All of the prmtop actions used in PARMED. Each class is a separate action.
"""
from __future__ import absolute_import, division, print_function
import copy
import math
import os
import sys
import warnings
from collections import Counter, OrderedDict
import numpy as np
from .. import gromacs, unit as u
from ..amber import (AmberAsciiRestart, AmberMask, AmberMdcrd, AmberParm,
AmoebaParm, ChamberParm, NetCDFRestart, NetCDFTraj)
from ..amber._chamberparm import ConvertFromPSF
from ..charmm import CharmmParameterSet, CharmmPsfFile
from ..exceptions import ParmedError
from ..formats import CIFFile, Mol2File, PDBFile
from ..formats.registry import load_file
from ..modeller import AmberOFFLibrary, ResidueTemplateContainer
from ..periodic_table import Element as _Element
from ..residue import (ANION_NAMES, CATION_NAMES, SOLVENT_NAMES,
AminoAcidResidue, DNAResidue, RNAResidue)
from ..structure import Structure
from ..topologyobjects import (Angle, AngleType, Bond, BondType, Dihedral,
DihedralType, DihedralTypeList)
from ..utils.six import add_metaclass, iteritems, string_types
from ..utils.six.moves import range, zip
from .argumentlist import ArgumentList
from .exceptions import (AddPDBError, AddPDBWarning, # CoarseGrainError,
AmbiguousParmError, ArgumentError, ChamberError,
ChangeLJPairError, ChangeStateError,
DeleteDihedralError, FileDoesNotExist, FileExists,
HMassRepartitionError, IncompatibleParmsError,
InputError, LJ12_6_4Error, NoArgument,
NonexistentParm, ParmedChangeError, ParmError,
ParmFileNotFound, ParmWarning, SeriousParmWarning,
SetParamError, SimulationError, TiMergeError,
UnhandledArgumentWarning, WriteOFFError)
from .parmlist import ParmList
GB_RADII = ['amber6', 'bondi', 'mbondi', 'mbondi2', 'mbondi3']
__all__ = [] # This gets populated by the ActionType metaclass
COMMANDMAP = dict()
Usages = dict(go='go', quit='quit', help='help [<action>]', history='history')
class ActionType(type):
"""
Metaclass for defining Action types, to add it to the list of available
actions and provide a hash lookup to the interpreter to support
case-insensitive command access
Parameters
----------
cls : class type
The class that is being generated by this metaclass
name : str
The name of the class being created
bases : tuple of types
Tuple of all base class types for this class
dct: dict
The list of options and attributes currently present in this class
"""
def __init__(cls, name, bases, dct):
global COMMANDMAP, __all__, Usages
__all__.append(name)
if name != 'Action':
COMMANDMAP[name.lower()] = cls
if 'usage' in dct:
Usages[name.lower()] = '%s %s' % (name, dct['usage'])
elif name != 'Action': # Skip base action class
Usages[name.lower()] = name
super(ActionType, cls).__init__(name, bases, dct)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
lawsuit = object
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]@add_metaclass(ActionType)
class Action(lawsuit):
"""
The base class for all ParmEd actions. The metaclass for Action adds the
action and its usage statement to a map that is imported and used by the
ParmEd interpreter. So just adding the subclass, wherever it happens, and
giving it a 'usage' attribute is sufficient to have it added to the
interpreter map. There are a number of attributes that can refine the
behavior of individual Actions in the interpreter:
Parameters
----------
input_parm : Structure
arg_list : str or ArgumentList
Attributes
----------
stderr : file-like, optional
Having a ``write`` attribute, this is where errors and warnings will be
written. Default is sys.stderr
needs_parm : bool, optional
If this Action needs a ``parm`` instance to operate on. If your Action
either acts on the interpreter or *creates* a parm, this can be False.
Default is True
supported_subclasses : tuple
A tuple of all types whose subclasses can be acted upon by this Action.
This is useful for permissive Action classes
strictly_supported : tuple
This Action can *only* act on a particular set of class instances. If
the parm does not have one of the listed types, it is not supported
(subclasses do *not* count here)
not_supported : tuple
These classes are *not* supported (supplements supported_subclasses)
"""
stderr = sys.stderr
# Does this action need a populated parm_list? If yes, bail out if it's
# unpopulated
needs_parm = True
# Do we allow any action to overwrite existing files? Static variable.
overwrite = True # Default for API. parmed sets this every time
# We do one of two kinds of filtering for each action. Each action can
# support either a strict subset of classes, in which case the type *must*
# be an exact member (not a subclass) of one of that list. If that list is
# empty, then we loop through the supported_subclasses static tuple and see
# if our parm is a subclass of any of those classes
supported_subclasses = (Structure,)
strictly_supported = ()
not_supported = ()
def __init__(self, input_parm, arg_list=None, *args, **kwargs):
# Is this action 'valid' (i.e., can we run 'execute')?
self.valid = False
# Accept both an AmberParm or ParmList instance to simplify the API
if isinstance(input_parm, ParmList):
self.parm_list = input_parm
elif isinstance(input_parm, Structure):
self.parm_list = ParmList()
self.parm_list.add_parm(input_parm)
elif input_parm is None:
self.parm_list = ParmList()
else:
raise TypeError('input_parm expected to be a ParmList or Structure instance')
# Set active parm
self.parm = self.parm_list.parm # Current active parm
if self.needs_parm and self.parm_list.empty():
raise ParmError('Action requires a loaded topology file')
# Fill the argument list
if args or kwargs:
if arg_list is None:
arg_list = ''
elif isinstance(arg_list, string_types) and ' ' in arg_list.strip():
# arg_list has a space, so enclose it in quotes
arg_list = '"%s" ' % arg_list
else:
arg_list = '%s ' % arg_list
arg_list += ' '.join([self._format_arg(a) for a in args])
for kw, item in iteritems(kwargs):
arg_list += ' %s %s ' % (kw, self._format_arg(item))
elif arg_list is None:
arg_list = ArgumentList('')
elif isinstance(arg_list, string_types):
arg_list = self._format_arg(arg_list)
# If our arg_list is a string, convert it to an ArgumentList
if isinstance(arg_list, string_types):
arg_list = ArgumentList(arg_list)
else:
arg_list = ArgumentList(str(arg_list))
# Now that we have the argument list, see if we have requested a
# specific parm. If so, set self.parm equal to that object, instead
parm = arg_list.get_key_string('parm', None)
# If it is an integer, see if it is a valid index. If it's not a valid
# index, try it as a string. Otherwise just try it as a string.
if parm is not None:
try:
parm = int(parm)
except ValueError:
if parm in self.parm_list:
print('Using parm %s' % parm)
self.parm = self.parm_list[parm]
else:
warnings.warn('Cannot find parm %s. Skipping this action' % parm,
SeriousParmWarning)
return # pragma: no cover
else:
if parm >= 0 and parm < len(self.parm_list):
print('Using parm %s' % self.parm_list[parm])
self.parm = self.parm_list[parm]
else:
warnings.warn('Cannot find parm %s. Skipping this action' % parm,
SeriousParmWarning)
return # pragma: no cover
if self.needs_parm:
if (self.strictly_supported and
type(self.parm) not in self.strictly_supported):
raise ParmError('%s objects are not supported by this action' %
type(self.parm).__name__)
elif not self.strictly_supported:
for cls in self.not_supported:
if type(self.parm) is cls:
raise ParmError('%s objects are not supported by this action' %
type(self.parm).__name__)
for cls in self.supported_subclasses:
if isinstance(self.parm, cls):
break
else:
raise ParmError('%s objects are not supported by this action' %
type(self.parm).__name__)
self._arg_list = arg_list
try:
self.init(arg_list)
except NoArgument:
cmdname = type(self).__name__
try:
usage = '%s %s' % (cmdname, self.usage)
except AttributeError:
usage = cmdname
Action.stderr.write("Bad command %s:\n\t%s\n" % (cmdname, usage))
return
# Check any unmarked commands
unmarked_cmds = arg_list.unmarked()
if len(unmarked_cmds) > 0:
warnings.warn(' '.join(unmarked_cmds), UnhandledArgumentWarning)
self.valid = True
[docs] def init(self, arg_list):
""" This should be overridden if anything needs to be done """
raise NotImplementedError('init must be overridden by Action subclass')
[docs] def execute(self):
""" Commands involved in executing the action """
pass
def __str__(self):
return ''
@staticmethod
def _format_arg(arg):
"""
Formats an argument so that if it's a string with a space in it, the argument will be
encapsulated with quotes
"""
arg = str(arg)
if ' ' in arg or '\t' in arg:
return ' "%s" ' % arg
return arg
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class parmout(Action):
"""
Final prmtop written after all actions are complete
"""
usage = '<prmtop_name> [<inpcrd_name>] [netcdf]'
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
self.filename = arg_list.get_next_string()
self.rst_name = arg_list.get_next_string(optional=True)
if arg_list.has_key('netcdf'):
self.rst7_format = 'NCRST'
else:
self.rst7_format = 'RST7'
def __str__(self):
if self.rst_name is not None:
return 'Outputting Amber topology file %s and restart %s' % (self.filename,
self.rst_name)
return 'Outputting Amber topology file %s' % self.filename
[docs] def execute(self):
if not Action.overwrite and os.path.exists(self.filename):
raise FileExists('%s exists; not overwriting.' % self.filename)
if self.rst_name is not None:
if not Action.overwrite and os.path.exists(self.rst_name):
raise FileExists('%s exists; not overwriting.' % self.rst_name)
self.parm.write_parm(self.filename)
if self.rst_name is not None:
self.parm.save(self.rst_name, format=self.rst7_format, overwrite=Action.overwrite)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class setOverwrite(Action):
"""
Necessary to overwrite original topology file name.
"""
usage = '[True|False]'
needs_parm = False
[docs] def init(self, arg_list):
# see if we got an argument
arg = arg_list.get_next_string(optional=True)
if arg is not None:
if not arg.lower() in ('false', 'true'):
warnings.warn("setOverwrite: unrecognized argument. Assuming False",
SeriousParmWarning)
self._overwrite = (arg is None or arg.lower() == "true")
def __str__(self):
if self._overwrite:
return 'Files are overwritable'
else:
return 'Files are NOT overwritable'
[docs] def execute(self):
Action.overwrite = self._overwrite
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class writeFrcmod(Action):
"""
Writes an frcmod file from all of the parameters in the topology file.
"""
usage = '<frcmod_name>'
strictly_supported = (AmberParm,)
[docs] def init(self, arg_list):
self.frcmod_name = arg_list.get_next_string(optional=True)
if self.frcmod_name is None: self.frcmod_name = 'frcmod'
# Emit warnings about 10-12 prmtops if we detect any 10-12 parameters
hbond_indexes = set()
for idx in self.parm.parm_data['NONBONDED_PARM_INDEX']:
if idx < 0: hbond_indexes.add(abs(idx))
try:
for idx in hbond_indexes:
if (self.parm.parm_data['HBOND_ACOEF'][idx-1] > 0 or
self.parm.parm_data['HBOND_BCOEF'][idx-1] > 0 or
self.parm.parm_data['HBCUT'][idx-1] > 0):
warnings.warn('Frcmod dumping does not work with 10-12 prmtops',
SeriousParmWarning)
break # pragma: no cover
except IndexError:
pass # pragma: no cover
def __str__(self):
return 'Dumping FRCMOD file %s with parameters from %s' % (
self.frcmod_name, self.parm)
[docs] def execute(self):
""" Writes the frcmod file """
from parmed.amber.parameters import AmberParameterSet
if not Action.overwrite and os.path.exists(self.frcmod_name):
raise FileExists('%s exists; not overwriting' % self.frcmod_name)
parmset = AmberParameterSet.from_structure(self.parm)
title = 'Force field parameters from %s' % os.path.split(str(self.parm))[1]
parmset.write(self.frcmod_name, title=title)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class loadRestrt(Action):
"""
Loads a restart file so we have coordinates. Necessary for distance-based
mask criteria and writeOFF
"""
usage = '<restrt_filename>'
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
self.rst_name = arg_list.get_next_string()
def __str__(self):
return 'Loading restart file %s' % self.rst_name
[docs] def execute(self):
self.parm.load_rst7(self.rst_name)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class loadCoordinates(Action):
"""
Reads a coordinate file and loads the first set of coordinates found into
the active structure. File type is auto-detected. Supported file formats
include:
- Amber restart file
- Amber NetCDF restart file
- CHARMM coordinate file
- CHARMM restart file
- Amber mdcrd trajectory
- Amber NetCDF trajectory
- PDB file
- PDBx/mmCIF file
- Gromacs GRO file
- Mol2 file
"""
usage = '<filename>'
[docs] def init(self, arg_list):
self.filename = arg_list.get_next_string()
def __str__(self):
return 'Adding coordinates to %s from %s' % (self.parm.name, self.filename)
[docs] def execute(self):
crd = load_file(self.filename, natom=len(self.parm.atoms), hasbox=self.parm.box is not None)
try:
self.parm.coordinates = crd.coordinates.copy()
except AttributeError:
raise ParmError('Cannot get coordinates from %s' % self.filename)
self.parm.box = copy.copy(crd.box)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class writeCoordinates(Action):
"""
Writes the coordinates of the active structure to a coordinate file. The
format of the file is detected from filename extension, with the following
extensions being recognized:
- ``.nc``: Amber NetCDF trajectory with a single frame
- ``.ncrst``: Amber NetCDF restart file
- ``.pdb``: PDB file
- ``.cif``: PDBx/mmCIF file
- ``.rst7``, ``.restrt``, ``.inpcrd``: Amber restart file
- ``.mdcrd``: Amber mdcrd file
- ``.mol2``: Sybyl Mol2 file
- Default is Amber restart file
Alternatively, the following keywords can be used to override the filename
extension:
- ``netcdftraj``: Amber NetCDF trajectory with a single frame
- ``netcdf``: Amber NetCDF restart file
- ``pdb``: PDB file
- ``cif``: PDBx/mmCIF file
- ``restart``: Amber restart/inpcrd file
- ``mdcrd``: Amber mdcrd file
- ``mol2``: Sybyl mol2 file
Note, NetCDF files require scipy or netCDF4 to be installed.
"""
usage = ('<filename> [netcdftraj | netcdf | pdb | cif | restart | mdcrd | '
'mol2]')
[docs] def init(self, arg_list):
self.filename = filename = arg_list.get_next_string()
self.filetype = arg_list.get_next_string(optional=True, default='').lower()
if not self.filetype:
if filename.endswith('.nc'):
self.filetype = 'NCTRAJ'
elif filename.endswith('.ncrst'):
self.filetype = 'NCRESTART'
elif filename.endswith('.pdb'):
self.filetype = 'PDB'
elif filename.endswith('.cif'):
self.filetype = 'CIF'
elif (filename.endswith('.rst7') or filename.endswith('.restrt') or
filename.endswith('.inpcrd')):
self.filetype = 'RESTART'
elif filename.endswith('.mdcrd'):
self.filetype = 'MDCRD'
elif filename.endswith('.mol2'):
self.filetype = 'MOL2'
else:
self.filetype = 'RESTART'
else:
if self.filetype == 'netcdftraj':
self.filetype = 'NCTRAJ'
elif self.filetype == 'netcdf':
self.filetype = 'NCRESTART'
elif self.filetype == 'pdb':
self.filetype = 'PDB'
elif self.filetype == 'cif':
self.filetype = 'CIF'
elif self.filetype == 'restart':
self.filetype = 'RESTART'
elif self.filetype == 'mdcrd':
self.filetype = 'MDCRD'
elif self.filetype == 'mol2':
self.filetype = 'MOL2'
else:
raise InputError('Unrecognized file format %s' % self.filetype)
def __str__(self):
return 'Writing coordinates to %s as type %s' % (self.filename, self.filetype)
[docs] def execute(self):
if not Action.overwrite and os.path.exists(self.filename):
raise FileExists('%s exists; not overwriting' % self.filename)
coordinates = self.parm.coordinates
velocities = self.parm.velocities
if self.filetype == 'NCTRAJ':
traj = NetCDFTraj.open_new(self.filename, natom=len(self.parm.atoms),
box=self.parm.box is not None, vels=velocities is not None)
traj.add_time(0)
traj.add_coordinates(coordinates)
if velocities is not None:
traj.add_velocities(velocities)
if self.parm.box is not None:
traj.add_box(self.parm.box)
traj.close()
elif self.filetype == 'NCRESTART':
rst = NetCDFRestart.open_new(self.filename, natom=len(self.parm.atoms),
box=self.parm.box is not None, vels=velocities is not None)
rst.coordinates = coordinates
if velocities is not None:
rst.velocities = velocities
if self.parm.box is not None:
rst.box = self.parm.box
rst.time = 0
rst.close()
elif self.filetype == 'PDB':
PDBFile.write(self.parm, self.filename, renumber=True)
elif self.filetype == 'CIF':
CIFFile.write(self.parm, self.filename, renumber=True)
elif self.filetype == 'MOL2':
Mol2File.write(self.parm, self.filename)
elif self.filetype == 'MDCRD':
traj = AmberMdcrd(self.filename, natom=len(self.parm.atoms),
hasbox=self.parm.box is not None, mode='w')
traj.add_coordinates(coordinates)
if self.parm.box is not None:
traj.add_box(self.parm.box)
traj.close()
elif self.filetype == 'RESTART':
rst = AmberAsciiRestart(self.filename, natom=len(self.parm.atoms), mode='w',
hasbox=self.parm.box is not None)
rst.coordinates = coordinates
if velocities is not None:
rst.velocities = velocities
if self.parm.box is not None:
rst.box = self.parm.box
rst.close()
else:
raise AssertionError('Should not be here. Unrecognized coordinate file format type.')
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class writeOFF(Action):
"""
Writes an Amber OFF Library with all of the residues found in the topology
"""
usage = '<OFF_filename>'
strictly_supported = (AmberParm,)
[docs] def init(self, arg_list):
self.off_file = arg_list.get_next_string()
def __str__(self):
return 'Writing Amber OFF file %s' % self.off_file
[docs] def execute(self):
if not Action.overwrite and os.path.exists(self.off_file):
raise FileExists('%s exists; not overwriting' % self.off_file)
if self.parm.coordinates is None:
raise WriteOFFError('You must load a restart for WriteOFF!')
lib = ResidueTemplateContainer.from_structure(self.parm).to_library()
AmberOFFLibrary.write(lib, self.off_file)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class changeRadii(Action):
"""
Changes intrinsic GB radii to the specified set: Allowed values are
amber6, bondi, mbondi, mbondi2, mbondi3
"""
usage = '<radii_set>'
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.radii = arg_list.get_next_string()
def __str__(self):
return 'Changing PB/GB radii to %s' % self.radii
[docs] def execute(self):
from parmed.tools.changeradii import ChRad
# Add RADIUS_SET to prmtop if it's not there already, and a blank
# description, since it's about to be set here
if isinstance(self.parm, AmberParm):
if not 'RADIUS_SET' in self.parm.flag_list:
self.parm.add_flag('RADIUS_SET', '1a80', num_items=1)
if not 'RADII' in self.parm.flag_list:
self.parm.add_flag('RADII', '5E16.8', num_items=len(self.parm.atoms))
if not 'SCREEN' in self.parm.flag_list:
self.parm.add_flag('SCREEN', '5E16.8', num_items=len(self.parm.atoms))
ChRad(self.parm, self.radii)
# Load the data into the parm arrays
for i, atom in enumerate(self.parm.atoms):
self.parm.parm_data['RADII'][i] = atom.solvent_radius
self.parm.parm_data['SCREEN'][i] = atom.screen
else:
# Otherwise, just set the radii
ChRad(self.parm, self.radii)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class changeLJPair(Action):
"""
Changes a particular Lennard Jones pair based on a given (pre-combined)
epsilon/Rmin
"""
usage = '<mask1> <mask2> <Rmin> <epsilon>'
strictly_supported = (AmberParm, ChamberParm)
[docs] def init(self, arg_list):
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
self.rmin = arg_list.get_next_float()
self.eps = arg_list.get_next_float()
def __str__(self):
return ('Setting LJ %s-%s pairwise interaction to have Rmin = %16.5f and Epsilon = '
'%16.5f') % (self.mask1, self.mask2, self.rmin, self.eps)
[docs] def execute(self):
selection1 = self.mask1.Selection()
selection2 = self.mask2.Selection()
if sum(selection1) == 0 or sum(selection2) == 0:
Action.stderr.write('Skipping empty masks in changeLJPair\n')
return 0
# Make sure we only selected 1 atom type in each mask
attype1 = None
attype2 = None
for i, atom in enumerate(self.parm.atoms):
if selection1[i] == 1:
if attype1 is None:
attype1 = atom.nb_idx
else:
if attype1 != atom.nb_idx:
raise ChangeLJPairError('First mask matches multiple atom types!')
if selection2[i] == 1:
if attype2 is None:
attype2 = atom.nb_idx
else:
if attype2 != atom.nb_idx:
raise ChangeLJPairError('Second mask matches multiple atom types!')
_change_lj_pair(self.parm, attype1, attype2, self.rmin, self.eps)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class changeLJ14Pair(Action):
"""
Changes a particular 1-4 Lennard Jones pair based on a given (pre-combined)
epsilon/Rmin. Only valid for CHAMBER prmtops
"""
usage = '<mask1> <mask2> <Rmin> <epsilon>'
strictly_supported = (ChamberParm,)
[docs] def init(self, arg_list):
# Make sure this is a chamber prmtop
assert self.parm.chamber, 'Chamber-style prmtop required!'
# If not, initiate instance data
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
self.rmin = arg_list.get_next_float()
self.eps = arg_list.get_next_float()
def __str__(self):
return ('Setting LJ 1-4 %s-%s pairwise interaction to have 1-4 Rmin = %16.5f and 1-4 '
'Epsilon = %16.5f' % (self.mask1, self.mask2, self.rmin, self.eps))
[docs] def execute(self):
selection1 = self.mask1.Selection()
selection2 = self.mask2.Selection()
if sum(selection1) == 0 or sum(selection2) == 0:
Action.stderr.write('Skipping empty masks in changeLJ14Pair\n')
return None
# Make sure we only selected 1 atom type, and figure out what it is
attype1 = None
attype2 = None
for i, atom in enumerate(self.parm.atoms):
if selection1[i] == 1:
if attype1 is None:
attype1 = atom.nb_idx
else:
if attype1 != atom.nb_idx:
raise ChangeLJPairError('First mask matches multiple atom types!')
if selection2[i] == 1:
if attype2 is None:
attype2 = atom.nb_idx
else:
if attype2 != atom.nb_idx:
raise ChangeLJPairError('Second mask matches multiple atom types!')
# Adjust 1-4 non-bonded terms as well if we're using a chamber-prmtop
_change_lj_pair(self.parm, attype1, attype2, self.rmin, self.eps, True)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class checkValidity(Action):
"""
Basic checks for prmtop validity.
"""
output = sys.stdout
strictly_supported = (AmberParm, ChamberParm)
[docs] def init(self, arg_list):
pass
def __str__(self):
return 'Determining validity of prmtop'
[docs] def execute(self):
from parmed.tools.checkvalidity import check_validity
from parmed.tools.exceptions import WarningList
# Clear our warnings and start logging them, since check_validity
# reports concerns about the prmtop through the warning system.
warning_log = WarningList(empty_msg=('%s looks OK to me!' % self.parm))
check_validity(self.parm, warning_log)
warning_log.dump(self.output)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class change(Action):
"""
Changes the property of given atoms to a new value. <property> can be
CHARGE, MASS, RADII, SCREEN, ATOM_NAME, ATOM_TYPE, ATOM_TYPE_INDEX,
or ATOMIC_NUMBER (note, changing elements with this command will NOT
change their assignment for SHAKE!).
If given, the [quiet] keyword will prevent ParmEd from printing out a
summary with every property changed for every atom that was changed
useful for suppressing overwhelming output if you are zeroing every
charge, for instance)
"""
usage = '<property> <mask> <new_value> [quiet]'
[docs] def init(self, arg_list):
self.quiet = arg_list.has_key('quiet')
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
self.prop = self.flag_name = prop = arg_list.get_next_string().upper()
self.add_flag = False
if prop in ('CHARGE', 'RADII', 'SCREEN', 'MASS'):
self.new_val = arg_list.get_next_float()
self.new_val_str = '%.4f' % self.new_val
elif prop in ('ATOM_TYPE_INDEX', 'ATOMIC_NUMBER'):
self.new_val = arg_list.get_next_int()
self.new_val_str = '%4i' % self.new_val
elif prop in ('ATOM_NAME', 'AMBER_ATOM_TYPE', 'TREE_CHAIN_CLASSIFICATION', 'ATOM_TYPE'):
self.new_val = arg_list.get_next_string()
if len(self.new_val) > 4:
warnings.warn('Only 4 letters allowed for %s entries! Truncating remaining '
'letters.' % prop, ParmWarning)
self.new_val = self.new_val[:4]
self.new_val_str = '%-4s' % self.new_val
else:
raise ParmedChangeError('You may only use "change" with CHARGE, MASS, RADII, SCREEN, '
'ATOM_NAME, ATOM_TYPE, ATOM_TYPE_INDEX, ATOMIC_NUMBER, or '
'TREE_CHAIN_CLASSIFICATION')
if type(self.parm) is AmoebaParm:
# Catch illegal values for Amoeba topologies
if prop in ('CHARGE', 'RADII', 'SCREEN', 'ATOM_TYPE_INDEX'):
raise ParmedChangeError('You cannot change %s in Amoeba topologies' % prop)
if prop == 'ATOMIC_NUMBER':
# AmoebaParm can have an AMOEBA_ATOMIC_NUMBER section instead of
# ATOMIC_NUMBER. So see which of them is available and change
# that one
if self.parm.amoeba:
self.flag_name = 'AMOEBA_ATOMIC_NUMBER'
if prop not in self.parm.parm_data:
# Add it
self.add_flag = True
def __str__(self):
atnums = self.mask.Selection()
if sum(atnums) == 0:
return "change %s: Nothing to do" % self.prop
if self.quiet:
return "Changing %s of %s to %s" % (self.prop, self.mask, self.new_val)
string = ''
for i, atom in enumerate(self.parm.atoms):
if atnums[i] == 1:
string += "Changing %s of atom # %d (%s) from %s to %s\n" % (
self.prop, i+1, atom.name, self.parm.parm_data[self.prop][i], self.new_val_str
)
return string
[docs] def execute(self):
atnums = self.mask.Selection()
if sum(atnums) == 0:
warnings.warn('change %s: %s matches no atoms' % (self.prop,self.mask), ParmWarning)
return
if self.prop == 'ATOM_TYPE_INDEX':
prop = 'nb_idx'
elif self.prop == 'ATOM_NAME':
prop = 'name'
elif self.prop in ('AMBER_ATOM_TYPE', 'ATOM_TYPE'):
prop = 'type'
elif self.prop == 'TREE_CHAIN_CLASSIFICATION':
prop = 'tree'
elif self.prop == 'RADII':
prop = 'solvent_radius'
else:
prop = self.prop.lower()
for i, atom in enumerate(self.parm.atoms):
if atnums[i] == 1:
setattr(atom, prop, self.new_val)
# Update the raw data for any Amber topologies
if isinstance(self.parm, AmberParm):
if self.add_flag:
addAtomicNumber(self.parm).execute()
for i, atom in enumerate(self.parm.atoms):
self.parm.parm_data[self.flag_name][i] = getattr(atom, prop)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printInfo(Action):
"""
Prints all prmtop data corresponding to the given %FLAG
"""
usage = '<flag>'
outfile = sys.stdout
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
self.flag = arg_list.get_next_string().upper()
if not self.flag in self.parm.flag_list:
warnings.warn('%%FLAG %s not found!' % self.flag,
SeriousParmWarning)
self.found = False # pragma: no cover
else:
if self.parm.formats[self.flag].type is float:
self.format = '%16.5f '
else:
self.format = '%-16s '
self.found = True
def __repr__(self):
ret_str = []
if self.found:
for i, item in enumerate(self.parm.parm_data[self.flag]):
ret_str.append(self.format % item)
if i % 5 == 4:
ret_str.append('\n')
return ''.join(ret_str)
def __str__(self):
return self.__repr__()
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class addLJType(Action):
"""
Turns given mask into a new LJ atom type. It uses the radius and Rmin from
the first atom type in <mask> if new_radius or new_epsilon aren't provided
"""
usage = ('<mask> [radius <new_radius>] [epsilon <new_epsilon>] '
'[radius_14 <new_radius14>] [epsilon_14 <new_epsilon14>]')
strictly_supported = (AmberParm, ChamberParm)
[docs] def init(self, arg_list):
self.new_radius_14 = arg_list.get_key_float('radius_14', None)
self.new_epsilon_14 = arg_list.get_key_float('epsilon_14', None)
self.new_radius = arg_list.get_key_float('radius', None)
self.new_epsilon = arg_list.get_key_float('epsilon', None)
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
def __str__(self):
return 'Making atoms %s into a new LJ atom type' % self.mask
[docs] def execute(self):
from parmed.tools.addljtype import AddLJType
# Find the first atom that's selected in this selection. We've
# already made sure that at least one atom was selected
sel_atms = self.mask.Selection()
for i, atom in enumerate(self.parm.atoms):
if sel_atms[i] == 1:
break
# If either the radius or epsilon were not specified, then pull it
# from the *first* atom in the selection
if self.new_radius is None:
self.new_radius = self.parm.LJ_radius[atom.nb_idx-1]
else:
self.new_radius = self.new_radius
if self.new_epsilon is None:
self.new_epsilon = self.parm.LJ_depth[atom.nb_idx-1]
else:
self.new_epsilon = self.new_epsilon
# Now do the same for chamber prmtops
if self.new_radius_14 is None and self.parm.chamber:
self.new_radius_14 = self.parm.LJ_14_radius[atom.nb_idx-1]
elif not self.parm.chamber:
self.new_radius_14 = None
if self.new_epsilon_14 is None and self.parm.chamber:
self.new_epsilon_14 = self.parm.LJ_14_depth[atom.nb_idx-1]
elif not self.parm.chamber:
self.new_epsilon_14 = None
AddLJType(self.parm, sel_atms, self.new_radius, self.new_epsilon,
self.new_radius_14, self.new_epsilon_14)
self.parm.load_atom_info()
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class outparm(parmout):
"""
Prints a new prmtop like parmout, but keeps its place in the action stack so
several can be written out in 1 parmed session
"""
usage = '<prmtop_name> [<inpcrd_name>] [netcdf]'
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printLJTypes(Action):
"""
Prints the Lennard Jones type index for the given atom mask or, if no value
is given, the LJ type index for each atom.
"""
usage = '[<mask>|<type idx>]'
strictly_supported = (AmberParm, ChamberParm)
[docs] def init(self, arg_list):
# Compile the list of indices
try:
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
self.type_list = None
try_type = False
except NoArgument:
try_type = True
if try_type:
try:
self.type_list = arg_list.get_next_string()
self.mask = None
except NoArgument:
# Our default is print all indices...
self.mask = AmberMask(self.parm, ':*')
if self.mask is None:
type_fields = self.type_list.strip().split(',')
self.type_list = []
for field in type_fields:
if len(field.strip()) == 0: continue
if '-' in field:
begin = int(field.split('-')[0])
end = min(int(field.split('-')[1]), self.parm.ptr('ntypes'))
if (begin <= 0 or end < begin or
begin > self.parm.ptr('ntypes')):
raise ParmError('printLJTypes: Bad atom type range')
self.type_list.extend([i for i in range(begin, end+1)])
else:
idx = int(field)
if idx <= 0 or idx > self.parm.ptr('ntypes'):
raise ParmError('printLJTypes: Bad atom type index')
self.type_list.append(int(field))
def __str__(self):
return self.__repr__()
def __repr__(self):
# Construct the atom selections and related atom types lists
if self.mask:
selection = self.mask.Selection()
elif self.type_list:
selection = [0 for atom in self.parm.atoms]
for i, atom in enumerate(self.parm.atoms):
if atom.nb_idx in self.type_list:
selection[i] = 1
if sum(selection) == 0:
return 'Nothing to do for printLJTypes'
indices = set()
for i, atom in enumerate(self.parm.atoms):
if selection[i] == 1:
indices.add(atom.nb_idx)
string = '\n%15s %4s %4s\n' % (" ATOM NUMBER ", 'NAME', 'TYPE')
string += '---------------------------------------------\n'
for i, atom in enumerate(self.parm.atoms):
if atom.nb_idx in indices:
string += 'ATOM %-10d %-4s %-4s: Type index: %d\n' % (i+1, atom.name, atom.type,
atom.nb_idx)
return string
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class scee(Action):
"""
Sets the 1-4 EEL scaling factor in the prmtop
"""
usage = '<scee_value>'
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.scee_value = arg_list.get_next_float()
def __str__(self):
return ("Setting default SCEE electrostatic scaling value to %.4f" % self.scee_value)
[docs] def execute(self):
for dt in self.parm.dihedral_types:
dt.scee = self.scee_value
if isinstance(self.parm, AmberParm):
nptra = self.parm.ptr('nptra')
if not 'SCEE_SCALE_FACTOR' in self.parm.flag_list:
self.parm.add_flag('SCEE_SCALE_FACTOR', '5E16.8',
data=[self.scee_value for i in range(nptra)])
else:
self.parm.parm_data['SCEE_SCALE_FACTOR'] = [self.scee_value for i in range(nptra)]
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class scnb(Action):
"""
Sets the 1-4 VDW scaling factor in the prmtop
"""
usage = '<scnb_value>'
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.scnb_value = arg_list.get_next_float()
def __str__(self):
return "Setting default SCNB van der Waals scaling value to %.4f" % self.scnb_value
[docs] def execute(self):
for dt in self.parm.dihedral_types:
dt.scnb = self.scnb_value
if isinstance(self.parm, AmberParm):
nptra = self.parm.ptr('nptra')
if not 'SCNB_SCALE_FACTOR' in self.parm.flag_list:
self.parm.add_flag('SCNB_SCALE_FACTOR','5E16.8',
data=[self.scnb_value for i in range(nptra)])
else:
self.parm.parm_data['SCNB_SCALE_FACTOR'] = [self.scnb_value for i in range(nptra)]
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class changeLJSingleType(Action):
"""
Allows you to change the radius/well depth of a single LJ type specified by
<mask>. Note, this may change more than the atoms selected in just the mask!
To find out what else will be changed, look at the output of printLJTypes.
Use addLJType to change the Lennard-Jones parameters on a set of specific
atoms.
- <mask> : The selection of atoms to change the LJ type for. All atoms
selected must have the same LJ type
- <radius> : Rmin/2 (van der Waals radius of the new type)
- <depth> : Well-depth (a.k.a., epsilon)
"""
usage = '<mask> <radius> <depth>'
strictly_supported = (AmberParm, ChamberParm)
[docs] def init(self, arg_list):
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
self.radius = arg_list.get_next_float()
self.depth = arg_list.get_next_float()
self.orig_radius = self.orig_depth = None
selection = self.mask.Selection()
if 1 in selection:
first_loc = selection.index(1)
attype = self.parm.atoms[first_loc].nb_idx
self.orig_radius = self.parm.LJ_radius[attype-1]
self.orig_depth = self.parm.LJ_depth[attype-1]
def __str__(self):
if sum(self.mask.Selection()) == 0:
return "No atoms selected in %s. Nothing to do." % self.mask
return ("Changing %s Lennard-Jones well depth from %.4f to %.4f (kal/mol) and radius from "
"%.4f to %.4f (Angstroms)" % (self.mask, self.orig_depth, self.depth,
self.orig_radius, self.radius)
)
[docs] def execute(self):
from math import sqrt
from parmed.tools.exceptions import LJ_TypeError
# If this is an empty mask do nothing
if self.orig_radius is None: return
# Make sure we've only selected a single atom type with our mask
attype = None
selection = self.mask.Selection()
for i, atom in enumerate(self.parm.atoms):
if selection[i] == 1:
if attype is None:
attype = atom.nb_idx
else:
if attype != atom.nb_idx:
raise LJ_TypeError('changeLJSingleType: Mask has multiple atom types!')
# Fill the Lennard-Jones radius and depth arrays to make sure they're
# up-to-date
self.parm.fill_LJ()
self.parm.LJ_radius[attype-1] = self.radius
self.parm.LJ_depth[attype-1] = self.depth
ntypes = self.parm.ptr('NTYPES')
for i in range(ntypes):
lj_index = self.parm.parm_data['NONBONDED_PARM_INDEX'][ntypes*i+attype-1] - 1
rij = self.parm.LJ_radius[i] + self.radius
wij = sqrt(self.parm.LJ_depth[i] * self.depth)
acoef = wij * rij ** 12
bcoef = 2 * wij * rij ** 6
self.parm.parm_data['LENNARD_JONES_ACOEF'][lj_index] = acoef
self.parm.parm_data['LENNARD_JONES_BCOEF'][lj_index] = bcoef
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printDetails(Action):
"""
Returns information about all atoms in a given mask
"""
usage = '<mask>'
[docs] def init(self, arg_list):
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
def __str__(self):
return self.__repr__()
def __repr__(self):
selection = self.mask.Selection()
retstr = ["\nThe mask %s matches %d atoms:\n\n" % (self.mask, sum(selection))]
# Separate printout for Amoeba-style prmtop files
if isinstance(self.parm, AmoebaParm):
retstr.append('%7s%7s%9s%6s%6s%7s%10s\n' % ('ATOM', 'RES', 'RESNAME', 'NAME', 'TYPE',
'At.#', 'Mass')
)
for i, atm in enumerate(self.parm.atoms):
if not selection[i]: continue
retstr.append('%7d%7d%9s%6s%6s%7d%10.4f\n' % (
i+1, atm.residue.idx+1, atm.residue.name, atm.name, atm.type,
atm.atomic_number, atm.mass)
)
else:
retstr.append("%7s%7s%9s%6s%6s%7s%12s%12s%10s%10s%10s%10s\n" %
('ATOM', 'RES', 'RESNAME', 'NAME', 'TYPE', 'At.#', 'LJ Radius',
'LJ Depth', 'Mass', 'Charge', 'GB Radius', 'GB Screen')
)
for i, atm in enumerate(self.parm.atoms):
if not selection[i]: continue
retstr.append("%7d%7d%9s%6s%6s%7d%12.4f%12.4f%10.4f%10.4f%10.4f%10.4f\n" %
(i+1, atm.residue.idx+1, atm.residue.name, atm.name, atm.type,
atm.atomic_number, atm.rmin, atm.epsilon, atm.mass, atm.charge,
atm.solvent_radius, atm.screen)
)
return ''.join(retstr)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printFlags(Action):
"""
Prints all %FLAGs found in the topology file
"""
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
pass
def __str__(self):
return self.__repr__()
def __repr__(self):
string = ['\n']
string.extend('%%FLAG %s\n' % flag for flag in self.parm.flag_list)
return ''.join(string)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printPointers(Action):
"""
Prints a list of all the POINTERS and their values
"""
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
pass
def __str__(self):
return self.__repr__()
def __repr__(self):
ptrs = self.parm.parm_data['POINTERS']
ret_str = """
NATOM (number of atoms in system)................= %d
NTYPES (number of atom type names)...............= %d
NBONH (number of bonds containing H).............= %d
MBONA (number of bonds without H)................= %d
NTHETH (number of angles containing H)...........= %d
MTHETA (number of angles without H)..............= %d
NPHIH (number of dihedrals containing H).........= %d
MPHIA (number of dihedrals without H)............= %d
NHPARM (currently unused)........................= %d
NPARM (1 if made with addles, 0 if not)..........= %d
NNB (number of excluded atoms)...................= %d
NRES (number of residues in system)..............= %d
NBONA (MBONA + constraint bonds).................= %d
NTHETA (MTHETA + constraint angles)..............= %d
NPHIA (MPHIA + constraint dihedrals).............= %d
NUMBND (number of unique bond types).............= %d
NUMANG (number of unique angle types)............= %d
NPTRA (number of unique dihedral types)..........= %d
NATYP (number of nonbonded atom types)...........= %d
NPHB (number of distinct 10-12 H-bond pairs).....= %d
IFPERT (1 if prmtop is perturbed; not used)......= %d
NBPER (perturbed bonds; not used)................= %d
NGPER (perturbed angles; not used)...............= %d
NDPER (perturbed dihedrals; not used)............= %d
MBPER (bonds in perturbed group; not used).......= %d
MGPER (angles in perturbed group; not used)......= %d
MDPER (diheds in perturbed group; not used)......= %d
IFBOX (Type of box: 1=orthogonal, 2=not, 0=none).= %d
NMXRS (number of atoms in largest residue).......= %d
IFCAP (1 if solvent cap exists)..................= %d
NUMEXTRA (number of extra points in topology)....= %d
""" % tuple(ptrs[:31])
if len(ptrs) == 32:
ret_str += "NCOPY (number of PIMD slices/number of beads)....= %d\n" % ptrs[31]
if self.parm.ptr('IFBOX'):
ret_str += "\nSOLVENT POINTERS\n" + """
IPTRES (Final solute residue)....................= %d
NSPM (Total number of molecules).................= %d
NSPSOL (The first solvent "molecule")............= %d
""" % (self.parm.parm_data['SOLVENT_POINTERS'][0],
self.parm.parm_data['SOLVENT_POINTERS'][1],
self.parm.parm_data['SOLVENT_POINTERS'][2])
return ret_str
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class setMolecules(Action):
"""
Determines the molecularity of the system based on the bonding network and
correctly determines the SOLVENT_POINTERS and ATOMS_PER_MOLECULE sections of
the topology file. It will consider the ions to be part of the solute if
True is passed or not if False is passed. Defaults to True.
"""
usage = '[solute_ions True|False]'
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
# Allow solute_ions to be a keyword argument or the only argument.
# Default is True
solute_ions = arg_list.get_key_string('solute_ions', None)
if solute_ions is None:
self.solute_ions = arg_list.get_next_string(optional=True)
elif solute_ions.lower() == 'true':
self.solute_ions = True
elif solute_ions.lower() == 'false':
self.solute_ions = False
else:
warnings.warn("Value of solute_ions is unrecognized [%s]! Assuming True" % solute_ions,
SeriousParmWarning)
self.solute_ions = True # pragma: no cover
def __str__(self):
return ("Setting MOLECULE properties of the prmtop (SOLVENT_POINTERS "
"and ATOMS_PER_MOLECULE)")
[docs] def execute(self):
owner = self.parm.rediscover_molecules(self.solute_ions)
if owner is not None:
if self.parm.coordinates is None:
warnings.warn('The atoms in %s were reordered to correct molecule ordering. Any '
'topology printed from now on will *not* work with the original '
'inpcrd or trajectory files created with this prmtop! Consider '
'quitting and loading a restart prior to using setMolecules' %
self.parm, ParmWarning)
# If we had to reorder our atoms, we need to remake our parm
self.parm.remake_parm()
self.parm.load_pointers()
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
#
#class addCoarseGrain(Action):
# """
# Adds coarse graining information to the Amber topology file according to
# a given coarse graining parameter file
# """
# strictly_supported = (AmberParm,)
# usage = '<parameter_file>
# def init(self, arg_list):
# self.cg_param_file = arg_list.get_next_string()
# if not os.path.exists(self.cg_param_file):
# raise CoarseGrainError('Cannot find parameter file %s' %
# self.cg_param_file)
# # Check to see if we've already made this a coarsegrained file...
# if 'ANGLE_COEF_A' in self.parm.flag_list:
# warnings.warn('Prmtop already has coarse grained sections',
# ParmWarning)
#
# def __str__(self):
# return ("Setting up coarse graining for topology file using parameter "
# "file " + self.cg_param_file)
#
# def execute(self):
# from parmed.tools.coarsegrain import addCoarseGrain as addCG
# if 'ANGLE_COEF_A' in self.parm.flag_list: return
# addCG(self.parm, self.cg_param_file)
#
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class changeProtState(Action):
"""
Changes the protonation state of a given titratable residue that can be
treated via constant pH MD in Amber.
"""
usage = '<mask> <state #>'
strictly_supported = (AmberParm,)
[docs] def init(self, arg_list):
self.state = arg_list.get_next_int()
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
def __str__(self):
sel = self.mask.Selection()
if sum(sel) == 0:
return "No residues selected for state change"
res = self.parm.atoms[sel.index(1)].residue
return 'Changing protonation state of residue %d (%s) to %d' % (res.idx+1, res.name,
self.state)
@staticmethod
def _add_ash_glh(residues):
"""
Adds ASH and GLH to the titratable residue list unless it's already
there
"""
if 'ASH' in residues.titratable_residues: return None
dummyrefene1 = residues._ReferenceEnergy()
dummyrefene1_old = residues._ReferenceEnergy()
dummyrefene1_old.set_pKa(1.0)
dummyrefene2 = residues._ReferenceEnergy()
atomnames = ['N', 'H', 'CA', 'HA', 'CB', 'HB2', 'HB3', 'CG', 'OD1', 'OD2', 'HD21', 'C', 'O']
ash = residues.TitratableResidue('ASH', atomnames, pka=4.0, typ="ph")
ash.add_state(protcnt=0, refene=dummyrefene1, refene_old=dummyrefene1_old, pka_corr=0.0,
charges=[-0.4157, 0.2719, 0.0341, 0.0864, -0.1783, -0.0122, -0.0122, 0.7994,
-0.8014, -0.8014, 0.0, 0.5973, -0.5679]
)
ash.add_state(protcnt=1, refene=dummyrefene2, refene_old=dummyrefene2, pka_corr=4.0,
charges=[-0.4157, 0.2719, 0.0341, 0.0864, -0.0316, 0.0488, 0.0488, 0.6462,
-0.5554, -0.6376, 0.4747, 0.5973, -0.5679]
)
ash.check()
atomnames = ['N', 'H', 'CA', 'HA', 'CB', 'HB2', 'HB3', 'CG', 'HG2', 'HG3', 'CD', 'OE1',
'OE2', 'HE21', 'C', 'O']
glh = residues.TitratableResidue('GLH', atomnames, pka=4.4, typ="ph")
glh.add_state(protcnt=0, refene=dummyrefene1, refene_old=dummyrefene1_old, pka_corr=0.0,
charges=[-0.4157, 0.2719, 0.0145, 0.0779, -0.0398, -0.0173, -0.0173, 0.0136,
-0.0425, -0.0425, 0.8054, -0.8188, -0.8188, 0.0, 0.5973, -0.5679]
)
glh.add_state(protcnt=1, refene=dummyrefene2, refene_old=dummyrefene2, pka_corr=4.4,
charges=[-0.4157, 0.2719, 0.0145, 0.0779, -0.0071, 0.0256, 0.0256, -0.0174,
0.0430, 0.0430, 0.6801, -0.5838, -0.6511, 0.4641, 0.5973, -0.5679]
)
glh.check()
residues.ASH, residues.GLH = ash, glh
residues.titratable_residues.extend(['ASH', 'GLH'])
[docs] def execute(self):
from parmed.amber import titratable_residues as residues
changeProtState._add_ash_glh(residues)
sel = self.mask.Selection()
# If we didn't select any residues, just return
if sum(sel) == 0: return
residue = self.parm.atoms[sel.index(1)].residue
resname = residue.name
# Get the charges from cpin_data. The first 2 elements are energy and
# proton count so the charges are chgs[2:]
if not resname in residues.titratable_residues:
raise ChangeStateError("Residue %s isn't defined as a titratable "
"residue in titratable_residues.py" % resname)
if not getattr(residues, resname).typ == "ph":
raise ChangeStateError('Redidue %s is not a pH titratable residue' % resname)
res = getattr(residues, resname)
if self.state >= len(res.states):
raise ChangeStateError('Residue %s only has titratable states 0--%d. You chose state %d'
% (resname, len(res.states)-1, self.state))
if sum(sel) != len(res.states[self.state].charges):
raise ChangeStateError('You must select one and only one entire titratable residue')
charges = res.states[self.state].charges
for i, atom in enumerate(residue.atoms):
if sel[atom.idx] != 1:
raise ChangeStateError('You must select 1 and only 1 entire residue of which to '
'change the protonation state')
# Actually make the change
self.parm.parm_data['CHARGE'][atom.idx] = atom.charge = charges[i]
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class changeRedoxState(Action):
"""
Changes the reduction state of a given titratable residue that can be
treated via constant redox potential MD in Amber.
"""
usage = '<mask> <state #>'
strictly_supported = (AmberParm,)
[docs] def init(self, arg_list):
self.state = arg_list.get_next_int()
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
def __str__(self):
sel = self.mask.Selection()
if sum(sel) == 0:
return "No residues selected for state change"
res = self.parm.atoms[sel.index(1)].residue
return 'Changing reduction state of residue %d (%s) to %d' % (res.idx+1, res.name,
self.state)
[docs] def execute(self):
from parmed.amber import titratable_residues as residues
sel = self.mask.Selection()
# If we didn't select any residues, just return
if sum(sel) == 0: return
residue = self.parm.atoms[sel.index(1)].residue
resname = residue.name
# Get the charges from cein_data. The first 2 elements are energy and
# electron count so the charges are chgs[2:]
if not resname in residues.titratable_residues:
raise ChangeStateError("Residue %s isn't defined as a titratable "
"residue in titratable_residues.py" % resname)
if not getattr(residues, resname).typ == "redox":
raise ChangeStateError('Redidue %s is not a redox potential titratable residue' % resname)
res = getattr(residues, resname)
if self.state >= len(res.states):
raise ChangeStateError('Residue %s only has titratable states 0--%d. You chose state %d'
% (resname, len(res.states)-1, self.state))
if sum(sel) != len(res.states[self.state].charges):
raise ChangeStateError('You must select one and only one entire titratable residue')
charges = res.states[self.state].charges
for i, atom in enumerate(residue.atoms):
if sel[atom.idx] != 1:
raise ChangeStateError('You must select 1 and only 1 entire residue of which to '
'change the reduction state')
# Actually make the change
self.parm.parm_data['CHARGE'][atom.idx] = atom.charge = charges[i]
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class netCharge(Action):
"""
Prints the total charge of all of the atoms given by the mask. Defaults to all atoms
"""
usage = '[<mask>]'
outfile = sys.stdout
[docs] def init(self, arg_list):
mask = arg_list.get_next_mask(optional=True)
if mask is None: mask = ':*'
self.mask = AmberMask(self.parm, mask)
def __str__(self):
return 'The net charge of %s is %.4f' % (self.mask,
sum([self.parm.atoms[i].charge for i in self.mask.Selected()]))
[docs] def execute(self):
""" Calculates the charge of all atoms selected in mask """
return sum([self.parm.atoms[i].charge for i in self.mask.Selected()])
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class strip(Action):
"""
Deletes the atoms specified by <mask> from the topology file and rebuilds
the topology file according to the parameters that remain. If nobox is
provided, the unit cell information is discarded (useful when stripping
solvent to run an aperiodic implicit solvent calculation).
"""
usage = '<mask> [nobox]'
[docs] def init(self, arg_list):
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
self.nobox = arg_list.has_key('nobox')
self.num_atoms = sum(self.mask.Selection())
def __str__(self):
retstr = ["Removing mask '%s' (%d atoms) from the topology file." %
(self.mask, self.num_atoms)]
if self.nobox:
retstr.append('Deleting box info.')
return ' '.join(retstr)
[docs] def execute(self):
self.parm.strip(self.mask)
if self.nobox:
self.parm.box = None
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class defineSolvent(Action):
"""
Allows you to change what parmed will consider to be "solvent".
<residue list> must be a comma-separated set of residue names with no
spaces between them.
"""
usage = '<residue_list>'
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
from parmed import residue
res_list = arg_list.get_next_string()
res_list.replace(' ', '')
if res_list.endswith(','):
self.res_list = res_list[:len(res_list)-1]
else:
self.res_list = res_list
residue.SOLVENT_NAMES = self.res_list.split(',')
def __str__(self):
return "Residues %s are now considered to be solvent" % self.res_list
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class addExclusions(Action):
"""
Allows you to add arbitrary exclusions to the exclusion list. Every atom in
<mask2> is added to the exclusion list for each atom in <mask1> so that
non-bonded interactions between those atom pairs will not be computed. NOTE
that this ONLY applies to direct-space (short-range) non-bonded potentials.
For PME simulations, long-range electrostatics between these atom pairs are
still computed (in different unit cells).
"""
usage = '<mask1> <mask2>'
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
def __str__(self):
return 'Adding atoms from %s to exclusion lists of atoms in %s' % (
self.mask2, self.mask1)
[docs] def execute(self):
# Loop through both selections and add each selected atom in sel2 to
# the exclusion list for selected atoms in sel1 (and vice-versa).
for i in self.mask1.Selected():
atm1 = self.parm.atoms[i]
for j in self.mask2.Selected():
atm2 = self.parm.atoms[j]
# Skip over atm1 == atm2
if atm1 is atm2: continue
# Add each other to each other's exclusion lists.
atm1.exclude(atm2)
self.parm.atoms.changed = True
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printBonds(Action):
"""
Prints all of the bonds (with their details) for the given atoms in the
mask. If a second mask is given, only bonds in which one atom appears in
*each* list will be printed. If coordinates and parameter types are present,
also print the actual distance (in Angstroms) and energy (in kcal/mol) for
each printed bond.
"""
usage = '[<mask> [<mask>] ]'
[docs] def init(self, arg_list):
mask = arg_list.get_next_mask(optional=True, default='*')
self.mask = AmberMask(self.parm, mask)
mask = arg_list.get_next_mask(optional=True, default='*')
self.mask2 = AmberMask(self.parm, mask)
def __str__(self):
return self.__repr__()
def __repr__(self):
retstr = ['%19s %19s %10s %10s' %
('Atom 1', 'Atom 2', 'R eq', 'Frc Cnst')]
# Loop through all of the bonds without and inc hydrogen
atomsel = set(self.mask.Selected())
atomsel2 = set(self.mask2.Selected())
do_measured = self.parm.coordinates is not None
do_energy = all(b.type is not None for b in self.parm.bonds)
if do_measured:
retstr.append(' %10s' % 'Distance')
if do_energy:
retstr.append(' %10s' % 'Energy')
retstr.append('\n')
for bond in self.parm.bonds:
atom1, atom2 = bond.atom1, bond.atom2
found = False
if atom1.idx in atomsel and atom2.idx in atomsel2:
found = True
elif atom2.idx in atomsel and atom1.idx in atomsel2:
found = True
if not found: continue
if bond.type is not None:
retstr.append('%7d %4s (%4s) %7d %4s (%4s) %10.4f %10.4f' % (
atom1.idx+1, atom1.name, atom1.type, atom2.idx+1,
atom2.name, atom2.type, bond.type.req, bond.type.k)
)
else:
retstr.append('%7d %4s (%4s) %7d %4s (%4s) %-10s %-10s' % (
atom1.idx+1, atom1.name, atom1.type, atom2.idx+1,
atom2.name, atom2.type, 'N/A', 'N/A'))
if do_measured:
retstr.append(' %10.4f' % bond.measure())
if do_energy:
retstr.append(' %10.4f' % bond.energy())
retstr.append('\n')
return ''.join(retstr)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printAngles(Action):
"""
Prints all of the angles (with their details) for the given atoms in the
mask. If a second mask is given, only atoms whose central atom is in the
second mask and another atom is in the first mask is printed. If a third
mask is given, the central atom must be in the second mask and the other two
atoms must appear in the first *and* third masks (in any order).
If coordinates and parameter types are present, the value of the angle (in
degrees) and its energy (in kcal/mol) are reported for each printed angle.
"""
usage = '[<mask> [<mask> [<mask>] ] ]'
[docs] def init(self, arg_list):
mask = arg_list.get_next_mask(optional=True, default=':*')
self.mask = AmberMask(self.parm, mask)
arg2 = arg_list.get_next_mask(optional=True)
arg3 = arg_list.get_next_mask(optional=True)
if arg2 is None:
self.one_arg = True
else:
self.one_arg = False
self.mask2 = AmberMask(self.parm, arg2)
if arg3 is None: arg3 = '*'
self.mask3 = AmberMask(self.parm, arg3)
def __str__(self):
return self.__repr__()
def __repr__(self):
retstr = ['%19s %19s %19s %10s %10s' %
('Atom 1', 'Atom 2', 'Atom 3', 'Frc Cnst', 'Theta eq')]
do_measured = self.parm.coordinates is not None
do_energy = all(a.type is not None for a in self.parm.angles)
if do_measured:
retstr.append(' %10s' % 'Angle')
if do_energy:
retstr.append(' %10s' % 'Energy')
retstr.append('\n')
if self.one_arg:
atomsel = self.mask.Selection()
for angle in self.parm.angles:
atom1, atom2, atom3 = angle.atom1, angle.atom2, angle.atom3
if not (atomsel[atom1.idx] or atomsel[atom2.idx] or
atomsel[atom3.idx]):
continue
if angle.type is not None:
retstr.append('%7d %4s (%4s) %7d %4s (%4s) %7d %4s (%4s) '
'%10.4f %10.4f' % (atom1.idx+1, atom1.name,
atom1.type, atom2.idx+1, atom2.name, atom2.type,
atom3.idx+1, atom3.name, atom3.type, angle.type.k,
angle.type.theteq)
)
else:
retstr.append('%7d %4s (%4s) %7d %4s (%4s) %7d %4s (%4s) '
'%-10s %-10s' % (atom1.idx+1, atom1.name,
atom1.type, atom2.idx+1, atom2.name, atom2.type,
atom3.idx+1, atom3.name, atom3.type, 'N/A', 'N/A')
)
if do_measured:
retstr.append(' %10.4f' % angle.measure())
if do_energy:
retstr.append(' %10.4f' % angle.energy())
retstr.append('\n')
else:
atomsel = set(self.mask.Selected())
atomsel2 = set(self.mask2.Selected())
atomsel3 = set(self.mask3.Selected())
for angle in self.parm.angles:
atom1, atom2, atom3 = angle.atom1, angle.atom2, angle.atom3
if atom2.idx not in atomsel2:
continue
if atom1.idx not in atomsel and atom1.idx not in atomsel3:
continue
found = False
if atom1.idx in atomsel and atom3.idx in atomsel3:
found = True
elif atom3.idx in atomsel and atom1.idx in atomsel3:
found = True
if not found: continue
if angle.type is not None:
retstr.append('%7d %4s (%4s) %7d %4s (%4s) %7d %4s (%4s) '
'%10.4f %10.4f' % (atom1.idx+1, atom1.name,
atom1.type, atom2.idx+1, atom2.name, atom2.type,
atom3.idx+1, atom3.name, atom3.type, angle.type.k,
angle.type.theteq)
)
else:
retstr.append('%7d %4s (%4s) %7d %4s (%4s) %7d %4s (%4s) '
'%-10s %-10s' % (atom1.idx+1, atom1.name,
atom1.type, atom2.idx+1, atom2.name, atom2.type,
atom3.idx+1, atom3.name, atom3.type, 'N/A', 'N/A')
)
if do_measured:
retstr.append(' %10.4f' % angle.measure())
if do_energy:
retstr.append(' %10.4f' % angle.energy())
retstr.append('\n')
return ''.join(retstr)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printDihedrals(Action):
"""
Prints all of the dihedrals (with their details) for the given atoms in the
mask. If multiple masks are given, only dihedrals that have one atom in each
mask are printed. Ordering is important here, so the first atom must be in
the first mask, the second atom in the second, etc. The order can be
precisely reversed, but no other ordering is recognized.
If coordinates and parameter types are present, the value of the torsion
angle (in degrees) and the energy of each dihedral (in kcal/mol) are
reported for each printed dihedral.
"""
usage = '[<mask> [<mask> [<mask> [<mask>] ] ] ]'
[docs] def init(self, arg_list):
mask = arg_list.get_next_mask(optional=True, default=':*')
self.mask = AmberMask(self.parm, mask)
arg2 = arg_list.get_next_mask(optional=True)
arg3 = arg_list.get_next_mask(optional=True)
arg4 = arg_list.get_next_mask(optional=True)
if arg2 is None:
self.one_mask = True
else:
self.one_mask = False
self.mask2 = AmberMask(self.parm, arg2)
if arg3 is None: arg3 = '*'
self.mask3 = AmberMask(self.parm, arg3)
if arg4 is None: arg4 = '*'
self.mask4 = AmberMask(self.parm, arg4)
def __str__(self):
return self.__repr__()
def __repr__(self):
retstr = ['%21s %19s %19s %19s %10s %10s %10s %10s %10s' % (
'Atom 1', 'Atom 2', 'Atom 3', 'Atom 4', 'Height', 'Periodic.',
'Phase', 'EEL Scale', 'VDW Scale')]
do_measured = self.parm.coordinates is not None
do_energy = all(isinstance(d.type, (DihedralType, DihedralTypeList))
for d in self.parm.dihedrals)
if do_measured:
retstr.append(' %10s' % 'Dihedral')
if do_energy:
retstr.append(' %10s' % 'Energy')
retstr.append('\n')
# Loop through all of the bonds without and inc hydrogen
if self.one_mask:
atomsel = self.mask.Selection()
for dihedral in self.parm.dihedrals:
atom1 = dihedral.atom1
atom2 = dihedral.atom2
atom3 = dihedral.atom3
atom4 = dihedral.atom4
if not (atomsel[atom1.idx] or atomsel[atom2.idx] or
atomsel[atom3.idx] or atomsel[atom4.idx]):
continue
# Determine if it's an Improper, Multiterm, or neither
if isinstance(self.parm, AmoebaParm):
char = ' '
scee = scnb = 'N/A'
k = '%10.4f' % dihedral.type.phi_k
per = '%10.4f' % dihedral.type.per
phase = '%10.4f' % dihedral.type.phase
elif dihedral.improper:
char = 'I'
if dihedral.type is not None:
scee = '%10.4f' % dihedral.type.scee
scnb = '%10.4f' % dihedral.type.scnb
k = '%10.4f' % dihedral.type.phi_k
per = '%10.4f' % dihedral.type.per
phase = '%10.4f' % dihedral.type.phase
else:
scee = scnb = k = per = phase = 'N/A'
elif dihedral.ignore_end:
char = 'M'
if dihedral.type is not None:
scee = '%10.4f' % dihedral.type.scee
scnb = '%10.4f' % dihedral.type.scnb
k = '%10.4f' % dihedral.type.phi_k
per = '%10.4f' % dihedral.type.per
phase = '%10.4f' % dihedral.type.phase
else:
scee = scnb = k = per = phase = 'N/A'
else:
char = ' '
if dihedral.type is not None:
scee = '%10.4f' % dihedral.type.scee
scnb = '%10.4f' % dihedral.type.scnb
k = '%10.4f' % dihedral.type.phi_k
per = '%10.4f' % dihedral.type.per
phase = '%10.4f' % dihedral.type.phase
else:
scee = scnb = k = per = phase = 'N/A'
retstr.append('%1s %7d %4s (%4s) %7d %4s (%4s) %7d %4s (%4s) '
' %7d %4s (%4s) %10s %10s %10s %10s %10s' %
(char, atom1.idx+1, atom1.name, atom1.type, atom2.idx+1,
atom2.name, atom2.type, atom3.idx+1, atom3.name,
atom3.type, atom4.idx+1, atom4.name, atom4.type,
k, per, phase, scee, scnb)
)
if do_measured:
retstr.append(' %10.4f' % dihedral.measure())
if do_energy:
retstr.append(' %10.4f' % dihedral.energy())
retstr.append('\n')
else:
atomsel = set(self.mask.Selected())
atomsel2 = set(self.mask2.Selected())
atomsel3 = set(self.mask3.Selected())
atomsel4 = set(self.mask4.Selected())
for dihedral in self.parm.dihedrals:
atom1 = dihedral.atom1
atom2 = dihedral.atom2
atom3 = dihedral.atom3
atom4 = dihedral.atom4
if atom1.idx not in atomsel and atom1.idx not in atomsel4:
continue
found = False
if (atom1.idx in atomsel and atom2.idx in atomsel2 and
atom3.idx in atomsel3 and atom4.idx in atomsel4):
found = True
elif (atom1.idx in atomsel4 and atom2.idx in atomsel3 and
atom3.idx in atomsel2 and atom4.idx in atomsel):
found = True
if not found: continue
if isinstance(self.parm, AmoebaParm):
char = ' '
scee = scnb = 'N/A'
k = '%10.4f' % dihedral.type.phi_k
per = '%10.4f' % dihedral.type.per
phase = '%10.4f' % dihedral.type.phase
elif dihedral.improper:
char = 'I'
if dihedral.type is not None:
scee = '%10.4f' % dihedral.type.scee
scnb = '%10.4f' % dihedral.type.scnb
k = '%10.4f' % dihedral.type.phi_k
per = '%10.4f' % dihedral.type.per
phase = '%10.4f' % dihedral.type.phase
else:
scee = scnb = k = per = phase = 'N/A'
elif dihedral.ignore_end:
char = 'M'
if dihedral.type is not None:
scee = '%10.4f' % dihedral.type.scee
scnb = '%10.4f' % dihedral.type.scnb
k = '%10.4f' % dihedral.type.phi_k
per = '%10.4f' % dihedral.type.per
phase = '%10.4f' % dihedral.type.phase
else:
scee = scnb = k = per = phase = 'N/A'
else:
char = ' '
if dihedral.type is not None:
scee = '%10.4f' % dihedral.type.scee
scnb = '%10.4f' % dihedral.type.scnb
k = '%10.4f' % dihedral.type.phi_k
per = '%10.4f' % dihedral.type.per
phase = '%10.4f' % dihedral.type.phase
else:
scee = scnb = k = per = phase = 'N/A'
retstr.append('%1s %7d %4s (%4s) %7d %4s (%4s) %7d %4s '
'(%4s) %7d %4s (%4s) %10s %10s %10s %10s '
'%10s' % (char, atom1.idx+1, atom1.name,
atom1.type, atom2.idx+1, atom2.name, atom2.type,
atom3.idx+1, atom3.name, atom3.type, atom4.idx+1,
atom4.name, atom4.type, k, per, phase, scee, scnb)
)
if do_measured:
retstr.append(' %10.4f' % dihedral.measure())
if do_energy:
retstr.append(' %10.4f' % dihedral.energy())
retstr.append('\n')
return ''.join(retstr)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class setBond(Action):
"""
Changes (or adds a non-existent) bond in the topology file. Each mask must
select the same number of atoms, and a bond will be placed between the
atoms in mask1 and mask2 (one bond between atom1 from mask1 and atom1
from mask2 and another bond between atom2 from mask1 and atom2 from mask2,
etc.)
- <mask1> : Selection of first atoms in each bond
- <mask2> : Selection of second atoms in each bond
- <k> : Force constant (kcal/mol/A^2) in energy expression k(R-Req)^2
- <Req> : Equilibrium distance (A)
"""
usage = '<mask1> <mask2> <k> <Req>'
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.k = arg_list.get_next_float()
self.req = arg_list.get_next_float()
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
def __str__(self):
return ('Set a bond between %s and %s with k = %f kcal/(mol Angstrom'
'**2) and Req = %f Angstroms') % (self.mask1, self.mask2,
self.k, self.req)
[docs] def execute(self):
sel1 = self.mask1.Selection()
sel2 = self.mask2.Selection()
if sum(sel1) != sum(sel2):
raise SetParamError('setBond: Each mask must select the same '
'number of atoms!')
# If no atoms, nothing to do
if sum(sel1) == 0: return
# Create the new bond type
new_bnd_typ = BondType(self.k, self.req)
# Does that bond type exist in the list already? If it does, re-bind
# new_bnd to that bond type reference
exists = False
for bnd_typ in self.parm.bond_types:
if new_bnd_typ == bnd_typ:
new_bnd_typ = bnd_typ
exists = True
break
# If the bond is new, add it to the type list
if not exists:
self.parm.bond_types.append(new_bnd_typ)
new_bnd_typ.list = self.parm.bond_types
atnum1, atnum2 = -1, -1
# Loop through all of the selected atoms
for it in range(sum(sel1)):
# Collect the atoms involved
atnum1 = sel1.index(1, atnum1+1)
atnum2 = sel2.index(1, atnum2+1)
atm1 = self.parm.atoms[atnum1]
atm2 = self.parm.atoms[atnum2]
# See if the bond exists in the first place, and if so, replace its
# bond type with our new bond type (new_bnd)
if atm2 in atm1.bond_partners and atm1 in atm2.bond_partners:
for bond in atm1.bonds:
if atm2 in bond:
bond.type = new_bnd_typ
self.parm.bonds.changed = True
break
# Otherwise, it doesn't exist, so we just create a new one
else:
self.parm.bonds.append(Bond(atm1, atm2, new_bnd_typ))
# Make sure we update 1-4 exception handling if we created any rings
self.parm.update_dihedral_exclusions()
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class setAngle(Action):
"""
Changes (or adds a non-existent) angle in the topology file. Each mask must
select the same number of atoms, and an angle will be placed between the
atoms in mask1, mask2, and mask3 (one angle between atom1 from mask1, atom1
from mask2, and atom1 from mask3, another angle between atom2 from mask1,
atom2 from mask2, and atom2 from mask3, etc.)
- <mask1> : The selection of one of the end-atoms in each angle
- <mask2> : The selection of central atoms in each angle
- <mask3> : The selection of other end-atoms in each angle
- <k> : Force constant in kcal/mol/radians^2 in energy expression
k(THET - THETeq)^2
- <THETeq> : Equilibrium angle (in *degrees*)
"""
usage = '<mask1> <mask2> <mask3> <k> <THETeq>'
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.k = arg_list.get_next_float()
self.theteq = arg_list.get_next_float()
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask3 = AmberMask(self.parm, arg_list.get_next_mask())
def __str__(self):
return ('Set an angle between %s, %s and %s with k = %f kcal/(mol '
'rad**2) and THETeq = %f degrees') % (self.mask1, self.mask2,
self.mask3, self.k, self.theteq)
[docs] def execute(self):
sel1 = self.mask1.Selection()
sel2 = self.mask2.Selection()
sel3 = self.mask3.Selection()
if sum(sel1) != sum(sel2) or sum(sel1) != sum(sel3):
raise SetParamError('Each mask in setAngle must select the same '
'number of atoms!')
if sum(sel1) == 0: return
# Create the new angle type
new_ang_typ = AngleType(self.k, self.theteq)
# Does that angle type exist in the list already? If it does, re-bind
# new_ang to that angle type reference
exists = False
for ang_typ in self.parm.angle_types:
if new_ang_typ == ang_typ:
new_ang_typ = ang_typ
exists = True
break
# If the angle is new, add it to the type list
if not exists:
self.parm.angle_types.append(new_ang_typ)
new_ang_typ.list = self.parm.angle_types
atnum1, atnum2, atnum3 = -1, -1, -1
# Loop through all of the selections
for it in range(sum(sel1)):
# Collect the atoms involved
atnum1 = sel1.index(1, atnum1+1)
atnum2 = sel2.index(1, atnum2+1)
atnum3 = sel3.index(1, atnum3+1)
atm1 = self.parm.atoms[atnum1]
atm2 = self.parm.atoms[atnum2]
atm3 = self.parm.atoms[atnum3]
# See if the angle exists in the first place, and if so, replace its
# angle type with our new angle type (new_ang)
found = False
if atm1 in atm3.angle_partners:
for ang in atm1.angles:
if atm2 in ang and atm3 in ang:
ang.type = new_ang_typ
self.parm.angles.changed = True
found = True
break
# If not found, create a new angle
if not found:
self.parm.angles.append(Angle(atm1, atm2, atm3, new_ang_typ))
# Make sure we update 1-4 exception handling if we created any rings
self.parm.update_dihedral_exclusions()
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class addDihedral(Action):
"""
Adds a dihedral between mask1, mask2, mask3, and mask4. Each mask must
specify the same number of atoms, and the dihedral is defined around the
bond between atoms in mask 2 and 3. If each mask selects 2 atoms, for
instance, a dihedral will be placed around atom1 in mask 1, atom1 in mask 2,
atom1 in mask 3, and atom1 in mask 4. A second dihedral will be placed
around atom2 in mask 1, atom2 in mask 2, atom2 in mask 3, and atom2 in
mask4.
- <mask1> : Selection of one of the end-atoms for each dihedral
- <mask2> : Selection of the middle atom bonded to <mask1> and
<mask3> in each dihedral
- <mask3> : Selection of the other middle atom bonded to <mask2>
and <mask4>
- <mask4> : Selection of the other end-atom in each dihedral
- <phi_k> : Force constant in kcal/mol
- <per> : Periodicity
- <phase> : Torsion phase shift
- <scnb> : 1-4 Lennard-Jones scaling constant (default 2.0)
- <scee> : 1-4 electrostatic scaling constant (default 1.2)
- <type> : Type of dihedral, either "improper" or "normal".
Default is "normal"
"""
usage = ('<mask1> <mask2> <mask3> <mask4> <phi_k> <per> <phase> [<scee>] '
'[<scnb>] [type <type>]')
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.phi_k = arg_list.get_next_float()
self.per = arg_list.get_next_float()
self.phase = arg_list.get_next_float()
self.scee = arg_list.get_next_float(optional=True, default=1.2)
self.scnb = arg_list.get_next_float(optional=True, default=2.0)
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask3 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask4 = AmberMask(self.parm, arg_list.get_next_mask())
dihed_type = arg_list.get_key_string('type', 'normal')
if dihed_type.lower() == 'normal'[:len(dihed_type)]:
self.improper = False
self.type = 'a normal'
elif dihed_type.lower() == 'improper'[:len(dihed_type)]:
self.improper = True
self.type = 'an improper'
else:
raise InputError('addDihedral: type must be "normal" or "improper"')
def __str__(self):
return ('Set %s dihedral between %s, %s, %s, and %s with phi_k = %f '
'kcal/mol periodicity = %f phase = %f degrees scee = %f '
'scnb = %f' %
(self.type, self.mask1, self.mask2, self.mask3, self.mask4,
self.phi_k, self.per, self.phase, self.scee, self.scnb)
)
[docs] def execute(self):
sel1 = self.mask1.Selection()
sel2 = self.mask2.Selection()
sel3 = self.mask3.Selection()
sel4 = self.mask4.Selection()
if (sum(sel1) != sum(sel2) or sum(sel1) != sum(sel3) or
sum(sel1) != sum(sel4)):
raise SetParamError('addDihedral: Each mask must select the same '
'number of atoms!')
# If we do nothing, just return
if sum(sel1) == 0: return
# Create the new dihedral type
new_dih_typ = DihedralType(self.phi_k, self.per, self.phase, self.scee,
self.scnb)
exists = False
# Do not add a duplicate dihedral type
for dih_typ in self.parm.dihedral_types:
if new_dih_typ == dih_typ:
new_dih_typ = dih_typ
exists = True
break
if not exists:
self.parm.dihedral_types.append(new_dih_typ)
new_dih_typ.list = self.parm.dihedral_types
# Loop through all of the atoms
atnum1, atnum2, atnum3, atnum4 = -1, -1, -1, -1
for it in range(sum(sel1)):
# Collect the atoms involved
atnum1 = sel1.index(1, atnum1+1)
atnum2 = sel2.index(1, atnum2+1)
atnum3 = sel3.index(1, atnum3+1)
atnum4 = sel4.index(1, atnum4+1)
atm1 = self.parm.atoms[atnum1]
atm2 = self.parm.atoms[atnum2]
atm3 = self.parm.atoms[atnum3]
atm4 = self.parm.atoms[atnum4]
if (atm1 is atm2 or atm1 is atm3 or atm1 is atm4 or
atm2 is atm3 or atm2 is atm4 or atm3 is atm4):
raise SetParamError('addDihedral: Duplicate atoms found!')
# Determine if end-group interactions need to be ignored
ignore_end = (atm1 in atm4.bond_partners or
atm1 in atm4.angle_partners or
atm1 in atm4.dihedral_partners)
self.parm.dihedrals.append(
Dihedral(atm1, atm2, atm3, atm4, improper=self.improper,
ignore_end=ignore_end, type=new_dih_typ)
)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class addAtomicNumber(Action):
"""
Adds the atomic number of each atom to a new section titled "ATOMIC_NUMBER"
in the topology file. Elements are identified by the atomic mass found in
the MASS section of the topology files. Elements are matched by picking the
element whose average atomic mass in the periodic table is closest to each
atom, which should work appropriately for all isotopes of all atoms, except
possibly Tritium
"""
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
if self.parm.amoeba:
self.present = 'AMOEBA_ATOMIC_NUMBER' in self.parm.flag_list
else:
self.present = 'ATOMIC_NUMBER' in self.parm.flag_list
def __str__(self):
if self.present:
return 'ATOMIC_NUMBER already in [%s] -- Doing nothing.' % self.parm
return "Adding ATOMIC_NUMBER to [%s]" % self.parm
[docs] def execute(self):
if self.present: return
if self.parm.amoeba:
self.parm.add_flag('AMOEBA_ATOMIC_NUMBER', '10I8',
num_items=len(self.parm.atoms))
flag = 'AMOEBA_ATOMIC_NUMBER'
else:
self.parm.add_flag('ATOMIC_NUMBER', '10I8',
num_items=len(self.parm.atoms))
flag = 'ATOMIC_NUMBER'
for i, atm in enumerate(self.parm.atoms):
self.parm.parm_data[flag][i] = atm.atomic_number
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class deleteDihedral(Action):
"""
Deletes the dihedral around <mask2> and <mask3> in which the end-groups are
<mask1> and <mask4>. For multi-term dihedrals, it removes each term.
"""
usage = '<mask1> <mask2> <mask3> <mask4>'
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask3 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask4 = AmberMask(self.parm, arg_list.get_next_mask())
if (sum(self.mask1.Selection()) != sum(self.mask2.Selection()) or
sum(self.mask1.Selection()) != sum(self.mask3.Selection()) or
sum(self.mask1.Selection()) != sum(self.mask4.Selection())):
raise DeleteDihedralError('All masks must select the same number '
'of atoms!. They selected %d, %d, %d, '
'and %d, respectively' % (
sum(self.mask1.Selection()), sum(self.mask2.Selection()),
sum(self.mask3.Selection()), sum(self.mask4.Selection()))
)
def __str__(self):
if sum(self.mask1.Selection()) == 0:
return 'No specified dihedrals to delete'
return ('Deleting dihedral terms involving [%s]-[%s]-[%s]-[%s]'
' (At most %d total, distinct, dihedrals)' %
(self.mask1, self.mask2, self.mask3, self.mask4,
sum(self.mask1.Selection()))
)
[docs] def execute(self):
""" Returns the total number of dihedrals deleted """
sel1, sel2 = self.mask1.Selection(), self.mask2.Selection()
sel3, sel4 = self.mask3.Selection(), self.mask4.Selection()
# Bail out if we're deleting nothing
if sum(sel1) == 0: return
# Keep track of the dihedrals we want to delete from each
# dihedral list (dihedrals_inc_h, dihedrals_without_h)
deleting_dihedrals = []
# We have already checked that they are the same number of atoms
# Now, loop through the atoms and see if any dihedrals match that spec
atnum1 = atnum2 = atnum3 = atnum4 = -1
total_diheds = 0
for i in range(sum(sel1)):
# Collect the atoms involved
atnum1 = sel1.index(1, atnum1+1)
atnum2 = sel2.index(1, atnum2+1)
atnum3 = sel3.index(1, atnum3+1)
atnum4 = sel4.index(1, atnum4+1)
atm1 = self.parm.atoms[atnum1]
atm2 = self.parm.atoms[atnum2]
atm3 = self.parm.atoms[atnum3]
atm4 = self.parm.atoms[atnum4]
# Make sure none of the indices are the same
if (atm1 is atm2 or atm1 is atm3 or atm1 is atm4 or
atm2 is atm3 or atm2 is atm4 or atm3 is atm4):
warnings.warn('Skipping %d-%d-%d-%d dihedral deletion -- '
'duplicate atoms!' %
(atnum1, atnum2, atnum3, atnum4),
SeriousParmWarning)
continue # pragma: no cover
# Now search through our dihedral list to see which indexes (if any)
# we have to remove. Keep tabs of them so we can pop them in reverse
# order (so we don't have to re-figure indices) afterwards
proposed_dihedral = (atnum1, atnum2, atnum3, atnum4)
for j, dihed in enumerate(self.parm.dihedrals):
if dihed.same_atoms(proposed_dihedral):
total_diheds += 1
deleting_dihedrals.append(j)
if not deleting_dihedrals:
return 0
# At this point, we've collected all of our dihedrals, now sort them
deleting_dihedrals.sort()
# deleting_dihedrals now contains all of our dihedral indexes
while deleting_dihedrals:
idx = deleting_dihedrals.pop()
self.parm.dihedrals[idx].delete()
del self.parm.dihedrals[idx]
return total_diheds
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class printLJMatrix(Action):
"""
This function prints out how every atom type interacts with the atom type(s)
in <mask>. The atom types are printed as all type names that have at least
one atom with the Lennard Jones index given in square brackets at the end.
Alternatively, you can request a particular atom type index
"""
usage = '<mask>|<index>'
strictly_supported = (AmberParm, ChamberParm)
[docs] def init(self, arg_list):
self.idx = arg_list.get_next_int(optional=True)
self.mask = None
if self.idx is None:
self.mask = AmberMask(self.parm, arg_list.get_next_mask())
def __str__(self):
return self.__repr__()
def __repr__(self):
has_1264 = (hasattr(self.parm, 'parm_data') and
'LENNARD_JONES_CCOEF' in self.parm.parm_data)
ntypes = self.parm.ptr('NTYPES')
ret_str = []
if self.idx is not None:
sel = [0 for i in self.parm.atoms]
for i, atom in enumerate(self.parm.atoms):
if atom.nb_idx == self.idx: sel[i] = 1
else:
sel = self.mask.Selection()
# If we selected no atoms, bail out
if sum(sel) == 0: return 'No atom types selected'
# Figure out which types correspond to which names
typenames = [set() for i in range(self.parm.ptr('NTYPES'))]
for i, atom in enumerate(self.parm.atoms):
typenames[atom.nb_idx-1].add(atom.type)
# Otherwise, collect our list of atom types that we selected
sel_types = set()
for i, atom in enumerate(self.parm.atoms):
if sel[i] == 0: continue
sel_types.add(atom.nb_idx)
sel_types = sorted(list(sel_types)) # sort the atom types
# Convert all of the typenames into strings, then find the longest one so
# we can properly format the string
maxlen = 0
for i, names in enumerate(typenames):
typenames[i] = ','.join(sorted(list(names))) + ' [%d]' % (i+1)
maxlen = max(maxlen, len(typenames[i]))
if has_1264:
fmt = '\n%%%ds %%%ds %%15s %%15s %%15s %%10s %%10s' % (maxlen, maxlen)
args = ('Atom Type 1', 'Atom Type 2', 'A coefficient',
'B coefficient', 'C coefficient', 'R i,j', 'Eps i,j')
else:
fmt = '\n%%%ds %%%ds %%15s %%15s %%10s %%10s' % (maxlen, maxlen)
args = ('Atom Type 1', 'Atom Type 2', 'A coefficient',
'B coefficient', 'R i,j', 'Eps i,j')
ret_str.append(fmt % args)
ret_str.extend(['\n','-'*len(ret_str[-1]), '\n'])
for ty in sel_types:
for ty2 in range(1,ntypes+1):
type1, type2 = min(ty, ty2), max(ty, ty2)
idx = self.parm.parm_data['NONBONDED_PARM_INDEX'][
ntypes*(type1-1)+type2-1]
acoef = self.parm.parm_data['LENNARD_JONES_ACOEF'][idx-1]
bcoef = self.parm.parm_data['LENNARD_JONES_BCOEF'][idx-1]
if has_1264:
ccoef = self.parm.parm_data['LENNARD_JONES_CCOEF'][idx-1]
if bcoef == 0 or acoef == 0:
rij = eij = 0.0
else:
rij = (2 * acoef / bcoef) ** (1 / 6)
eij = (bcoef * bcoef / (4 * acoef))
if has_1264:
ret_str.append('%%%ds %%%ds %%15.6f %%15.6f %%15.6f %%10.6f %%10.6f\n' %
(maxlen, maxlen) %
(typenames[type1-1], typenames[type2-1],
acoef, bcoef, ccoef, rij, eij)
)
else:
ret_str.append('%%%ds %%%ds %%15.6f %%15.6f %%10.6f %%10.6f\n' %
(maxlen, maxlen) %
(typenames[type1-1], typenames[type2-1],
acoef, bcoef, rij, eij)
)
return ''.join(ret_str)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class tiMerge(Action):
"""
Merges molecules removing redundant bonding terms. Input amber masks
corresponding to molecules 1/2 <mol1mask>/<mol2mask>, and the soft core
atoms in each molecule as <scmask1>/<scmask2>. The input topology can be
created using leap, with the two molecules to be merged adjacent to each
other in residue number. This improves the efficiency for pmemd TI when only
part of a molecule is being perturbed.
<scmask1/2N> are for softcore molecules that are not going to be merged.
These options will just add these atoms to the timask output, correcting for
any changes in atom number.
This can also be used for non-softcore simulations, where
<scmask1>/<scmask2> represent the perturbed atoms. The output will give the
scmask1/scmask2 flags, which can just be ignored.
<tol> is the tolerence to use when matching coordinates (default 0.01).
This is used when the atoms in molecules 1/2 are not in the same order and
for checking the input coordinates.
"""
usage = '<mol1mask> <mol2mask> <scmask1> <scmask2> [<scmask1N>] [<scmask2N>] [tol <tol>]'
strictly_supported = (AmberParm, ChamberParm)
output = sys.stdout
[docs] def init(self, arg_list):
self.tol = arg_list.get_key_float('tol', 0.01)
self.molmask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.molmask2 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
self.sc_mask1 = ''
self.sc_mask2 = ''
molmask1N = arg_list.get_next_mask(optional=True)
if molmask1N is not None:
self.molmask1N = AmberMask(self.parm, molmask1N)
else:
self.molmask1N = None
molmask2N = arg_list.get_next_mask(optional=True)
if molmask2N is not None:
self.molmask2N = AmberMask(self.parm, molmask2N)
else:
self.molmask2N = None
def __str__(self):
return ('Merging molecules [%s] [%s] with sc mask [%s] [%s]' % (
self.molmask1, self.molmask2, self.mask1, self.mask2 )
)
[docs] def execute(self):
sel1 = self.mask1.Selection()
sel2 = self.mask2.Selection()
molsel1 = self.molmask1.Selection()
molsel2 = self.molmask2.Selection()
natom = len(self.parm.atoms)
if self.molmask1N is not None:
molsel1N = self.molmask1N.Selection()
else:
molsel1N = [0 for i in range(natom)]
if self.molmask2N is not None:
molsel2N = self.molmask2N.Selection()
else:
molsel2N = [0 for i in range(natom)]
for i in range(natom):
if sel1[i] and not molsel1[i]:
raise TiMergeError('scmask1 must be a subset of mol1mask.')
if sel2[i] and not molsel2[i]:
raise TiMergeError('scmask2 must be a subset of mol2mask.')
if molsel1[i] and molsel2[i]:
raise TiMergeError('mol1mask can not overlap with mol2mask.')
if i < natom-1:
if molsel1[i] and not (molsel1[i+1] or molsel2[i+1]):
raise TiMergeError('mol1mask and mol2mask must be adjacent in topology. '
'Recreate topology with the molecules to be merged adjacent '
'in the PDB file.')
# self.parm.coordinates is a property that manipulates a numpy array, so get a copy of it
# here so we do not wind up accessing that property in an inner loop... EXPENSIVE!
coordinates = self.parm.coordinates
if coordinates is None:
raise TiMergeError('Load coordinates before merging topology.')
# we now have enough info to remap the atom indicies if an atom in
# molsel2 has no overlap (dihedrals, angles, bonds) with sel2 then we
# can just delete it (it is redundant).
keep_mask = [0 for i in range(natom)]
for i in range(natom):
if molsel2[i]:
for j in range(natom):
if sel2[j]:
atm1 = self.parm.atoms[i]
atm2 = self.parm.atoms[j]
if (atm1 in atm2.bond_partners or atm1 in atm2.angle_partners or
atm1 in atm2.dihedral_partners):
keep_mask[i] = 1
nremove = sum(molsel2) - sum(sel2)
# We use an ambmask here to minimize changes to the code and
# not introduce extra maintenance issues
remove_mask = []
# remove_map[old_atm_idx] = new_atm_idx
remove_map = [0 for i in range(natom)]
new_atm_idx = 0
for i in range(natom):
if molsel2[i] == 1 and sel2[i] == 0:
remove_mask.append('%d' % (i+1))
else:
remove_map[i] = new_atm_idx
new_atm_idx += 1
remove_str = '@' + ','.join(remove_mask)
# However, if there is overlap, we need to re-index the atoms involved
# Create a map from molsel2 to molsel1 excluding sel2/sel1 respectively
mol1common = []
mol2common = []
for i in range(natom):
if molsel1[i] == 1 and sel1[i] == 0:
mol1common.append(i)
for i in range(natom):
if molsel2[i] == 1 and sel2[i] == 0:
mol2common.append(i)
if len(mol1common) != len(mol2common):
raise TiMergeError('The number of nonsoftcore atoms in mol1mask and mol2mask must be '
'the same.')
mol2common_sort = []
# reorder mol2common so that it matches mol1common
for i in range(len(mol1common)):
atm_i = mol1common[i]
for j in range(len(mol2common)):
atm_j = mol2common[j]
diff = coordinates[atm_i] - coordinates[atm_j]
if (np.abs(diff) < self.tol).sum() == 3:
mol2common_sort.append(atm_j)
mol2common = mol2common_sort
# check again if we didn't match all coords
if len(mol1common) != len(mol2common):
raise TiMergeError('The number of nonsoftcore atoms in mol1mask and mol2mask must be '
'the same. Check the masks. If these look correct try using a '
'larger tolerance.')
for i in range(len(mol1common)):
atm_i = mol1common[i]
atm_j = mol2common[i]
diff = coordinates[atm_i] - coordinates[atm_j]
if (np.abs(diff) > self.tol).any():
raise TiMergeError('Common (nonsoftcore) atoms must have the ' # pragma: no cover
'same coordinates.')
for j in range(natom):
if keep_mask[j] == 1 and sel2[j] == 0:
atm = self.parm.atoms[j]
idx = mol1common[mol2common.index(j)]
atm_new = self.parm.atoms[idx]
# What is happening here???
for k in range(natom):
if sel2[k]:
atm2 = self.parm.atoms[k]
# update partners -- the exclusion list will be updated
# when the file is written out
if atm in atm2.bond_partners:
atm._bond_partners.remove(atm2)
atm2._bond_partners.remove(atm)
atm2.bond_to(atm_new)
if atm in atm2.angle_partners:
atm._angle_partners.remove(atm2)
atm2._angle_partners.remove(atm)
atm2.angle_to(atm_new)
if atm in atm2.dihedral_partners:
atm._dihedral_partners.remove(atm2)
atm2._dihedral_partners.remove(atm)
atm2.dihedral_to(atm_new)
# Now go through each array re-indexing the atoms
# Check to make sure that this is a bond/angle/dihed
# involving the common atom j and the softcore atom k
for bond in self.parm.bonds:
if (bond.atom1.idx == j and bond.atom2.idx == k):
bond.atom1 = atm_new
elif (bond.atom2.idx == j and bond.atom1.idx == k):
bond.atom2 = atm_new
for angle in self.parm.angles:
if angle.atom1.idx == j:
if angle.atom2.idx == k or angle.atom3.idx == k:
angle.atom1 = atm_new
elif angle.atom2.idx == j:
if angle.atom1.idx == k or angle.atom3.idx == k:
angle.atom2 = atm_new
elif angle.atom3.idx == j:
if angle.atom1.idx == k or angle.atom2.idx == k:
angle.atom3 = atm_new
for dihed in self.parm.dihedrals:
if dihed.atom1.idx == j:
if (dihed.atom2.idx == k or dihed.atom3.idx == k
or dihed.atom4.idx == k):
dihed.atom1 = atm_new
elif dihed.atom2.idx == j:
if (dihed.atom1.idx == k or dihed.atom3.idx == k
or dihed.atom4.idx == k):
dihed.atom2 = atm_new
elif dihed.atom3.idx == j:
if (dihed.atom1.idx == k or dihed.atom2.idx == k
or dihed.atom4.idx == k):
dihed.atom3 = atm_new
elif dihed.atom4.idx == j:
if (dihed.atom1.idx == k or dihed.atom2.idx == k
or dihed.atom3.idx == k):
dihed.atom4 = atm_new
for imp in self.parm.impropers:
if imp.atom1.idx == j:
if (imp.atom2.idx == k or imp.atom3.idx == k
or imp.atom4.idx == k):
imp.atom1 = atm_new
elif imp.atom2.idx == j:
if (imp.atom1.idx == k or imp.atom3.idx == k
or imp.atom4.idx == k):
imp.atom2 = atm_new
elif imp.atom3.idx == j:
if (imp.atom1.idx == k or imp.atom2.idx == k
or imp.atom4.idx == k):
imp.atom3 = atm_new
elif imp.atom4.idx == j:
if (imp.atom1.idx == k or imp.atom2.idx == k
or imp.atom3.idx == k):
imp.atom4 = atm_new
for cmap in self.parm.cmaps:
if cmap.atom1.idx == j:
if (cmap.atom2.idx == k or cmap.atom3.idx == k
or cmap.atom4.idx == k or cmap.atom5.idx == k):
cmap.atom1 = atm_new
elif cmap.atom2.idx == j:
if (cmap.atom1.idx == k or cmap.atom3.idx == k
or cmap.atom4.idx == k or cmap.atom5.idx == k):
cmap.atom2 = atm_new
elif cmap.atom3.idx == j:
if (cmap.atom1.idx == k or cmap.atom2.idx == k
or cmap.atom4.idx == k or cmap.atom5.idx == k):
cmap.atom3 = atm_new
elif cmap.atom4.idx == j:
if (cmap.atom1.idx == k or cmap.atom2.idx == k
or cmap.atom3.idx == k or cmap.atom5.idx == k):
cmap.atom4 = atm_new
elif cmap.atom5.idx == j:
if (cmap.atom1.idx == k or cmap.atom2.idx == k
or cmap.atom3.idx == k or cmap.atom4.idx == k):
cmap.atom5 = atm_new
self.parm.atoms.changed = True
if nremove > 0:
self.parm.strip(remove_str)
new_sc_atm1 = []
new_sc_atm2 = []
new_sc_atm1_int = []
new_sc_atm2_int = []
for i in range(natom):
if sel1[i] or molsel1N[i]:
new_sc_atm1_int.append(remove_map[i])
new_sc_atm1.append('%d' % (remove_map[i]+1))
elif sel2[i] or molsel2N[i]:
new_sc_atm2_int.append(remove_map[i])
new_sc_atm2.append('%d' % (remove_map[i]+1))
# Do not allow definition where dihedrals cross through the softcore
# region. This is generally breaking a ring (and possibly other cases),
# and can cause problems with the 1-4 nonbonded calculations.
# This can be worked-around:
# Define your softcore region so that it includes the ring.
for dihed in self.parm.dihedrals:
# skip impropers, these are not used to define 1-4 interactions
# so these can cross through the softcore region
if dihed.improper: continue
atmi = dihed.atom1.idx
atmj = dihed.atom2.idx
atmk = dihed.atom3.idx
atml = dihed.atom4.idx
if (atmj in new_sc_atm1_int or atmk in new_sc_atm1_int or atmj in new_sc_atm2_int or
atmk in new_sc_atm2_int): #dihedral includes sc atoms
#endpoint atoms are not softcore
#we are crossing through the softcore region
if (atmi not in new_sc_atm1_int and atmi not in new_sc_atm2_int and
atml not in new_sc_atm1_int and atml not in new_sc_atm2_int):
raise TiMergeError( # pragma: no cover
'Cannot have dihedral cross through softcore '
'region. (DIHED : %d %d %d %d). Usually this means '
'you have defined the softcore region in a way '
'that breaks a ring. Try redefining your softcore '
'region to include the ring or at least three '
'consecutive atoms.' % (atmi+1, atmj+1, abs(atmk)+1, abs(atml)+1)
)
self.sc_mask1 = '@' + ','.join(new_sc_atm1)
self.sc_mask2 = '@' + ','.join(new_sc_atm2)
ret_str = ("Merging molecules %s and %s into the same molecule.\n"
% (self.molmask1, self.molmask2))
ret_str2 = ("Use softcore mask:\ntimask1=\'%s\',\ntimask2=\'%s\',"
% (self.sc_mask1, self.sc_mask2))
ret_str3 = ("\nscmask1=\'%s\',\nscmask2=\'%s\',"
% (self.sc_mask1, self.sc_mask2))
self.output.write('%s%s%s\n' % (ret_str, ret_str2, ret_str3))
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class source(Action):
"""
Sources a file with a list of parmed commands
"""
usage = '<file>'
needs_parm = False
[docs] def init(self, arg_list):
self.filename = arg_list.get_next_string()
def __str__(self):
return 'Sourcing %s' % self.filename
[docs] def execute(self):
"""
This is a no-op, since a separate command interpreter for this file is
launched inside parmed_cmd.py
"""
pass # pragma: no cover
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class parm(Action):
"""
Either adds a new parm to the list of available parms to edit in ParmEd, or
it sets a new 'active' parm to edit by default with new commands
"""
usage = ('<filename> [<filename> [<filename> ...]] || parm copy <filename>|'
'<index> || parm select <filename>|<index>')
needs_parm = False
[docs] def init(self, arg_list):
from glob import glob
self.new_active_parm = arg_list.get_key_string('select', None)
self.copied_parm = arg_list.get_key_string('copy', None)
self.new_parm = None
# Add as many parms as we want with one command, and support globbing
new_parm = len(arg_list.unmarked()) or None
# Make sure we specified only one operating mode
if (self.new_active_parm is None and self.copied_parm is None and
new_parm is None):
raise ParmError('Improper usage of `parm\' command')
opts = (self.new_active_parm, self.copied_parm, new_parm)
nnone = 0
for opt in opts:
if opt is None:
nnone += 1
if nnone != 2:
raise ParmError('Improper usage of `parm\' -- choose one behavior')
# Now get our indexes for copied and new parms if they are integers,
# otherwise we index based on filename
if self.new_active_parm is not None:
try:
self.new_active_parm = int(self.new_active_parm)
except ValueError:
pass
elif self.copied_parm is not None:
try:
self.copied_parm = int(self.copied_parm)
except ValueError:
pass
else:
self.new_parm = []
new_parm = arg_list.get_next_string(optional=True)
while new_parm is not None:
listparms = glob(new_parm)
if not listparms:
raise ParmFileNotFound('No files matching %s' % new_parm)
self.new_parm.extend(glob(new_parm))
new_parm = arg_list.get_next_string(optional=True)
assert self.new_parm, 'No matching parm files? should not happen'
def __str__(self):
if self.new_active_parm is not None:
try:
idx = self.parm_list.index(self.new_active_parm)
except IndexError:
return '%s not in parm list. Doing nothing' % self.new_active_parm
return 'Setting new active parm [%s]' % self.parm_list[idx]
elif self.new_parm is not None:
return ('Adding prmtop %s to parm list. %s is the active parm.' %
(', '.join(self.new_parm), self.new_parm[-1]))
elif self.copied_parm is not None:
try:
idx = self.parm_list.index(self.copied_parm)
except IndexError:
return '%s not in parm list. Doing nothing' % self.new_active_parm
return ("Copying prmtop %s to parm list. %s's copy is the active parm." %
(self.parm_list[idx], self.parm_list[idx]))
raise AssertionError('Should not be here')
[docs] def execute(self):
""" Either set the new active parm or add the new parm """
from copy import copy
if self.new_parm is not None:
# Add the new parm
for new_parm in self.new_parm:
try:
self.parm_list.add_parm(new_parm)
except IOError:
warnings.warn('Could not open %s for reading' % new_parm, SeriousParmWarning)
elif self.new_active_parm is not None:
try:
self.parm_list.set_new_active(self.new_active_parm)
except IndexError:
warnings.warn('%s is not in the parm list!' % self.new_active_parm,
SeriousParmWarning)
elif self.copied_parm is not None:
try:
self.parm_list.add_parm(copy(self.parm_list[self.copied_parm]))
except IndexError:
warnings.warn('%s is not in the parm list!' % self.copied_parm, SeriousParmWarning)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class ls(Action):
"""
Lists directory contents. Like UNIX 'ls'
"""
usage = '[Unix ls options]'
needs_parm = False
[docs] def init(self, arg_list):
from glob import glob
self.args = []
# Process the argument list to mimic the real ls as much as possible
while True:
try:
arg = arg_list.get_next_string()
if not arg.startswith('-'):
# Glob this argument
globarg = glob(arg)
if len(globarg) > 0:
self.args.extend(globarg)
else:
self.args.append(arg)
else:
self.args.append(arg)
except NoArgument:
break
def __str__(self):
from subprocess import Popen, PIPE
process = Popen(['/bin/ls', '-C'] + self.args, stdout=PIPE, stderr=PIPE)
out, err = process.communicate('')
process.wait()
return (out + err).decode('UTF-8')
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class cd(Action):
"""
Changes to a new directory like UNIX 'cd'
"""
usage = '<directory>'
needs_parm = False
[docs] def init(self, arg_list):
from glob import glob
from os.path import expanduser, expandvars
mydir = expanduser(expandvars(arg_list.get_next_string()))
self.directory = glob(mydir)
def __str__(self):
if len(self.directory) != 1:
return 'Change directory failed'
if not os.path.isdir(self.directory[0]):
return '%s does not exist. cd failed.' % self.directory[0]
return 'New working directory: %s' % self.directory[0]
[docs] def execute(self):
if len(self.directory) < 1:
warnings.warn('No recognized directories given to cd',
SeriousParmWarning)
return # pragma: no cover
elif len(self.directory) > 1:
warnings.warn('More than one file/directory given to cd',
SeriousParmWarning)
return # pragma: no cover
if not os.path.isdir(self.directory[0]):
warnings.warn('%s is not a directory' % self.directory[0],
SeriousParmWarning)
return # pragma: no cover
# If we've gotten this far, go ahead and change to the directory
os.chdir(self.directory[0])
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class listParms(Action):
"""
Lists all of the loaded topology files
"""
needs_parm = False
[docs] def init(self, arg_list):
pass
def __repr__(self):
if self.parm_list.empty():
return "No topology files are loaded"
retstr = 'Loaded topology files:'
for i, parm in enumerate(self.parm_list):
retstr += '\n[%d]\t%s' % (i, parm)
if parm is self.parm_list.parm:
retstr += ' (active)'
return retstr
__str__ = __repr__
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class interpolate(Action):
"""
Interpolates between two topology files (VDW and electrostatic terms only).
If [eleconly] is present, only the charges will be interpolated. <nparm> is
the number of 'interpolated' topology files you want (in addition to the two
end-points). The second prmtop must be specified if there are more than 2
parms currently loaded in the ParmEd parm list. startnum has no effect on
the generated prmtops, but it can be used to control the names of the
outputted topology files.
- <nparm> : Number of topology files that will be generated
- <other_parm> : The other parm object used in interpolation if more
than 2 parms are present (first parm is active one)
- eleconly : Only do charge interpolation
- <prefix> : Generated parm objects will be written as <prefix>.#, where
# starts from <num> and increases by 1
- <num> : Starting number for file names (see <prefix> above)
"""
usage = ('<nparm> [parm2 <other_parm>] [eleconly] [prefix <prefix>] '
'[startnum <num>]')
strictly_supported = (AmberParm, ChamberParm)
[docs] def init(self, arg_list):
# Make sure we have at least 2 prmtops
if len(self.parm_list) < 2:
raise NonexistentParm('Must have 2 topology files to interpolate!')
parm2 = arg_list.get_key_string('parm2', None)
if parm2 is None and len(self.parm_list) == 2:
self.parm2 = self.parm_list[1]
if self.parm_list[0] is not self.parm:
self.parm2 = self.parm_list[0]
elif parm2 is None:
raise AmbiguousParmError('You must identify parm2 if more than 2 '
'parm instances exist!')
else:
try:
parm2 = int(parm2)
except ValueError:
pass
self.parm2 = self.parm_list[parm2]
self.startnum = arg_list.get_key_int('startnum', 1)
self.prefix = arg_list.get_key_string('prefix', str(self.parm))
self.eleconly = arg_list.has_key('eleconly')
self.nparm = arg_list.get_next_int()
if self.nparm <= 0:
raise ArgumentError('Must have >= 1 prmtop')
self.diff_vdw = False
self._check_parms()
def __str__(self):
extra = ''
if self.eleconly and self.diff_vdw:
extra = ' [only interpolating charges]'
return 'Creating %d interpolated prmtops between %s and %s' % (
self.nparm, self.parm, self.parm2) + extra
def _check_parms(self):
""" Makes sure that the atoms in both parms are all the same """
parm1, parm2 = self.parm, self.parm2
if len(parm1.atoms) != len(parm2.atoms):
raise IncompatibleParmsError('%s and %s have different #s of '
'atoms!' % (parm1, parm2))
ndiff = 0
for atom1, atom2 in zip(parm1.atoms, parm2.atoms):
if atom1.name != atom2.name:
ndiff += 1
if ndiff > 0:
warnings.warn('%d atoms have different names b/w %s and %s' %
(ndiff, parm1, parm2), SeriousParmWarning)
for atm1, atm2 in zip(parm1.atoms, parm2.atoms):
i1, i2 = atm1.nb_idx-1, atm2.nb_idx-1
if ((parm1.LJ_radius[i1] != parm2.LJ_radius[i2]) or
(parm1.LJ_depth[i1] != parm2.LJ_radius[i2])):
self.diff_vdw = True
[docs] def execute(self):
""" Interpolates the prmtops """
if self.diff_vdw and not self.eleconly:
raise NotImplementedError('No support for scaling vdW '
'parameters yet!')
parm1, parm2 = self.parm, self.parm2
# Original charges for parm 1
orig_chg1 = parm1.parm_data['CHARGE']
chg1 = np.array(parm1.parm_data['CHARGE'])
chg2 = np.array(parm2.parm_data['CHARGE'])
diff = chg2 - chg1
diff *= 1 / (self.nparm + 1)
for i in range(self.nparm):
new_chg = chg1 + diff * (i + 1)
parm1.parm_data['CHARGE'] = new_chg.tolist()
parm1.load_atom_info()
newname = '%s.%d' % (self.prefix, i+self.startnum)
print('Printing %s' % newname)
parm1.write_parm(newname)
# Restore the original charges
parm1.parm_data['CHARGE'] = orig_chg1
parm1.load_atom_info()
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
def _split_range(chunksize, start, stop):
'''split a given range to n_chunks. taken from pytraj.
Examples
--------
[(0, 3), (3, 6), (6, 10)]
'''
n_chunks = (stop - start)//chunksize
if ((stop - start) % chunksize ) != 0:
n_chunks += 1
for i in range(n_chunks):
if i < n_chunks - 1:
_stop = start + (i + 1) * chunksize
else:
_stop = stop
yield start + i * chunksize, _stop
def _reformat_long_sentence(long_sentence, title, n_words=6):
words = long_sentence.split(', ')
empty = "\n" + " " * len(title)
lines = [words[slice(*idx)] for idx in _split_range(n_words, 0, len(words))]
sentences = [', '.join(line) for line in lines]
sentences[0] = title[:] + sentences[0]
return empty.join(sentences)
[docs]class summary(Action):
"""
Prints out a summary of prmtop contents
"""
[docs] def init(self, arg_list):
pass
def __str__(self):
return self.__repr__()
def __repr__(self):
""" Collect statistics """
nnuc = namin = ncion = naion = nwat = nunk = 0
for res in self.parm.residues:
if RNAResidue.has(res.name) or DNAResidue.has(res.name):
nnuc += 1
elif res.name in AminoAcidResidue._all_residues_by_abbr:
namin += 1
elif res.name in SOLVENT_NAMES:
nwat += 1
elif res.name in ANION_NAMES:
naion += 1
elif res.name in CATION_NAMES:
ncion += 1
else:
nunk += 1
tmass = sum(atom.mass for atom in self.parm.atoms)
tchg = sum(atom.charge for atom in self.parm.atoms)
retval = [('Amino Acid Residues: %d\n'
'Nucleic Acid Residues: %d\n'
'Number of cations: %d\n'
'Number of anions: %d\n'
'Num. of solvent mols: %d\n'
'Num. of unknown res: %d\n'
'Total charge (e-): %.4f\n'
'Total mass (amu): %.4f\n'
'Number of atoms: %d\n'
'Number of residues: %d' %
(namin, nnuc, ncion, naion, nwat, nunk, tchg, tmass,
len(self.parm.atoms), len(self.parm.residues))
)]
rset = ", ".join(sorted(set(res.name for res in self.parm.residues)))
retval.append(
_reformat_long_sentence(rset, 'Residue set: ',
n_words=7)
)
rcount = ','.join('%s: %d' % (x, y)
for x, y in iteritems(
Counter(res.name for res in self.parm.residues)
)
)
retval.append(
_reformat_long_sentence(', '.join((sorted(rcount.split(',')))),
'Residue count: ', n_words=7)
)
if self.parm.box is not None and set(self.parm.box[3:]) == set([90]):
a, b, c = self.parm.box[:3]
v = a * b * c
# Get the total volume (and density) of orthorhombic box
retval.append('System volume (ang^3): %.2f\n'
'System density (g/mL): %f' %
(v, tmass / (v * 0.602204))
)
elif self.parm.box is not None:
# General triclinic cell
a, b, c, alpha, beta, gamma = self.parm.box[:]
# Convert to radians
cosa = math.cos(alpha * math.pi / 180)
cosb = math.cos(beta * math.pi / 180)
cosg = math.cos(gamma * math.pi / 180)
v = a * b * c * math.sqrt(1 - cosa*cosa - cosb*cosb - cosg*cosg +
2 * cosa*cosb*cosg)
retval.append('System volume (ang^3): %.2f\n'
'System density (g/mL): %f' %
(v, tmass / (v * 0.602204))
)
return '%s\n' % '\n'.join(retval)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class scale(Action):
"""
Multiplies all values in a particular section of an Amber prmtop by a scalar
factor
"""
usage = '<FLAG> <factor>'
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
self.flag = arg_list.get_next_string().upper()
self.factor = arg_list.get_next_float()
if not self.flag in self.parm.flag_list:
raise ArgumentError('%s is not a valid prmtop flag!' % self.flag)
def __str__(self):
return 'Scaling data in %s by %f' % (self.flag, self.factor)
[docs] def execute(self):
try:
for i in range(len(self.parm.parm_data[self.flag])):
self.parm.parm_data[self.flag][i] *= self.factor
except TypeError:
raise ArgumentError('Cannot scale data in %%FLAG %s' % self.flag)
else:
self.parm.load_structure()
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class lmod(Action):
"""
Adjusts Lennard Jones parameters to work with the LMOD code in Amber
(changes LJ A coefficients that are 0 to 1000).
"""
strictly_supported = (AmberParm, ChamberParm)
[docs] def init(self, arg_list):
pass
def __str__(self):
return ('Making adjustments for LMOD (LJ A-coef. for H atoms bonded '
'to O)')
[docs] def execute(self):
for i, val in enumerate(self.parm.parm_data['LENNARD_JONES_ACOEF']):
self.parm.parm_data['LENNARD_JONES_ACOEF'][i] = val or 1000.0
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class addPDB(Action):
"""
Adds PDB information to new flags in an Amber topology file to enable
analyses based on the original residue information in the PDB file,
<filename>. It adds the flags:
Residue Properties
------------------
- RESIDUE_CHAINID: The chain ID of each residue (* if LEaP added it)
- RESIDUE_ICODE: Insertion code, if it exists
- RESIDUE_NUMBER: Original residue serial number in the PDB
Atom Properties
---------------
- ATOM_ELEMENT: Atomic element (redundant now, not printed by default)
- ATOM_OCCUPANCY: The occupancy of each atom
- ATOM_BFACTOR: The temperature factor of each atom
- ATOM_NUMBER: The original atom serial number in the PDB
The 'strict' keyword turns residue mismatches (NOT solvent) into errors
The 'elem' keyword will force printing of the element names.
The 'allicodes' keyword forces insertion codes to be printed, even if every
one will be blank (so that parsers can count on that section existing).
Residues *not* in the PDB will be assigned a CHAINID of '*' and
RESIDUE_NUMBER of 0. Any occupancy or temperature (B) factor that is not
present in the input PDB file is assigned a number of 0
Historical info:
This action is based on, and reproduces the key results of, the
historical 'add_pdb' program. It is a bit more flexible, though.
"""
usage = '<filename> [elem] [strict] [allicodes]'
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
self.pdbfile = arg_list.get_next_string()
self.elements = arg_list.has_key('elem')
self.printicodes = arg_list.has_key('allicodes')
if arg_list.has_key('strict'):
warnings.filterwarnings('error', category=AddPDBWarning)
else:
warnings.filterwarnings('always', category=AddPDBWarning)
self.pdbpresent = ('RESIDUE_NUMBER' in self.parm.flag_list or
'RESIDUE_CHAINID' in self.parm.flag_list or
'RESIDUE_ICODE' in self.parm.flag_list or
'ATOM_ELEMENT' in self.parm.flag_list or
'ATOM_OCCUPANCY' in self.parm.flag_list or
'ATOM_BFACTOR' in self.parm.flag_list or
'ATOM_NUMBER' in self.parm.flag_list
)
def __str__(self):
if self.pdbpresent:
return 'PDB information already present in %s. Doing nothing'
retstr = 'Adding PDB information from %s' % self.pdbfile
if self.elements: retstr += '\n\t[printing elements from prmtop]'
return retstr
[docs] def execute(self):
from parmed import read_PDB
if self.pdbpresent: return
pdb = read_PDB(self.pdbfile)
resnums = [0 for res in self.parm.residues]
chainids = ['Z' for res in self.parm.residues]
icodes = ['' for res in self.parm.residues]
tempfac = [0.0 for atom in self.parm.atoms]
occupancies = [0.0 for atom in self.parm.atoms]
atomnums = [-1 for atom in self.parm.atoms]
for res in self.parm.residues:
res.number = 0
res.chain = '*'
res.insertion_code = ''
for i, res in enumerate(pdb.residues):
try:
parmres = self.parm.residues[i]
reslab = parmres.name
resname = res.name.strip()
if resname != reslab:
if (not resname in ('WAT', 'HOH') or
not reslab in ('WAT', 'HOH')):
# Allow for 3's and 5's in terminal nucleic acid
# residues
if resname in ('A', 'C', 'G', 'U', 'DA',
'DG', 'DC', 'DT'):
if reslab[-1] in '35':
reslab = reslab[:-1]
if reslab != resname:
needs_warn = True
# Support other Amber residue name replacements
if reslab in ('ASP', 'ASH', 'AS4') and \
resname in ('ASP', 'ASH', 'AS4'):
needs_warn = False
elif reslab in ('GLU', 'GLH', 'GL4') and \
resname in ('GLU', 'GLH', 'GL4'):
needs_warn = False
elif reslab in ('HIP', 'HIS', 'HIE', 'HID') and \
resname in ('HIP', 'HIS', 'HIE', 'HID'):
needs_warn = False
elif reslab in ('LYS', 'LYN') and \
resname in ('LYS', 'LYN'):
needs_warn = False
elif reslab in ('CYS', 'CYX', 'CYM') and \
resname in ('CYS', 'CYX', 'CYM'):
needs_warn = False
if needs_warn:
warnings.warn('Residue name mismatch [#%d] %s '
'vs. %s' % (i+1, resname, reslab),
AddPDBWarning)
resnums[i] = parmres.number = res.number
chainids[i] = parmres.chain = res.chain
icodes[i] = parmres.insertion_code = res.insertion_code
except IndexError:
raise AddPDBError('PDB %s has more residues than prmtop %s' %
(self.pdbfile, self.parm))
# Now loop through all of the atoms in the parm residue, look for
# the atom with the same name in the PDB residue
for atom in parmres:
for pdbatom in res:
if atom.name == pdbatom.name:
tempfac[atom.idx] = pdbatom.bfactor
occupancies[atom.idx] = pdbatom.occupancy
atomnums[atom.idx] = pdbatom.number
break
ncmts = ['Residue number (resSeq) read from PDB file; DIMENSION(NRES)']
if self.printicodes or any(icodes):
haveicodes = True
ncmts += ['Residue insertion codes (iCode) present in %FLAG '
'RESIDUE_ICODE']
else:
haveicodes = False
ncmts += ['Residue insertion code (iCode) not present in PDB file',
'If present: %FLAG RESIDUE_ICODE, %FORMAT(20a4)']
self.parm.add_flag('RESIDUE_NUMBER', '20I4', data=resnums,
comments=ncmts)
self.parm.add_flag('RESIDUE_CHAINID', '20a4', data=chainids,
comments=['Residue chain ID (chainId) read from PDB '
'file; DIMENSION(NRES)']
)
if haveicodes:
self.parm.add_flag('RESIDUE_ICODE', '20a4', data=icodes,
comments=['Residue insertion code (iCode) read from PDB file; '
'DIMENSION(NRES)']
)
if self.elements:
self.parm.add_flag('ATOM_ELEMENT', '20a4',
data=['%2s' % (_Element[atm.atomic_number].upper())
for atm in self.parm.atoms
], comments=['Atom element name read from topology file']
)
self.parm.add_flag('ATOM_OCCUPANCY', '10F8.2', data=occupancies,
comments=['Atom occupancies read from the PDB file'])
self.parm.add_flag('ATOM_BFACTOR', '10F8.2', data=tempfac,
comments=['Atom temperature factor from the PDB file'])
self.parm.add_flag('ATOM_NUMBER', '10I8', data=atomnums,
comments=['Atom serial number from the PDB file'])
self.parm.load_atom_info() # Get that information saved
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class deletePDB(Action):
"""
This action deletes the flags added by addPDB if they are present.
"""
supported_subclasses = (AmberParm,)
[docs] def init(self, arg_list):
self.pdbpresent = ('RESIDUE_NUMBER' in self.parm.flag_list or
'RESIDUE_CHAINID' in self.parm.flag_list or
'RESIDUE_ICODE' in self.parm.flag_list or
'ATOM_ELEMENT' in self.parm.flag_list or
'ATOM_BFACTOR' in self.parm.flag_list or
'ATOM_OCCUPANCY' in self.parm.flag_list or
'ATOM_NUMBER' in self.parm.flag_list
)
def __str__(self):
if self.pdbpresent:
return 'Deleting PDB info from %s' % self.parm
return 'No PDB information present in %s. Doing nothing' % self.parm
[docs] def execute(self):
if not self.pdbpresent: return
self.parm.delete_flag('RESIDUE_NUMBER')
self.parm.delete_flag('RESIDUE_CHAINID')
self.parm.delete_flag('RESIDUE_ICODE')
self.parm.delete_flag('ATOM_ELEMENT')
self.parm.delete_flag('ATOM_BFACTOR')
self.parm.delete_flag('ATOM_OCCUPANCY')
self.parm.delete_flag('ATOM_NUMBER')
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class add12_6_4(Action):
"""
Adds the LENNARD_JONES_CCOEF term for the new divalent metal ion 12-6-4
Lennard-Jones potential term. If provided, the mask will allow you to
specify which ions are divalent. The C4 parameter between the metal ion and
water can either be taken from Ref. [1] for the requested [watermodel]
(TIP3P, TIP4PEW, or SPCE) or provided in the file specified by the c4file
keyword. The polarizabilities must be present in the the polfile file. The
chemical symbol of the element will be used to determine the atom type.
Parameters are expected in a file with 2 columns:
<atom type> <parameter>
All defaults come from Ref. [1], [2] and [3]
[1] Pengfei Li and Kenneth M. Merz, J. Chem. Theory Comput., 2014, 10,
289-297
[2] Pengfei Li, Lin F. Song and Kenneth M. Merz, J. Phys. Chem. B, 2015,
119, 883-895
[3] Pengfei Li, Lin F. Song and Kenneth M. Merz, J. Chem. Theory Comput.,
2015, 11, 1645-1657.
"""
usage = ('[<divalent ion mask>] [c4file <C4 Param. File> | watermodel '
'<water model>] [polfile <Pol. Param File> [tunfactor <tunfactor>]')
strictly_supported = (AmberParm, ChamberParm)
_supported_wms = ('TIP3P', 'TIP4PEW', 'SPCE')
[docs] def init(self, arg_list):
self.mask = AmberMask(self.parm,
arg_list.get_next_mask(optional=True, default=':ZN'))
self.c4file = arg_list.get_key_string('c4file', None)
self.watermodel = arg_list.get_key_string('watermodel', None)
self.polfile = arg_list.get_key_string('polfile',
os.path.join(os.getenv('AMBERHOME') or '', 'dat',
'leap', 'parm', 'lj_1264_pol.dat'))
self.tunfactor = arg_list.get_key_float('tunfactor', 1.0)
if self.c4file is None:
if self.watermodel is None:
self.watermodel = 'TIP3P'
else:
self.watermodel = self.watermodel.upper()
if not self.watermodel in self._supported_wms:
raise LJ12_6_4Error('Unrecognized water model %s; pick one '
'of the following: %s' %
(self.watermodel,
', '.join(self._supported_wms))
)
else:
if self.watermodel is not None:
raise LJ12_6_4Error('c4file and watermodel are mutually '
'exclusive')
def __str__(self):
retstr = ('Adding Lennard-Jones C-coefficient for 12-6-4 potential. '
'Polarizabilities read from [%s]. ' % self.polfile)
if self.c4file is None:
retstr += ('Using default C4 parameters for water model [%s].' %
self.watermodel)
else:
retstr += 'Reading C4 parameters from [%s].' % self.c4file
return retstr
[docs] def execute(self):
from parmed.tools.add1264 import params1264 as params
self.parm.delete_flag('LENNARD_JONES_CCOEF')
self.parm.add_flag('LENNARD_JONES_CCOEF', '5E16.8',
num_items=len(self.parm.parm_data['LENNARD_JONES_ACOEF']),
comments=['For 12-6-4 potential used for divalent metal ions'])
for i, param in enumerate(params(self.parm, self.mask, self.c4file,
self.watermodel, self.polfile,
self.tunfactor)):
self.parm.parm_data['LENNARD_JONES_CCOEF'][i] = param
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class HMassRepartition(Action):
"""
This action implements hydrogen mass repartitioning in the system by
changing the mass of each hydrogen to the desired value (the default new
hydrogen mass is 3.024 daltons) and adjusting the mass of the atom to which
it is bonded by the amount required to leave the total mass unchanged. By
default, water hydrogen masses are unchanged (the SETTLE algorithm for
implementing constraints on water molecules is analytical). Water masses can
be repartitioned as well with the 'dowater' keyword.
"""
usage = '[<mass>] [dowater]'
[docs] def init(self, arg_list):
self.changewater = arg_list.has_key('dowater')
self.new_h_mass = arg_list.get_next_float(optional=True, default=3.024)
def __str__(self):
retstr = ('Repartitioning hydrogen masses to %s daltons. ' %
self.new_h_mass)
if self.changewater:
return retstr + 'Also changing water hydrogen masses.'
return retstr + 'Not changing water hydrogen masses.'
[docs] def execute(self):
# Back up the masses in case something goes wrong
original_masses = [atom.mass for atom in self.parm.atoms]
water = SOLVENT_NAMES
for i, atom in enumerate(self.parm.atoms):
if atom.atomic_number != 1: continue
if not self.changewater and atom.residue.name in water:
continue
heteroatom = None
heteroidx = 0
bondeds = list(atom.bond_partners)
while heteroidx < len(bondeds):
if bondeds[heteroidx].atomic_number != 1:
heteroatom = bondeds[heteroidx]
break
heteroidx += 1
if heteroatom is None:
# Only bonded to other hydrogens. Weird, but do not repartition
warnings.warn('H atom detected not bound to heteroatom. '
'Ignoring.', ParmWarning)
continue
transfermass = self.new_h_mass - atom.mass
atom.mass = self.new_h_mass
heteroatom.mass -= transfermass
if isinstance(self.parm, AmberParm):
self.parm.parm_data['MASS'][i] = self.new_h_mass
self.parm.parm_data['MASS'][heteroatom.idx] -= transfermass
# Now make sure that all masses are positive, or revert masses and
# raise an exception
for atom in self.parm.atoms:
if atom.mass <= 0 and atom.atomic_number > 0:
for i, atom in enumerate(self.parm.atoms):
atom.mass = original_masses[i]
if isinstance(self.parm, AmberParm):
self.parm.parm_data['MASS'][i] = original_masses[i]
raise HMassRepartitionError('Too much mass removed from atom '
'%d. Hydrogen masses must be '
'smaller.' % i)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class OpenMM(Action):
"""
This class will read a sander/pmemd input file and run an equivalent
simulation using the OpenMM Python API. It uses a topology file that has
been defined (and/or modified) here. The command-line flags are identical to
those used for sander and pmemd (described in the Amber manual). You can
additionally specify a computational platform available through OpenMM such
as CUDA, OpenCL, Reference, and CPU.
The default prmtop will be the 'active' parm. The prmtop can be changed
using either the 'parm' keyword or the '-p' flag, but it must resolve to one
of the stored topology files. Any parm given with the '-p' flag has
precedence. If present, the 'dcd' keyword triggers writting DCD-formatted
trajectory files instead of NetCDF when ioutfm=1. The DCD trajectory writing
offers binary trajectory support without requiring a NetCDF-Python library.
The 'progress' keyword triggers printing of ParmEd's progress in setting up
and starting the calculation.
The 'script' keyword provides a file to write an OpenMM script that performs
the same calculation requested by this ParmEd command. The 'norun' command
will prevent anything from actually being run and can only be used when a
script file is provided.
"""
usage = ('[-p <parm>|<parm index>] [sander/pmemd options] [-platform '
'<platform>] [-precision <precision model>] [dcd] [progress] '
'[script <script_file.py>] [norun]')
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
parm = arg_list.get_key_string('-p', default=None)
if parm is not None:
self.parm = self.parm_list[parm]
# This consumes the remaining arguments
self.arg_list = ArgumentList(arg_list)
def __str__(self):
return ("Running Amber-style simulation with OpenMM using the command "
"string:\n\t%s\nThis may take awhile..." % self.arg_list)
[docs] def execute(self):
""" Runs the OpenMM simulation """
from parmed.tools.simulations.openmm import simulate, HAS_OPENMM
if not HAS_OPENMM:
raise SimulationError('OpenMM could not be imported. Skipping.') # pragma: no cover
# First try to load a restart file if it was supplied
inptraj = self.arg_list.has_key('-y', mark=False)
has_inpcrd = self.arg_list.has_key('-c', mark=False)
if self.parm.coordinates is None and not inptraj and not has_inpcrd:
raise SimulationError('No input coordinates provided.')
# Eliminate some incompatibilities that are easy to catch now
simulate(self.parm, self.arg_list)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class energy(Action):
"""
This action will compute a single-point energy for the loaded structure (you
must use 'loadRestart' prior to this command to load coordinates). The
following options and keywords are supported:
Options
-------
- cutoff <cut> : The size of the non-bonded cutoff, in Angstroms.
Default 8 A for periodic systems or infinite for
nonperiodic systems
For systems with no periodic box information
--------------------------------------------
- igb <IGB> : An integer corresponding to the desired GB model. May be
1, 2, 5, 7, or 8 as described by the sander and pmemd
manual. Default 5.
- saltcon <conc> : The salt concentration, in mol/L, modeled by a Debye
screening parameter. Default 0.0
For periodic systems
--------------------
- Ewald : Use an Ewald sum to compute long-range electrostatics instead
of PME
- nodisper : Do not use a long-range vdW dispersion correction
OpenMM-specific options
-----------------------
- omm : If present, this keyword will instruct ParmEd to use OpenMM to
compute the energies (and optionally forces) using OpenMM
instead of sander.
- platform <platform> : OpenMM compute platform to use. Options are
CUDA, OpenCL, Reference, and CPU. Consult the OpenMM manual for
details (only used if 'omm' is provided)
- precision <precision model> : Precision model to use. Options are
single, double, and mixed. Reference platform is always double
and CPU platform is always single. Mixed (default) uses single
precision for calculations and double for accumulation (only
used if 'omm' is provided)
- decompose : Print bond, angle, dihedral, and nonbonded energies
separately
- applayer : Use OpenMM's class to compute the energy
Examples
--------
# Using AMBER
import parmed as pmd
parm = pmd.load_file('prmtop', xyz='rst7')
pmd.tools.energy(parm, 'igb 8').execute()
# Using Openmm
pmd.tools.energy(parm, 'igb 5 omm').execute()
"""
usage = ('[cutoff <cut>] [[igb <IGB>] [saltcon <conc>] | [Ewald]] '
'[nodisper] [omm] [applayer] [platform <platform>] [precision '
'<precision model>] [decompose]')
output = sys.stdout
[docs] def init(self, arg_list):
self.use_openmm = (arg_list.has_key('omm') or
not isinstance(self.parm, AmberParm))
self.arg_list = ArgumentList(arg_list)
if self.use_openmm and isinstance(self.parm, AmoebaParm):
raise NotImplementedError('Amoeba prmtops can only get energies '
'from sander')
def __str__(self):
return 'Computing a single-point energy for %s' % self.parm
[docs] def execute(self):
if self.use_openmm:
from parmed.tools.simulations.openmm import energy, HAS_OPENMM
if not HAS_OPENMM:
raise SimulationError('OpenMM could not be imported. Skipping.') # pragma: no cover
energy(self.parm, self.arg_list, self.output)
else:
from parmed.tools.simulations.sanderapi import energy, HAS_SANDER
if not HAS_SANDER:
raise SimulationError('sander could not be imported. Skipping.') # pragma: no cover
# Consume the OMM-specific arguments so we don't have any apparently
# unused arguments
self.arg_list.get_key_string('platform', None)
self.arg_list.get_key_string('precision', None)
self.arg_list.has_key('decompose')
self.arg_list.has_key('applayer')
energy(self.parm, self.arg_list, self.output)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class deleteBond(Action):
"""
This action deletes any bonds that occur between the atoms in two masks.
- <mask1> : Amber mask defining one of the atoms in a bond
- <mask2> : Amber mask defining the other atom in the bond
- [verbose] : Print out every bond that is deleted as well as the
number of other valence terms that were eliminated.
All bonds will be matched in which one atom comes from <mask1> and the other
atom comes from <mask2>. This action will also delete any other valence term
(like angles and dihedrals) that would get severed by the deletion of one of
the bonds.
"""
usage = '<mask1> <mask2>'
[docs] def init(self, arg_list):
self.mask1 = AmberMask(self.parm, arg_list.get_next_mask())
self.mask2 = AmberMask(self.parm, arg_list.get_next_mask())
self.verbose = arg_list.has_key('verbose')
# Go through each atom in mask1 and see if it forms a bond with any atom
# in mask2.
bonds_to_delete = set()
self.del_bonds = set()
self.del_angles = set()
self.del_dihedrals = set()
self.del_rbtorsions = set()
self.del_urey_bradleys = set()
self.del_impropers = set()
self.del_cmaps = set()
self.del_trigonal_angles = set()
self.del_oopbends = set()
self.del_pi_torsions = set()
self.del_strbnds = set()
self.del_tortors = set()
for i in self.mask1.Selected():
ai = self.parm.atoms[i]
for j in self.mask2.Selected():
aj = self.parm.atoms[j]
# Skip if these two atoms are identical
if ai is aj: continue
for bond in ai.bonds:
if aj not in bond: continue
bonds_to_delete.add(bond)
# Find other valence terms we need to delete
for i, bond in enumerate(self.parm.bonds):
if bond in bonds_to_delete:
self.del_bonds.add(i)
for bond in bonds_to_delete:
for i, angle in enumerate(self.parm.angles):
if bond in angle:
self.del_angles.add(i)
for i, dihed in enumerate(self.parm.dihedrals):
if bond in dihed:
self.del_dihedrals.add(i)
for i, dihed in enumerate(self.parm.rb_torsions):
if bond in dihed:
self.del_rbtorsions.add(i)
for i, urey in enumerate(self.parm.urey_bradleys):
if bond in urey:
self.del_urey_bradleys.add(i)
for i, improper in enumerate(self.parm.impropers):
if bond in improper:
self.del_impropers.add(i)
for i, cmap in enumerate(self.parm.cmaps):
if bond in cmap:
self.del_cmaps.add(i)
for i, trigonal_angle in enumerate(self.parm.trigonal_angles):
if bond in trigonal_angle:
self.del_trigonal_angles.add(i)
for i, oopbend in enumerate(self.parm.out_of_plane_bends):
if bond in oopbend:
self.del_oopbends.add(i)
for i, pitor in enumerate(self.parm.pi_torsions):
if bond in pitor:
self.del_pi_torsions.add(i)
for i, strbnd in enumerate(self.parm.stretch_bends):
if bond in strbnd:
self.del_strbnds.add(i)
for i, tortor in enumerate(self.parm.torsion_torsions):
if bond in tortor:
self.del_tortors.add(i)
def __str__(self):
if not self.del_bonds:
return 'No bonds to delete'
if not self.verbose:
return 'Deleting the %d bonds found between %s and %s' % (
len(self.del_bonds), self.mask1, self.mask2)
# Now we want full statistics
retstr = 'Deleting %d bonds between %s and %s:\n' % (
len(self.del_bonds), self.mask1, self.mask2)
for bond in self.del_bonds:
a1, a2 = self.parm.bonds[bond].atom1, self.parm.bonds[bond].atom2
retstr += '\t%d [%s %d] %s --- %d [%s %d] %s\n' % (a1.idx+1,
a1.residue.name, a1.residue.idx+1, a1.name, a2.idx+1,
a2.residue.name, a2.residue.idx+1, a2.name)
retstr += 'Deleting %d angles, ' % (len(self.del_angles))
if self.parm.trigonal_angles or self.parm.out_of_plane_bends or \
self.parm.stretch_bends or self.parm.torsion_torsions:
retstr += ('%d Urey-Bradleys, %d trigonal angles,\n '
'%d dihedrals, %d out-of-plane bends, %d stretch-bends\n'
' and %d torsion-torsions' % (
len(self.del_urey_bradleys),
len(self.del_trigonal_angles),
len(self.del_dihedrals)+len(self.del_rbtorsions),
len(self.del_oopbends), len(self.del_strbnds),
len(self.del_tortors))
)
elif self.parm.urey_bradleys or self.parm.impropers or self.parm.cmaps:
retstr += ('%d Urey-Bradleys, %d impropers,\n %d dihedrals '
'and %d CMAPs' % (
len(self.del_urey_bradleys), len(self.del_impropers),
len(self.del_dihedrals)+len(self.del_rbtorsions),
len(self.del_cmaps))
)
else:
retstr += 'and %d dihedrals' % (len(self.del_dihedrals) +
len(self.del_rbtorsions))
return retstr
[docs] def execute(self):
if not self.del_bonds: return
for i in sorted(self.del_bonds, reverse=True):
self.parm.bonds[i].delete()
del self.parm.bonds[i]
for i in sorted(self.del_angles, reverse=True):
self.parm.angles[i].delete()
del self.parm.angles[i]
for i in sorted(self.del_dihedrals, reverse=True):
self.parm.dihedrals[i].delete()
del self.parm.dihedrals[i]
for i in sorted(self.del_rbtorsions, reverse=True):
self.parm.rb_torsions[i].delete()
del self.parm.rb_torsions[i]
for i in sorted(self.del_urey_bradleys, reverse=True):
self.parm.urey_bradleys[i].delete()
del self.parm.urey_bradleys[i]
for i in sorted(self.del_impropers, reverse=True):
self.parm.impropers[i].delete()
del self.parm.impropers[i]
for i in sorted(self.del_cmaps, reverse=True):
self.parm.cmaps[i].delete()
del self.parm.cmaps[i]
for i in sorted(self.del_trigonal_angles, reverse=True):
del self.parm.trigonal_angles[i]
for i in sorted(self.del_oopbends, reverse=True):
del self.parm.out_of_plane_bends[i]
for i in sorted(self.del_tortors, reverse=True):
self.parm.torsion_torsions[i].delete()
del self.parm.torsion_torsions[i]
for i in sorted(self.del_strbnds, reverse=True):
del self.parm.stretch_bends[i]
try:
self.parm.remake_parm()
except AttributeError:
self.parm.prune_empty_terms()
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class chamber(Action):
"""
This action will read CHARMM parameter, topology (RTF), and stream (STR)
files and apply those parameters to a structure defined in a CHARMM PSF
(protein structure file). XPLOR, CHARMM, and VMD-generated PSF files are all
supported. You may specify -top, -param, and -str as many times as you want
to provide multiple parameter files. All topology files are read first (in
the order they are specified), followed by all parameter files, and finally
all stream files are read in. Topology files are only necessary if the
parameter files do not define the atom types (newer CHARMM force fields
define atom types directly in the parameter file. Environment variables and
shell wild-cards (* and ?), as well as the home shortcut character ~ are
properly expanded when looking for files.
Options
-------
- -top CHARMM topology, or Residue Topology File (RTF) file
- -param CHARMM parameter file
- -str CHARMM stream file. Only RTF and parameter sections read
- -toppar CHARMM topology, parameter, and/or stream file(s). File
type is detected from extension (or in the case of .inp
files, the presence of 'top', 'par' in the name)
- -psf CHARMM PSF file
- -crd Coordinate file (PDB, CHARMM CRD, Amber restart, etc.)
- -nocmap Do not use any CMAP parameters
- -box Box dimensions. If no angles are defined, they are assumed
to be 90 degrees (orthorhombic box). Alternatively, you
can use the word 'bounding' to define a box that encloses
the centers of all atoms.
- -radii Implicit solvent solvation radii. <radiusset> can be
amber6, bondi, mbondi, mbondi2, mbondi3
Same effect as the changeRadii command. Default is mbondi.
- nosettle This avoids changing the names of the water residue and
oxygen atoms from TIP3 and OH2 to WAT and O
If the PDB file has a CRYST1 record, the box information will be set from
there. Any box info given on the command-line will override the box found in
the crd file.
"""
usage = ('-top <RTF> -param <PAR> [-str <STR>] -psf <PSF> [-crd <CRD>] '
'[-nocmap] [-box a,b,c[,alpha,beta,gamma]|bounding] [-radii '
'<radiusset>] [nosettle]')
needs_parm = False
[docs] def init(self, arg_list):
from os.path import expandvars, expanduser
import glob
# First get all of the topology files
self.topfiles, self.paramfiles, self.streamfiles = [], [], []
while arg_list.has_key('-top', mark=False):
topfile = expanduser(arg_list.get_key_string('-top', None))
topfile = expandvars(topfile)
topfiles = glob.glob(topfile)
if not topfiles:
raise FileDoesNotExist('CHARMM RTF file %s does not exist' %
topfile)
self.topfiles.extend(topfiles)
while arg_list.has_key('-param', mark=False):
param = expanduser(arg_list.get_key_string('-param', None))
param = expandvars(param)
params = glob.glob(param)
if not params:
raise FileDoesNotExist('CHARMM parameter file %s does not exist'
% param)
self.paramfiles.extend(params)
while arg_list.has_key('-str', mark=False):
stream = expanduser(arg_list.get_key_string('-str', None))
stream = expandvars(stream)
streams = glob.glob(stream)
if not streams:
raise FileDoesNotExist('CHARMM stream file %s does not exist' %
stream)
self.streamfiles.extend(streams)
while arg_list.has_key('-toppar', mark=False):
toppar = expanduser(arg_list.get_key_string('-toppar', None))
toppar = expandvars(toppar)
toppars = glob.glob(toppar)
if not toppars:
raise FileDoesNotExist('CHARMM toppar file %s does not exist' %
toppar)
for fname in toppars:
base = os.path.split(fname)[1] # ignore the directory name
if base.endswith('.rtf') or base.endswith('.top'):
self.topfiles.append(fname)
elif base.endswith('.par') or base.endswith('.prm'):
self.paramfiles.append(fname)
elif base.endswith('.str'):
self.streamfiles.append(fname)
elif base.endswith('.inp'):
if 'par' in base:
self.paramfiles.append(fname)
elif 'top' in base:
self.topfiles.append(fname)
else:
raise InputError('Cannot detect filetype of %s' % fname)
else:
raise InputError('Cannot detect filetype of %s' % fname)
crdfile = arg_list.get_key_string('-crd', None)
if crdfile is not None:
crdfiles = glob.glob(expanduser(expandvars(crdfile)))
if not crdfiles:
raise FileDoesNotExist('Coordinate file %s does not exist' %
crdfile)
if len(crdfiles) > 1:
raise InputError('Too many coordinate files selected through '
'globbing')
self.crdfile = crdfiles[0]
else:
self.crdfile = None
self.cmap = not arg_list.has_key('-nocmap')
psf = arg_list.get_key_string('-psf', None)
if psf is None:
raise InputError('chamber requires a PSF file')
self.psf = glob.glob(expanduser(expandvars(psf)))
if not self.psf:
raise FileDoesNotExist('chamber PSF file %s cannot be found' % psf)
if len(self.psf) > 1:
raise InputError('Too many PSF files selected through globbing')
self.psf = self.psf[0]
box = arg_list.get_key_string('-box', None)
if box is None:
self.box = None
elif box.lower() != 'bounding':
try:
self.box = [float(w) for w in box.replace(',', ' ').split()]
except ValueError:
raise InputError('Box info must be comma-delimited floats')
if len(self.box) == 3:
self.box += [90.0, 90.0, 90.0]
elif len(self.box) != 6:
raise InputError('Box must be 3 lengths or 3 lengths and 3 '
'angles!')
else:
self.box = 'bounding'
self.radii = arg_list.get_key_string('-radii', 'mbondi')
# Make sure we have legal input
if not self.paramfiles and not self.streamfiles:
raise InputError('No parameter files were provided!')
if not self.radii in GB_RADII:
raise InputError('Illegal radius set %s -- must be one of '
'%s' % (self.radii, ', '.join(GB_RADII)))
arg_list.has_key('nocondense')
self.nosettle = arg_list.has_key('nosettle')
def __str__(self):
retstr = 'Creating chamber topology file from PSF %s, ' % self.psf
if self.topfiles:
retstr += 'RTF files [%s], ' % (', '.join(self.topfiles))
if not self.streamfiles or not self.paramfiles:
retstr += 'and '
if self.paramfiles:
retstr += 'PAR files [%s]' % (', '.join(self.paramfiles))
if self.streamfiles:
retstr += ', and '
if self.streamfiles:
retstr += 'STR files [%s].' % (', '.join(self.streamfiles))
if self.crdfile is not None:
retstr += ' Coords from %s.' % self.crdfile
if self.cmap:
retstr += ' Using CMAP.'
else:
retstr += ' NO CMAP.'
if self.box == 'bounding':
retstr += ' Defining a bounding box.'
elif self.box is not None:
retstr += ' Box info %s.' % (self.box)
retstr += ' GB Radius set %s.' % self.radii
if self.nosettle:
retstr += ' Not changing water names.'
return retstr
[docs] def execute(self):
parmset = CharmmParameterSet()
for tfile in self.topfiles:
parmset.read_topology_file(tfile)
for pfile in self.paramfiles:
parmset.read_parameter_file(pfile)
for sfile in self.streamfiles:
parmset.read_stream_file(sfile)
# Now read the PSF
psf = CharmmPsfFile(self.psf)
# Read the PDB and set the box information
if self.crdfile is not None:
crdbox = None
# Automatic format determination
crd = load_file(self.crdfile)
try:
if len(crd.coordinates.shape) == 3:
coords = crd.coordinates[0]
else:
coords = crd.coordinates
except (AttributeError, TypeError):
raise ChamberError('No coordinates in %s' % self.crdfile)
if hasattr(crd, 'box') and crd.box is not None:
if len(crd.box.shape) == 1:
crdbox = crd.box
else:
# Trajectory
crdbox = crd.box[0]
if coords.shape != (len(psf.atoms), 3):
raise ChamberError('Mismatch in number of coordinates (%d) and '
'3*number of atoms (%d)' % (len(coords),
3*len(psf.atoms)))
# Set the coordinates now, since creating the parm may re-order the
# atoms in order to maintain contiguous molecules
psf.coordinates = coords
# Set the box info from self.box if set
if self.box is None and crdbox is not None:
if len(crdbox) == 3:
psf.box = list(crdbox) + [90.0, 90.0, 90.0]
elif len(crdbox) == 6:
psf.box = list(crdbox)
else:
raise ValueError('Unexpected box array shape')
elif self.box == 'bounding':
# Define the bounding box
extent = coords.max(axis=0) - coords.min(axis=0)
psf.box = list(extent) + [90.0, 90.0, 90.0]
elif self.box is not None:
if len(self.box) == 3:
psf.box = list(self.box) + [90.0, 90.0, 90.0]
elif len(self.box) == 6:
psf.box = copy.copy(self.box)
else:
raise ValueError('Unexpected box array shape')
else:
psf.box = None
else:
# Set the box information
if self.box is None:
psf.box = None
else:
if len(self.box) == 3:
psf.box = list(self.box) + [90.0, 90.0, 90.0]
elif len(self.box) == 6:
psf.box = self.box[:]
else:
raise ValueError('Unexpected box array shape')
nsets = len(parmset.parametersets)
if nsets > 0:
frcfield = []
for pset in parmset.parametersets:
frcfield.extend([nsets, pset])
else:
frcfield = [1, 'No FF information parsed...']
# Delete the CMAP list if requested
if not self.cmap:
psf.clear_cmap()
# Now load the parameters
try:
psf.load_parameters(parmset)
except ParmedError as e:
raise ChamberError('Problem assigning parameters to PSF: %s' % e)
parm = ConvertFromPSF(psf, parmset)
parm.name = self.psf
changeRadii(parm, self.radii).execute()
self.parm_list.add_parm(parm)
self.parm = parm
if not self.nosettle:
# Change water residue names from TIP3 to WAT and OH2 to O
for res, nam in zip(parm.residues, parm.parm_data['RESIDUE_LABEL']):
if nam != 'TIP3':
continue
res.name = parm.parm_data['RESIDUE_LABEL'][res.idx] = 'WAT'
for atom in res:
if atom.name == 'OH2':
atom.name = parm.parm_data['ATOM_NAME'][atom.idx] = 'O'
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class minimize(Action):
"""
This action takes a structure and minimizes the energy using either
scipy.optimize.minimize (with BFGS) and the sander API *or* OpenMM.
Following this action, the coordinates stored in the topology will be the
minimized structure
General Options
---------------
- omm Use OpenMM instead of sander+scipy to minimize
- cutoff <cutoff> Nonbonded cutoff in Angstroms
- tol <tolerance> The largest energy difference between successive
steps in the minimization allow the minimization
to stop (Default 0.001)
- maxcyc <cycles> The maximum number of minimization cycles
permitted. No limit by default (minimization
stops when tolerance is satisfied)
Implicit Solvent options
------------------------
- igb <IGB> GB model to use (0=No GB, 1,2,5,7,8 correspond to
Amber models)
- saltcon <conc> Salt concentration for GB in Molarity
OpenMM-specific options
-----------------------
- restrain <mask> Selection of atoms to restrain
- weight <k> Weight of positional restraints (kcal/mol/A^2)
- norun Do not run the calculation
- script <script> Name of the Python script to write to run this
calculation
- platform <platform> Which compute platform to use. Options are CUDA,
OpenCL, CPU, and Reference. Consult OpenMM docs
for more information
- precision <precision_model> Which precision model to use. Can be
"single", "double", or "mixed", and only has an
effect on the CUDA and OpenCL platforms.
The implicit solvent options will be ignored for systems with periodic
boundary conditions. The precision option will be ignored for platforms that
do not support user-specified precision. The platform, precision, and script
arguments will be ignored unless ``omm`` is specified.
"""
usage = ('[cutoff <cut>] [[igb <IGB>] [saltcon <conc>]] [[restrain <mask>] '
'[weight <k>]] [norun] [script <script_file.py>] [platform '
'<platform>] [precision <precision model>] [tol <tolerance>] '
'[maxcyc <cycles>] [omm]')
not_supported = (AmoebaParm,)
[docs] def init(self, arg_list):
self.use_openmm = arg_list.has_key('omm')
self.cutoff = arg_list.get_key_float('cutoff', None)
mask = arg_list.get_key_mask('restrain', None)
self.igb = arg_list.get_key_int('igb', 0)
self.saltcon = arg_list.get_key_float('saltcon', 0.0)
self.weight = arg_list.get_key_float('weight', 0.0)
self.norun = arg_list.has_key('norun')
self.script = arg_list.get_key_string('script', None)
self.platform = arg_list.get_key_string('platform', None)
self.precision = arg_list.get_key_string('precision', 'mixed')
self.tol = arg_list.get_key_float('tol', 0.001)
self.maxcyc = arg_list.get_key_int('maxcyc', None)
# Check for legal values
if self.parm.ptr('ifbox') == 0:
if self.cutoff is None or self.cutoff > 500:
self.cutoff = None # no cutoff
elif self.cutoff is None:
self.cutoff = 8.0
elif self.cutoff < 7:
raise InputError('Cutoff unreasonably small.')
if self.parm.ptr('ifbox') == 0 and self.saltcon < 0:
raise InputError('Salt concentration must be non-negative')
if mask is not None:
self.restrain = AmberMask(self.parm, mask)
if not self.use_openmm:
raise InputError('Restraints only permitted with omm option')
if self.weight <= 0:
raise InputError('Restraints require a restraint stiffness '
'larger than 0 kcal/mol/A^2')
else:
self.restrain = None
if self.platform not in ('CUDA', 'OpenCL', 'CPU', 'Reference', None):
raise InputError('platform must be CUDA, OpenCL, CPU or Reference '
'(NOT %s)' % self.platform)
if self.precision not in ('mixed', 'single', 'double'):
raise InputError('precision must be single, double, or mixed.')
if self.parm.ptr('ifbox') == 0 and not self.igb in (0,1,2,5,6,7,8):
raise InputError('Illegal igb value (%d). Must be 0, 1, 2, 5, 6, '
'7, or 8' % self.igb)
if self.maxcyc is not None and self.maxcyc <= 0:
raise InputError('maxcyc must be a positive integer')
if self.tol < 0:
raise InputError('tolerance must be positive')
if self.tol == 0 and self.maxcyc is None:
raise InputError('tolerance must be > 0 if maxcyc is not set.')
if self.norun and not self.use_openmm:
raise InputError('norun only makes sense with the script and omm '
'options')
def __str__(self):
retstr = 'Minimizing %s ' % self.parm
if self.parm.ptr('ifbox'):
retstr += 'with PME '
elif self.igb == 0 or self.igb == 6:
retstr += 'in gas phase '
elif self.igb == 1:
retstr += 'with GB(HCT) '
elif self.igb == 2:
retstr += 'with GB(OBC1) '
elif self.igb == 5:
retstr += 'with GB(OBC2) '
elif self.igb == 7:
retstr += 'with GB(GBneck) '
elif self.igb == 8:
retstr += 'with GB(GBneck2) '
if self.cutoff is None:
retstr += 'and no cutoff '
else:
retstr += 'and a cutoff of %.2f Angstroms ' % self.cutoff
retstr += 'to a tolerance of %s. ' % self.tol
if self.restrain is not None:
retstr += 'Restraining %s with weights %f. ' % (self.restrain,
self.weight)
if self.maxcyc is not None:
retstr += 'Using at most %d minimization steps.' % self.maxcyc
return retstr.strip()
[docs] def execute(self):
from parmed.tools.simulations.openmm import minimize as omm_min
from parmed.tools.simulations.sanderapi import minimize as san_min
if self.parm.coordinates is None:
raise SimulationError('You must load coordinates before "minimize"')
if self.use_openmm:
omm_min(self.parm, self.igb, self.saltcon, self.cutoff,
self.restrain, self.weight, self.script, self.platform,
self.precision, self.norun, self.tol, self.maxcyc)
else:
san_min(self.parm, self.igb, self.saltcon, self.cutoff, self.tol,
self.maxcyc)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class outPDB(Action):
"""
Write a PDB file from the currently active system to <file>
- <file>: The PDB file to write
- [norenumber]: Use the original atom and residue numbering if available
- [charmm]: Put the SEGID, if available, in columns 72 to 76
- [anisou]: Write anisotropic B-factors if available
"""
usage = "<file> [norenumber] [charmm] [anisou]"
[docs] def init(self, arg_list):
self.renumber = not arg_list.has_key('norenumber')
self.charmm = arg_list.has_key('charmm')
self.anisou = arg_list.has_key('anisou')
self.filename = arg_list.get_next_string()
if self.parm.coordinates is None:
raise InputError('Parm %s does not have loaded coordinates' %
self.parm)
def __str__(self):
retstr = 'Writing PDB file %s' % self.filename
if self.renumber:
retstr += ' renumbering atoms and residues'
else:
retstr += ' not renumbering atoms and residues'
if self.charmm:
retstr += ' and adding CHARMM SEGIDs'
if self.anisou:
retstr += ' and adding anisotropic B-factors'
return retstr
[docs] def execute(self):
self.parm.write_pdb(self.filename, renumber=self.renumber,
charmm=self.charmm, write_anisou=self.anisou)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class outCIF(Action):
"""
Write a PDBx/mmCIF file from the currently active system to <file>
- <file>: The PDBx/mmCIF file to write
- [norenumber]: Use the original atom and residue numbering if available
- [anisou]: Write anisotropic B-factors if available
"""
usage = "<file> [norenumber] [charmm] [anisou]"
[docs] def init(self, arg_list):
self.renumber = not arg_list.has_key('norenumber')
self.anisou = arg_list.has_key('anisou')
self.filename = arg_list.get_next_string()
if self.parm.coordinates is None:
raise InputError('Parm %s does not have loaded coordinates' %
self.parm)
def __str__(self):
retstr = 'Writing PDB file %s' % self.filename
if self.renumber:
retstr += ' renumbering atoms and residues'
else:
retstr += ' not renumbering atoms and residues'
if self.anisou:
retstr += ' and adding anisotropic B-factors'
return retstr
[docs] def execute(self):
self.parm.write_cif(self.filename, renumber=self.renumber,
write_anisou=self.anisou)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class gromber(Action):
"""
Load a Gromacs topology file with parameters as an Amber-formatted system.
Note, if your Gromacs topology file requires any include topology files, you
will need to have Gromacs installed for this to work.
- <top_file>: The Gromacs topology file to load
- <coord_file>: The coordinate file to load into the system. Can be any
recognized format (GRO, PDB, mmCIF, inpcrd, etc.)
- define <DEFINE[=VAR]>: Preprocessor defines that control the
processing of the Gromacs topology file.
- topdir <directory>: The directory containing all Gromacs include
topology files. This is only necessary if Gromacs is not
installed in a location that ParmEd can find.
- radii <radiusset>: The GB radius set to use. Can be 'mbondi', 'bondi',
'mbondi2', or 'amber6'. Default is mbondi
Gromacs topology files do not store the unit cell information. Therefore, in
order to make sure that unit cell information is properly assigned to the
resulting system, the provided ``<coord_file>`` should contain unit cell
information (e.g., GRO, PDB, PDBx/mmCIF, and inpcrd files can all store box
information).
ParmEd will try to locate the Gromacs topology directory using either the
GMXDATA or GMXBIN environment variables (which should point to the
$PREFIX/share/gromacs or $PREFIX/bin directories, respectively, where
$PREFIX is the install prefix). If neither is set, the topology directory is
located relative to the location of the ``gmx`` (Gromacs 5+) or ``pdb2gmx``
(Gromacs 4 or older) in the user's PATH. If none of the above is true, the
default installation location (/usr/local/gromacs/share/gromacs/top) is
used. Any provided ``topdir`` will override default choices (but only for
this particular command -- future ``gromber`` actions will use the default
location again).
You can provide as many defines as you wish, and the ordering you specify
them is preserved. The default value assigned to each define is "1". To
provide multiple defines, use the ``define`` keyword multiple times, for
example:
define MYVAR=something define MYVAR2=something_else ...
"""
usage = ("<top_file> [<coord_file>] [define <DEFINE[=VAR]>] "
"[topdir <directory>] [radii <radiusset>]")
needs_parm = False
[docs] def init(self, arg_list):
self.topfile = arg_list.get_next_string()
defines = OrderedDict()
current_define = arg_list.get_key_string('define', None)
while current_define is not None:
if '=' in current_define:
parts = current_define.split('=')
if len(parts) != 2:
raise InputError('Illegal define: %s' % current_define)
key, val = parts
defines[key] = val
else:
defines[current_define] = '1'
current_define = arg_list.get_key_string('define', None)
self.defines = defines or None
self.orig_topdir = gromacs.GROMACS_TOPDIR
topdir = arg_list.get_key_string('topdir', gromacs.GROMACS_TOPDIR)
gromacs.GROMACS_TOPDIR = topdir
self.coordinate_file = arg_list.get_next_string(optional=True)
self.radii = arg_list.get_key_string('radii', 'mbondi')
def __str__(self):
retstr = ['Converting Gromacs topology %s to Amber. ' % self.topfile]
if self.defines:
retstr.append('Using the following defines:\n')
for key, val in iteritems(self.defines):
retstr.append('\t%s=%s\n' % (key, val))
retstr.append('Using topology directory [%s]. ' %
gromacs.GROMACS_TOPDIR)
if self.coordinate_file is None:
retstr.append('No coordinates provided.')
else:
retstr.append('Getting coordinates (and box) from %s.' %
self.coordinate_file)
return ''.join(retstr)
[docs] def execute(self):
try:
top = gromacs.GromacsTopologyFile(self.topfile,defines=self.defines)
finally:
gromacs.GROMACS_TOPDIR = self.orig_topdir
if self.coordinate_file is not None:
crd = load_file(self.coordinate_file)
if isinstance(crd, Structure):
if len(top.atoms) != len(crd.atoms):
raise InputError('Coordinate/topology file size mismatch')
for a1, a2 in zip(top.atoms, crd.atoms):
try:
a1.xx = a2.xx
a1.xy = a2.xy
a1.xz = a2.xz
except AttributeError:
raise InputError('%s does not contain coordinates' %
self.coordinate_file)
try:
a1.vx = a2.vx
a1.vy = a2.vy
a1.vz = a2.vz
except AttributeError:
pass
elif hasattr(crd, 'coordinates') and not callable(crd.coordinates):
it = iter(crd.coordinates)
for a, x, y, z in zip(top.atoms, it, it, it):
a.xx, a.xy, a.xz = x, y, z
elif hasattr(crd, 'coordinates') and callable(crd.coordinates):
it = iter(crd.coordinates(0)) # first frame
for a, x, y, z in zip(top.atoms, it, it, it):
a.xx, a.xy, a.xz = x, y, z
elif hasattr(crd, 'positions'):
pos = crd.positions.value_in_unit(u.angstroms)
for a, xyz in zip(top.atoms, pos):
a.xx, a.xy, a.xz = pos
else:
raise InputError('Cannot find coordinates in %s' %
self.coordinate_file)
if hasattr(crd, 'box') and callable(crd.box):
top.box = copy.copy(crd.box(0))
else:
top.box = copy.copy(crd.box)
try:
if top.impropers or top.urey_bradleys or top.cmaps:
self.parm = ChamberParm.from_structure(top)
else:
self.parm = AmberParm.from_structure(top)
except TypeError as err:
raise InputError(str(err))
self.parm.name = self.topfile
changeRadii(self.parm, self.radii).execute()
self.parm_list.add_parm(self.parm)
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
# Private helper methods
def _change_lj_pair(parm, atom_1, atom_2, rmin, eps, one_4=False):
if one_4:
key = 'LENNARD_JONES_14'
else:
key = 'LENNARD_JONES'
# Make sure that atom type 1 comes first
a1, a2 = sorted([atom_1, atom_2])
ntypes = parm.pointers['NTYPES']
# Find the atom1 - atom2 interaction (adjusting for indexing from 0)
term_idx = parm.parm_data['NONBONDED_PARM_INDEX'][ntypes*(a1-1)+a2-1] - 1
# Now change the ACOEF and BCOEF arrays, assuming pre-combined values
parm.parm_data['%s_ACOEF' % key][term_idx] = eps * rmin**12
parm.parm_data['%s_BCOEF' % key][term_idx] = 2 * eps * rmin**6
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+