Skip to content

GitLab

"wrappers/python/vscode:/vscode.git/clone" did not exist on "311cc01ee59eed662cdcda793a56cb58df6e9ca4"
Commit 74415dd9 authored by Lee-Ping Wang's avatar Lee-Ping Wang
Browse files

Merge branch 'master' of https://github.com/SimTk/openmm

parents 25c22045 6701dacc
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 823 additions and 823 deletions
wrappers/python/simtk/openmm/__init__.py
View file @ 74415dd9
...@@ -35,4 +35,4 @@ else: ...@@ -35,4 +35,4 @@ else:
from simtk.openmm.openmm import * from simtk.openmm.openmm import *
from simtk.openmm.vec3 import Vec3 from simtk.openmm.vec3 import Vec3
pluginLoadedLibNames = Platform.loadPluginsFromDirectory(Platform.getDefaultPluginsDirectory()) pluginLoadedLibNames = Platform.loadPluginsFromDirectory(Platform.getDefaultPluginsDirectory())
\ No newline at end of file
wrappers/python/simtk/openmm/app/amberinpcrdfile.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -40,15 +40,15 @@ except: ...@@ -40,15 +40,15 @@ except:
class AmberInpcrdFile(object): class AmberInpcrdFile(object):
"""AmberInpcrdFile parses an AMBER inpcrd file and loads the data stored in it.""" """AmberInpcrdFile parses an AMBER inpcrd file and loads the data stored in it."""
def __init__(self, file, loadVelocities=False, loadBoxVectors=False): def __init__(self, file, loadVelocities=False, loadBoxVectors=False):
"""Load an inpcrd file. """Load an inpcrd file.
An inpcrd file contains atom positions and, optionally, velocities and periodic box dimensions. An inpcrd file contains atom positions and, optionally, velocities and periodic box dimensions.
Unfortunately, it is sometimes impossible to determine from the file itself exactly what data Unfortunately, it is sometimes impossible to determine from the file itself exactly what data
it contains. You therefore must specify in advance what data to load. It is stored into this it contains. You therefore must specify in advance what data to load. It is stored into this
object's "positions", "velocities", and "boxVectors" fields. object's "positions", "velocities", and "boxVectors" fields.
Parameters: Parameters:
- file (string) the name of the file to load - file (string) the name of the file to load
- loadVelocities (boolean=False) whether to load velocities from the file - loadVelocities (boolean=False) whether to load velocities from the file
...@@ -78,7 +78,7 @@ class AmberInpcrdFile(object): ...@@ -78,7 +78,7 @@ class AmberInpcrdFile(object):
def getPositions(self, asNumpy=False): def getPositions(self, asNumpy=False):
"""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 - asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s
""" """
...@@ -87,10 +87,10 @@ class AmberInpcrdFile(object): ...@@ -87,10 +87,10 @@ class AmberInpcrdFile(object):
self._numpyPositions = Quantity(numpy.array(self.positions.value_in_unit(nanometers)), nanometers) self._numpyPositions = Quantity(numpy.array(self.positions.value_in_unit(nanometers)), nanometers)
return self._numpyPositions return self._numpyPositions
return self.positions return self.positions
def getVelocities(self, asNumpy=False): def getVelocities(self, asNumpy=False):
"""Get the atomic velocities. """Get the atomic velocities.
Parameters: Parameters:
- asNumpy (boolean=False) if true, the vectors are returned as numpy arrays instead of Vec3s - asNumpy (boolean=False) if true, the vectors are returned as numpy arrays instead of Vec3s
""" """
...@@ -99,10 +99,10 @@ class AmberInpcrdFile(object): ...@@ -99,10 +99,10 @@ class AmberInpcrdFile(object):
self._numpyVelocities = Quantity(numpy.array(self.velocities.value_in_unit(nanometers/picoseconds)), nanometers/picoseconds) self._numpyVelocities = Quantity(numpy.array(self.velocities.value_in_unit(nanometers/picoseconds)), nanometers/picoseconds)
return self._numpyVelocities return self._numpyVelocities
return self.velocities return self.velocities
def getBoxVectors(self, asNumpy=False): def getBoxVectors(self, asNumpy=False):
"""Get the periodic box vectors. """Get the periodic box vectors.
Parameters: Parameters:
- asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s - asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s
""" """
...@@ -114,4 +114,4 @@ class AmberInpcrdFile(object): ...@@ -114,4 +114,4 @@ class AmberInpcrdFile(object):
self._numpyBoxVectors.append(Quantity(numpy.array(self.boxVectors[2].value_in_unit(nanometers)), nanometers)) self._numpyBoxVectors.append(Quantity(numpy.array(self.boxVectors[2].value_in_unit(nanometers)), nanometers))
return self._numpyBoxVectors return self._numpyBoxVectors
return self.boxVectors return self.boxVectors
wrappers/python/simtk/openmm/app/amberprmtopfile.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -48,15 +48,15 @@ GBn = object() ...@@ -48,15 +48,15 @@ GBn = object()
class AmberPrmtopFile(object): class AmberPrmtopFile(object):
"""AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it.""" """AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it."""
def __init__(self, file): def __init__(self, file):
"""Load a prmtop file.""" """Load a prmtop file."""
top = Topology() top = Topology()
## The Topology read from the prmtop file ## The Topology read from the prmtop file
self.topology = top self.topology = top
# Load the prmtop file # Load the prmtop file
prmtop = amber_file_parser.PrmtopLoader(file) prmtop = amber_file_parser.PrmtopLoader(file)
self._prmtop = prmtop self._prmtop = prmtop
...@@ -82,7 +82,7 @@ class AmberPrmtopFile(object): ...@@ -82,7 +82,7 @@ class AmberPrmtopFile(object):
atomName = atomReplacements[atomName] atomName = atomReplacements[atomName]
# Try to guess the element. # Try to guess the element.
upper = atomName.upper() upper = atomName.upper()
if upper.startswith('CL'): if upper.startswith('CL'):
element = elem.chlorine element = elem.chlorine
...@@ -96,17 +96,17 @@ class AmberPrmtopFile(object): ...@@ -96,17 +96,17 @@ class AmberPrmtopFile(object):
except KeyError: except KeyError:
element = None element = None
top.addAtom(atomName, element, r) top.addAtom(atomName, element, r)
# Add bonds to the topology # Add bonds to the topology
atoms = list(top.atoms()) atoms = list(top.atoms())
for bond in prmtop.getBondsWithH(): for bond in prmtop.getBondsWithH():
top.addBond(atoms[bond[0]], atoms[bond[1]]) top.addBond(atoms[bond[0]], atoms[bond[1]])
for bond in prmtop.getBondsNoH(): for bond in prmtop.getBondsNoH():
top.addBond(atoms[bond[0]], atoms[bond[1]]) top.addBond(atoms[bond[0]], atoms[bond[1]])
# Set the periodic box size. # Set the periodic box size.
if prmtop.getIfBox(): if prmtop.getIfBox():
top.setUnitCellDimensions(tuple(x.value_in_unit(unit.nanometer) for x in prmtop.getBoxBetaAndDimensions()[1:4])*unit.nanometer) top.setUnitCellDimensions(tuple(x.value_in_unit(unit.nanometer) for x in prmtop.getBoxBetaAndDimensions()[1:4])*unit.nanometer)
...@@ -114,7 +114,7 @@ class AmberPrmtopFile(object): ...@@ -114,7 +114,7 @@ class AmberPrmtopFile(object):
constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, removeCMMotion=True, constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, removeCMMotion=True,
ewaldErrorTolerance=0.0005): ewaldErrorTolerance=0.0005):
"""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: Parameters:
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are - nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
...@@ -168,4 +168,4 @@ class AmberPrmtopFile(object): ...@@ -168,4 +168,4 @@ class AmberPrmtopFile(object):
force.setEwaldErrorTolerance(ewaldErrorTolerance) force.setEwaldErrorTolerance(ewaldErrorTolerance)
if removeCMMotion: if removeCMMotion:
sys.addForce(mm.CMMotionRemover()) sys.addForce(mm.CMMotionRemover())
return sys return sys
\ No newline at end of file
wrappers/python/simtk/openmm/app/dcdfile.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -40,18 +40,18 @@ from simtk.unit import picoseconds, nanometers, angstroms, is_quantity, norm ...@@ -40,18 +40,18 @@ from simtk.unit import picoseconds, nanometers, angstroms, is_quantity, norm
class DCDFile(object): class DCDFile(object):
"""DCDFile provides methods for creating DCD files. """DCDFile provides methods for creating DCD files.
DCD is a file format for storing simulation trajectories. It is supported by many programs, such DCD is a file format for storing simulation trajectories. It is supported by many programs, such
as CHARMM, NAMD, and X-PLOR. Note, however, that different programs produce subtly different as CHARMM, NAMD, and X-PLOR. Note, however, that different programs produce subtly different
versions of the format. This class generates the CHARMM version. Also note that there is no versions of the format. This class generates the CHARMM version. Also note that there is no
standard byte ordering (big-endian or little-endian) for this format. This class always generates standard byte ordering (big-endian or little-endian) for this format. This class always generates
files with little-endian ordering. files with little-endian ordering.
To use this class, create a DCDFile object, then call writeModel() once for each model in the file.""" To use this class, create a DCDFile object, then call writeModel() once for each model in the file."""
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 - file (file) A file to write to
- topology (Topology) The Topology defining the molecular system being written - topology (Topology) The Topology defining the molecular system being written
...@@ -76,10 +76,10 @@ class DCDFile(object): ...@@ -76,10 +76,10 @@ class DCDFile(object):
header += struct.pack('<80s', 'Created '+time.asctime(time.localtime(time.time()))) header += struct.pack('<80s', 'Created '+time.asctime(time.localtime(time.time())))
header += struct.pack('<4i', 164, 4, len(list(topology.atoms())), 4) header += struct.pack('<4i', 164, 4, len(list(topology.atoms())), 4)
file.write(header) file.write(header)
def writeModel(self, positions, unitCellDimensions=None): def writeModel(self, positions, unitCellDimensions=None):
"""Write out a model to the DCD file. """Write out a model to the DCD file.
Parameters: Parameters:
- positions (list) The list of atomic positions to write - positions (list) The list of atomic positions to write
- unitCellDimensions (Vec3=None) The dimensions of the crystallographic unit cell. If None, the dimensions specified in - unitCellDimensions (Vec3=None) The dimensions of the crystallographic unit cell. If None, the dimensions specified in
...@@ -87,7 +87,7 @@ class DCDFile(object): ...@@ -87,7 +87,7 @@ class DCDFile(object):
represent a periodic system. represent a periodic system.
""" """
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')
if is_quantity(positions): if is_quantity(positions):
positions = positions.value_in_unit(nanometers) positions = positions.value_in_unit(nanometers)
if any(math.isnan(norm(pos)) for pos in positions): if any(math.isnan(norm(pos)) for pos in positions):
...@@ -95,17 +95,17 @@ class DCDFile(object): ...@@ -95,17 +95,17 @@ class DCDFile(object):
if any(math.isinf(norm(pos)) for pos in positions): if any(math.isinf(norm(pos)) for pos in positions):
raise ValueError('Particle position is infinite') raise ValueError('Particle position is infinite')
file = self._file file = self._file
# Update the header. # Update the header.
self._modelCount += 1 self._modelCount += 1
file.seek(8, os.SEEK_SET) file.seek(8, os.SEEK_SET)
file.write(struct.pack('<i', self._modelCount)) file.write(struct.pack('<i', self._modelCount))
file.seek(20, os.SEEK_SET) file.seek(20, os.SEEK_SET)
file.write(struct.pack('<i', self._firstStep+self._modelCount*self._interval)) file.write(struct.pack('<i', self._firstStep+self._modelCount*self._interval))
# Write the data. # Write the data.
file.seek(0, os.SEEK_END) file.seek(0, os.SEEK_END)
boxSize = self._topology.getUnitCellDimensions() boxSize = self._topology.getUnitCellDimensions()
if boxSize is not None: if boxSize is not None:
......
wrappers/python/simtk/openmm/app/dcdreporter.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -34,16 +34,16 @@ __version__ = "1.0" ...@@ -34,16 +34,16 @@ __version__ = "1.0"
import simtk.openmm as mm import simtk.openmm as mm
from simtk.openmm.app import DCDFile from simtk.openmm.app import DCDFile
from simtk.unit import nanometer from simtk.unit import nanometer
class DCDReporter(object): class DCDReporter(object):
"""DCDReporter outputs a series of frames from a Simulation to a DCD file. """DCDReporter outputs a series of frames from a Simulation to a DCD file.
To use it, create a DCDReporter, then add it to the Simulation's list of reporters. To use it, create a DCDReporter, then add it to the Simulation's list of reporters.
""" """
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 - file (string) The file to write to
- reportInterval (int) The interval (in time steps) at which to write frames - reportInterval (int) The interval (in time steps) at which to write frames
...@@ -51,10 +51,10 @@ class DCDReporter(object): ...@@ -51,10 +51,10 @@ class DCDReporter(object):
self._reportInterval = reportInterval self._reportInterval = reportInterval
self._out = open(file, 'wb') self._out = open(file, 'wb')
self._dcd = None self._dcd = None
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) The Simulation to generate a report for
Returns: A five element tuple. The first element is the number of steps until the Returns: A five element tuple. The first element is the number of steps until the
...@@ -63,10 +63,10 @@ class DCDReporter(object): ...@@ -63,10 +63,10 @@ class DCDReporter(object):
""" """
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)
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 - simulation (Simulation) The Simulation to generate a report for
- state (State) The current state of the simulation - state (State) The current state of the simulation
...@@ -75,6 +75,6 @@ class DCDReporter(object): ...@@ -75,6 +75,6 @@ class DCDReporter(object):
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()
self._dcd.writeModel(state.getPositions(), mm.Vec3(a[0].value_in_unit(nanometer), b[1].value_in_unit(nanometer), c[2].value_in_unit(nanometer))*nanometer) self._dcd.writeModel(state.getPositions(), mm.Vec3(a[0].value_in_unit(nanometer), b[1].value_in_unit(nanometer), c[2].value_in_unit(nanometer))*nanometer)
def __del__(self): def __del__(self):
self._out.close() self._out.close()
wrappers/python/simtk/openmm/app/element.py
View file @ 74415dd9
...@@ -12,7 +12,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -12,7 +12,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Christopher M. Bruns Authors: Christopher M. Bruns
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -38,13 +38,13 @@ from simtk.unit import daltons ...@@ -38,13 +38,13 @@ from simtk.unit import daltons
class Element: class Element:
"""An Element represents a chemical element. """An Element represents a chemical element.
The simtk.openmm.app.element module contains objects for all the standard chemical elements, The simtk.openmm.app.element module contains objects for all the standard chemical elements,
such as element.hydrogen or element.carbon. You can also call the static method Element.getBySymbol() to such as element.hydrogen or element.carbon. You can also call the static method Element.getBySymbol() to
look up the Element with a particular chemical symbol.""" look up the Element with a particular chemical symbol."""
_elements_by_symbol = {} _elements_by_symbol = {}
def __init__(self, number, name, symbol, mass): def __init__(self, number, name, symbol, mass):
## The atomic number of the element ## The atomic number of the element
self.atomic_number = number self.atomic_number = number
......
wrappers/python/simtk/openmm/app/gromacsgrofile.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Lee-Ping Wang, Peter Eastman Authors: Lee-Ping Wang, Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -47,7 +47,7 @@ def _isint(word): ...@@ -47,7 +47,7 @@ def _isint(word):
@param[in] word String (for instance, '123', '153.0', '2.', '-354') @param[in] word String (for instance, '123', '153.0', '2.', '-354')
@return answer Boolean which specifies whether the string is an integer (only +/- sign followed by digits) @return answer Boolean which specifies whether the string is an integer (only +/- sign followed by digits)
""" """
return match('^[-+]?[0-9]+$',word) return match('^[-+]?[0-9]+$',word)
...@@ -57,7 +57,7 @@ def _isfloat(word): ...@@ -57,7 +57,7 @@ def _isfloat(word):
@param[in] word String (for instance, '123', '153.0', '2.', '-354') @param[in] word String (for instance, '123', '153.0', '2.', '-354')
@return answer Boolean which specifies whether the string is any number @return answer Boolean which specifies whether the string is any number
""" """
return match('^[-+]?[0-9]*\.?[0-9]*([eEdD][-+]?[0-9]+)?$',word) return match('^[-+]?[0-9]*\.?[0-9]*([eEdD][-+]?[0-9]+)?$',word)
...@@ -65,7 +65,7 @@ def _is_gro_coord(line): ...@@ -65,7 +65,7 @@ def _is_gro_coord(line):
""" Determines whether a line contains GROMACS data or not """ Determines whether a line contains GROMACS data or not
@param[in] line The line to be tested @param[in] line The line to be tested
""" """
sline = line.split() sline = line.split()
if len(sline) == 6 or len(sline) == 9: if len(sline) == 6 or len(sline) == 9:
...@@ -79,7 +79,7 @@ def _is_gro_box(line): ...@@ -79,7 +79,7 @@ def _is_gro_box(line):
""" Determines whether a line contains a GROMACS box vector or not """ Determines whether a line contains a GROMACS box vector or not
@param[in] line The line to be tested @param[in] line The line to be tested
""" """
sline = line.split() sline = line.split()
if len(sline) == 9 and all([_isfloat(i) for i in sline]): if len(sline) == 9 and all([_isfloat(i) for i in sline]):
...@@ -91,16 +91,16 @@ def _is_gro_box(line): ...@@ -91,16 +91,16 @@ def _is_gro_box(line):
class GromacsGroFile(object): class GromacsGroFile(object):
"""GromacsGroFile parses a Gromacs .gro file and constructs a set of atom positions from it. """GromacsGroFile parses a Gromacs .gro file and constructs a set of atom positions from it.
A .gro file also contains some topological information, such as elements and residue names, A .gro file also contains some topological information, such as elements and residue names,
but not enough to construct a full Topology object. This information is recorded and stored but not enough to construct a full Topology object. This information is recorded and stored
in the object's public fields.""" in the object's public fields."""
def __init__(self, file): def __init__(self, file):
"""Load a .gro file. """Load a .gro file.
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
""" """
...@@ -145,7 +145,7 @@ class GromacsGroFile(object): ...@@ -145,7 +145,7 @@ class GromacsGroFile(object):
else: else:
raise Exception("Unexpected line in .gro file: "+line) raise Exception("Unexpected line in .gro file: "+line)
ln += 1 ln += 1
## The atom positions read from the file. If the file contains multiple frames, these are the positions in the first frame. ## The atom positions read from the file. If the file contains multiple frames, these are the positions in the first frame.
self.positions = xyzs[0] self.positions = xyzs[0]
## A list containing the element of each atom stored in the file ## A list containing the element of each atom stored in the file
...@@ -159,14 +159,14 @@ class GromacsGroFile(object): ...@@ -159,14 +159,14 @@ class GromacsGroFile(object):
self._positions = xyzs self._positions = xyzs
self._unitCellDimensions = boxes self._unitCellDimensions = boxes
self._numpyPositions = None self._numpyPositions = None
def getNumFrames(self): def getNumFrames(self):
"""Get the number of frames stored in the file.""" """Get the number of frames stored in the file."""
return len(self._positions) return len(self._positions)
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 - 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 - frame (int=0) the index of the frame for which to get positions
...@@ -178,10 +178,10 @@ class GromacsGroFile(object): ...@@ -178,10 +178,10 @@ class GromacsGroFile(object):
self._numpyPositions[frame] = Quantity(numpy.array(self._positions[frame].value_in_unit(nanometers)), nanometers) self._numpyPositions[frame] = Quantity(numpy.array(self._positions[frame].value_in_unit(nanometers)), nanometers)
return self._numpyPositions[frame] return self._numpyPositions[frame]
return self._positions[frame] return self._positions[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
""" """
......
wrappers/python/simtk/openmm/app/gromacstopfile.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -49,7 +49,7 @@ OBC2 = prmtop.OBC2 ...@@ -49,7 +49,7 @@ OBC2 = prmtop.OBC2
class GromacsTopFile(object): class GromacsTopFile(object):
"""GromacsTopFile parses a Gromacs top file and constructs a Topology and (optionally) an OpenMM System from it.""" """GromacsTopFile parses a Gromacs top file and constructs a Topology and (optionally) an OpenMM System from it."""
class _MoleculeType(object): class _MoleculeType(object):
"""Inner class to store information about a molecule type.""" """Inner class to store information about a molecule type."""
def __init__(self): def __init__(self):
...@@ -60,7 +60,7 @@ class GromacsTopFile(object): ...@@ -60,7 +60,7 @@ class GromacsTopFile(object):
self.exclusions = [] self.exclusions = []
self.pairs = [] self.pairs = []
self.cmaps = [] self.cmaps = []
def _processFile(self, file): def _processFile(self, file):
append = '' append = ''
for line in open(file): for line in open(file):
...@@ -69,7 +69,7 @@ class GromacsTopFile(object): ...@@ -69,7 +69,7 @@ class GromacsTopFile(object):
else: else:
self._processLine(append+' '+line, file) self._processLine(append+' '+line, file)
append = '' append = ''
def _processLine(self, line, file): def _processLine(self, line, file):
"""Process one line from a file.""" """Process one line from a file."""
if ';' in line: if ';' in line:
...@@ -79,13 +79,13 @@ class GromacsTopFile(object): ...@@ -79,13 +79,13 @@ class GromacsTopFile(object):
if stripped.startswith('*') or len(stripped) == 0: if stripped.startswith('*') or len(stripped) == 0:
# A comment or empty line. # A comment or empty line.
return return
elif stripped.startswith('[') and not ignore: elif stripped.startswith('[') and not ignore:
# The start of a category. # The start of a category.
if not stripped.endswith(']'): if not stripped.endswith(']'):
raise ValueError('Illegal line in .top file: '+line) raise ValueError('Illegal line in .top file: '+line)
self._currentCategory = stripped[1:-1].strip() self._currentCategory = stripped[1:-1].strip()
elif stripped.startswith('#'): elif stripped.startswith('#'):
# A preprocessor command. # A preprocessor command.
fields = stripped.split() fields = stripped.split()
...@@ -127,7 +127,7 @@ class GromacsTopFile(object): ...@@ -127,7 +127,7 @@ class GromacsTopFile(object):
if len(self._ifStack) == 0: if len(self._ifStack) == 0:
raise ValueError('Unexpected line in .top file: '+line) raise ValueError('Unexpected line in .top file: '+line)
del(self._ifStack[-1]) del(self._ifStack[-1])
elif not ignore: elif not ignore:
# A line of data for the current category # A line of data for the current category
if self._currentCategory is None: if self._currentCategory is None:
...@@ -166,7 +166,7 @@ class GromacsTopFile(object): ...@@ -166,7 +166,7 @@ class GromacsTopFile(object):
self._processPairType(line) self._processPairType(line)
elif self._currentCategory == 'cmaptypes': elif self._currentCategory == 'cmaptypes':
self._processCmapType(line) self._processCmapType(line)
def _processDefaults(self, line): def _processDefaults(self, line):
"""Process the [ defaults ] line.""" """Process the [ defaults ] line."""
fields = line.split() fields = line.split()
...@@ -179,7 +179,7 @@ class GromacsTopFile(object): ...@@ -179,7 +179,7 @@ class GromacsTopFile(object):
if fields[2].lower() == 'no': if fields[2].lower() == 'no':
raise ValueError('gen_pairs=no is not supported') raise ValueError('gen_pairs=no is not supported')
self._defaults = fields self._defaults = fields
def _processMoleculeType(self, line): def _processMoleculeType(self, line):
"""Process a line in the [ moleculetypes ] category.""" """Process a line in the [ moleculetypes ] category."""
fields = line.split() fields = line.split()
...@@ -188,14 +188,14 @@ class GromacsTopFile(object): ...@@ -188,14 +188,14 @@ class GromacsTopFile(object):
type = GromacsTopFile._MoleculeType() type = GromacsTopFile._MoleculeType()
self._moleculeTypes[fields[0]] = type self._moleculeTypes[fields[0]] = type
self._currentMoleculeType = type self._currentMoleculeType = type
def _processMolecule(self, line): def _processMolecule(self, line):
"""Process a line in the [ molecules ] category.""" """Process a line in the [ molecules ] category."""
fields = line.split() fields = line.split()
if len(fields) < 2: if len(fields) < 2:
raise ValueError('Too few fields in [ molecules ] line: '+line); raise ValueError('Too few fields in [ molecules ] line: '+line);
self._molecules.append((fields[0], int(fields[1]))) self._molecules.append((fields[0], int(fields[1])))
def _processAtom(self, line): def _processAtom(self, line):
"""Process a line in the [ atoms ] category.""" """Process a line in the [ atoms ] category."""
if self._currentMoleculeType is None: if self._currentMoleculeType is None:
...@@ -204,7 +204,7 @@ class GromacsTopFile(object): ...@@ -204,7 +204,7 @@ class GromacsTopFile(object):
if len(fields) < 5: if len(fields) < 5:
raise ValueError('Too few fields in [ atoms ] line: '+line); raise ValueError('Too few fields in [ atoms ] line: '+line);
self._currentMoleculeType.atoms.append(fields) self._currentMoleculeType.atoms.append(fields)
def _processBond(self, line): def _processBond(self, line):
"""Process a line in the [ bonds ] category.""" """Process a line in the [ bonds ] category."""
if self._currentMoleculeType is None: if self._currentMoleculeType is None:
...@@ -215,7 +215,7 @@ class GromacsTopFile(object): ...@@ -215,7 +215,7 @@ class GromacsTopFile(object):
if fields[2] != '1': if fields[2] != '1':
raise ValueError('Unsupported function type in [ bonds ] line: '+line); raise ValueError('Unsupported function type in [ bonds ] line: '+line);
self._currentMoleculeType.bonds.append(fields) self._currentMoleculeType.bonds.append(fields)
def _processAngle(self, line): def _processAngle(self, line):
"""Process a line in the [ angles ] category.""" """Process a line in the [ angles ] category."""
if self._currentMoleculeType is None: if self._currentMoleculeType is None:
...@@ -226,7 +226,7 @@ class GromacsTopFile(object): ...@@ -226,7 +226,7 @@ class GromacsTopFile(object):
if fields[3] not in ('1', '5'): if fields[3] not in ('1', '5'):
raise ValueError('Unsupported function type in [ angles ] line: '+line); raise ValueError('Unsupported function type in [ angles ] line: '+line);
self._currentMoleculeType.angles.append(fields) self._currentMoleculeType.angles.append(fields)
def _processDihedral(self, line): def _processDihedral(self, line):
"""Process a line in the [ dihedrals ] category.""" """Process a line in the [ dihedrals ] category."""
if self._currentMoleculeType is None: if self._currentMoleculeType is None:
...@@ -237,7 +237,7 @@ class GromacsTopFile(object): ...@@ -237,7 +237,7 @@ class GromacsTopFile(object):
if fields[4] not in ('1', '2', '3', '4', '9'): if fields[4] not in ('1', '2', '3', '4', '9'):
raise ValueError('Unsupported function type in [ dihedrals ] line: '+line); raise ValueError('Unsupported function type in [ dihedrals ] line: '+line);
self._currentMoleculeType.dihedrals.append(fields) self._currentMoleculeType.dihedrals.append(fields)
def _processExclusion(self, line): def _processExclusion(self, line):
"""Process a line in the [ exclusions ] category.""" """Process a line in the [ exclusions ] category."""
if self._currentMoleculeType is None: if self._currentMoleculeType is None:
...@@ -246,7 +246,7 @@ class GromacsTopFile(object): ...@@ -246,7 +246,7 @@ class GromacsTopFile(object):
if len(fields) < 2: if len(fields) < 2:
raise ValueError('Too few fields in [ exclusions ] line: '+line); raise ValueError('Too few fields in [ exclusions ] line: '+line);
self._currentMoleculeType.exclusions.append(fields) self._currentMoleculeType.exclusions.append(fields)
def _processPair(self, line): def _processPair(self, line):
"""Process a line in the [ pairs ] category.""" """Process a line in the [ pairs ] category."""
if self._currentMoleculeType is None: if self._currentMoleculeType is None:
...@@ -257,7 +257,7 @@ class GromacsTopFile(object): ...@@ -257,7 +257,7 @@ class GromacsTopFile(object):
if fields[2] != '1': if fields[2] != '1':
raise ValueError('Unsupported function type in [ pairs ] line: '+line); raise ValueError('Unsupported function type in [ pairs ] line: '+line);
self._currentMoleculeType.pairs.append(fields) self._currentMoleculeType.pairs.append(fields)
def _processCmap(self, line): def _processCmap(self, line):
"""Process a line in the [ cmaps ] category.""" """Process a line in the [ cmaps ] category."""
if self._currentMoleculeType is None: if self._currentMoleculeType is None:
...@@ -266,14 +266,14 @@ class GromacsTopFile(object): ...@@ -266,14 +266,14 @@ class GromacsTopFile(object):
if len(fields) < 6: if len(fields) < 6:
raise ValueError('Too few fields in [ pairs ] line: '+line); raise ValueError('Too few fields in [ pairs ] line: '+line);
self._currentMoleculeType.cmaps.append(fields) self._currentMoleculeType.cmaps.append(fields)
def _processAtomType(self, line): def _processAtomType(self, line):
"""Process a line in the [ atomtypes ] category.""" """Process a line in the [ atomtypes ] category."""
fields = line.split() fields = line.split()
if len(fields) < 7: if len(fields) < 7:
raise ValueError('Too few fields in [ atomtypes ] line: '+line); raise ValueError('Too few fields in [ atomtypes ] line: '+line);
self._atomTypes[fields[0]] = fields self._atomTypes[fields[0]] = fields
def _processBondType(self, line): def _processBondType(self, line):
"""Process a line in the [ bondtypes ] category.""" """Process a line in the [ bondtypes ] category."""
fields = line.split() fields = line.split()
...@@ -282,7 +282,7 @@ class GromacsTopFile(object): ...@@ -282,7 +282,7 @@ class GromacsTopFile(object):
if fields[2] != '1': if fields[2] != '1':
raise ValueError('Unsupported function type in [ bondtypes ] line: '+line); raise ValueError('Unsupported function type in [ bondtypes ] line: '+line);
self._bondTypes[tuple(fields[:2])] = fields self._bondTypes[tuple(fields[:2])] = fields
def _processAngleType(self, line): def _processAngleType(self, line):
"""Process a line in the [ angletypes ] category.""" """Process a line in the [ angletypes ] category."""
fields = line.split() fields = line.split()
...@@ -291,7 +291,7 @@ class GromacsTopFile(object): ...@@ -291,7 +291,7 @@ class GromacsTopFile(object):
if fields[3] not in ('1', '5'): if fields[3] not in ('1', '5'):
raise ValueError('Unsupported function type in [ angletypes ] line: '+line); raise ValueError('Unsupported function type in [ angletypes ] line: '+line);
self._angleTypes[tuple(fields[:3])] = fields self._angleTypes[tuple(fields[:3])] = fields
def _processDihedralType(self, line): def _processDihedralType(self, line):
"""Process a line in the [ dihedraltypes ] category.""" """Process a line in the [ dihedraltypes ] category."""
fields = line.split() fields = line.split()
...@@ -305,14 +305,14 @@ class GromacsTopFile(object): ...@@ -305,14 +305,14 @@ class GromacsTopFile(object):
self._dihedralTypes[key].append(fields) self._dihedralTypes[key].append(fields)
else: else:
self._dihedralTypes[key] = [fields] self._dihedralTypes[key] = [fields]
def _processImplicitType(self, line): def _processImplicitType(self, line):
"""Process a line in the [ implicit_genborn_params ] category.""" """Process a line in the [ implicit_genborn_params ] category."""
fields = line.split() fields = line.split()
if len(fields) < 6: if len(fields) < 6:
raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line); raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line);
self._implicitTypes[fields[0]] = fields self._implicitTypes[fields[0]] = fields
def _processPairType(self, line): def _processPairType(self, line):
"""Process a line in the [ pairtypes ] category.""" """Process a line in the [ pairtypes ] category."""
fields = line.split() fields = line.split()
...@@ -321,7 +321,7 @@ class GromacsTopFile(object): ...@@ -321,7 +321,7 @@ class GromacsTopFile(object):
if fields[2] != '1': if fields[2] != '1':
raise ValueError('Unsupported function type in [ pairtypes ] line: '+line); raise ValueError('Unsupported function type in [ pairtypes ] line: '+line);
self._pairTypes[tuple(fields[:2])] = fields self._pairTypes[tuple(fields[:2])] = fields
def _processCmapType(self, line): def _processCmapType(self, line):
"""Process a line in the [ cmaptypes ] category.""" """Process a line in the [ cmaptypes ] category."""
fields = line.split() fields = line.split()
...@@ -330,10 +330,10 @@ class GromacsTopFile(object): ...@@ -330,10 +330,10 @@ class GromacsTopFile(object):
if fields[5] != '1': if fields[5] != '1':
raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line); raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line);
self._cmapTypes[tuple(fields[:5])] = fields self._cmapTypes[tuple(fields[:5])] = fields
def __init__(self, file, unitCellDimensions=None, includeDir='/usr/local/gromacs/share/gromacs/top', defines={}): def __init__(self, file, unitCellDimensions=None, includeDir='/usr/local/gromacs/share/gromacs/top', defines={}):
"""Load a top file. """Load a top file.
Parameters: Parameters:
- file (string) the name of the file to load - file (string) the name of the file to load
- unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell - unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell
...@@ -343,9 +343,9 @@ class GromacsTopFile(object): ...@@ -343,9 +343,9 @@ class GromacsTopFile(object):
""" """
self._includeDirs = (os.path.dirname(file), includeDir) self._includeDirs = (os.path.dirname(file), includeDir)
self._defines = defines self._defines = defines
# Parse the file. # Parse the file.
self._currentCategory = None self._currentCategory = None
self._ifStack = [] self._ifStack = []
self._moleculeTypes = {} self._moleculeTypes = {}
...@@ -359,9 +359,9 @@ class GromacsTopFile(object): ...@@ -359,9 +359,9 @@ class GromacsTopFile(object):
self._pairTypes = {} self._pairTypes = {}
self._cmapTypes = {} self._cmapTypes = {}
self._processFile(file) self._processFile(file)
# Create the Topology from it. # Create the Topology from it.
top = Topology() top = Topology()
## The Topology read from the prmtop file ## The Topology read from the prmtop file
self.topology = top self.topology = top
...@@ -371,9 +371,9 @@ class GromacsTopFile(object): ...@@ -371,9 +371,9 @@ class GromacsTopFile(object):
if moleculeName not in self._moleculeTypes: if moleculeName not in self._moleculeTypes:
raise ValueError("Unknown molecule type: "+moleculeName) raise ValueError("Unknown molecule type: "+moleculeName)
moleculeType = self._moleculeTypes[moleculeName] moleculeType = self._moleculeTypes[moleculeName]
# Create the specified number of molecules of this type. # Create the specified number of molecules of this type.
for i in range(moleculeCount): for i in range(moleculeCount):
atoms = [] atoms = []
lastResidue = None lastResidue = None
...@@ -393,9 +393,9 @@ class GromacsTopFile(object): ...@@ -393,9 +393,9 @@ class GromacsTopFile(object):
atomName = fields[4] atomName = fields[4]
if atomName in atomReplacements: if atomName in atomReplacements:
atomName = atomReplacements[atomName] atomName = atomReplacements[atomName]
# Try to guess the element. # Try to guess the element.
upper = atomName.upper() upper = atomName.upper()
if upper.startswith('CL'): if upper.startswith('CL'):
element = elem.chlorine element = elem.chlorine
...@@ -409,16 +409,16 @@ class GromacsTopFile(object): ...@@ -409,16 +409,16 @@ class GromacsTopFile(object):
except KeyError: except KeyError:
element = None element = None
atoms.append(top.addAtom(atomName, element, r)) atoms.append(top.addAtom(atomName, element, r))
# Add bonds to the topology # Add bonds to the topology
for fields in moleculeType.bonds: for fields in moleculeType.bonds:
top.addBond(atoms[int(fields[0])-1], atoms[int(fields[1])-1]) top.addBond(atoms[int(fields[0])-1], atoms[int(fields[1])-1])
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): constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, ewaldErrorTolerance=0.0005, removeCMMotion=True):
"""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: Parameters:
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are - nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME. NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
...@@ -434,7 +434,7 @@ class GromacsTopFile(object): ...@@ -434,7 +434,7 @@ class GromacsTopFile(object):
Returns: the newly created System Returns: the newly created System
""" """
# Create the System. # Create the System.
sys = mm.System() sys = mm.System()
boxSize = self.topology.getUnitCellDimensions() boxSize = self.topology.getUnitCellDimensions()
if boxSize is not None: if boxSize is not None:
...@@ -462,33 +462,33 @@ class GromacsTopFile(object): ...@@ -462,33 +462,33 @@ 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])
# Loop over molecules and create the specified number of each type. # Loop over molecules and create the specified number of each type.
for moleculeName, moleculeCount in self._molecules: for moleculeName, moleculeCount in self._molecules:
moleculeType = self._moleculeTypes[moleculeName] moleculeType = self._moleculeTypes[moleculeName]
for i in range(moleculeCount): for i in range(moleculeCount):
# Record the types of all atoms. # Record the types of all atoms.
baseAtomIndex = sys.getNumParticles() baseAtomIndex = sys.getNumParticles()
atomTypes = [atom[1] for atom in moleculeType.atoms] atomTypes = [atom[1] for atom in moleculeType.atoms]
try: try:
[self._atomTypes[t][1] for t in atomTypes] [self._atomTypes[t][1] for t in atomTypes]
except KeyError as e: except KeyError as e:
raise ValueError('Unknown atom type: '+e.message) raise ValueError('Unknown atom type: '+e.message)
# Add atoms. # Add atoms.
for fields in moleculeType.atoms: for fields in moleculeType.atoms:
if len(fields) >= 8: if len(fields) >= 8:
mass = float(fields[7]) mass = float(fields[7])
else: else:
mass = float(self._atomTypes[fields[1]][2]) mass = float(self._atomTypes[fields[1]][2])
sys.addParticle(mass) sys.addParticle(mass)
# Add bonds. # Add bonds.
atomBonds = [{} for x in range(len(moleculeType.atoms))] atomBonds = [{} for x in range(len(moleculeType.atoms))]
for fields in moleculeType.bonds: for fields in moleculeType.bonds:
atoms = [int(x)-1 for x in fields[:2]] atoms = [int(x)-1 for x in fields[:2]]
...@@ -523,9 +523,9 @@ class GromacsTopFile(object): ...@@ -523,9 +523,9 @@ class GromacsTopFile(object):
# Record information that will be needed for constraining angles. # Record information that will be needed for constraining angles.
atomBonds[atoms[0]][atoms[1]] = length atomBonds[atoms[0]][atoms[1]] = length
atomBonds[atoms[1]][atoms[0]] = length atomBonds[atoms[1]][atoms[0]] = length
# Add angles. # Add angles.
degToRad = math.pi/180 degToRad = math.pi/180
for fields in moleculeType.angles: for fields in moleculeType.angles:
atoms = [int(x)-1 for x in fields[:3]] atoms = [int(x)-1 for x in fields[:3]]
...@@ -572,7 +572,7 @@ class GromacsTopFile(object): ...@@ -572,7 +572,7 @@ class GromacsTopFile(object):
bonds.addBond(baseAtomIndex+atoms[0], baseAtomIndex+atoms[2], float(params[2]), k) bonds.addBond(baseAtomIndex+atoms[0], baseAtomIndex+atoms[2], float(params[2]), k)
# Add torsions. # Add torsions.
for fields in moleculeType.dihedrals: for fields in moleculeType.dihedrals:
atoms = [int(x)-1 for x in fields[:4]] atoms = [int(x)-1 for x in fields[:4]]
types = tuple(atomTypes[i] for i in atoms) types = tuple(atomTypes[i] for i in atoms)
...@@ -618,7 +618,7 @@ class GromacsTopFile(object): ...@@ -618,7 +618,7 @@ class GromacsTopFile(object):
rb = mm.RBTorsionForce() rb = mm.RBTorsionForce()
sys.addForce(rb) sys.addForce(rb)
rb.addTorsion(baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], c[0], c[1], c[2], c[3], c[4], c[5]) rb.addTorsion(baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], c[0], c[1], c[2], c[3], c[4], c[5])
# Add CMAP terms. # Add CMAP terms.
for fields in moleculeType.cmaps: for fields in moleculeType.cmaps:
...@@ -649,7 +649,7 @@ class GromacsTopFile(object): ...@@ -649,7 +649,7 @@ class GromacsTopFile(object):
baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], baseAtomIndex+atoms[4]) baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], baseAtomIndex+atoms[4])
# Set nonbonded parameters for particles. # Set nonbonded parameters for particles.
for fields in moleculeType.atoms: for fields in moleculeType.atoms:
params = self._atomTypes[fields[1]] params = self._atomTypes[fields[1]]
if len(fields) > 6: if len(fields) > 6:
...@@ -665,9 +665,9 @@ class GromacsTopFile(object): ...@@ -665,9 +665,9 @@ class GromacsTopFile(object):
for fields in moleculeType.bonds: for fields in moleculeType.bonds:
atoms = [int(x)-1 for x in fields[:2]] atoms = [int(x)-1 for x in fields[:2]]
bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1])) bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1]))
# Record nonbonded exceptions. # Record nonbonded exceptions.
for fields in moleculeType.pairs: for fields in moleculeType.pairs:
atoms = [int(x)-1 for x in fields[:2]] atoms = [int(x)-1 for x in fields[:2]]
types = tuple(atomTypes[i] for i in atoms) types = tuple(atomTypes[i] for i in atoms)
...@@ -682,15 +682,15 @@ class GromacsTopFile(object): ...@@ -682,15 +682,15 @@ class GromacsTopFile(object):
atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0]) atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0])
atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1]) atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1])
exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1])) exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
# Create nonbonded exceptions. # Create nonbonded exceptions.
nb.createExceptionsFromBonds(bondIndices, fudgeQQ, float(self._defaults[3])) nb.createExceptionsFromBonds(bondIndices, fudgeQQ, float(self._defaults[3]))
for exception in exceptions: for exception in exceptions:
nb.addException(exception[0], exception[1], exception[2], float(exception[3]), float(exception[4]), True) nb.addException(exception[0], exception[1], exception[2], float(exception[3]), float(exception[4]), True)
# Finish configuring the NonbondedForce. # Finish configuring the NonbondedForce.
methodMap = {ff.NoCutoff:mm.NonbondedForce.NoCutoff, methodMap = {ff.NoCutoff:mm.NonbondedForce.NoCutoff,
ff.CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic, ff.CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
ff.CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic, ff.CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic,
......
wrappers/python/simtk/openmm/app/internal/amber_file_parser.py
View file @ 74415dd9
...@@ -16,7 +16,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors. ...@@ -16,7 +16,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors.
Authors: Randall J. Radmer, John D. Chodera, Peter Eastman Authors: Randall J. Radmer, John D. Chodera, Peter Eastman
Contributors: Christoph Klein, Michael R. Shirts Contributors: Christoph Klein, Michael R. Shirts
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -92,7 +92,7 @@ class PrmtopLoader(object): ...@@ -92,7 +92,7 @@ class PrmtopLoader(object):
>>> import os, os.path >>> import os, os.path
>>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit') >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
>>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop') >>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
>>> prmtop = PrmtopLoader(prmtop_filename) >>> prmtop = PrmtopLoader(prmtop_filename)
""" """
def __init__(self, inFilename): def __init__(self, inFilename):
...@@ -101,7 +101,7 @@ class PrmtopLoader(object): ...@@ -101,7 +101,7 @@ class PrmtopLoader(object):
ARGUMENTS ARGUMENTS
inFilename (string) - AMBER 'new-style' prmtop file, probably generated with one of the AMBER tleap/xleap/sleap inFilename (string) - AMBER 'new-style' prmtop file, probably generated with one of the AMBER tleap/xleap/sleap
""" """
...@@ -152,7 +152,7 @@ class PrmtopLoader(object): ...@@ -152,7 +152,7 @@ class PrmtopLoader(object):
Parameter: Parameter:
- pointerLabel: a string matching one of the following: - pointerLabel: a string matching one of the following:
NATOM : total number of atoms NATOM : total number of atoms
NTYPES : total number of distinct atom types NTYPES : total number of distinct atom types
NBONH : number of bonds containing hydrogen NBONH : number of bonds containing hydrogen
MBONA : number of bonds not containing hydrogen MBONA : number of bonds not containing hydrogen
...@@ -183,7 +183,7 @@ class PrmtopLoader(object): ...@@ -183,7 +183,7 @@ class PrmtopLoader(object):
NMXRS : number of atoms in the largest residue NMXRS : number of atoms in the largest residue
IFCAP : set to 1 if the CAP option from edit was specified IFCAP : set to 1 if the CAP option from edit was specified
""" """
index = POINTER_LABEL_LIST.index(pointerLabel) index = POINTER_LABEL_LIST.index(pointerLabel)
return float(self._raw_data['POINTERS'][index]) return float(self._raw_data['POINTERS'][index])
def getNumAtoms(self): def getNumAtoms(self):
...@@ -350,7 +350,7 @@ class PrmtopLoader(object): ...@@ -350,7 +350,7 @@ class PrmtopLoader(object):
bondPointers=self._raw_data["BONDS_INC_HYDROGEN"] bondPointers=self._raw_data["BONDS_INC_HYDROGEN"]
self._bondListWithH = self._getBonds(bondPointers) self._bondListWithH = self._getBonds(bondPointers)
return self._bondListWithH return self._bondListWithH
def getBondsNoH(self): def getBondsNoH(self):
"""Return list of bonded atom pairs, K, and Rmin for each bond with no hydrogen""" """Return list of bonded atom pairs, K, and Rmin for each bond with no hydrogen"""
...@@ -509,11 +509,11 @@ class PrmtopLoader(object): ...@@ -509,11 +509,11 @@ class PrmtopLoader(object):
def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5, nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=1.2, scnb=2.0, mm=None, verbose=False, EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True): def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5, nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=1.2, scnb=2.0, mm=None, verbose=False, EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True):
""" """
Create an OpenMM System from an Amber prmtop file. Create an OpenMM System from an Amber prmtop file.
ARGUMENTS (specify one or the other, but not both) ARGUMENTS (specify one or the other, but not both)
prmtop_filename (String) - name of Amber prmtop file (new-style only) prmtop_filename (String) - name of Amber prmtop file (new-style only)
prmtop_loader (PrmtopLoader) - the loaded prmtop file prmtop_loader (PrmtopLoader) - the loaded prmtop file
OPTIONAL ARGUMENTS OPTIONAL ARGUMENTS
shake (String) - if 'h-bonds', will SHAKE all bonds to hydrogen and water; if 'all-bonds', will SHAKE all bonds and water (default: None) shake (String) - if 'h-bonds', will SHAKE all bonds to hydrogen and water; if 'all-bonds', will SHAKE all bonds and water (default: None)
gbmodel (String) - if 'OBC', OBC GBSA will be used; if 'GBVI', GB/VI will be used (default: None) gbmodel (String) - if 'OBC', OBC GBSA will be used; if 'GBVI', GB/VI will be used (default: None)
...@@ -547,10 +547,10 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -547,10 +547,10 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
>>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit') >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
>>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop') >>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
>>> system = readAmberSystem(prmtop_filename) >>> system = readAmberSystem(prmtop_filename)
""" """
if prmtop_filename is None and prmtop_loader is None: if prmtop_filename is None and prmtop_loader is None:
raise Exception("Must specify a filename or loader") raise Exception("Must specify a filename or loader")
if prmtop_filename is not None and prmtop_loader is not None: if prmtop_filename is not None and prmtop_loader is not None:
...@@ -567,7 +567,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -567,7 +567,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
if prmtop.getIfPert()>0: if prmtop.getIfPert()>0:
raise Exception("perturbation not currently supported") raise Exception("perturbation not currently supported")
if prmtop.getIfBox()>1: if prmtop.getIfBox()>1:
raise Exception("only standard periodic boxes are currently supported") raise Exception("only standard periodic boxes are currently supported")
...@@ -596,9 +596,9 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -596,9 +596,9 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH(): for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
if isWater[iAtom] and isWater[jAtom]: if isWater[iAtom] and isWater[jAtom]:
system.addConstraint(iAtom, jAtom, rMin) system.addConstraint(iAtom, jAtom, rMin)
# Add harmonic bonds. # Add harmonic bonds.
if verbose: print "Adding bonds..." if verbose: print "Adding bonds..."
force = mm.HarmonicBondForce() force = mm.HarmonicBondForce()
if flexibleConstraints or (shake not in ('h-bonds', 'all-bonds', 'h-angles')): if flexibleConstraints or (shake not in ('h-bonds', 'all-bonds', 'h-angles')):
for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH(): for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
...@@ -610,7 +610,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -610,7 +610,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
system.addForce(force) system.addForce(force)
# Add harmonic angles. # Add harmonic angles.
if verbose: print "Adding angles..." if verbose: print "Adding angles..."
force = mm.HarmonicAngleForce() force = mm.HarmonicAngleForce()
if shake == 'h-angles': if shake == 'h-angles':
numConstrainedBonds = system.getNumConstraints() numConstrainedBonds = system.getNumConstraints()
...@@ -638,7 +638,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -638,7 +638,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
l1 = bond[1] l1 = bond[1]
elif bond[0] == kAtom: elif bond[0] == kAtom:
l2 = bond[1] l2 = bond[1]
# Compute the distance between atoms and add a constraint # Compute the distance between atoms and add a constraint
length = math.sqrt(l1*l1 + l2*l2 - 2*l1*l2*math.cos(aMin)) length = math.sqrt(l1*l1 + l2*l2 - 2*l1*l2*math.cos(aMin))
system.addConstraint(iAtom, kAtom, length) system.addConstraint(iAtom, kAtom, length)
...@@ -647,14 +647,14 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -647,14 +647,14 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
system.addForce(force) system.addForce(force)
# Add torsions. # Add torsions.
if verbose: print "Adding torsions..." if verbose: print "Adding torsions..."
force = mm.PeriodicTorsionForce() force = mm.PeriodicTorsionForce()
for (iAtom, jAtom, kAtom, lAtom, forceConstant, phase, periodicity) in prmtop.getDihedrals(): for (iAtom, jAtom, kAtom, lAtom, forceConstant, phase, periodicity) in prmtop.getDihedrals():
force.addTorsion(iAtom, jAtom, kAtom, lAtom, periodicity, phase, forceConstant) force.addTorsion(iAtom, jAtom, kAtom, lAtom, periodicity, phase, forceConstant)
system.addForce(force) system.addForce(force)
# Add nonbonded interactions. # Add nonbonded interactions.
if verbose: print "Adding nonbonded interactions..." if verbose: print "Adding nonbonded interactions..."
force = mm.NonbondedForce() force = mm.NonbondedForce()
if (prmtop.getIfBox() == 0): if (prmtop.getIfBox() == 0):
# System is non-periodic. # System is non-periodic.
...@@ -663,12 +663,12 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -663,12 +663,12 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
elif nonbondedMethod == 'CutoffNonPeriodic': elif nonbondedMethod == 'CutoffNonPeriodic':
if nonbondedCutoff is None: if nonbondedCutoff is None:
raise Exception("No cutoff value specified") raise Exception("No cutoff value specified")
force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic) force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
force.setCutoffDistance(nonbondedCutoff) force.setCutoffDistance(nonbondedCutoff)
else: else:
raise Exception("Illegal nonbonded method for a non-periodic system") raise Exception("Illegal nonbonded method for a non-periodic system")
else: else:
# System is periodic. # System is periodic.
# Set periodic box vectors for periodic system # Set periodic box vectors for periodic system
(boxBeta, boxX, boxY, boxZ) = prmtop.getBoxBetaAndDimensions() (boxBeta, boxX, boxY, boxZ) = prmtop.getBoxBetaAndDimensions()
d0 = units.Quantity(0.0, units.angstroms) d0 = units.Quantity(0.0, units.angstroms)
...@@ -676,16 +676,16 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -676,16 +676,16 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
yVec = units.Quantity((d0, boxY, d0)) yVec = units.Quantity((d0, boxY, d0))
zVec = units.Quantity((d0, d0, boxZ)) zVec = units.Quantity((d0, d0, boxZ))
system.setDefaultPeriodicBoxVectors(xVec, yVec, zVec) system.setDefaultPeriodicBoxVectors(xVec, yVec, zVec)
# Set cutoff. # Set cutoff.
if nonbondedCutoff is None: if nonbondedCutoff is None:
# Compute cutoff automatically. # Compute cutoff automatically.
min_box_width = min([boxX / units.nanometers, boxY / units.nanometers, boxZ / units.nanometers]) min_box_width = min([boxX / units.nanometers, boxY / units.nanometers, boxZ / units.nanometers])
CLEARANCE_FACTOR = 0.97 # reduce the cutoff to be a bit smaller than 1/2 smallest box length CLEARANCE_FACTOR = 0.97 # reduce the cutoff to be a bit smaller than 1/2 smallest box length
nonbondedCutoff = units.Quantity((min_box_width * CLEARANCE_FACTOR) / 2.0, units.nanometers) nonbondedCutoff = units.Quantity((min_box_width * CLEARANCE_FACTOR) / 2.0, units.nanometers)
if nonbondedMethod != 'NoCutoff': if nonbondedMethod != 'NoCutoff':
force.setCutoffDistance(nonbondedCutoff) force.setCutoffDistance(nonbondedCutoff)
# Set nonbonded method. # Set nonbonded method.
if nonbondedMethod == 'NoCutoff': if nonbondedMethod == 'NoCutoff':
force.setNonbondedMethod(mm.NonbondedForce.NoCutoff) force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
...@@ -724,7 +724,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -724,7 +724,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
excludeParams = (0.0, 0.1, 0.0) excludeParams = (0.0, 0.1, 0.0)
for iAtom in range(prmtop.getNumAtoms()): for iAtom in range(prmtop.getNumAtoms()):
for jAtom in excludedAtoms[iAtom]: for jAtom in excludedAtoms[iAtom]:
if min((iAtom, jAtom), (jAtom, iAtom)) in excludedAtomPairs: continue if min((iAtom, jAtom), (jAtom, iAtom)) in excludedAtomPairs: continue
force.addException(iAtom, jAtom, excludeParams[0], excludeParams[1], excludeParams[2]) force.addException(iAtom, jAtom, excludeParams[0], excludeParams[1], excludeParams[2])
system.addForce(force) system.addForce(force)
...@@ -812,7 +812,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -812,7 +812,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
if gbmodel == 'OBC2': if gbmodel == 'OBC2':
gb.addParticle(charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]) gb.addParticle(charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1])
else: else:
gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]]) gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]])
system.addForce(gb) system.addForce(gb)
if nonbondedMethod == 'NoCutoff': if nonbondedMethod == 'NoCutoff':
gb.setNonbondedMethod(mm.NonbondedForce.NoCutoff) gb.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
...@@ -836,7 +836,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode ...@@ -836,7 +836,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbose=False, asNumpy=False): def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbose=False, asNumpy=False):
""" """
Read atomic coordinates (and optionally, box vectors) from Amber formatted coordinate file. Read atomic coordinates (and optionally, box vectors) from Amber formatted coordinate file.
ARGUMENTS ARGUMENTS
...@@ -853,13 +853,13 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos ...@@ -853,13 +853,13 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
Read coordinates in vacuum. Read coordinates in vacuum.
>>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-gbsa') >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-gbsa')
>>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd') >>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')
>>> coordinates = readAmberCoordinates(crd_filename) >>> coordinates = readAmberCoordinates(crd_filename)
Read coordinates in solvent. Read coordinates in solvent.
>>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit') >>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
>>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd') >>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')
>>> [coordinates, box_vectors] = readAmberCoordinates(crd_filename, read_box=True) >>> [coordinates, box_vectors] = readAmberCoordinates(crd_filename, read_box=True)
""" """
...@@ -922,7 +922,7 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos ...@@ -922,7 +922,7 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
velocities = newvel velocities = newvel
# Assign units. # Assign units.
velocities = units.Quantity(velocities, units.angstroms/units.picoseconds) velocities = units.Quantity(velocities, units.angstroms/units.picoseconds)
# Read box size if present # Read box size if present
box_vectors = None box_vectors = None
if (read_box): if (read_box):
...@@ -944,20 +944,20 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos ...@@ -944,20 +944,20 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
else: else:
a = units.Quantity(mm.Vec3(box_dimensions[0], 0.0, 0.0), units.angstroms) a = units.Quantity(mm.Vec3(box_dimensions[0], 0.0, 0.0), units.angstroms)
b = units.Quantity(mm.Vec3(0.0, box_dimensions[1], 0.0), units.angstroms) b = units.Quantity(mm.Vec3(0.0, box_dimensions[1], 0.0), units.angstroms)
c = units.Quantity(mm.Vec3(0.0, 0.0, box_dimensions[2]), units.angstroms) c = units.Quantity(mm.Vec3(0.0, 0.0, box_dimensions[2]), units.angstroms)
box_vectors = [a,b,c] box_vectors = [a,b,c]
else: else:
raise Exception("Don't know what to do with box vectors: %s" % line) raise Exception("Don't know what to do with box vectors: %s" % line)
# Close file # Close file
infile.close() infile.close()
if box_vectors and velocities: if box_vectors and velocities:
return (coordinates, box_vectors, velocities) return (coordinates, box_vectors, velocities)
if box_vectors: if box_vectors:
return (coordinates, box_vectors) return (coordinates, box_vectors)
if velocities: if velocities:
return (coordinates, velocities) return (coordinates, velocities)
return coordinates return coordinates
#============================================================================================= #=============================================================================================
......
wrappers/python/simtk/openmm/app/internal/customgbforces.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 University of Virginia and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 University of Virginia and the Authors.
Authors: Christoph Klein, Michael R. Shirts Authors: Christoph Klein, Michael R. Shirts
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -35,15 +35,15 @@ import sys, pdb, pickle ...@@ -35,15 +35,15 @@ import sys, pdb, pickle
d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444, d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
2.73212,2.78965,2.84705,2.9043,2.96141,3.0184,3.07524,3.13196, 2.73212,2.78965,2.84705,2.9043,2.96141,3.0184,3.07524,3.13196,
3.18854,3.24498,3.30132,3.35752,3.4136, 3.18854,3.24498,3.30132,3.35752,3.4136,
2.31191,2.37017,2.4283,2.48632,2.5442,2.60197,2.65961,2.71711, 2.31191,2.37017,2.4283,2.48632,2.5442,2.60197,2.65961,2.71711,
2.77449,2.83175,2.88887,2.94586,3.00273,3.05948,3.1161,3.1726, 2.77449,2.83175,2.88887,2.94586,3.00273,3.05948,3.1161,3.1726,
3.22897,3.28522,3.34136,3.39738,3.45072, 3.22897,3.28522,3.34136,3.39738,3.45072,
2.35759,2.41549,2.47329,2.53097,2.58854,2.646,2.70333,2.76056, 2.35759,2.41549,2.47329,2.53097,2.58854,2.646,2.70333,2.76056,
2.81766,2.87465,2.93152,2.98827,3.0449,3.10142,3.15782,3.21411, 2.81766,2.87465,2.93152,2.98827,3.0449,3.10142,3.15782,3.21411,
3.27028,3.32634,3.3823,3.43813,3.49387, 3.27028,3.32634,3.3823,3.43813,3.49387,
2.4038,2.46138,2.51885,2.57623,2.63351,2.69067,2.74773,2.80469, 2.4038,2.46138,2.51885,2.57623,2.63351,2.69067,2.74773,2.80469,
2.86152,2.91826,2.97489,3.0314,3.08781,3.1441,3.20031,3.25638, 2.86152,2.91826,2.97489,3.0314,3.08781,3.1441,3.20031,3.25638,
3.31237,3.36825,3.42402,3.4797,3.53527, 3.31237,3.36825,3.42402,3.4797,3.53527,
2.45045,2.50773,2.56492,2.62201,2.679,2.7359,2.7927,2.8494,2.90599, 2.45045,2.50773,2.56492,2.62201,2.679,2.7359,2.7927,2.8494,2.90599,
2.9625,3.0189,3.07518,3.13138,3.18748,3.24347,3.29937,3.35515, 2.9625,3.0189,3.07518,3.13138,3.18748,3.24347,3.29937,3.35515,
3.41085,3.46646,3.52196,3.57738, 3.41085,3.46646,3.52196,3.57738,
...@@ -53,134 +53,134 @@ d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444, ...@@ -53,134 +53,134 @@ d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
2.54489,2.60164,2.6583,2.71488,2.77134,2.8278,2.88412,2.94034, 2.54489,2.60164,2.6583,2.71488,2.77134,2.8278,2.88412,2.94034,
2.9965,3.05256,3.10853,3.16442,3.22021,3.27592,3.33154,3.38707, 2.9965,3.05256,3.10853,3.16442,3.22021,3.27592,3.33154,3.38707,
3.44253,3.49789,3.55316,3.60836,3.66348, 3.44253,3.49789,3.55316,3.60836,3.66348,
2.59259,2.6491,2.70553,2.76188,2.81815,2.87434,2.93044,2.98646, 2.59259,2.6491,2.70553,2.76188,2.81815,2.87434,2.93044,2.98646,
3.04241,3.09827,3.15404,3.20974,3.26536,3.32089,3.37633,3.4317, 3.04241,3.09827,3.15404,3.20974,3.26536,3.32089,3.37633,3.4317,
3.48699,3.54219,3.59731,3.65237,3.70734, 3.48699,3.54219,3.59731,3.65237,3.70734,
2.64054,2.69684,2.75305,2.80918,2.86523,2.92122,2.97712,3.03295, 2.64054,2.69684,2.75305,2.80918,2.86523,2.92122,2.97712,3.03295,
3.0887,3.14437,3.19996,3.25548,3.31091,3.36627,3.42156,3.47677, 3.0887,3.14437,3.19996,3.25548,3.31091,3.36627,3.42156,3.47677,
3.5319,3.58695,3.64193,3.69684,3.75167, 3.5319,3.58695,3.64193,3.69684,3.75167,
2.68873,2.74482,2.80083,2.85676,2.91262,2.96841,3.02412,3.07976, 2.68873,2.74482,2.80083,2.85676,2.91262,2.96841,3.02412,3.07976,
3.13533,3.19082,3.24623,3.30157,3.35685,3.41205,3.46718,3.52223, 3.13533,3.19082,3.24623,3.30157,3.35685,3.41205,3.46718,3.52223,
3.57721,3.63213,3.68696,3.74174,3.79644, 3.57721,3.63213,3.68696,3.74174,3.79644,
2.73713,2.79302,2.84884,2.90459,2.96027,3.01587,3.0714,3.12686, 2.73713,2.79302,2.84884,2.90459,2.96027,3.01587,3.0714,3.12686,
3.18225,3.23757,3.29282,3.34801,3.40313,3.45815,3.51315,3.56805, 3.18225,3.23757,3.29282,3.34801,3.40313,3.45815,3.51315,3.56805,
3.6229,3.67767,3.73237,3.78701,3.84159, 3.6229,3.67767,3.73237,3.78701,3.84159,
2.78572,2.84143,2.89707,2.95264,3.00813,3.06356,3.11892,3.17422, 2.78572,2.84143,2.89707,2.95264,3.00813,3.06356,3.11892,3.17422,
3.22946,3.28462,3.33971,3.39474,3.44971,3.5046,3.55944,3.61421, 3.22946,3.28462,3.33971,3.39474,3.44971,3.5046,3.55944,3.61421,
3.66891,3.72356,3.77814,3.83264,3.8871, 3.66891,3.72356,3.77814,3.83264,3.8871,
2.83446,2.89,2.94547,3.00088,3.05621,3.11147,3.16669,3.22183, 2.83446,2.89,2.94547,3.00088,3.05621,3.11147,3.16669,3.22183,
3.27689,3.33191,3.38685,3.44174,3.49656,3.55132,3.60602,3.66066, 3.27689,3.33191,3.38685,3.44174,3.49656,3.55132,3.60602,3.66066,
3.71523,3.76975,3.82421,3.8786,3.93293, 3.71523,3.76975,3.82421,3.8786,3.93293,
2.88335,2.93873,2.99404,3.04929,3.10447,3.15959,3.21464,3.26963, 2.88335,2.93873,2.99404,3.04929,3.10447,3.15959,3.21464,3.26963,
3.32456,3.37943,3.43424,3.48898,3.54366,3.5983,3.65287,3.70737, 3.32456,3.37943,3.43424,3.48898,3.54366,3.5983,3.65287,3.70737,
3.76183,3.81622,3.87056,3.92484,3.97905, 3.76183,3.81622,3.87056,3.92484,3.97905,
2.93234,2.9876,3.04277,3.09786,3.15291,3.20787,3.26278,3.31764, 2.93234,2.9876,3.04277,3.09786,3.15291,3.20787,3.26278,3.31764,
3.37242,3.42716,3.48184,3.53662,3.591,3.64551,3.69995,3.75435, 3.37242,3.42716,3.48184,3.53662,3.591,3.64551,3.69995,3.75435,
3.80867,3.86295,3.91718,3.97134,4.02545, 3.80867,3.86295,3.91718,3.97134,4.02545,
2.98151,3.0366,3.09163,3.14659,3.20149,3.25632,3.3111,3.36581, 2.98151,3.0366,3.09163,3.14659,3.20149,3.25632,3.3111,3.36581,
3.42047,3.47507,3.52963,3.58411,3.63855,3.69293,3.74725,3.80153, 3.42047,3.47507,3.52963,3.58411,3.63855,3.69293,3.74725,3.80153,
3.85575,3.90991,3.96403,4.01809,4.07211, 3.85575,3.90991,3.96403,4.01809,4.07211,
3.03074,3.08571,3.14061,3.19543,3.25021,3.30491,3.35956,3.41415, 3.03074,3.08571,3.14061,3.19543,3.25021,3.30491,3.35956,3.41415,
3.46869,3.52317,3.57759,3.63196,3.68628,3.74054,3.79476,3.84893, 3.46869,3.52317,3.57759,3.63196,3.68628,3.74054,3.79476,3.84893,
3.90303,3.95709,4.01111,4.06506,4.11897, 3.90303,3.95709,4.01111,4.06506,4.11897,
3.08008,3.13492,3.1897,3.2444,3.29905,3.35363,3.40815,3.46263, 3.08008,3.13492,3.1897,3.2444,3.29905,3.35363,3.40815,3.46263,
3.51704,3.57141,3.62572,3.67998,3.73418,3.78834,3.84244,3.8965, 3.51704,3.57141,3.62572,3.67998,3.73418,3.78834,3.84244,3.8965,
3.95051,4.00447,4.05837,4.11224,4.16605, 3.95051,4.00447,4.05837,4.11224,4.16605,
3.12949,3.18422,3.23888,3.29347,3.348,3.40247,3.45688,3.51124, 3.12949,3.18422,3.23888,3.29347,3.348,3.40247,3.45688,3.51124,
3.56554,3.6198,3.674,3.72815,3.78225,3.83629,3.8903,3.94425, 3.56554,3.6198,3.674,3.72815,3.78225,3.83629,3.8903,3.94425,
3.99816,4.05203,4.10583,4.15961,4.21333, 3.99816,4.05203,4.10583,4.15961,4.21333,
3.17899,3.23361,3.28815,3.34264,3.39706,3.45142,3.50571,3.55997, 3.17899,3.23361,3.28815,3.34264,3.39706,3.45142,3.50571,3.55997,
3.61416,3.66831,3.72241,3.77645,3.83046,3.8844,3.93831,3.99216, 3.61416,3.66831,3.72241,3.77645,3.83046,3.8844,3.93831,3.99216,
4.04598,4.09974,4.15347,4.20715,4.26078, 4.04598,4.09974,4.15347,4.20715,4.26078,
3.22855,3.28307,3.33751,3.39188,3.4462,3.50046,3.55466,3.6088, 3.22855,3.28307,3.33751,3.39188,3.4462,3.50046,3.55466,3.6088,
3.6629,3.71694,3.77095,3.82489,3.8788,3.93265,3.98646,4.04022, 3.6629,3.71694,3.77095,3.82489,3.8788,3.93265,3.98646,4.04022,
4.09395,4.14762,4.20126,4.25485,4.3084] 4.09395,4.14762,4.20126,4.25485,4.3084]
m0=[0.0381511,0.0338587,0.0301776,0.027003,0.0242506,0.0218529, m0=[0.0381511,0.0338587,0.0301776,0.027003,0.0242506,0.0218529,
0.0197547,0.0179109,0.0162844,0.0148442,0.0135647,0.0124243, 0.0197547,0.0179109,0.0162844,0.0148442,0.0135647,0.0124243,
0.0114047,0.0104906,0.00966876,0.008928,0.0082587,0.00765255, 0.0114047,0.0104906,0.00966876,0.008928,0.0082587,0.00765255,
0.00710237,0.00660196,0.00614589, 0.00710237,0.00660196,0.00614589,
0.0396198,0.0351837,0.0313767,0.0280911,0.0252409,0.0227563, 0.0396198,0.0351837,0.0313767,0.0280911,0.0252409,0.0227563,
0.0205808,0.0186681,0.0169799,0.0154843,0.014155,0.0129696, 0.0205808,0.0186681,0.0169799,0.0154843,0.014155,0.0129696,
0.0119094,0.0109584,0.0101031,0.00933189,0.0086348,0.00800326, 0.0119094,0.0109584,0.0101031,0.00933189,0.0086348,0.00800326,
0.00742986,0.00690814,0.00643255, 0.00742986,0.00690814,0.00643255,
0.041048,0.0364738,0.0325456,0.0291532,0.0262084,0.0236399, 0.041048,0.0364738,0.0325456,0.0291532,0.0262084,0.0236399,
0.0213897,0.0194102,0.0176622,0.0161129,0.0147351,0.0135059, 0.0213897,0.0194102,0.0176622,0.0161129,0.0147351,0.0135059,
0.0124061,0.0114192,0.0105312,0.00973027,0.00900602,0.00834965, 0.0124061,0.0114192,0.0105312,0.00973027,0.00900602,0.00834965,
0.0077535,0.00721091,0.00671609, 0.0077535,0.00721091,0.00671609,
0.0424365,0.0377295,0.0336846,0.0301893,0.0271533,0.0245038, 0.0424365,0.0377295,0.0336846,0.0301893,0.0271533,0.0245038,
0.0221813,0.0201371,0.018331,0.0167295,0.0153047,0.014033, 0.0221813,0.0201371,0.018331,0.0167295,0.0153047,0.014033,
0.0128946,0.0118727,0.0109529,0.0101229,0.00937212,0.00869147, 0.0128946,0.0118727,0.0109529,0.0101229,0.00937212,0.00869147,
0.00807306,0.00751003,0.00699641, 0.00807306,0.00751003,0.00699641,
0.0437861,0.0389516,0.0347944,0.0311998,0.0280758,0.0253479, 0.0437861,0.0389516,0.0347944,0.0311998,0.0280758,0.0253479,
0.0229555,0.0208487,0.0189864,0.0173343,0.0158637,0.0145507, 0.0229555,0.0208487,0.0189864,0.0173343,0.0158637,0.0145507,
0.0133748,0.0123188,0.0113679,0.0105096,0.0097329,0.00902853, 0.0133748,0.0123188,0.0113679,0.0105096,0.0097329,0.00902853,
0.00838835,0.00780533,0.0072733, 0.00838835,0.00780533,0.0072733,
0.0450979,0.0401406,0.0358753,0.0321851,0.0289761,0.0261726, 0.0450979,0.0401406,0.0358753,0.0321851,0.0289761,0.0261726,
0.0237125,0.0215451,0.0196282,0.017927,0.0164121,0.0150588, 0.0237125,0.0215451,0.0196282,0.017927,0.0164121,0.0150588,
0.0138465,0.0127573,0.0117761,0.0108902,0.0100882,0.00936068, 0.0138465,0.0127573,0.0117761,0.0108902,0.0100882,0.00936068,
0.00869923,0.00809665,0.00754661, 0.00869923,0.00809665,0.00754661,
0.0463729,0.0412976,0.0369281,0.0331456,0.0298547,0.026978, 0.0463729,0.0412976,0.0369281,0.0331456,0.0298547,0.026978,
0.0244525,0.0222264,0.0202567,0.0185078,0.0169498,0.0155575, 0.0244525,0.0222264,0.0202567,0.0185078,0.0169498,0.0155575,
0.0143096,0.0131881,0.0121775,0.0112646,0.010438,0.00968781, 0.0143096,0.0131881,0.0121775,0.0112646,0.010438,0.00968781,
0.00900559,0.00838388,0.00781622, 0.00900559,0.00838388,0.00781622,
0.0476123,0.0424233,0.0379534,0.034082,0.0307118,0.0277645, 0.0476123,0.0424233,0.0379534,0.034082,0.0307118,0.0277645,
0.0251757,0.0228927,0.0208718,0.0190767,0.0174768,0.0160466, 0.0251757,0.0228927,0.0208718,0.0190767,0.0174768,0.0160466,
0.0147642,0.0136112,0.0125719,0.0116328,0.0107821,0.0100099, 0.0147642,0.0136112,0.0125719,0.0116328,0.0107821,0.0100099,
0.00930735,0.00866695,0.00808206, 0.00930735,0.00866695,0.00808206,
0.0488171,0.0435186,0.038952,0.0349947,0.0315481,0.0285324, 0.0488171,0.0435186,0.038952,0.0349947,0.0315481,0.0285324,
0.0258824,0.0235443,0.0214738,0.0196339,0.0179934,0.0165262, 0.0258824,0.0235443,0.0214738,0.0196339,0.0179934,0.0165262,
0.0152103,0.0140267,0.0129595,0.0119947,0.0111206,0.0103268, 0.0152103,0.0140267,0.0129595,0.0119947,0.0111206,0.0103268,
0.00960445,0.00894579,0.00834405, 0.00960445,0.00894579,0.00834405,
0.0499883,0.0445845,0.0399246,0.0358844,0.032364,0.0292822, 0.0499883,0.0445845,0.0399246,0.0358844,0.032364,0.0292822,
0.0265729,0.0241815,0.0220629,0.0201794,0.0184994,0.0169964, 0.0265729,0.0241815,0.0220629,0.0201794,0.0184994,0.0169964,
0.0156479,0.0144345,0.0133401,0.0123504,0.0114534,0.0106386, 0.0156479,0.0144345,0.0133401,0.0123504,0.0114534,0.0106386,
0.00989687,0.00922037,0.00860216, 0.00989687,0.00922037,0.00860216,
0.0511272,0.0456219,0.040872,0.0367518,0.0331599,0.0300142, 0.0511272,0.0456219,0.040872,0.0367518,0.0331599,0.0300142,
0.0272475,0.0248045,0.0226392,0.0207135,0.0189952,0.0174574, 0.0272475,0.0248045,0.0226392,0.0207135,0.0189952,0.0174574,
0.0160771,0.0148348,0.0137138,0.0126998,0.0117805,0.0109452, 0.0160771,0.0148348,0.0137138,0.0126998,0.0117805,0.0109452,
0.0101846,0.00949067,0.00885636, 0.0101846,0.00949067,0.00885636,
0.0522348,0.0466315,0.0417948,0.0375973,0.0339365,0.030729, 0.0522348,0.0466315,0.0417948,0.0375973,0.0339365,0.030729,
0.0279067,0.0254136,0.023203,0.0212363,0.0194809,0.0179092, 0.0279067,0.0254136,0.023203,0.0212363,0.0194809,0.0179092,
0.016498,0.0152275,0.0140807,0.013043,0.012102,0.0112466, 0.016498,0.0152275,0.0140807,0.013043,0.012102,0.0112466,
0.0104676,0.00975668,0.00910664, 0.0104676,0.00975668,0.00910664,
0.0533123,0.0476145,0.042694,0.0384218,0.0346942,0.0314268, 0.0533123,0.0476145,0.042694,0.0384218,0.0346942,0.0314268,
0.0285507,0.026009,0.0237547,0.0217482,0.0199566,0.018352, 0.0285507,0.026009,0.0237547,0.0217482,0.0199566,0.018352,
0.0169108,0.0156128,0.0144408,0.0133801,0.0124179,0.011543, 0.0169108,0.0156128,0.0144408,0.0133801,0.0124179,0.011543,
0.010746,0.0100184,0.00935302, 0.010746,0.0100184,0.00935302,
0.0543606,0.0485716,0.04357,0.0392257,0.0354335,0.0321082, 0.0543606,0.0485716,0.04357,0.0392257,0.0354335,0.0321082,
0.02918,0.0265913,0.0242943,0.0222492,0.0204225,0.0187859, 0.02918,0.0265913,0.0242943,0.0222492,0.0204225,0.0187859,
0.0173155,0.0159908,0.0147943,0.0137111,0.0127282,0.0118343, 0.0173155,0.0159908,0.0147943,0.0137111,0.0127282,0.0118343,
0.0110197,0.0102759,0.00959549, 0.0110197,0.0102759,0.00959549,
0.0553807,0.0495037,0.0444239,0.0400097,0.0361551,0.0327736, 0.0553807,0.0495037,0.0444239,0.0400097,0.0361551,0.0327736,
0.0297949,0.0271605,0.0248222,0.0227396,0.0208788,0.0192111, 0.0297949,0.0271605,0.0248222,0.0227396,0.0208788,0.0192111,
0.0177122,0.0163615,0.0151413,0.0140361,0.013033,0.0121206, 0.0177122,0.0163615,0.0151413,0.0140361,0.013033,0.0121206,
0.0112888,0.0105292,0.00983409, 0.0112888,0.0105292,0.00983409,
0.0563738,0.0504116,0.0452562,0.0407745,0.0368593,0.0334235, 0.0563738,0.0504116,0.0452562,0.0407745,0.0368593,0.0334235,
0.0303958,0.0277171,0.0253387,0.0232197,0.0213257,0.0196277, 0.0303958,0.0277171,0.0253387,0.0232197,0.0213257,0.0196277,
0.0181013,0.0167252,0.0154817,0.0143552,0.0133325,0.0124019, 0.0181013,0.0167252,0.0154817,0.0143552,0.0133325,0.0124019,
0.0115534,0.0107783,0.0100688, 0.0115534,0.0107783,0.0100688,
0.0573406,0.0512963,0.0460676,0.0415206,0.0375468,0.0340583, 0.0573406,0.0512963,0.0460676,0.0415206,0.0375468,0.0340583,
0.030983,0.0282614,0.0258441,0.0236896,0.0217634,0.020036, 0.030983,0.0282614,0.0258441,0.0236896,0.0217634,0.020036,
0.0184826,0.017082,0.0158158,0.0146685,0.0136266,0.0126783, 0.0184826,0.017082,0.0158158,0.0146685,0.0136266,0.0126783,
0.0118135,0.0110232,0.0102998, 0.0118135,0.0110232,0.0102998,
0.0582822,0.0521584,0.0468589,0.0422486,0.038218,0.0346784, 0.0582822,0.0521584,0.0468589,0.0422486,0.038218,0.0346784,
0.0315571,0.0287938,0.0263386,0.0241497,0.0221922,0.0204362, 0.0315571,0.0287938,0.0263386,0.0241497,0.0221922,0.0204362,
0.0188566,0.0174319,0.0161437,0.0149761,0.0139154,0.0129499, 0.0188566,0.0174319,0.0161437,0.0149761,0.0139154,0.0129499,
0.0120691,0.0112641,0.0105269, 0.0120691,0.0112641,0.0105269,
0.0591994,0.0529987,0.0476307,0.042959,0.0388734,0.0352843, 0.0591994,0.0529987,0.0476307,0.042959,0.0388734,0.0352843,
0.0321182,0.0293144,0.0268225,0.0246002,0.0226121,0.0208283, 0.0321182,0.0293144,0.0268225,0.0246002,0.0226121,0.0208283,
0.0192232,0.0177751,0.0164654,0.015278,0.0141991,0.0132167, 0.0192232,0.0177751,0.0164654,0.015278,0.0141991,0.0132167,
0.0123204,0.0115009,0.0107504, 0.0123204,0.0115009,0.0107504,
0.0600932,0.053818,0.0483836,0.0436525,0.0395136,0.0358764, 0.0600932,0.053818,0.0483836,0.0436525,0.0395136,0.0358764,
0.0326669,0.0298237,0.0272961,0.0250413,0.0230236,0.0212126, 0.0326669,0.0298237,0.0272961,0.0250413,0.0230236,0.0212126,
0.0195826,0.0181118,0.0167811,0.0155744,0.0144778,0.0134789, 0.0195826,0.0181118,0.0167811,0.0155744,0.0144778,0.0134789,
0.0125673,0.0117338,0.0109702, 0.0125673,0.0117338,0.0109702,
0.0609642,0.0546169,0.0491183,0.0443295,0.0401388,0.036455, 0.0609642,0.0546169,0.0491183,0.0443295,0.0401388,0.036455,
0.0332033,0.030322,0.0277596,0.0254732,0.0234266,0.0215892, 0.0332033,0.030322,0.0277596,0.0254732,0.0234266,0.0215892,
0.0199351,0.018442,0.0170909,0.0158654,0.0147514,0.0137365, 0.0199351,0.018442,0.0170909,0.0158654,0.0147514,0.0137365,
0.0128101,0.0119627,0.0111863] 0.0128101,0.0119627,0.0111863]
# Rescale to nm # Rescale to nm
for i in range (len(d0)): for i in range (len(d0)):
...@@ -196,7 +196,7 @@ Amber Equivalent: igb = 1 ...@@ -196,7 +196,7 @@ Amber Equivalent: igb = 1
def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None): def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
custom = CustomGBForce() custom = CustomGBForce()
custom.addPerParticleParameter("q"); custom.addPerParticleParameter("q");
custom.addPerParticleParameter("radius"); custom.addPerParticleParameter("radius");
custom.addPerParticleParameter("scale"); custom.addPerParticleParameter("scale");
...@@ -213,7 +213,7 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -213,7 +213,7 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
custom.addComputedValue("B", "1/(1/or-I);" custom.addComputedValue("B", "1/(1/or-I);"
"or=radius-offset", CustomGBForce.SingleParticle) "or=radius-offset", CustomGBForce.SingleParticle)
custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle) custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
if SA=='ACE': if SA=='ACE':
custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle) custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
...@@ -225,7 +225,7 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -225,7 +225,7 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
return custom return custom
""" """
Amber Equivalents: igb = 2 Amber Equivalents: igb = 2
""" """
def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None): def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None):
...@@ -258,7 +258,7 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -258,7 +258,7 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None):
return custom return custom
""" """
Amber Equivalents: igb = 5 Amber Equivalents: igb = 5
""" """
def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None): def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None):
...@@ -295,27 +295,27 @@ Amber Equivalents: igb = 7 ...@@ -295,27 +295,27 @@ Amber Equivalents: igb = 7
""" """
def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None): def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
""" """
Indexing for tables: Indexing for tables:
input: radius1, radius2 input: radius1, radius2
index = (radius2*200-20)*21 + (radius1*200-20) index = (radius2*200-20)*21 + (radius1*200-20)
output: index of desired value in row-by-row, 1D version of Tables 3 & 4 output: index of desired value in row-by-row, 1D version of Tables 3 & 4
""" """
custom = CustomGBForce() custom = CustomGBForce()
custom.addPerParticleParameter("q"); custom.addPerParticleParameter("q");
custom.addPerParticleParameter("radius"); custom.addPerParticleParameter("radius");
custom.addPerParticleParameter("scale"); custom.addPerParticleParameter("scale");
custom.addGlobalParameter("solventDielectric", solventDielectric); custom.addGlobalParameter("solventDielectric", solventDielectric);
custom.addGlobalParameter("soluteDielectric", soluteDielectric); custom.addGlobalParameter("soluteDielectric", soluteDielectric);
custom.addGlobalParameter("offset", 0.009) custom.addGlobalParameter("offset", 0.009)
custom.addGlobalParameter("neckScale", 0.361825) custom.addGlobalParameter("neckScale", 0.361825)
custom.addGlobalParameter("neckCut", 0.68) custom.addGlobalParameter("neckCut", 0.68)
custom.addFunction("getd0", d0, 0, 440) custom.addFunction("getd0", d0, 0, 440)
custom.addFunction("getm0", m0, 0, 440) custom.addFunction("getm0", m0, 0, 440)
...@@ -328,10 +328,10 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -328,10 +328,10 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
"D=abs(r-sr2);" "D=abs(r-sr2);"
"sr2 = scale2*or2;" "sr2 = scale2*or2;"
"or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions) "or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions)
custom.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);" custom.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);"
"psi=I*or; or=radius-offset", CustomGBForce.SingleParticle) "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle) custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
if SA=='ACE': if SA=='ACE':
custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle) custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
...@@ -340,4 +340,4 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None): ...@@ -340,4 +340,4 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;" custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
"f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions) "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
return custom return custom
\ No newline at end of file
wrappers/python/simtk/openmm/app/internal/pdbstructure.py
View file @ 74415dd9
...@@ -12,7 +12,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -12,7 +12,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Christopher M. Bruns Authors: Christopher M. Bruns
Contributors: Peter Eastman Contributors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -260,7 +260,7 @@ class PdbStructure(object): ...@@ -260,7 +260,7 @@ class PdbStructure(object):
"""Establish first and last residues, atoms, etc.""" """Establish first and last residues, atoms, etc."""
for model in self.models: for model in self.models:
model._finalize() model._finalize()
def get_unit_cell_dimensions(self): def get_unit_cell_dimensions(self):
"""Get the dimensions of the crystallographic unit cell (may be None).""" """Get the dimensions of the crystallographic unit cell (may be None)."""
return self._unit_cell_dimensions return self._unit_cell_dimensions
...@@ -547,12 +547,12 @@ class Residue(object): ...@@ -547,12 +547,12 @@ class Residue(object):
... ...
>>> for atom in res: >>> for atom in res:
... print atom ... print atom
ATOM 188 N CYS A 42 40.714 -5.292 12.123 1.00 11.29 N ATOM 188 N CYS A 42 40.714 -5.292 12.123 1.00 11.29 N
ATOM 189 CA CYS A 42 39.736 -5.883 12.911 1.00 10.01 C ATOM 189 CA CYS A 42 39.736 -5.883 12.911 1.00 10.01 C
ATOM 190 C CYS A 42 40.339 -6.654 14.087 1.00 22.28 C ATOM 190 C CYS A 42 40.339 -6.654 14.087 1.00 22.28 C
ATOM 191 O CYS A 42 41.181 -7.530 13.859 1.00 13.70 O ATOM 191 O CYS A 42 41.181 -7.530 13.859 1.00 13.70 O
ATOM 192 CB CYS A 42 38.949 -6.825 12.002 1.00 9.67 C ATOM 192 CB CYS A 42 38.949 -6.825 12.002 1.00 9.67 C
ATOM 193 SG CYS A 42 37.557 -7.514 12.922 1.00 20.12 S ATOM 193 SG CYS A 42 37.557 -7.514 12.922 1.00 20.12 S
""" """
for atom in self.iter_atoms(): for atom in self.iter_atoms():
yield atom yield atom
......
wrappers/python/simtk/openmm/app/pdbfile.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -48,28 +48,28 @@ except: ...@@ -48,28 +48,28 @@ except:
class PDBFile(object): class PDBFile(object):
"""PDBFile parses a Protein Data Bank (PDB) file and constructs a Topology and a set of atom positions from it. """PDBFile parses a Protein Data Bank (PDB) file and constructs a Topology and a set of atom positions from it.
This class also provides methods for creating PDB files. To write a file containing a single model, call This class also provides methods for creating PDB files. To write a file containing a single model, call
writeFile(). You also can create files that contain multiple models. To do this, first call writeHeader(), writeFile(). You also can create files that contain multiple models. To do this, first call writeHeader(),
then writeModel() once for each model in the file, and finally writeFooter() to complete the file.""" then writeModel() once for each model in the file, and finally writeFooter() to complete the file."""
_residueNameReplacements = {} _residueNameReplacements = {}
_atomNameReplacements = {} _atomNameReplacements = {}
def __init__(self, file): def __init__(self, file):
"""Load a PDB file. """Load a PDB file.
The atom positions and Topology can be retrieved by calling getPositions() and getTopology(). The atom positions and Topology can be retrieved by calling getPositions() and getTopology().
Parameters: Parameters:
- file (string) the name of the file to load - file (string) the name of the file to load
""" """
top = Topology() top = Topology()
## The Topology read from the PDB file ## The Topology read from the PDB file
self.topology = top self.topology = top
# Load the PDB file # Load the PDB file
pdb = PdbStructure(open(file), load_all_models=True) pdb = PdbStructure(open(file), load_all_models=True)
PDBFile._loadNameReplacementTables() PDBFile._loadNameReplacementTables()
...@@ -95,7 +95,7 @@ class PDBFile(object): ...@@ -95,7 +95,7 @@ class PDBFile(object):
element = atom.element element = atom.element
if element is None: if element is None:
# Try to guess the element. # Try to guess the element.
upper = atomName.upper() upper = atomName.upper()
if upper.startswith('CL'): if upper.startswith('CL'):
element = elem.chlorine element = elem.chlorine
...@@ -133,9 +133,9 @@ class PDBFile(object): ...@@ -133,9 +133,9 @@ class PDBFile(object):
self.topology.createStandardBonds() self.topology.createStandardBonds()
self.topology.createDisulfideBonds(self.positions) self.topology.createDisulfideBonds(self.positions)
self._numpyPositions = None self._numpyPositions = None
# Add bonds based on CONECT records. # Add bonds based on CONECT records.
connectBonds = [] connectBonds = []
for connect in pdb.models[0].connects: for connect in pdb.models[0].connects:
i = connect[0] i = connect[0]
...@@ -152,14 +152,14 @@ class PDBFile(object): ...@@ -152,14 +152,14 @@ class PDBFile(object):
def getTopology(self): def getTopology(self):
"""Get the Topology of the model.""" """Get the Topology of the model."""
return self.topology return self.topology
def getNumFrames(self): def getNumFrames(self):
"""Get the number of frames stored in the file.""" """Get the number of frames stored in the file."""
return len(self._positions) return len(self._positions)
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 - 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 - frame (int=0) the index of the frame for which to get positions
...@@ -213,11 +213,11 @@ class PDBFile(object): ...@@ -213,11 +213,11 @@ class PDBFile(object):
name = atom.attrib['name'] name = atom.attrib['name']
for id in atom.attrib: for id in atom.attrib:
map[atom.attrib[id]] = name map[atom.attrib[id]] = name
@staticmethod @staticmethod
def writeFile(topology, positions, file=sys.stdout, modelIndex=None): def writeFile(topology, positions, file=sys.stdout, modelIndex=None):
"""Write a PDB file containing a single model. """Write a PDB file containing a single model.
Parameters: Parameters:
- topology (Topology) The Topology defining the model to write - topology (Topology) The Topology defining the model to write
- positions (list) The list of atomic positions to write - positions (list) The list of atomic positions to write
...@@ -230,7 +230,7 @@ class PDBFile(object): ...@@ -230,7 +230,7 @@ class PDBFile(object):
@staticmethod @staticmethod
def writeHeader(topology, file=sys.stdout): def writeHeader(topology, file=sys.stdout):
"""Write out the header for a PDB file. """Write out the header for a PDB file.
Parameters: Parameters:
- topology (Topology) The Topology defining the molecular system being written - topology (Topology) The Topology defining the molecular system being written
- file (file=stdout) A file to write the file to - file (file=stdout) A file to write the file to
...@@ -239,11 +239,11 @@ class PDBFile(object): ...@@ -239,11 +239,11 @@ class PDBFile(object):
if boxSize is not None: if boxSize is not None:
size = boxSize.value_in_unit(angstroms) size = boxSize.value_in_unit(angstroms)
print >>file, "CRYST1%9.3f%9.3f%9.3f 90.00 90.00 90.00 P 1 1 " % size print >>file, "CRYST1%9.3f%9.3f%9.3f 90.00 90.00 90.00 P 1 1 " % size
@staticmethod @staticmethod
def writeModel(topology, positions, file=sys.stdout, modelIndex=None): def writeModel(topology, positions, file=sys.stdout, modelIndex=None):
"""Write out a model to a PDB file. """Write out a model to a PDB file.
Parameters: Parameters:
- topology (Topology) The Topology defining the model to write - topology (Topology) The Topology defining the model to write
- positions (list) The list of atomic positions to write - positions (list) The list of atomic positions to write
...@@ -251,7 +251,7 @@ class PDBFile(object): ...@@ -251,7 +251,7 @@ class PDBFile(object):
- modelIndex (int=None) If not None, the model will be surrounded by MODEL/ENDMDL records with this index - modelIndex (int=None) If not None, the model will be surrounded by MODEL/ENDMDL records with this index
""" """
if len(list(topology.atoms())) != len(positions): if len(list(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')
if is_quantity(positions): if is_quantity(positions):
positions = positions.value_in_unit(angstroms) positions = positions.value_in_unit(angstroms)
if any(math.isnan(norm(pos)) for pos in positions): if any(math.isnan(norm(pos)) for pos in positions):
...@@ -295,7 +295,7 @@ class PDBFile(object): ...@@ -295,7 +295,7 @@ class PDBFile(object):
@staticmethod @staticmethod
def writeFooter(topology, file=sys.stdout): def writeFooter(topology, file=sys.stdout):
"""Write out the footer for a PDB file. """Write out the footer for a PDB file.
Parameters: Parameters:
- topology (Topology) The Topology defining the molecular system being written - topology (Topology) The Topology defining the molecular system being written
- file (file=stdout) A file to write the file to - file (file=stdout) A file to write the file to
......
wrappers/python/simtk/openmm/app/pdbreporter.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -33,16 +33,16 @@ __version__ = "1.0" ...@@ -33,16 +33,16 @@ __version__ = "1.0"
import simtk.openmm as mm import simtk.openmm as mm
from simtk.openmm.app import PDBFile from simtk.openmm.app import PDBFile
class PDBReporter(object): class PDBReporter(object):
"""PDBReporter outputs a series of frames from a Simulation to a PDB file. """PDBReporter outputs a series of frames from a Simulation to a PDB file.
To use it, create a PDBReporter, then add it to the Simulation's list of reporters. To use it, create a PDBReporter, then add it to the Simulation's list of reporters.
""" """
def __init__(self, file, reportInterval): def __init__(self, file, reportInterval):
"""Create a PDBReporter. """Create a PDBReporter.
Parameters: Parameters:
- file (string) The file to write to - file (string) The file to write to
- reportInterval (int) The interval (in time steps) at which to write frames - reportInterval (int) The interval (in time steps) at which to write frames
...@@ -51,10 +51,10 @@ class PDBReporter(object): ...@@ -51,10 +51,10 @@ class PDBReporter(object):
self._out = open(file, 'w') self._out = open(file, 'w')
self._topology = None self._topology = None
self._nextModel = 0 self._nextModel = 0
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) The Simulation to generate a report for
Returns: A five element tuple. The first element is the number of steps until the Returns: A five element tuple. The first element is the number of steps until the
...@@ -63,10 +63,10 @@ class PDBReporter(object): ...@@ -63,10 +63,10 @@ class PDBReporter(object):
""" """
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)
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 - simulation (Simulation) The Simulation to generate a report for
- state (State) The current state of the simulation - state (State) The current state of the simulation
...@@ -77,7 +77,7 @@ class PDBReporter(object): ...@@ -77,7 +77,7 @@ class PDBReporter(object):
self._nextModel += 1 self._nextModel += 1
PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel) PDBFile.writeModel(simulation.topology, state.getPositions(), self._out, self._nextModel)
self._nextModel += 1 self._nextModel += 1
def __del__(self): def __del__(self):
PDBFile.writeFooter(self._topology, self._out) PDBFile.writeFooter(self._topology, self._out)
self._out.close() self._out.close()
wrappers/python/simtk/openmm/app/simulation.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -33,25 +33,25 @@ __version__ = "1.0" ...@@ -33,25 +33,25 @@ __version__ = "1.0"
import simtk.openmm as mm import simtk.openmm as mm
import simtk.unit as unit import simtk.unit as unit
class Simulation(object): class Simulation(object):
"""Simulation provides a simplified API for running simulations with OpenMM and reporting results. """Simulation provides a simplified API for running simulations with OpenMM and reporting results.
A Simulation ties together various objects used for running a simulation: a Topology, System, A Simulation ties together various objects used for running a simulation: a Topology, System,
Integrator, and Context. To use it, you provide the Topology, System, and Integrator, and it Integrator, and Context. To use it, you provide the Topology, System, and Integrator, and it
creates the Context automatically. creates the Context automatically.
Simulation also maintains a list of "reporter" objects that record or analyze data as the simulation Simulation also maintains a list of "reporter" objects that record or analyze data as the simulation
runs, such as writing coordinates to files or displaying structures on the screen. For example, runs, such as writing coordinates to files or displaying structures on the screen. For example,
the following line will cause a file called "output.pdb" to be created, and a structure written to the following line will cause a file called "output.pdb" to be created, and a structure written to
it every 1000 time steps: it every 1000 time steps:
simulation.reporters.append(PDBReporter('output.pdb', 1000)) simulation.reporters.append(PDBReporter('output.pdb', 1000))
""" """
def __init__(self, topology, system, integrator, platform=None, platformProperties=None): def __init__(self, topology, system, integrator, platform=None, platformProperties=None):
"""Create a Simulation. """Create a Simulation.
Parameters: Parameters:
- topology (Topology) A Topology describing the the system to simulate - topology (Topology) A Topology describing the the system to simulate
- system (System) The OpenMM System object to simulate - system (System) The OpenMM System object to simulate
...@@ -77,17 +77,17 @@ class Simulation(object): ...@@ -77,17 +77,17 @@ class Simulation(object):
self.context = mm.Context(system, integrator, platform) self.context = mm.Context(system, integrator, platform)
else: else:
self.context = mm.Context(system, integrator, platform, platformProperties) self.context = mm.Context(system, integrator, platform, platformProperties)
def minimizeEnergy(self, tolerance=1*unit.kilojoule/unit.mole, maxIterations=0): def minimizeEnergy(self, tolerance=1*unit.kilojoule/unit.mole, maxIterations=0):
"""Perform a local energy minimization on the system. """Perform a local energy minimization on the system.
Parameters: Parameters:
- tolerance (energy=1*kilojoule/mole) The energy tolerance to which the system should be minimized - tolerance (energy=1*kilojoule/mole) The energy tolerance to which the system should be minimized
- maxIterations (int=0) The maximum number of iterations to perform. If this is 0, minimization is continued - maxIterations (int=0) The maximum number of iterations to perform. If this is 0, minimization is continued
until the results converge without regard to how many iterations it takes. until the results converge without regard to how many iterations it takes.
""" """
mm.LocalEnergyMinimizer.minimize(self.context, tolerance, maxIterations) mm.LocalEnergyMinimizer.minimize(self.context, tolerance, maxIterations)
def step(self, steps): def step(self, steps):
"""Advance the simulation by integrating a specified number of time steps.""" """Advance the simulation by integrating a specified number of time steps."""
stepTo = self.currentStep+steps stepTo = self.currentStep+steps
......
wrappers/python/simtk/openmm/app/statedatareporter.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
Contributors: Robert McGibbon Contributors: Robert McGibbon
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -35,7 +35,7 @@ import simtk.openmm as mm ...@@ -35,7 +35,7 @@ import simtk.openmm as mm
import simtk.unit as unit import simtk.unit as unit
from simtk.openmm.app import PDBFile from simtk.openmm.app import PDBFile
import math import math
class StateDataReporter(object): class StateDataReporter(object):
"""StateDataReporter outputs information about a simulation, such as energy and temperature, to a file. """StateDataReporter outputs information about a simulation, such as energy and temperature, to a file.
......
wrappers/python/simtk/openmm/app/topology.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -37,18 +37,18 @@ from simtk.unit import nanometers, sqrt ...@@ -37,18 +37,18 @@ from simtk.unit import nanometers, sqrt
class Topology(object): class Topology(object):
"""Topology stores the topological information about a system. """Topology stores the topological information about a system.
The structure of a Topology object is similar to that of a PDB file. It consists of a set of Chains The structure of a Topology object is similar to that of a PDB file. It consists of a set of Chains
(often but not always corresponding to polymer chains). Each Chain contains a set of Residues, (often but not always corresponding to polymer chains). Each Chain contains a set of Residues,
and each Residue contains a set of Atoms. In addition, the Topology stores a list of which atom and each Residue contains a set of Atoms. In addition, the Topology stores a list of which atom
pairs are bonded to each other, and the dimensions of the crystallographic unit cell. pairs are bonded to each other, and the dimensions of the crystallographic unit cell.
Atom and residue names should follow the PDB 3.0 nomenclature for all molecules for which one exists. Atom and residue names should follow the PDB 3.0 nomenclature for all molecules for which one exists.
""" """
_standardBonds = {} _standardBonds = {}
_hasLoadedStandardBonds = False _hasLoadedStandardBonds = False
def __init__(self): def __init__(self):
"""Create a new Topology object""" """Create a new Topology object"""
self._chains = [] self._chains = []
...@@ -56,19 +56,19 @@ class Topology(object): ...@@ -56,19 +56,19 @@ class Topology(object):
self._numAtoms = 0 self._numAtoms = 0
self._bonds = [] self._bonds = []
self._unitCellDimensions = None self._unitCellDimensions = None
def addChain(self): def addChain(self):
"""Create a new Chain and add it to the Topology. """Create a new Chain and add it to the Topology.
Returns: the newly created Chain Returns: the newly created Chain
""" """
chain = Chain(len(self._chains), self) chain = Chain(len(self._chains), self)
self._chains.append(chain) self._chains.append(chain)
return chain return chain
def addResidue(self, name, chain): def addResidue(self, name, chain):
"""Create a new Residue and add it to the Topology. """Create a new Residue and add it to the Topology.
Parameters: Parameters:
- name (string) The name of the residue to add - name (string) The name of the residue to add
- chain (Chain) The Chain to add it to - chain (Chain) The Chain to add it to
...@@ -78,10 +78,10 @@ class Topology(object): ...@@ -78,10 +78,10 @@ class Topology(object):
self._numResidues += 1 self._numResidues += 1
chain._residues.append(residue) chain._residues.append(residue)
return residue return residue
def addAtom(self, name, element, residue): def addAtom(self, name, element, residue):
"""Create a new Atom and add it to the Topology. """Create a new Atom and add it to the Topology.
Parameters: Parameters:
- name (string) The name of the atom to add - name (string) The name of the atom to add
- element (Element) The element of the atom to add - element (Element) The element of the atom to add
...@@ -92,52 +92,52 @@ class Topology(object): ...@@ -92,52 +92,52 @@ class Topology(object):
self._numAtoms += 1 self._numAtoms += 1
residue._atoms.append(atom) residue._atoms.append(atom)
return atom return atom
def addBond(self, atom1, atom2): def addBond(self, atom1, atom2):
"""Create a new bond and add it to the Topology. """Create a new bond and add it to the Topology.
Parameters: Parameters:
- atom1 (Atom) The first Atom connected by the bond - atom1 (Atom) The first Atom connected by the bond
- atom2 (Atom) The second Atom connected by the bond - atom2 (Atom) The second Atom connected by the bond
""" """
self._bonds.append((atom1, atom2)) self._bonds.append((atom1, atom2))
def chains(self): def chains(self):
"""Iterate over all Chains in the Topology.""" """Iterate over all Chains in the Topology."""
return iter(self._chains) return iter(self._chains)
def residues(self): def residues(self):
"""Iterate over all Residues in the Topology.""" """Iterate over all Residues in the Topology."""
for chain in self._chains: for chain in self._chains:
for residue in chain._residues: for residue in chain._residues:
yield residue yield residue
def atoms(self): def atoms(self):
"""Iterate over all Atoms in the Topology.""" """Iterate over all Atoms in the Topology."""
for chain in self._chains: for chain in self._chains:
for residue in chain._residues: for residue in chain._residues:
for atom in residue._atoms: for atom in residue._atoms:
yield atom yield atom
def bonds(self): def bonds(self):
"""Iterate over all bonds (each represented as a tuple of two Atoms) in the Topology.""" """Iterate over all bonds (each represented as a tuple of two Atoms) in the Topology."""
return iter(self._bonds) return iter(self._bonds)
def getUnitCellDimensions(self): def getUnitCellDimensions(self):
"""Get the dimensions of the crystallographic unit cell. """Get the dimensions of the crystallographic unit cell.
The return value may be None if this Topology does not represent a periodic structure. The return value may be None if this Topology does not represent a periodic structure.
""" """
return self._unitCellDimensions return self._unitCellDimensions
def setUnitCellDimensions(self, dimensions): def setUnitCellDimensions(self, dimensions):
"""Set the dimensions of the crystallographic unit cell.""" """Set the dimensions of the crystallographic unit cell."""
self._unitCellDimensions = dimensions self._unitCellDimensions = dimensions
@staticmethod @staticmethod
def loadBondDefinitions(file): def loadBondDefinitions(file):
"""Load an XML file containing definitions of bonds that should be used by createStandardBonds(). """Load an XML file containing definitions of bonds that should be used by createStandardBonds().
The built in residues.xml file containing definitions for standard amino acids and nucleotides is loaded automatically. The built in residues.xml file containing definitions for standard amino acids and nucleotides is loaded automatically.
This method can be used to load additional definitions for other residue types. They will then be used in subsequent This method can be used to load additional definitions for other residue types. They will then be used in subsequent
calls to createStandardBonds(). This is a static method, so it affects subsequent calls on all Topology objects. calls to createStandardBonds(). This is a static method, so it affects subsequent calls on all Topology objects.
...@@ -149,31 +149,31 @@ class Topology(object): ...@@ -149,31 +149,31 @@ class Topology(object):
bonds = [] bonds = []
Topology._standardBonds[residue.attrib['name']] = bonds Topology._standardBonds[residue.attrib['name']] = bonds
for bond in residue.findall('Bond'): for bond in residue.findall('Bond'):
bonds.append((bond.attrib['from'], bond.attrib['to'])) bonds.append((bond.attrib['from'], bond.attrib['to']))
def createStandardBonds(self): def createStandardBonds(self):
"""Create bonds based on the atom and residue names for all standard residue types. """Create bonds based on the atom and residue names for all standard residue types.
Definitions for standard amino acids and nucleotides are built in. You can call loadBondDefinitions() to load Definitions for standard amino acids and nucleotides are built in. You can call loadBondDefinitions() to load
additional definitions for other residue types. additional definitions for other residue types.
""" """
if not Topology._hasLoadedStandardBonds: if not Topology._hasLoadedStandardBonds:
# Load the standard bond definitions. # Load the standard bond definitions.
Topology.loadBondDefinitions(os.path.join(os.path.dirname(__file__), 'data', 'residues.xml')) Topology.loadBondDefinitions(os.path.join(os.path.dirname(__file__), 'data', 'residues.xml'))
Topology._hasLoadedStandardBonds = True Topology._hasLoadedStandardBonds = True
for chain in self._chains: for chain in self._chains:
# First build a map of atom names to atoms. # First build a map of atom names to atoms.
atomMaps = [] atomMaps = []
for residue in chain._residues: for residue in chain._residues:
atomMap = {} atomMap = {}
atomMaps.append(atomMap) atomMaps.append(atomMap)
for atom in residue._atoms: for atom in residue._atoms:
atomMap[atom.name] = atom atomMap[atom.name] = atom
# Loop over residues and construct bonds. # Loop over residues and construct bonds.
for i in range(len(chain._residues)): for i in range(len(chain._residues)):
name = chain._residues[i].name name = chain._residues[i].name
if name in Topology._standardBonds: if name in Topology._standardBonds:
...@@ -198,17 +198,17 @@ class Topology(object): ...@@ -198,17 +198,17 @@ class Topology(object):
toAtom = bond[1] toAtom = bond[1]
if fromAtom in atomMaps[fromResidue] and toAtom in atomMaps[toResidue]: if fromAtom in atomMaps[fromResidue] and toAtom in atomMaps[toResidue]:
self.addBond(atomMaps[fromResidue][fromAtom], atomMaps[toResidue][toAtom]) self.addBond(atomMaps[fromResidue][fromAtom], atomMaps[toResidue][toAtom])
def createDisulfideBonds(self, positions): def createDisulfideBonds(self, positions):
"""Identify disulfide bonds based on proximity and add them to the Topology. """Identify disulfide bonds based on proximity and add them to the Topology.
Parameters: Parameters:
- positions (list) The list of atomic positions based on which to identify bonded atoms - positions (list) The list of atomic positions based on which to identify bonded atoms
""" """
def isCyx(res): def isCyx(res):
names = [atom.name for atom in res._atoms] names = [atom.name for atom in res._atoms]
return 'SG' in names and 'HG' not in names return 'SG' in names and 'HG' not in names
cyx = [res for res in self.residues() if res.name == 'CYS' and isCyx(res)] cyx = [res for res in self.residues() if res.name == 'CYS' and isCyx(res)]
atomNames = [[atom.name for atom in res._atoms] for res in cyx] atomNames = [[atom.name for atom in res._atoms] for res in cyx]
for i in range(len(cyx)): for i in range(len(cyx)):
...@@ -231,7 +231,7 @@ class Chain(object): ...@@ -231,7 +231,7 @@ class Chain(object):
## The Topology this Chain belongs to ## The Topology this Chain belongs to
self.topology = topology self.topology = topology
self._residues = [] self._residues = []
def residues(self): def residues(self):
"""Iterate over all Residues in the Chain.""" """Iterate over all Residues in the Chain."""
return iter(self._residues) return iter(self._residues)
...@@ -253,14 +253,14 @@ class Residue(object): ...@@ -253,14 +253,14 @@ class Residue(object):
## The Chain this Residue belongs to ## The Chain this Residue belongs to
self.chain = chain self.chain = chain
self._atoms = [] self._atoms = []
def atoms(self): def atoms(self):
"""Iterate over all Atoms in the Residue.""" """Iterate over all Atoms in the Residue."""
return iter(self._atoms) return iter(self._atoms)
class Atom(object): class Atom(object):
"""An Atom object represents a residue within a Topology.""" """An Atom object represents a residue within a Topology."""
def __init__(self, name, element, index, residue): def __init__(self, name, element, index, residue):
"""Construct a new Atom. You should call addAtom() on the Topology instead of calling this directly.""" """Construct a new Atom. You should call addAtom() on the Topology instead of calling this directly."""
## The name of the Atom ## The name of the Atom
......
wrappers/python/simtk/openmm/vec3.py
View file @ 74415dd9
...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -35,43 +35,43 @@ import simtk.unit as unit ...@@ -35,43 +35,43 @@ import simtk.unit as unit
class Vec3(tuple): class Vec3(tuple):
"""Vec3 is a 3-element tuple that supports many math operations.""" """Vec3 is a 3-element tuple that supports many math operations."""
def __new__(cls, x, y, z): def __new__(cls, x, y, z):
"""Create a new Vec3.""" """Create a new Vec3."""
return tuple.__new__(cls, (x, y, z)) return tuple.__new__(cls, (x, y, z))
def __add__(self, other): def __add__(self, other):
"""Add two Vec3s.""" """Add two Vec3s."""
return Vec3(self[0]+other[0], self[1]+other[1], self[2]+other[2]) return Vec3(self[0]+other[0], self[1]+other[1], self[2]+other[2])
def __radd__(self, other): def __radd__(self, other):
"""Add two Vec3s.""" """Add two Vec3s."""
return Vec3(self[0]+other[0], self[1]+other[1], self[2]+other[2]) return Vec3(self[0]+other[0], self[1]+other[1], self[2]+other[2])
def __sub__(self, other): def __sub__(self, other):
"""Add two Vec3s.""" """Add two Vec3s."""
return Vec3(self[0]-other[0], self[1]-other[1], self[2]-other[2]) return Vec3(self[0]-other[0], self[1]-other[1], self[2]-other[2])
def __rsub__(self, other): def __rsub__(self, other):
"""Add two Vec3s.""" """Add two Vec3s."""
return Vec3(other[0]-self[0], other[1]-self[1], other[2]-self[2]) return Vec3(other[0]-self[0], other[1]-self[1], other[2]-self[2])
def __mul__(self, other): def __mul__(self, other):
"""Multiply a Vec3 by a constant.""" """Multiply a Vec3 by a constant."""
if unit.is_unit(other): if unit.is_unit(other):
return unit.Quantity(self, other) return unit.Quantity(self, other)
return Vec3(other*self[0], other*self[1], other*self[2]) return Vec3(other*self[0], other*self[1], other*self[2])
def __rmul__(self, other): def __rmul__(self, other):
"""Multiply a Vec3 by a constant.""" """Multiply a Vec3 by a constant."""
if unit.is_unit(other): if unit.is_unit(other):
return unit.Quantity(self, other) return unit.Quantity(self, other)
return Vec3(other*self[0], other*self[1], other*self[2]) return Vec3(other*self[0], other*self[1], other*self[2])
def __div__(self, other): def __div__(self, other):
"""Divide a Vec3 by a constant.""" """Divide a Vec3 by a constant."""
return Vec3(self[0]/other, self[1]/other, self[2]/other) return Vec3(self[0]/other, self[1]/other, self[2]/other)
__truediv__ = __div__ __truediv__ = __div__
def __deepcopy__(self, memo): def __deepcopy__(self, memo):
return Vec3(self[0], self[1], self[2]) return Vec3(self[0], self[1], self[2])
wrappers/python/simtk/unit/basedimension.py
View file @ 74415dd9
...@@ -16,7 +16,7 @@ Portions copyright (c) 2012 Stanford University and the Authors. ...@@ -16,7 +16,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Christopher M. Bruns Authors: Christopher M. Bruns
Contributors: Peter Eastman Contributors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense, the rights to use, copy, modify, merge, publish, distribute, sublicense,
...@@ -42,7 +42,7 @@ __version__ = "0.6" ...@@ -42,7 +42,7 @@ __version__ = "0.6"
class BaseDimension(object): class BaseDimension(object):
''' '''
A physical dimension such as length, mass, or temperature. A physical dimension such as length, mass, or temperature.
It is unlikely the user will need to create new ones. It is unlikely the user will need to create new ones.
''' '''
# Keep deterministic order of dimensions # Keep deterministic order of dimensions
...@@ -57,7 +57,7 @@ class BaseDimension(object): ...@@ -57,7 +57,7 @@ class BaseDimension(object):
'angle': 8, 'angle': 8,
} }
_next_unused_index = 9 _next_unused_index = 9
def __init__(self, name): def __init__(self, name):
"""Create a new BaseDimension. """Create a new BaseDimension.
...@@ -70,17 +70,17 @@ class BaseDimension(object): ...@@ -70,17 +70,17 @@ class BaseDimension(object):
BaseDimension._index_by_name[name] = BaseDimension._next_unused_index BaseDimension._index_by_name[name] = BaseDimension._next_unused_index
BaseDimension._next_unused_index += 1 BaseDimension._next_unused_index += 1
self._index = BaseDimension._index_by_name[name] self._index = BaseDimension._index_by_name[name]
def __lt__(self, other): def __lt__(self, other):
""" """
The implicit order of BaseDimensions is the order in which they were created. The implicit order of BaseDimensions is the order in which they were created.
This method is used for using BaseDimensions as hash keys, and also affects This method is used for using BaseDimensions as hash keys, and also affects
the order in which units appear in multi-dimensional Quantities. the order in which units appear in multi-dimensional Quantities.
Returns True if self < other, False otherwise. Returns True if self < other, False otherwise.
""" """
return self._index < other._index return self._index < other._index
def __hash__(self): def __hash__(self):
""" """
Needed for using BaseDimensions as hash keys. Needed for using BaseDimensions as hash keys.
......
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