Commit a1953312 authored by Robert McGibbon's avatar Robert McGibbon
Browse files

Merge branch 'docstrings' of github.com:rmcgibbo/openmm into docstrings

parents d137b536 d74a4cc1
...@@ -116,12 +116,8 @@ def _strip_optunit(thing, unit): ...@@ -116,12 +116,8 @@ def _strip_optunit(thing, unit):
_resre = re.compile(r'(\d+)([a-zA-Z]*)') _resre = re.compile(r'(\d+)([a-zA-Z]*)')
class CharmmPsfFile(object): class CharmmPsfFile(object):
""" """A chemical structure instantiated from CHARMM files.
A chemical structure instantiated from CHARMM files.
Example:
>>> cs = CharmmPsfFile("testfiles/test.psf")
This structure has numerous attributes that are lists of the elements of This structure has numerous attributes that are lists of the elements of
this structure, including atoms, bonds, torsions, etc. The attributes are this structure, including atoms, bonds, torsions, etc. The attributes are
- residue_list - residue_list
...@@ -138,13 +134,14 @@ class CharmmPsfFile(object): ...@@ -138,13 +134,14 @@ class CharmmPsfFile(object):
Additional attribute is available if a CharmmParameterSet is loaded into Additional attribute is available if a CharmmParameterSet is loaded into
this structure. this structure.
- urey_bradley_list - urey_bradley_list
The lengths of each of these lists gives the pointers (e.g., natom, nres, The lengths of each of these lists gives the pointers (e.g., natom, nres,
etc.) etc.)
Example: Examples
--------
>>> cs = CharmmPsfFile("testfiles/test.psf") >>> cs = CharmmPsfFile("testfiles/test.psf")
>>> len(cs.atom_list) >>> len(cs.atom_list)
33 33
...@@ -163,19 +160,21 @@ class CharmmPsfFile(object): ...@@ -163,19 +160,21 @@ class CharmmPsfFile(object):
CMAP_FORCE_GROUP = 5 CMAP_FORCE_GROUP = 5
NONBONDED_FORCE_GROUP = 6 NONBONDED_FORCE_GROUP = 6
GB_FORCE_GROUP = 6 GB_FORCE_GROUP = 6
@_catchindexerror @_catchindexerror
def __init__(self, psf_name): def __init__(self, psf_name):
""" """Opens and parses a PSF file, then instantiates a CharmmPsfFile
Opens and parses a PSF file, then instantiates a CharmmPsfFile
instance from the data. instance from the data.
Parameters: Parameters
psf_name (str) : Name of the PSF file (it must exist) ----------
psf_name : str
Exceptions Raised: Name of the PSF file (it must exist)
IOError : If file "psf_name" does not exist
CharmmPSFError: If any parsing errors are encountered Raises
------
IOError : If file "psf_name" does not exist
CharmmPSFError: If any parsing errors are encountered
""" """
conv = CharmmPsfFile._convert conv = CharmmPsfFile._convert
# Make sure the file exists # Make sure the file exists
...@@ -388,14 +387,17 @@ class CharmmPsfFile(object): ...@@ -388,14 +387,17 @@ class CharmmPsfFile(object):
@staticmethod @staticmethod
def _convert(string, type, message): def _convert(string, type, message):
""" """Converts a string to a specific type, making sure to raise
Converts a string to a specific type, making sure to raise
CharmmPSFError with the given message in the event of a failure. CharmmPSFError with the given message in the event of a failure.
Parameters: Parameters
- string (str) : Input string to process ----------
- type (type) : Type of data to convert to string : str
- message (str) : Error message to put in exception if failed Input string to process
type : type
Type of data to convert to
message : str
Error message to put in exception if failed
""" """
try: try:
return type(string) return type(string)
...@@ -405,23 +407,24 @@ class CharmmPsfFile(object): ...@@ -405,23 +407,24 @@ class CharmmPsfFile(object):
@staticmethod @staticmethod
def _parse_psf_section(psf): def _parse_psf_section(psf):
""" """This method parses a section of the PSF file
This method parses a section of the PSF file
Parameters
Parameters: ----------
- psf (CharmmFile) : Open file that is pointing to the first line psf : CharmmFile
of the section that is to be parsed Open file that is pointing to the first line of the section
that is to be parsed
Returns:
(title, pointers, data) Returns
--------
- title (str) : The label of the PSF section we are parsing str
- pointers (int/tuple of ints) : If one pointer is set, pointers is The label of the PSF section we are parsing
simply the integer that is value of that pointer. Otherwise int/tuple of ints
it is a tuple with every pointer value defined in the first If one pointer is set, pointers is simply the integer that is
line value of that pointer. Otherwise it is a tuple with every pointer
- data (list) : A list of all data in the parsed section converted value defined in the first line
to `dtype' list
A list of all data in the parsed section converted to `dtype'
""" """
conv = CharmmPsfFile._convert conv = CharmmPsfFile._convert
line = psf.readline() line = psf.readline()
...@@ -462,25 +465,25 @@ class CharmmPsfFile(object): ...@@ -462,25 +465,25 @@ class CharmmPsfFile(object):
return title, pointers, data return title, pointers, data
def loadParameters(self, parmset): def loadParameters(self, parmset):
""" """Loads parameters from a parameter set that was loaded via CHARMM RTF,
Loads parameters from a parameter set that was loaded via CHARMM RTF,
PAR, and STR files. PAR, and STR files.
Parameters: Parameters
- parmset (CharmmParameterSet) : List of all parameters ----------
parmset : CharmmParameterSet
Notes: List of all parameters
- If any parameters that are necessary cannot be found, a
MissingParameter exception is raised. Notes
-----
- If any dihedral or improper parameters cannot be found, I will try - If any parameters that are necessary cannot be found, a
inserting wildcards (at either end for dihedrals and as the two MissingParameter exception is raised.
central atoms in impropers) and see if that matches. Wild-cards - If any dihedral or improper parameters cannot be found, I will try
will apply ONLY if specific parameters cannot be found. inserting wildcards (at either end for dihedrals and as the two
central atoms in impropers) and see if that matches. Wild-cards
- This method will expand the dihedral_parameter_list attribute by will apply ONLY if specific parameters cannot be found.
adding a separate Dihedral object for each term for types that - This method will expand the dihedral_parameter_list attribute by
have a multi-term expansion adding a separate Dihedral object for each term for types that
have a multi-term expansion
""" """
# First load the atom types # First load the atom types
types_are_int = False types_are_int = False
...@@ -588,13 +591,22 @@ class CharmmPsfFile(object): ...@@ -588,13 +591,22 @@ class CharmmPsfFile(object):
def setBox(self, a, b, c, alpha=90.0*u.degrees, beta=90.0*u.degrees, def setBox(self, a, b, c, alpha=90.0*u.degrees, beta=90.0*u.degrees,
gamma=90.0*u.degrees): gamma=90.0*u.degrees):
""" """Sets the periodic box boundary conditions.
Sets the periodic box boundary conditions.
Parameters
Parameters: ----------
- a, b, c (floats) : Lengths of the periodic cell a : length
- alpha, beta, gamma (floats, optional) : Angles between the Lengths of the periodic cell
periodic cell vectors. b : length
Lengths of the periodic cell
c : length
Lengths of the periodic cell
alpha : floats, optional
Angles between the periodic cell vectors.
beta : floats, optional
Angles between the periodic cell vectors.
gamma : floats, optional
Angles between the periodic cell vectors.
""" """
try: try:
# Since we are setting the box, delete the cached box lengths if we # Since we are setting the box, delete the cached box lengths if we
...@@ -620,7 +632,7 @@ class CharmmPsfFile(object): ...@@ -620,7 +632,7 @@ class CharmmPsfFile(object):
pass pass
# Cache the topology for easy returning later # Cache the topology for easy returning later
self._topology = topology = Topology() self._topology = topology = Topology()
last_chain = None last_chain = None
last_residue = None last_residue = None
# Add each chain (separate 'system's) and residue # Add each chain (separate 'system's) and residue
...@@ -752,52 +764,60 @@ class CharmmPsfFile(object): ...@@ -752,52 +764,60 @@ class CharmmPsfFile(object):
ewaldErrorTolerance=0.0005, ewaldErrorTolerance=0.0005,
flexibleConstraints=True, flexibleConstraints=True,
verbose=False): verbose=False):
""" """Construct an OpenMM System representing the topology described by the
Construct an OpenMM System representing the topology described by the
prmtop file. You MUST have loaded a parameter set into this PSF before prmtop file. You MUST have loaded a parameter set into this PSF before
calling createSystem. If not, AttributeError will be raised. ValueError calling createSystem. If not, AttributeError will be raised. ValueError
is raised for illegal input. is raised for illegal input.
Parameters: Parameters
- params (CharmmParameterSet) The parameter set to use to parametrize ----------
this molecule params : CharmmParameterSet
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded The parameter set to use to parametrize this molecule
interactions. Allowed values are NoCutoff, CutoffNonPeriodic, nonbondedMethod : object=NoCutoff
CutoffPeriodic, Ewald, or PME. The method to use for nonbonded interactions. Allowed values are
- nonbondedCutoff (distance=1*nanometer) The cutoff distance to use NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
for nonbonded interactions. nonbondedCutoff : distance=1*nanometer
- switchDistance (distance=0*nanometer) The distance at which the The cutoff distance to use for nonbonded interactions.
switching function is active for nonbonded interactions. If the switchDistance : distance=0*nanometer
switchDistance evaluates to boolean False (if it is 0), no The distance at which the switching function is active for nonbonded
switching function will be used. Illegal values will raise a interactions. If the switchDistance evaluates to boolean False (if
ValueError it is 0), no switching function will be used. Illegal values will
- constraints (object=None) Specifies which bonds or angles should be raise a ValueError
implemented with constraints. Allowed values are None, HBonds, constraints : object=None
AllBonds, or HAngles. Specifies which bonds or angles should be implemented with
- rigidWater (boolean=True) If true, water molecules will be fully constraints. Allowed values are None, HBonds, AllBonds, or HAngles.
rigid regardless of the value passed for the constraints argument rigidWater : boolean=True
- implicitSolvent (object=None) If not None, the implicit solvent If true, water molecules will be fully rigid regardless of the value
model to use. Allowed values are HCT, OBC1, OBC2, or GBn passed for the constraints argument
- implicitSolventKappa (float=None): Debye screening parameter to implicitSolvent : object=None
model salt concentrations in GB solvent. If not None, the implicit solvent model to use. Allowed values are
- implicitSolventSaltConc (float=0.0*u.moles/u.liter): Salt HCT, OBC1, OBC2, or GBn
concentration for GB simulations. Converted to Debye length implicitSolventKappa : float=None
`kappa' Debye screening parameter to model salt concentrations in GB
- temperature (float=298.15*u.kelvin): Temperature used in the salt solvent.
concentration-to-kappa conversion for GB salt concentration term implicitSolventSaltConc : float=0.0*u.moles/u.liter
- soluteDielectric (float=1.0) The solute dielectric constant to use Salt concentration for GB simulations. Converted to Debye length
in the implicit solvent model. `kappa'
- solventDielectric (float=78.5) The solvent dielectric constant to temperature : float=298.15*u.kelvin
use in the implicit solvent model. Temperature used in the salt concentration-to-kappa conversion for
- removeCMMotion (boolean=True) If true, a CMMotionRemover will be GB salt concentration term
added to the System. soluteDielectric : float=1.0
- hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to The solute dielectric constant to use in the implicit solvent model.
heavy atoms. Any mass added to a hydrogen is subtracted from the solventDielectric : float=78.5
heavy atom to keep their total mass the same. The solvent dielectric constant to use in the implicit solvent
- ewaldErrorTolerance (float=0.0005) The error tolerance to use if the model.
nonbonded method is Ewald or PME. removeCMMotion : boolean=True
- flexibleConstraints (bool=True) Are our constraints flexible or not? If true, a CMMotionRemover will be added to the System.
- verbose (bool=False) Optionally prints out a running progress report hydrogenMass : mass=None
The mass to use for hydrogen atoms bound to heavy atoms. Any mass
added to a hydrogen is subtracted from the heavy atom to keep their
total mass the same.
ewaldErrorTolerance : float=0.0005
The error tolerance to use if the nonbonded method is Ewald or PME.
flexibleConstraints : bool=True
Are our constraints flexible or not?
verbose : bool=False
Optionally prints out a running progress report
""" """
# Load the parameter set # Load the parameter set
self.loadParameters(params.condense()) self.loadParameters(params.condense())
...@@ -821,7 +841,7 @@ class CharmmPsfFile(object): ...@@ -821,7 +841,7 @@ class CharmmPsfFile(object):
raise ValueError('Illegal implicit solvent model choice.') raise ValueError('Illegal implicit solvent model choice.')
if not constraints in (None, ff.HAngles, ff.HBonds, ff.AllBonds): if not constraints in (None, ff.HAngles, ff.HBonds, ff.AllBonds):
raise ValueError('Illegal constraints choice') raise ValueError('Illegal constraints choice')
# Define conversion factors # Define conversion factors
length_conv = u.angstrom.conversion_factor_to(u.nanometer) length_conv = u.angstrom.conversion_factor_to(u.nanometer)
_chmfrc = u.kilocalorie_per_mole/(u.angstrom*u.angstrom) _chmfrc = u.kilocalorie_per_mole/(u.angstrom*u.angstrom)
...@@ -833,7 +853,7 @@ class CharmmPsfFile(object): ...@@ -833,7 +853,7 @@ class CharmmPsfFile(object):
dihe_frc_conv = u.kilocalorie_per_mole.conversion_factor_to( dihe_frc_conv = u.kilocalorie_per_mole.conversion_factor_to(
u.kilojoule_per_mole) u.kilojoule_per_mole)
ene_conv = dihe_frc_conv ene_conv = dihe_frc_conv
# Create the system and determine if any of our atoms have NBFIX (and # Create the system and determine if any of our atoms have NBFIX (and
# therefore requires a CustomNonbondedForce instead) # therefore requires a CustomNonbondedForce instead)
typenames = set() typenames = set()
...@@ -1358,7 +1378,7 @@ class CharmmPsfFile(object): ...@@ -1358,7 +1378,7 @@ class CharmmPsfFile(object):
def boxLengths(self, stuff): def boxLengths(self, stuff):
raise RuntimeError('Use setBox to set a box with lengths and angles ' raise RuntimeError('Use setBox to set a box with lengths and angles '
'or set the boxVectors attribute with box vectors') 'or set the boxVectors attribute with box vectors')
@property @property
def boxVectors(self): def boxVectors(self):
""" Return the box vectors """ """ Return the box vectors """
...@@ -1403,12 +1423,12 @@ def set_molecules(atom_list): ...@@ -1403,12 +1423,12 @@ def set_molecules(atom_list):
owner = [] owner = []
# The way I do this is via a recursive algorithm, in which # The way I do this is via a recursive algorithm, in which
# the "set_owner" method is called for each bonded partner an atom # the "set_owner" method is called for each bonded partner an atom
# has, which in turn calls set_owner for each of its partners and # has, which in turn calls set_owner for each of its partners and
# so on until everything has been assigned. # so on until everything has been assigned.
molecule_number = 1 # which molecule number we are on molecule_number = 1 # which molecule number we are on
for i in range(len(atom_list)): for i in range(len(atom_list)):
# If this atom has not yet been "owned", make it the next molecule # If this atom has not yet been "owned", make it the next molecule
# However, we only increment which molecule number we're on if # However, we only increment which molecule number we're on if
# we actually assigned a new molecule (obviously) # we actually assigned a new molecule (obviously)
if not atom_list[i].marked: if not atom_list[i].marked:
tmp = [i] tmp = [i]
...@@ -1429,7 +1449,7 @@ def _set_owner(atom_list, owner_array, atm, mol_id): ...@@ -1429,7 +1449,7 @@ def _set_owner(atom_list, owner_array, atm, mol_id):
owner_array.append(partner.idx) owner_array.append(partner.idx)
_set_owner(atom_list, owner_array, partner.idx, mol_id) _set_owner(atom_list, owner_array, partner.idx, mol_id)
elif partner.marked != mol_id: elif partner.marked != mol_id:
raise MoleculeError('Atom %d in multiple molecules' % raise MoleculeError('Atom %d in multiple molecules' %
partner.idx) partner.idx)
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......
...@@ -8,7 +8,7 @@ Medical Research, grant U54 GM072970. See https://simtk.org. ...@@ -8,7 +8,7 @@ Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2014 Stanford University and the Authors. Portions copyright (c) 2014 Stanford University and the Authors.
Authors: Robert McGibbon Authors: Robert McGibbon
Contributors: Contributors:
Permission is hereby granted, free of charge, to any person obtaining a Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"), copy of this software and associated documentation files (the "Software"),
...@@ -71,10 +71,12 @@ class CheckpointReporter(object): ...@@ -71,10 +71,12 @@ class CheckpointReporter(object):
def __init__(self, file, reportInterval): def __init__(self, file, reportInterval):
"""Create a CheckpointReporter. """Create a CheckpointReporter.
Parameters: Parameters
- file (string or open file object) The file to write to. Any current ----------
contents will be overwritten. file : string or open file object
- reportInterval (int) The interval (in time steps) at which to write checkpoints The file to write to. Any current contents will be overwritten.
reportInterval : int
The interval (in time steps) at which to write checkpoints.
""" """
self._reportInterval = reportInterval self._reportInterval = reportInterval
...@@ -88,12 +90,23 @@ class CheckpointReporter(object): ...@@ -88,12 +90,23 @@ class CheckpointReporter(object):
def describeNextReport(self, simulation): def describeNextReport(self, simulation):
"""Get information about the next report this object will generate. """Get information about the next report this object will generate.
Parameters: Parameters
- simulation (Simulation) The Simulation to generate a report for ----------
simulation : Simulation
Returns: A five element tuple. The first element is the number of steps until the The Simulation to generate a report for
next report. The remaining elements specify whether that report will require
positions, velocities, forces, and energies respectively. Returns
-------
int
The number of steps until the
bool
Requires positions
bool
Requires velocities
bool
Requires forces
bool
Requires energies
""" """
steps = self._reportInterval - simulation.currentStep%self._reportInterval steps = self._reportInterval - simulation.currentStep%self._reportInterval
return (steps, False, False, False, False) return (steps, False, False, False, False)
...@@ -101,9 +114,12 @@ class CheckpointReporter(object): ...@@ -101,9 +114,12 @@ class CheckpointReporter(object):
def report(self, simulation, state): def report(self, simulation, state):
"""Generate a report. """Generate a report.
Parameters: Parameters
- simulation (Simulation) The Simulation to generate a report for ----------
- state (State) The current state of the simulation simulation : Simulation
The Simulation to generate a report for
state : State
The current state of the simulation
""" """
self._out.seek(0) self._out.seek(0)
chk = simulation.context.createCheckpoint() chk = simulation.context.createCheckpoint()
......
...@@ -55,12 +55,19 @@ class DCDFile(object): ...@@ -55,12 +55,19 @@ class DCDFile(object):
def __init__(self, file, topology, dt, firstStep=0, interval=1): def __init__(self, file, topology, dt, firstStep=0, interval=1):
"""Create a DCD file and write out the header. """Create a DCD file and write out the header.
Parameters: Parameters
- file (file) A file to write to ----------
- topology (Topology) The Topology defining the molecular system being written file : file
- dt (time) The time step used in the trajectory A file to write to
- firstStep (int=0) The index of the first step in the trajectory topology : Topology
- interval (int=1) The frequency (measured in time steps) at which states are written to the trajectory The Topology defining the molecular system being written
dt : time
The time step used in the trajectory
firstStep : int=0
The index of the first step in the trajectory
interval : int=1
The frequency (measured in time steps) at which states are written
to the trajectory
""" """
self._file = file self._file = file
self._topology = topology self._topology = topology
...@@ -83,14 +90,21 @@ class DCDFile(object): ...@@ -83,14 +90,21 @@ class DCDFile(object):
def writeModel(self, positions, unitCellDimensions=None, periodicBoxVectors=None): def writeModel(self, positions, unitCellDimensions=None, periodicBoxVectors=None):
"""Write out a model to the DCD file. """Write out a model to the DCD file.
The periodic box can be specified either by the unit cell dimensions (for a rectangular box), or the full set of box The periodic box can be specified either by the unit cell dimensions
vectors (for an arbitrary triclinic box). If neither is specified, the box vectors specified in the Topology will be (for a rectangular box), or the full set of box vectors (for an
used. Regardless of the value specified, no dimensions will be written if the Topology does not represent a periodic system. arbitrary triclinic box). If neither is specified, the box vectors
specified in the Topology will be used. Regardless of the value
Parameters: specified, no dimensions will be written if the Topology does not
- positions (list) The list of atomic positions to write represent a periodic system.
- unitCellDimensions (Vec3=None) The dimensions of the crystallographic unit cell.
- periodicBoxVectors (tuple of Vec3=None) The vectors defining the periodic box. Parameters
----------
positions : list
The list of atomic positions to write
unitCellDimensions : Vec3=None
The dimensions of the crystallographic unit cell.
periodicBoxVectors : tuple of Vec3=None
The vectors defining the periodic box.
""" """
if len(list(self._topology.atoms())) != len(positions): if len(list(self._topology.atoms())) != len(positions):
raise ValueError('The number of positions must match the number of atoms') raise ValueError('The number of positions must match the number of atoms')
......
...@@ -45,9 +45,12 @@ class DCDReporter(object): ...@@ -45,9 +45,12 @@ class DCDReporter(object):
def __init__(self, file, reportInterval): def __init__(self, file, reportInterval):
"""Create a DCDReporter. """Create a DCDReporter.
Parameters: Parameters
- file (string) The file to write to ----------
- reportInterval (int) The interval (in time steps) at which to write frames file : string
The file to write to
reportInterval : int
The interval (in time steps) at which to write frames
""" """
self._reportInterval = reportInterval self._reportInterval = reportInterval
self._out = open(file, 'wb') self._out = open(file, 'wb')
...@@ -56,11 +59,23 @@ class DCDReporter(object): ...@@ -56,11 +59,23 @@ class DCDReporter(object):
def describeNextReport(self, simulation): def describeNextReport(self, simulation):
"""Get information about the next report this object will generate. """Get information about the next report this object will generate.
Parameters: Parameters
- simulation (Simulation) The Simulation to generate a report for ----------
Returns: A five element tuple. The first element is the number of steps until the simulation : Simulation
next report. The remaining elements specify whether that report will require The Simulation to generate a report for
positions, velocities, forces, and energies respectively.
Returns
-------
int
The number of steps until the
bool
Requires positions
bool
Requires velocities
bool
Requires forces
bool
Requires energies
""" """
steps = self._reportInterval - simulation.currentStep%self._reportInterval steps = self._reportInterval - simulation.currentStep%self._reportInterval
return (steps, True, False, False, False) return (steps, True, False, False, False)
...@@ -68,10 +83,14 @@ class DCDReporter(object): ...@@ -68,10 +83,14 @@ class DCDReporter(object):
def report(self, simulation, state): def report(self, simulation, state):
"""Generate a report. """Generate a report.
Parameters: Parameters
- simulation (Simulation) The Simulation to generate a report for ----------
- state (State) The current state of the simulation simulation : Simulation
The Simulation to generate a report for
state : State
The current state of the simulation
""" """
if self._dcd is None: if self._dcd is None:
self._dcd = DCDFile(self._out, simulation.topology, simulation.integrator.getStepSize(), 0, self._reportInterval) self._dcd = DCDFile(self._out, simulation.topology, simulation.integrator.getStepSize(), 0, self._reportInterval)
a,b,c = state.getPeriodicBoxVectors() a,b,c = state.getPeriodicBoxVectors()
......
...@@ -46,10 +46,11 @@ class DesmondDMSFile(object): ...@@ -46,10 +46,11 @@ class DesmondDMSFile(object):
def __init__(self, file): def __init__(self, file):
"""Load a DMS file """Load a DMS file
Parameters: Parameters
- file (string) the name of the file to load ----------
file : string
the name of the file to load
""" """
# sqlite3 is included in the standard lib, but at python # sqlite3 is included in the standard lib, but at python
# compile time, you can disable support (I think), so it's # compile time, you can disable support (I think), so it's
# not *guarenteed* to be available. Doing the import here # not *guarenteed* to be available. Doing the import here
...@@ -157,16 +158,24 @@ class DesmondDMSFile(object): ...@@ -157,16 +158,24 @@ class DesmondDMSFile(object):
def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*nanometer, def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*nanometer,
ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None): ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None):
"""Construct an OpenMM System representing the topology described by this dms file """Construct an OpenMM System representing the topology described by this
DMS file
Parameters:
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are Parameters
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. ----------
- nonbondedCutoff (distance=1*nanometer) The cutoff distance to use for nonbonded interactions nonbondedMethod : object=NoCutoff
- ewaldErrorTolerance (float=0.0005) The error tolerance to use if nonbondedMethod is Ewald or PME. The method to use for nonbonded interactions. Allowed values are
- removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
- hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to heavy atoms. Any mass added to a hydrogen is nonbondedCutoff : distance=1*nanometer
subtracted from the heavy atom to keep their total mass the same. The cutoff distance to use for nonbonded interactions
ewaldErrorTolerance : float=0.0005
The error tolerance to use if nonbondedMethod is Ewald or PME.
removeCMMotion : boolean=True
If true, a CMMotionRemover will be added to the System
hydrogenMass : mass=None
The mass to use for hydrogen atoms bound to heavy atoms. Any mass
added to a hydrogen is subtracted from the heavy atom to keep their
total mass the same.
""" """
self._checkForUnsupportedTerms() self._checkForUnsupportedTerms()
sys = mm.System() sys = mm.System()
......
...@@ -58,11 +58,16 @@ class Element(object): ...@@ -58,11 +58,16 @@ class Element(object):
def __init__(self, number, name, symbol, mass): def __init__(self, number, name, symbol, mass):
"""Create a new element """Create a new element
Parameters: Parameters
number (int) The atomic number of the element ----------
name (string) The name of the element number : int
symbol (string) The chemical symbol of the element The atomic number of the element
mass (float) The atomic mass of the element name : string
The name of the element
symbol : string
The chemical symbol of the element
mass : float
The atomic mass of the element
""" """
## The atomic number of the element ## The atomic number of the element
self._atomic_number = number self._atomic_number = number
...@@ -115,7 +120,7 @@ class Element(object): ...@@ -115,7 +120,7 @@ class Element(object):
Returns Returns
------- -------
element : Element Element
The element whose atomic mass is closest to the input mass The element whose atomic mass is closest to the input mass
""" """
# Assume masses are in daltons if they are not units # Assume masses are in daltons if they are not units
......
...@@ -104,12 +104,14 @@ class ForceField(object): ...@@ -104,12 +104,14 @@ class ForceField(object):
def __init__(self, *files): def __init__(self, *files):
"""Load one or more XML files and create a ForceField object based on them. """Load one or more XML files and create a ForceField object based on them.
Parameters: Parameters
- files (list) A list of XML files defining the force field. Each entry may ----------
be an absolute file path, a path relative to the current working files : list
directory, a path relative to this module's data subdirectory A list of XML files defining the force field. Each entry may
(for built in force fields), or an open file-like object with a be an absolute file path, a path relative to the current working
read() method from which the forcefield XML data can be loaded. directory, a path relative to this module's data subdirectory
(for built in force fields), or an open file-like object with a
read() method from which the forcefield XML data can be loaded.
""" """
self._atomTypes = {} self._atomTypes = {}
self._templates = {} self._templates = {}
...@@ -123,12 +125,14 @@ class ForceField(object): ...@@ -123,12 +125,14 @@ class ForceField(object):
def loadFile(self, file): def loadFile(self, file):
"""Load an XML file and add the definitions from it to this FieldField. """Load an XML file and add the definitions from it to this FieldField.
Parameters: Parameters
- file (string or file) An XML file containing force field definitions. It may ----------
be either an absolute file path, a path relative to the current working file : string or file
directory, a path relative to this module's data subdirectory An XML file containing force field definitions. It may be either an
(for built in force fields), or an open file-like object with a absolute file path, a path relative to the current working
read() method from which the forcefield XML data can be loaded. directory, a path relative to this module's data subdirectory (for
built in force fields), or an open file-like object with a read()
method from which the forcefield XML data can be loaded.
""" """
try: try:
# this handles either filenames or open file-like objects # this handles either filenames or open file-like objects
...@@ -417,20 +421,36 @@ class ForceField(object): ...@@ -417,20 +421,36 @@ class ForceField(object):
constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args): constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args):
"""Construct an OpenMM System representing a Topology with this force field. """Construct an OpenMM System representing a Topology with this force field.
Parameters: Parameters
- topology (Topology) The Topology for which to create a System ----------
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are topology : Topology
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. The Topology for which to create a System
- nonbondedCutoff (distance=1*nanometer) The cutoff distance to use for nonbonded interactions nonbondedMethod : object=NoCutoff
- constraints (object=None) Specifies which bonds and angles should be implemented with constraints. The method to use for nonbonded interactions. Allowed values are
Allowed values are None, HBonds, AllBonds, or HAngles. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
- rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument nonbondedCutoff : distance=1*nanometer
- removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System The cutoff distance to use for nonbonded interactions
- hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to heavy atoms. Any mass added to a hydrogen is constraints : object=None
subtracted from the heavy atom to keep their total mass the same. Specifies which bonds and angles should be implemented with constraints.
- args Arbitrary additional keyword arguments may also be specified. This allows extra parameters to be specified that are specific to Allowed values are None, HBonds, AllBonds, or HAngles.
particular force fields. rigidWater : boolean=True
Returns: the newly created System If true, water molecules will be fully rigid regardless of the value
passed for the constraints argument
removeCMMotion : boolean=True
If true, a CMMotionRemover will be added to the System
hydrogenMass : mass=None
The mass to use for hydrogen atoms bound to heavy atoms. Any mass
added to a hydrogen is subtracted from the heavy atom to keep
their total mass the same.
args
Arbitrary additional keyword arguments may also be specified.
This allows extra parameters to be specified that are specific to
particular force fields.
Returns
-------
system
the newly created System
""" """
data = ForceField._SystemData() data = ForceField._SystemData()
data.atoms = list(topology.atoms()) data.atoms = list(topology.atoms())
...@@ -658,12 +678,20 @@ def _createResidueSignature(elements): ...@@ -658,12 +678,20 @@ def _createResidueSignature(elements):
def _matchResidue(res, template, bondedToAtom): def _matchResidue(res, template, bondedToAtom):
"""Determine whether a residue matches a template and return a list of corresponding atoms. """Determine whether a residue matches a template and return a list of corresponding atoms.
Parameters: Parameters
- res (Residue) The residue to check ----------
- template (_TemplateData) The template to compare it to res : Residue
- bondedToAtom (list) Enumerates which other atoms each atom is bonded to The residue to check
Returns: a list specifying which atom of the template each atom of the residue corresponds to, template : _TemplateData
or None if it does not match the template The template to compare it to
bondedToAtom : list
Enumerates which other atoms each atom is bonded to
Returns
-------
list
a list specifying which atom of the template each atom of the residue
corresponds to, or None if it does not match the template
""" """
atoms = list(res.atoms()) atoms = list(res.atoms())
if len(atoms) != len(template.atoms): if len(atoms) != len(template.atoms):
......
...@@ -114,10 +114,11 @@ class GromacsGroFile(object): ...@@ -114,10 +114,11 @@ class GromacsGroFile(object):
The atom positions can be retrieved by calling getPositions(). The atom positions can be retrieved by calling getPositions().
Parameters: Parameters
- file (string) the name of the file to load ----------
file : string
the name of the file to load
""" """
xyzs = [] xyzs = []
elements = [] # The element, most useful for quantum chemistry calculations elements = [] # The element, most useful for quantum chemistry calculations
atomname = [] # The atom name, for instance 'HW1' atomname = [] # The atom name, for instance 'HW1'
...@@ -183,10 +184,14 @@ class GromacsGroFile(object): ...@@ -183,10 +184,14 @@ class GromacsGroFile(object):
def getPositions(self, asNumpy=False, frame=0): def getPositions(self, asNumpy=False, frame=0):
"""Get the atomic positions. """Get the atomic positions.
Parameters: Parameters
- asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s ----------
- frame (int=0) the index of the frame for which to get positions asNumpy : boolean=False
""" if true, the values are returned as a numpy array instead of a list
of Vec3s
frame : int=0
the index of the frame for which to get positions
"""
if asNumpy: if asNumpy:
if self._numpyPositions is None: if self._numpyPositions is None:
self._numpyPositions = [None]*len(self._positions) self._numpyPositions = [None]*len(self._positions)
...@@ -198,16 +203,20 @@ class GromacsGroFile(object): ...@@ -198,16 +203,20 @@ class GromacsGroFile(object):
def getPeriodicBoxVectors(self, frame=0): def getPeriodicBoxVectors(self, frame=0):
"""Get the vectors defining the periodic box. """Get the vectors defining the periodic box.
Parameters: Parameters
- frame (int=0) the index of the frame for which to get the box vectors ----------
frame : int=0
the index of the frame for which to get the box vectors
""" """
return self._periodicBoxVectors[frame] return self._periodicBoxVectors[frame]
def getUnitCellDimensions(self, frame=0): def getUnitCellDimensions(self, frame=0):
"""Get the dimensions of the crystallographic unit cell. """Get the dimensions of the crystallographic unit cell.
Parameters: Parameters
- frame (int=0) the index of the frame for which to get the unit cell dimensions ----------
frame : int=0
the index of the frame for which to get the unit cell dimensions
""" """
xsize = self._periodicBoxVectors[frame][0][0].value_in_unit(nanometers) xsize = self._periodicBoxVectors[frame][0][0].value_in_unit(nanometers)
ysize = self._periodicBoxVectors[frame][1][1].value_in_unit(nanometers) ysize = self._periodicBoxVectors[frame][1][1].value_in_unit(nanometers)
......
...@@ -433,16 +433,22 @@ class GromacsTopFile(object): ...@@ -433,16 +433,22 @@ class GromacsTopFile(object):
def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None, includeDir=None, defines=None): def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None, includeDir=None, defines=None):
"""Load a top file. """Load a top file.
Parameters: Parameters
- file (string) the name of the file to load ----------
- periodicBoxVectors (tuple of Vec3=None) the vectors defining the periodic box file : str
- unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell. For the name of the file to load
non-rectangular unit cells, specify periodicBoxVectors instead. periodicBoxVectors : tuple of Vec3=None
- includeDir (string=None) A directory in which to look for other files the vectors defining the periodic box
included from the top file. If not specified, we will attempt to locate a gromacs unitCellDimensions : Vec3=None
installation on your system. When gromacs is installed in /usr/local, this will resolve the dimensions of the crystallographic unit cell. For
to /usr/local/gromacs/share/gromacs/top non-rectangular unit cells, specify periodicBoxVectors instead.
- defines (dict={}) preprocessor definitions that should be predefined when parsing the file includeDir : string=None
A directory in which to look for other files included from the
top file. If not specified, we will attempt to locate a gromacs
installation on your system. When gromacs is installed in
/usr/local, this will resolve to /usr/local/gromacs/share/gromacs/top
defines : dict={}
preprocessor definitions that should be predefined when parsing the file
""" """
if includeDir is None: if includeDir is None:
includeDir = _defaultGromacsIncludeDir() includeDir = _defaultGromacsIncludeDir()
...@@ -539,23 +545,43 @@ class GromacsTopFile(object): ...@@ -539,23 +545,43 @@ class GromacsTopFile(object):
def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer, def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None): constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None):
"""Construct an OpenMM System representing the topology described by this prmtop file. """Construct an OpenMM System representing the topology described by this
prmtop file.
Parameters:
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are Parameters
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. ----------
- nonbondedCutoff (distance=1*nanometer) The cutoff distance to use for nonbonded interactions nonbondedMethod : object=NoCutoff
- constraints (object=None) Specifies which bonds and angles should be implemented with constraints. The method to use for nonbonded interactions. Allowed values are
Allowed values are None, HBonds, AllBonds, or HAngles. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
- rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument nonbondedCutoff : distance=1*nanometer
- implicitSolvent (object=None) If not None, the implicit solvent model to use. The only allowed value is OBC2. The cutoff distance to use for nonbonded interactions
- soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model. constraints : object=None
- solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model. Specifies which bonds and angles should be implemented with
- ewaldErrorTolerance (float=0.0005) The error tolerance to use if nonbondedMethod is Ewald or PME. constraints. Allowed values are None, HBonds, AllBonds, or HAngles.
- removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System rigidWater : boolean=True
- hydrogenMass (mass=None) The mass to use for hydrogen atoms bound to heavy atoms. Any mass added to a hydrogen is If true, water molecules will be fully rigid regardless of the value
subtracted from the heavy atom to keep their total mass the same. passed for the constraints argument
Returns: the newly created System implicitSolvent : object=None
If not None, the implicit solvent model to use. The only allowed
value is OBC2.
soluteDielectric : float=1.0
The solute dielectric constant to use in the implicit solvent model.
solventDielectric : float=78.5
The solvent dielectric constant to use in the implicit solvent
model.
ewaldErrorTolerance : float=0.0005
The error tolerance to use if nonbondedMethod is Ewald or PME.
removeCMMotion : boolean=True
If true, a CMMotionRemover will be added to the System
hydrogenMass : mass=None
The mass to use for hydrogen atoms bound to heavy atoms. Any mass
added to a hydrogen is subtracted from the heavy atom to keep their
total mass the same.
Returns
-------
System
the newly created System
""" """
# Create the System. # Create the System.
...@@ -586,9 +612,9 @@ class GromacsTopFile(object): ...@@ -586,9 +612,9 @@ class GromacsTopFile(object):
topologyAtoms = list(self.topology.atoms()) topologyAtoms = list(self.topology.atoms())
exceptions = [] exceptions = []
fudgeQQ = float(self._defaults[4]) fudgeQQ = float(self._defaults[4])
# Build a lookup table to let us process dihedrals more quickly. # Build a lookup table to let us process dihedrals more quickly.
dihedralTypeTable = {} dihedralTypeTable = {}
for key in self._dihedralTypes: for key in self._dihedralTypes:
if key[1] != 'X' and key[2] != 'X': if key[1] != 'X' and key[2] != 'X':
...@@ -837,7 +863,7 @@ class GromacsTopFile(object): ...@@ -837,7 +863,7 @@ class GromacsTopFile(object):
for atom in atoms[1:]: for atom in atoms[1:]:
if atom > atoms[0]: if atom > atoms[0]:
exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom, 0, 0, 0)) exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atom, 0, 0, 0))
# Create nonbonded exceptions. # Create nonbonded exceptions.
...@@ -855,9 +881,9 @@ class GromacsTopFile(object): ...@@ -855,9 +881,9 @@ class GromacsTopFile(object):
nb.setNonbondedMethod(methodMap[nonbondedMethod]) nb.setNonbondedMethod(methodMap[nonbondedMethod])
nb.setCutoffDistance(nonbondedCutoff) nb.setCutoffDistance(nonbondedCutoff)
nb.setEwaldErrorTolerance(ewaldErrorTolerance) nb.setEwaldErrorTolerance(ewaldErrorTolerance)
# Adjust masses. # Adjust masses.
if hydrogenMass is not None: if hydrogenMass is not None:
for atom1, atom2 in self.topology.bonds(): for atom1, atom2 in self.topology.bonds():
if atom1.element == elem.hydrogen: if atom1.element == elem.hydrogen:
......
...@@ -123,15 +123,18 @@ class PdbStructure(object): ...@@ -123,15 +123,18 @@ class PdbStructure(object):
""" """
def __init__(self, input_stream, load_all_models = False): def __init__(self, input_stream, load_all_models=False):
"""Create a PDB model from a PDB file stream. """Create a PDB model from a PDB file stream.
Parameters: Parameters
- self (PdbStructure) The new object that is created. ----------
- input_stream (stream) An input file stream, probably created with self : PdbStructure
open(). The new object that is created.
- load_all_models (bool) Whether to load every model of an NMR input_stream : stream
structure or trajectory, or just load the first model, to save memory. An input file stream, probably created with open().
load_all_models : bool
Whether to load every model of an NMR structure or trajectory, or
just load the first model, to save memory.
""" """
# initialize models # initialize models
self.load_all_models = load_all_models self.load_all_models = load_all_models
...@@ -201,7 +204,7 @@ class PdbStructure(object): ...@@ -201,7 +204,7 @@ class PdbStructure(object):
def _reset_atom_numbers(self): def _reset_atom_numbers(self):
self._atom_numbers_are_hex = False self._atom_numbers_are_hex = False
self._next_atom_number = 1 self._next_atom_number = 1
def _reset_residue_numbers(self): def _reset_residue_numbers(self):
self._residue_numbers_are_hex = False self._residue_numbers_are_hex = False
self._next_residue_number = 1 self._next_residue_number = 1
...@@ -269,8 +272,11 @@ class PdbStructure(object): ...@@ -269,8 +272,11 @@ class PdbStructure(object):
Iterate over atomic positions. Iterate over atomic positions.
Parameters Parameters
- use_all_models (bool=False) Get positions from all models or just the first one. ----------
- include_alt_loc (bool=False) Get all positions for each atom, or just the first one. use_all_models : bool=False
Get positions from all models or just the first one.
include_alt_loc : bool=False
Get all positions for each atom, or just the first one.
""" """
for model in self.iter_models(use_all_models): for model in self.iter_models(use_all_models):
for loc in model.iter_positions(include_alt_loc): for loc in model.iter_positions(include_alt_loc):
...@@ -711,7 +717,7 @@ class Atom(object): ...@@ -711,7 +717,7 @@ class Atom(object):
self.is_first_atom_in_chain = False self.is_first_atom_in_chain = False
self.is_final_atom_in_chain = False self.is_final_atom_in_chain = False
self.is_first_residue_in_chain = False self.is_first_residue_in_chain = False
self.is_final_residue_in_chain = False self.is_final_residue_in_chain = False
# Start parsing fields from pdb line # Start parsing fields from pdb line
self.record_name = pdb_line[0:6].strip() self.record_name = pdb_line[0:6].strip()
if pdbstructure is not None and pdbstructure._atom_numbers_are_hex: if pdbstructure is not None and pdbstructure._atom_numbers_are_hex:
......
...@@ -63,9 +63,12 @@ class Modeller(object): ...@@ -63,9 +63,12 @@ class Modeller(object):
def __init__(self, topology, positions): def __init__(self, topology, positions):
"""Create a new Modeller object """Create a new Modeller object
Parameters: Parameters
- topology (Topology) the initial Topology of the model ----------
- positions (list) the initial atomic positions topology : Topology
the initial Topology of the model
positions : list
the initial atomic positions
""" """
## The Topology describing the structure of the system ## The Topology describing the structure of the system
self.topology = topology self.topology = topology
...@@ -85,12 +88,16 @@ class Modeller(object): ...@@ -85,12 +88,16 @@ class Modeller(object):
def add(self, addTopology, addPositions): def add(self, addTopology, addPositions):
"""Add chains, residues, atoms, and bonds to the model. """Add chains, residues, atoms, and bonds to the model.
Specify what to add by providing a new Topology object and the corresponding atomic positions. Specify what to add by providing a new Topology object and the
All chains, residues, atoms, and bonds contained in the Topology are added to the model. corresponding atomic positions. All chains, residues, atoms, and bonds
contained in the Topology are added to the model.
Parameters: Parameters
- addTopoology (Topology) a Topology whose contents should be added to the model ----------
- addPositions (list) the positions of the atoms to add addTopology : Topology
a Topology whose contents should be added to the model
addPositions : list
the positions of the atoms to add
""" """
# Copy over the existing model. # Copy over the existing model.
...@@ -137,8 +144,11 @@ class Modeller(object): ...@@ -137,8 +144,11 @@ class Modeller(object):
You also can specify a bond (as a tuple of Atom objects) to delete just that bond without You also can specify a bond (as a tuple of Atom objects) to delete just that bond without
deleting the Atoms it connects. deleting the Atoms it connects.
Parameters: Parameters
- toDelete (list) a list of Atoms, Residues, Chains, and bonds (specified as tuples of Atoms) to delete ----------
toDelete : list
a list of Atoms, Residues, Chains, and bonds (specified as tuples of
Atoms) to delete
""" """
newTopology = Topology() newTopology = Topology()
newTopology.setPeriodicBoxVectors(self.topology.getPeriodicBoxVectors()) newTopology.setPeriodicBoxVectors(self.topology.getPeriodicBoxVectors())
...@@ -176,10 +186,15 @@ class Modeller(object): ...@@ -176,10 +186,15 @@ class Modeller(object):
def convertWater(self, model='tip3p'): def convertWater(self, model='tip3p'):
"""Convert all water molecules to a different water model. """Convert all water molecules to a different water model.
Parameters: Parameters
- model (string='tip3p') the water model to convert to. Supported values are 'tip3p', 'spce', 'tip4pew', and 'tip5p'. ----------
model : string='tip3p'
@deprecated Use addExtraParticles() instead. It performs the same function but in a more general way. the water model to convert to. Supported values are 'tip3p',
'spce', 'tip4pew', and 'tip5p'.
@deprecated Use addExtraParticles() instead. It performs the same
function but in a more general way.
""" """
if model in ('tip3p', 'spce'): if model in ('tip3p', 'spce'):
sites = 3 sites = 3
...@@ -245,6 +260,7 @@ class Modeller(object): ...@@ -245,6 +260,7 @@ class Modeller(object):
"""Add solvent (both water and ions) to the model to fill a rectangular box. """Add solvent (both water and ions) to the model to fill a rectangular box.
The algorithm works as follows: The algorithm works as follows:
1. Water molecules are added to fill the box. 1. Water molecules are added to fill the box.
2. Water molecules are removed if their distance to any solute atom is less than the sum of their van der Waals radii. 2. Water molecules are removed if their distance to any solute atom is less than the sum of their van der Waals radii.
3. If the solute is charged and neutralize=True, enough positive or negative ions are added to neutralize it. Each ion is added by 3. If the solute is charged and neutralize=True, enough positive or negative ions are added to neutralize it. Each ion is added by
...@@ -252,28 +268,39 @@ class Modeller(object): ...@@ -252,28 +268,39 @@ class Modeller(object):
4. Ion pairs are added to give the requested total ionic strength. 4. Ion pairs are added to give the requested total ionic strength.
The box size can be specified in any of several ways: The box size can be specified in any of several ways:
1. You can explicitly give the vectors defining the periodic box to use. 1. You can explicitly give the vectors defining the periodic box to use.
2. Alternatively, for a rectangular box you can simply give the dimensions of the unit cell. 2. Alternatively, for a rectangular box you can simply give the dimensions of the unit cell.
3. You can give a padding distance. The largest dimension of the solute (along the x, y, or z axis) is determined, and a cubic 3. You can give a padding distance. The largest dimension of the solute (along the x, y, or z axis) is determined, and a cubic
box of size (largest dimension)+2*padding is used. box of size (largest dimension)+2*padding is used.
4. You can specify the total number of molecules (both waters and ions) to add. A cubic box is then created whose size is 4. You can specify the total number of molecules (both waters and ions) to add. A cubic box is then created whose size is
just large enough to hold the specified amount of solvent. just large enough hold the specified amount of solvent.
5. Finally, if none of the above options is specified, the existing Topology's box vectors are used. 5. Finally, if none of the above options is specified, the existing Topology's box vectors are used.
Parameters: Parameters
- forcefield (ForceField) the ForceField to use for determining van der Waals radii and atomic charges ----------
- model (string='tip3p') the water model to use. Supported values are 'tip3p', 'spce', 'tip4pew', and 'tip5p'. forcefield : ForceField
- boxSize (Vec3=None) the size of the box to fill with water the ForceField to use for determining van der Waals radii and atomic charges
- boxVectors (tuple of Vec3=None) the vectors defining the periodic box to fill with water model : str='tip3p'
- padding (distance=None) the padding distance to use the water model to use. Supported values are 'tip3p', 'spce', 'tip4pew', and 'tip5p'.
- numAdded (int=None) the total number of molecules (waters and ions) to add boxSize : Vec3=None
- positiveIon (string='Na+') the type of positive ion to add. Allowed values are 'Cs+', 'K+', 'Li+', 'Na+', and 'Rb+' the size of the box to fill with water
- negativeIon (string='Cl-') the type of negative ion to add. Allowed values are 'Cl-', 'Br-', 'F-', and 'I-'. Be aware boxVectors : tuple of Vec3=None
that not all force fields support all ion types. the vectors defining the periodic box to fill with water
- ionicStrength (concentration=0*molar) the total concentration of ions (both positive and negative) to add. This padding : distance=None
does not include ions that are added to neutralize the system. the padding distance to use
- neutralize (bool=True) whether to add ions to neutralize the system numAdded : int=None
the total number of molecules (waters and ions) to add
positiveIon : string='Na+'
the type of positive ion to add. Allowed values are 'Cs+', 'K+', 'Li+', 'Na+', and 'Rb+'
negativeIon : string='Cl-'
the type of negative ion to add. Allowed values are 'Cl-', 'Br-', 'F-', and 'I-'. Be aware
that not all force fields support all ion types.
ionicStrength : concentration=0*molar
the total concentration of ions (both positive and negative) to add. This
does not include ions that are added to neutralize the system.
neutralize : bool=True
whether to add ions to neutralize the system
""" """
if len([x for x in (boxSize, boxVectors, padding, numAdded) if x is not None]) > 1: if len([x for x in (boxSize, boxVectors, padding, numAdded) if x is not None]) > 1:
raise ValueError('At most one of the following arguments may be specified: boxSize, boxVectors, padding, numAdded') raise ValueError('At most one of the following arguments may be specified: boxSize, boxVectors, padding, numAdded')
...@@ -295,13 +322,13 @@ class Modeller(object): ...@@ -295,13 +322,13 @@ class Modeller(object):
pdbTopology = pdb.getTopology() pdbTopology = pdb.getTopology()
pdbPositions = pdb.getPositions().value_in_unit(nanometer) pdbPositions = pdb.getPositions().value_in_unit(nanometer)
pdbResidues = list(pdbTopology.residues()) pdbResidues = list(pdbTopology.residues())
pdbBoxSize = pdbTopology.getUnitCellDimensions().value_in_unit(nanometer) pdbBoxSize = pdbTopology.getUnitCellDimensions().value_in_unit(nanometer)
# Pick a unit cell size. # Pick a unit cell size.
if numAdded is not None: if numAdded is not None:
# Select a padding distance which is guaranteed to give more than the specified number of molecules. # Select a padding distance which is guaranteed to give more than the specified number of molecules.
padding = 1.1*(numAdded/((len(pdbResidues)/pdbBoxSize[0]**3)*8))**(1.0/3.0) padding = 1.1*(numAdded/((len(pdbResidues)/pdbBoxSize[0]**3)*8))**(1.0/3.0)
if padding < 0.5: if padding < 0.5:
padding = 0.5 # Ensure we have enough when adding very small numbers of molecules padding = 0.5 # Ensure we have enough when adding very small numbers of molecules
...@@ -443,20 +470,20 @@ class Modeller(object): ...@@ -443,20 +470,20 @@ class Modeller(object):
if numAdded is not None: if numAdded is not None:
# We added many more waters than we actually want. Sort them based on distance to the nearest box edge and # We added many more waters than we actually want. Sort them based on distance to the nearest box edge and
# only keep the ones in the middle. # only keep the ones in the middle.
lowerBound = center-box/2 lowerBound = center-box/2
upperBound = center+box/2 upperBound = center+box/2
distToEdge = (min(min(pos-lowerBound), min(upperBound-pos)) for index, pos in addedWaters) distToEdge = (min(min(pos-lowerBound), min(upperBound-pos)) for index, pos in addedWaters)
sortedIndex = [i[0] for i in sorted(enumerate(distToEdge), key=lambda x: -x[1])] sortedIndex = [i[0] for i in sorted(enumerate(distToEdge), key=lambda x: -x[1])]
addedWaters = [addedWaters[i] for i in sortedIndex[:numAdded]] addedWaters = [addedWaters[i] for i in sortedIndex[:numAdded]]
# Compute a new periodic box size. # Compute a new periodic box size.
maxSize = max(max((pos[i] for index, pos in addedWaters))-min((pos[i] for index, pos in addedWaters)) for i in range(3)) maxSize = max(max((pos[i] for index, pos in addedWaters))-min((pos[i] for index, pos in addedWaters)) for i in range(3))
newTopology.setUnitCellDimensions(Vec3(maxSize, maxSize, maxSize)) newTopology.setUnitCellDimensions(Vec3(maxSize, maxSize, maxSize))
else: else:
# There could be clashes between water molecules at the box edges. Find ones to remove. # There could be clashes between water molecules at the box edges. Find ones to remove.
upperCutoff = center+box/2-Vec3(waterCutoff, waterCutoff, waterCutoff) upperCutoff = center+box/2-Vec3(waterCutoff, waterCutoff, waterCutoff)
lowerCutoff = center-box/2+Vec3(waterCutoff, waterCutoff, waterCutoff) lowerCutoff = center-box/2+Vec3(waterCutoff, waterCutoff, waterCutoff)
lowerSkinPositions = [pos for index, pos in addedWaters if pos[0] < lowerCutoff[0] or pos[1] < lowerCutoff[1] or pos[2] < lowerCutoff[2]] lowerSkinPositions = [pos for index, pos in addedWaters if pos[0] < lowerCutoff[0] or pos[1] < lowerCutoff[1] or pos[2] < lowerCutoff[2]]
...@@ -606,16 +633,29 @@ class Modeller(object): ...@@ -606,16 +633,29 @@ class Modeller(object):
Definitions for standard amino acids and nucleotides are built in. You can call loadHydrogenDefinitions() to load Definitions for standard amino acids and nucleotides are built in. You can call loadHydrogenDefinitions() to load
additional definitions for other residue types. additional definitions for other residue types.
Parameters: Parameters
- forcefield (ForceField=None) the ForceField to use for determining the positions of hydrogens. If this is None, ----------
positions will be picked which are generally reasonable but not optimized for any particular ForceField. forcefield : ForceField=None
- pH (float=7.0) the pH based on which to select variants the ForceField to use for determining the positions of hydrogens.
- variants (list=None) an optional list of variants to use. If this is specified, its length must equal the number If this is None, positions will be picked which are generally
of residues in the model. variants[i] is the name of the variant to use for residue i (indexed starting at 0). reasonable but not optimized for any particular ForceField.
If an element is None, the standard rules will be followed to select a variant for that residue. pH : float=7.0
- platform (Platform=None) the Platform to use when computing the hydrogen atom positions. If this is None, the pH based on which to select variants
the default Platform will be used. variants : list=None
Returns: a list of what variant was actually selected for each residue, in the same format as the variants parameter an optional list of variants to use. If this is specified, its
length must equal the number of residues in the model. variants[i]
is the name of the variant to use for residue i (indexed starting at
0). If an element is None, the standard rules will be followed to
select a variant for that residue.
platform : Platform=None
the Platform to use when computing the hydrogen atom positions. If
this is None, the default Platform will be used.
Returns
-------
list
a list of what variant was actually selected for each residue,
in the same format as the variants parameter
""" """
# Check the list of variants. # Check the list of variants.
...@@ -804,7 +844,7 @@ class Modeller(object): ...@@ -804,7 +844,7 @@ class Modeller(object):
if forcefield is not None: if forcefield is not None:
# Use the ForceField the user specified. # Use the ForceField the user specified.
system = forcefield.createSystem(newTopology, rigidWater=False) system = forcefield.createSystem(newTopology, rigidWater=False)
atoms = list(newTopology.atoms()) atoms = list(newTopology.atoms())
for i in range(system.getNumParticles()): for i in range(system.getNumParticles()):
...@@ -814,7 +854,7 @@ class Modeller(object): ...@@ -814,7 +854,7 @@ class Modeller(object):
else: else:
# Create a System that restrains the distance of each hydrogen from its parent atom # Create a System that restrains the distance of each hydrogen from its parent atom
# and causes hydrogens to spread out evenly. # and causes hydrogens to spread out evenly.
system = System() system = System()
nonbonded = CustomNonbondedForce('100/((r/0.1)^4+1)') nonbonded = CustomNonbondedForce('100/((r/0.1)^4+1)')
bonds = HarmonicBondForce() bonds = HarmonicBondForce()
...@@ -838,7 +878,7 @@ class Modeller(object): ...@@ -838,7 +878,7 @@ class Modeller(object):
for residue in newTopology.residues(): for residue in newTopology.residues():
if residue.name == 'HOH': if residue.name == 'HOH':
# Add an angle term to make the water geometry correct. # Add an angle term to make the water geometry correct.
atoms = list(residue.atoms()) atoms = list(residue.atoms())
oindex = [i for i in range(len(atoms)) if atoms[i].element == elem.oxygen] oindex = [i for i in range(len(atoms)) if atoms[i].element == elem.oxygen]
if len(atoms) == 3 and len(oindex) == 1: if len(atoms) == 3 and len(oindex) == 1:
...@@ -846,12 +886,12 @@ class Modeller(object): ...@@ -846,12 +886,12 @@ class Modeller(object):
angles.addAngle(atoms[hindex[0]].index, atoms[oindex[0]].index, atoms[hindex[1]].index, 1.824, 836.8) angles.addAngle(atoms[hindex[0]].index, atoms[oindex[0]].index, atoms[hindex[1]].index, 1.824, 836.8)
else: else:
# Add angle terms for any hydroxyls. # Add angle terms for any hydroxyls.
for atom in residue.atoms(): for atom in residue.atoms():
index = atom.index index = atom.index
if atom.element == elem.oxygen and len(bondedTo[index]) == 2 and elem.hydrogen in (a.element for a in bondedTo[index]): if atom.element == elem.oxygen and len(bondedTo[index]) == 2 and elem.hydrogen in (a.element for a in bondedTo[index]):
angles.addAngle(bondedTo[index][0].index, index, bondedTo[index][1].index, 1.894, 460.24) angles.addAngle(bondedTo[index][0].index, index, bondedTo[index][1].index, 1.894, 460.24)
if platform is None: if platform is None:
context = Context(system, VerletIntegrator(0.0)) context = Context(system, VerletIntegrator(0.0))
else: else:
...@@ -864,18 +904,24 @@ class Modeller(object): ...@@ -864,18 +904,24 @@ class Modeller(object):
return actualVariants return actualVariants
def addExtraParticles(self, forcefield): def addExtraParticles(self, forcefield):
"""Add missing extra particles to the model that are required by a force field. """Add missing extra particles to the model that are required by a force
field.
Some force fields use "extra particles" that do not represent actual atoms, but still need to be included in
the System. Examples include lone pairs, Drude particles, and the virtual sites used in some water models Some force fields use "extra particles" that do not represent
to adjust the charge distribution. Extra particles can be recognized by the fact that their element is None. actual atoms, but still need to be included in the System. Examples
include lone pairs, Drude particles, and the virtual sites used in some
This method is primarily used to add extra particles, but it can also remove them. It tries to match every water models to adjust the charge distribution. Extra particles can be
residue in the Topology to a template in the force field. If there is no match, it will both add and remove recognized by the fact that their element is None.
extra particles as necessary to make it match.
This method is primarily used to add extra particles, but it can also
Parameters: remove them. It tries to match every residue in the Topology to a
- forcefield (ForceField) the ForceField defining what extra particles should be present template in the force field. If there is no match, it will both add
and remove extra particles as necessary to make it match.
Parameters
----------
forcefield : ForceField
the ForceField defining what extra particles should be present
""" """
# Create copies of all residue templates that have had all extra points removed. # Create copies of all residue templates that have had all extra points removed.
......
This diff is collapsed.
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment