Skip to content

GitLab

Commit 1c938ceb authored by Jason Swails's avatar Jason Swails
Browse files

Merge branch 'master' into amber-switching 

 Conflicts:
	wrappers/python/simtk/openmm/app/amberprmtopfile.py

In fixing the merge conflict, I went ahead and fixed up the switchDistance logic
to match what I did in CharmmPsfFile.
parents a1113e7b 167ae8a0
Hide whitespace changes
Inline Side-by-side
Showing with 1057 additions and 490 deletions
wrappers/python/simtk/openmm/app/pdbreporter.py
View file @ 1c938ceb
......@@ -44,9 +44,12 @@ class PDBReporter(object):
def __init__(self, file, reportInterval):
"""Create a PDBReporter.
Parameters:
- file (string) The file to write to
- reportInterval (int) The interval (in time steps) at which to write frames
Parameters
----------
file : string
The file to write to
reportInterval : int
The interval (in time steps) at which to write frames
"""
self._reportInterval = reportInterval
self._out = open(file, 'w')
......@@ -56,11 +59,18 @@ class PDBReporter(object):
def describeNextReport(self, simulation):
"""Get information about the next report this object will generate.
Parameters:
- simulation (Simulation) The Simulation to generate a report for
Returns: A five element tuple. The first element is the number of steps until the
next report. The remaining elements specify whether that report will require
positions, velocities, forces, and energies respectively.
Parameters
----------
simulation : Simulation
The Simulation to generate a report for
Returns
-------
tuple
A five element tuple. The first element is the number of steps
until the next report. The remaining elements specify whether
that report will require positions, velocities, forces, and
energies respectively.
"""
steps = self._reportInterval - simulation.currentStep%self._reportInterval
return (steps, True, False, False, False)
......@@ -68,9 +78,12 @@ class PDBReporter(object):
def report(self, simulation, state):
"""Generate a report.
Parameters:
- simulation (Simulation) The Simulation to generate a report for
- state (State) The current state of the simulation
Parameters
----------
simulation : Simulation
The Simulation to generate a report for
state : State
The current state of the simulation
"""
if self._nextModel == 0:
PDBFile.writeHeader(simulation.topology, self._out)
......@@ -95,9 +108,12 @@ class PDBxReporter(PDBReporter):
def report(self, simulation, state):
"""Generate a report.
Parameters:
- simulation (Simulation) The Simulation to generate a report for
- state (State) The current state of the simulation
Parameters
----------
simulation : Simulation
The Simulation to generate a report for
state : State
The current state of the simulation
"""
if self._nextModel == 0:
PDBxFile.writeHeader(simulation.topology, self._out)
......
wrappers/python/simtk/openmm/app/pdbxfile.py
View file @ 1c938ceb
......@@ -53,10 +53,14 @@ class PDBxFile(object):
def __init__(self, file):
"""Load a PDBx/mmCIF file.
The atom positions and Topology can be retrieved by calling getPositions() and getTopology().
Parameters:
- file (string) the name of the file to load. Alternatively you can pass an open file object.
The atom positions and Topology can be retrieved by calling
getPositions() and getTopology().
Parameters
----------
file : string
the name of the file to load. Alternatively you can pass an open
file object.
"""
top = Topology()
## The Topology read from the PDBx/mmCIF file
......@@ -193,10 +197,14 @@ class PDBxFile(object):
def getPositions(self, asNumpy=False, frame=0):
"""Get the atomic positions.
Parameters:
- asNumpy (boolean=False) if true, the values are returned as a numpy array instead of a list of Vec3s
- frame (int=0) the index of the frame for which to get positions
"""
Parameters
----------
asNumpy : bool=False
if true, the values are returned as a numpy array instead of a list
of Vec3s
frame : int=0
the index of the frame for which to get positions
"""
if asNumpy:
if self._numpyPositions is None:
self._numpyPositions = [None]*len(self._positions)
......@@ -210,14 +218,21 @@ class PDBxFile(object):
entry=None):
"""Write a PDBx/mmCIF file containing a single model.
Parameters:
- topology (Topology) The Topology defining the model to write
- positions (list) The list of atomic positions to write
- file (file=stdout) A file to write to
- keepIds (bool=False) If True, keep the residue and chain IDs specified in the Topology rather than generating
new ones. Warning: It is up to the caller to make sure these are valid IDs that satisfy the requirements of
the PDBx/mmCIF format. Otherwise, the output file will be invalid.
- entry (str=None) The entry ID to assign to the CIF file
Parameters
----------
topology : Topology
The Topology defining the model to write
positions : list
The list of atomic positions to write
file : file=stdout
A file to write to
keepIds : bool=False
If True, keep the residue and chain IDs specified in the Topology
rather than generating new ones. Warning: It is up to the caller to
make sure these are valid IDs that satisfy the requirements of the
PDBx/mmCIF format. Otherwise, the output file will be invalid.
entry : str=None
The entry ID to assign to the CIF file
"""
PDBxFile.writeHeader(topology, file, entry)
PDBxFile.writeModel(topology, positions, file, keepIds=keepIds)
......@@ -226,10 +241,14 @@ class PDBxFile(object):
def writeHeader(topology, file=sys.stdout, entry=None):
"""Write out the header for a PDBx/mmCIF file.
Parameters:
- topology (Topology) The Topology defining the molecular system being written
- file (file=stdout) A file to write the file to
- entry (str=None) The entry ID to assign to the CIF file
Parameters
----------
topology : Topology
The Topology defining the molecular system being written
file : file=stdout
A file to write the file to
entry : str=None
The entry ID to assign to the CIF file
"""
if entry is not None:
print('data_%s' % entry, file=file)
......@@ -280,14 +299,21 @@ class PDBxFile(object):
def writeModel(topology, positions, file=sys.stdout, modelIndex=1, keepIds=False):
"""Write out a model to a PDBx/mmCIF file.
Parameters:
- topology (Topology) The Topology defining the model to write
- positions (list) The list of atomic positions to write
- file (file=stdout) A file to write the model to
- modelIndex (int=1) The model number of this frame
- keepIds (bool=False) If True, keep the residue and chain IDs specified in the Topology rather than generating
new ones. Warning: It is up to the caller to make sure these are valid IDs that satisfy the requirements of
the PDBx/mmCIF format. Otherwise, the output file will be invalid.
Parameters
----------
topology : Topology
The Topology defining the model to write
positions : list
The list of atomic positions to write
file : file=stdout
A file to write the model to
modelIndex : int=1
The model number of this frame
keepIds : bool=False
If True, keep the residue and chain IDs specified in the Topology
rather than generating new ones. Warning: It is up to the caller to
make sure these are valid IDs that satisfy the requirements of the
PDBx/mmCIF format. Otherwise, the output file will be invalid.
"""
if len(list(topology.atoms())) != len(positions):
raise ValueError('The number of positions must match the number of atoms')
......
wrappers/python/simtk/openmm/app/simulation.py
View file @ 1c938ceb
......@@ -55,13 +55,19 @@ class Simulation(object):
def __init__(self, topology, system, integrator, platform=None, platformProperties=None):
"""Create a Simulation.
Parameters:
- topology (Topology) A Topology describing the the system to simulate
- system (System) The OpenMM System object to simulate
- integrator (Integrator) The OpenMM Integrator to use for simulating the System
- platform (Platform=None) If not None, the OpenMM Platform to use
- platformProperties (map=None) If not None, a set of platform-specific properties to pass
to the Context's constructor
Parameters
----------
topology : Topology
A Topology describing the the system to simulate
system : System
The OpenMM System object to simulate
integrator : Integrator
The OpenMM Integrator to use for simulating the System
platform : Platform=None
If not None, the OpenMM Platform to use
platformProperties : map=None
If not None, a set of platform-specific properties to pass to the
Context's constructor
"""
## The Topology describing the system being simulated
self.topology = topology
......@@ -84,35 +90,48 @@ class Simulation(object):
def minimizeEnergy(self, tolerance=10*unit.kilojoule/unit.mole, maxIterations=0):
"""Perform a local energy minimization on the system.
Parameters:
- tolerance (energy=10*kilojoules/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
until the results converge without regard to how many iterations it takes.
Parameters
----------
tolerance : energy=10*kilojoules/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 until the results converge without regard
to how many iterations it takes.
"""
mm.LocalEnergyMinimizer.minimize(self.context, tolerance, maxIterations)
def step(self, steps):
"""Advance the simulation by integrating a specified number of time steps."""
self._simulate(endStep=self.currentStep+steps)
def runForClockTime(self, time, checkpointFile=None, stateFile=None, checkpointInterval=None):
"""Advance the simulation by integrating time steps until a fixed amount of clock time has elapsed.
This is useful when you have a limited amount of computer time available, and want to run the longest simulation
possible in that time. This method will continue taking time steps until the specified clock time has elapsed,
then return. It also can automatically write out a checkpoint and/or state file before returning, so you can
later resume the simulation. Another option allows it to write checkpoints or states at regular intervals, so
you can resume even if the simulation is interrupted before the time limit is reached.
Parameters:
- time (time) the amount of time to run for. If no units are specified, it is assumed to be a number of hours.
- checkpointFile (string or file=None) if specified, a checkpoint file will be written at the end of the
simulation (and optionally at regular intervals before then) by passing this to saveCheckpoint().
- stateFile (string or file=None) if specified, a state file will be written at the end of the
simulation (and optionally at regular intervals before then) by passing this to saveState().
- checkpointInterval (time=None) if specified, checkpoints and/or states will be written at regular intervals
during the simulation, in addition to writing a final version at the end. If no units are specified, this is
assumed to be in hours.
Parameters
----------
time : time
the amount of time to run for. If no units are specified, it is
assumed to be a number of hours.
checkpointFile : string or file=None
if specified, a checkpoint file will be written at the end of the
simulation (and optionally at regular intervals before then) by
passing this to saveCheckpoint().
stateFile : string or file=None
if specified, a state file will be written at the end of the
simulation (and optionally at regular intervals before then) by
passing this to saveState().
checkpointInterval : time=None
if specified, checkpoints and/or states will be written at regular
intervals during the simulation, in addition to writing a final
version at the end. If no units are specified, this is assumed to
be in hours.
"""
if unit.is_quantity(time):
time = time.value_in_unit(unit.hours)
......@@ -131,7 +150,7 @@ class Simulation(object):
self.saveCheckpoint(checkpointFile)
if stateFile is not None:
self.saveState(stateFile)
def _simulate(self, endStep=None, endTime=None):
if endStep is None:
endStep = sys.maxsize
......@@ -174,30 +193,36 @@ class Simulation(object):
def saveCheckpoint(self, file):
"""Save a checkpoint of the simulation to a file.
The output is a binary file that contains a complete representation of the current state of the Simulation.
It includes both publicly visible data such as the particle positions and velocities, and also internal data
such as the states of random number generators. Reloading the checkpoint will put the Simulation back into
precisely the same state it had before, so it can be exactly continued.
A checkpoint file is highly specific to the Simulation it was created from. It can only be loaded into
another Simulation that has an identical System, uses the same Platform and OpenMM version, and is running on
identical hardware. If you need a more portable way to resume simulations, consider using saveState() instead.
Parameters:
- file (string or file) a File-like object to write the checkpoint to, or alternatively a filename
Parameters
----------
file : string or file
a File-like object to write the checkpoint to, or alternatively a
filename
"""
if isinstance(file, str):
with open(file, 'wb') as f:
f.write(self.context.createCheckpoint())
else:
file.write(self.context.createCheckpoint())
def loadCheckpoint(self, file):
"""Load a checkpoint file that was created with saveCheckpoint().
Parameters:
- file (string or file) a File-like object to load the checkpoint from, or alternatively a filename
Parameters
----------
file : string or file
a File-like object to load the checkpoint from, or alternatively a
filename
"""
if isinstance(file, str):
with open(file, 'rb') as f:
......@@ -207,18 +232,21 @@ class Simulation(object):
def saveState(self, file):
"""Save the current state of the simulation to a file.
The output is an XML file containing a serialized State object. It includes all publicly visible data,
including positions, velocities, and parameters. Reloading the State will put the Simulation back into
approximately the same state it had before.
Unlike saveCheckpoint(), this does not store internal data such as the states of random number generators.
Therefore, you should not expect the following trajectory to be identical to what would have been produced
with the original Simulation. On the other hand, this means it is portable across different Platforms or
hardware.
Parameters:
- file (string or file) a File-like object to write the state to, or alternatively a filename
Parameters
----------
file : string or file
a File-like object to write the state to, or alternatively a
filename
"""
state = self.context.getState(getPositions=True, getVelocities=True, getParameters=True)
xml = mm.XmlSerializer.serialize(state)
......@@ -227,12 +255,15 @@ class Simulation(object):
f.write(xml)
else:
file.write(xml)
def loadState(self, file):
"""Load a State file that was created with saveState().
Parameters:
- file (string or file) a File-like object to load the state from, or alternatively a filename
Parameters
----------
file : string or file
a File-like object to load the state from, or alternatively a
filename
"""
if isinstance(file, str):
with open(file, 'r') as f:
......
wrappers/python/simtk/openmm/app/statedatareporter.py
View file @ 1c938ceb
......@@ -60,31 +60,54 @@ class StateDataReporter(object):
progress=False, remainingTime=False, speed=False, elapsedTime=False, separator=',', systemMass=None, totalSteps=None):
"""Create a StateDataReporter.
Parameters:
- file (string or file) The file to write to, specified as a file name or file object
- reportInterval (int) The interval (in time steps) at which to write frames
- step (boolean=False) Whether to write the current step index to the file
- time (boolean=False) Whether to write the current time to the file
- potentialEnergy (boolean=False) Whether to write the potential energy to the file
- kineticEnergy (boolean=False) Whether to write the kinetic energy to the file
- totalEnergy (boolean=False) Whether to write the total energy to the file
- temperature (boolean=False) Whether to write the instantaneous temperature to the file
- volume (boolean=False) Whether to write the periodic box volume to the file
- density (boolean=False) Whether to write the system density to the file
- progress (boolean=False) Whether to write current progress (percent completion) to the file.
If this is True, you must also specify totalSteps.
- remainingTime (boolean=False) Whether to write an estimate of the remaining clock time until
completion to the file. If this is True, you must also specify totalSteps.
- speed (bool=False) Whether to write an estimate of the simulation speed in ns/day to the file
- elapsedTime (bool=False) Whether to write the elapsed time of the simulation in seconds to the file.
- separator (string=',') The separator to use between columns in the file
- systemMass (mass=None) The total mass to use for the system when reporting density. If this is
None (the default), the system mass is computed by summing the masses of all particles. This
parameter is useful when the particle masses do not reflect their actual physical mass, such as
when some particles have had their masses set to 0 to immobilize them.
- totalSteps (int=None) The total number of steps that will be included in the simulation. This
is required if either progress or remainingTime is set to True, and defines how many steps will
indicate 100% completion.
Parameters
----------
file : string or file
The file to write to, specified as a file name or file object
reportInterval : int
The interval (in time steps) at which to write frames
step : bool=False
Whether to write the current step index to the file
time : bool=False
Whether to write the current time to the file
potentialEnergy : bool=False
Whether to write the potential energy to the file
kineticEnergy : bool=False
Whether to write the kinetic energy to the file
totalEnergy : bool=False
Whether to write the total energy to the file
temperature : bool=False
Whether to write the instantaneous temperature to the file
volume : bool=False
Whether to write the periodic box volume to the file
density : bool=False
Whether to write the system density to the file
progress : bool=False
Whether to write current progress (percent completion) to the file.
If this is True, you must also specify totalSteps.
remainingTime : bool=False
Whether to write an estimate of the remaining clock time until
completion to the file. If this is True, you must also specify
totalSteps.
speed : bool=False
Whether to write an estimate of the simulation speed in ns/day to
the file
elapsedTime : bool=False
Whether to write the elapsed time of the simulation in seconds to
the file.
separator : string=','
The separator to use between columns in the file
systemMass : mass=None
The total mass to use for the system when reporting density. If
this is None (the default), the system mass is computed by summing
the masses of all particles. This parameter is useful when the
particle masses do not reflect their actual physical mass, such as
when some particles have had their masses set to 0 to immobilize
them.
totalSteps : int=None
The total number of steps that will be included in the simulation.
This is required if either progress or remainingTime is set to True,
and defines how many steps will indicate 100% completion.
"""
self._reportInterval = reportInterval
self._openedFile = isinstance(file, str)
......@@ -129,11 +152,18 @@ class StateDataReporter(object):
def describeNextReport(self, simulation):
"""Get information about the next report this object will generate.
Parameters:
- simulation (Simulation) The Simulation to generate a report for
Returns: A five element tuple. The first element is the number of steps until the
next report. The remaining elements specify whether that report will require
positions, velocities, forces, and energies respectively.
Parameters
----------
simulation : Simulation
The Simulation to generate a report for
Returns
-------
tuple
A five element tuple. The first element is the number of steps
until the next report. The remaining elements specify whether
that report will require positions, velocities, forces, and
energies respectively.
"""
steps = self._reportInterval - simulation.currentStep%self._reportInterval
return (steps, self._needsPositions, self._needsVelocities, self._needsForces, self._needEnergy)
......@@ -141,9 +171,12 @@ class StateDataReporter(object):
def report(self, simulation, state):
"""Generate a report.
Parameters:
- simulation (Simulation) The Simulation to generate a report for
- state (State) The current state of the simulation
Parameters
----------
simulation : Simulation
The Simulation to generate a report for
state : State
The current state of the simulation
"""
if not self._hasInitialized:
self._initializeConstants(simulation)
......@@ -174,11 +207,16 @@ class StateDataReporter(object):
def _constructReportValues(self, simulation, state):
"""Query the simulation for the current state of our observables of interest.
Parameters:
- simulation (Simulation) The Simulation to generate a report for
- state (State) The current state of the simulation
Parameters
----------
simulation : Simulation
The Simulation to generate a report for
state : State
The current state of the simulation
Returns: A list of values summarizing the current state of
Returns
-------
A list of values summarizing the current state of
the simulation, to be printed or saved. Each element in the list
corresponds to one of the columns in the resulting CSV file.
"""
......
wrappers/python/simtk/openmm/app/topology.py
View file @ 1c938ceb
......@@ -62,19 +62,43 @@ class Topology(object):
def __repr__(self):
nchains = len(self._chains)
nres = sum(1 for r in self.residues())
natom = sum(1 for a in self.atoms())
nres = self._numResidues
natom = self._numAtoms
nbond = len(self._bonds)
return '<%s; %d chains, %d residues, %d atoms, %d bonds>' % (
type(self).__name__, nchains, nres, natom, nbond)
def getNumAtoms(self):
"""Return the number of atoms in the Topology.
"""
natom = self._numAtoms
return natom
def getNumResidues(self):
"""Return the number of residues in the Topology.
"""
nres = self._numResidues
return nres
def getNumChains(self):
"""Return the number of chains in the Topology.
"""
nchain = len(self._chains)
return nchain
def addChain(self, id=None):
"""Create a new Chain and add it to the Topology.
Parameters:
- id (string=None) An optional identifier for the chain. If this is omitted, an id
is generated based on the chain index.
Returns: the newly created Chain
Parameters
----------
id : string=None
An optional identifier for the chain. If this is omitted, an id is
generated based on the chain index.
Returns
-------
Chain
the newly created Chain
"""
if id is None:
id = str(len(self._chains)+1)
......@@ -85,12 +109,20 @@ class Topology(object):
def addResidue(self, name, chain, id=None):
"""Create a new Residue and add it to the Topology.
Parameters:
- name (string) The name of the residue to add
- chain (Chain) The Chain to add it to
- id (string=None) An optional identifier for the residue. If this is omitted, an id
is generated based on the residue index.
Returns: the newly created Residue
Parameters
----------
name : string
The name of the residue to add
chain : Chain
The Chain to add it to
id : string=None
An optional identifier for the residue. If this is omitted, an id
is generated based on the residue index.
Returns
-------
Residue
the newly created Residue
"""
if id is None:
id = str(self._numResidues+1)
......@@ -102,13 +134,22 @@ class Topology(object):
def addAtom(self, name, element, residue, id=None):
"""Create a new Atom and add it to the Topology.
Parameters:
- name (string) The name of the atom to add
- element (Element) The element of the atom to add
- residue (Residue) The Residue to add it to
- id (string=None) An optional identifier for the atom. If this is omitted, an id
is generated based on the atom index.
Returns: the newly created Atom
Parameters
----------
name : string
The name of the atom to add
element : Element
The element of the atom to add
residue : Residue
The Residue to add it to
id : string=None
An optional identifier for the atom. If this is omitted, an id is
generated based on the atom index.
Returns
-------
Atom
the newly created Atom
"""
if id is None:
id = str(self._numAtoms+1)
......@@ -120,9 +161,12 @@ class Topology(object):
def addBond(self, atom1, atom2):
"""Create a new bond and add it to the Topology.
Parameters:
- atom1 (Atom) The first Atom connected by the bond
- atom2 (Atom) The second Atom connected by the bond
Parameters
----------
atom1 : Atom
The first Atom connected by the bond
atom2 : Atom
The second Atom connected by the bond
"""
self._bonds.append((atom1, atom2))
......@@ -256,10 +300,13 @@ class Topology(object):
self.addBond(atomMaps[fromResidue][fromAtom], atomMaps[toResidue][toAtom])
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:
- positions (list) The list of atomic positions based on which to identify bonded atoms
Parameters
----------
positions : list
The list of atomic positions based on which to identify bonded atoms
"""
def isCyx(res):
names = [atom.name for atom in res._atoms]
......
wrappers/python/simtk/openmm/mtsintegrator.py
View file @ 1c938ceb
......@@ -10,7 +10,7 @@ Portions copyright (c) 2013-2015 Stanford University and the Authors.
Authors: Peter Eastman
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"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
......@@ -36,49 +36,52 @@ from simtk.openmm import CustomIntegrator
class MTSIntegrator(CustomIntegrator):
"""MTSIntegrator implements the rRESPA multiple time step integration algorithm.
This integrator allows different forces to be evaluated at different frequencies,
for example to evaluate the expensive, slowly changing forces less frequently than
the inexpensive, quickly changing forces.
To use it, you must first divide your forces into two or more groups (by calling
setForceGroup() on them) that should be evaluated at different frequencies. When
you create the integrator, you provide a tuple for each group specifying the index
of the force group and the frequency (as a fraction of the outermost time step) at
which to evaluate it. For example:
<pre>
integrator = MTSIntegrator(4*femtoseconds, [(0,1), (1,2), (2,8)])
</pre>
This specifies that the outermost time step is 4 fs, so each step of the integrator
will advance time by that much. It also says that force group 0 should be evaluated
once per time step, force group 1 should be evaluated twice per time step (every 2 fs),
and force group 2 should be evaluated eight times per time step (every 0.5 fs).
A common use of this algorithm is to evaluate reciprocal space nonbonded interactions
less often than the bonded and direct space nonbonded interactions. The following
example looks up the NonbondedForce, sets the reciprocal space interactions to their
own force group, and then creates an integrator that evaluates them once every 4 fs,
but all other interactions every 2 fs.
<pre>
nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
nonbonded.setReciprocalSpaceForceGroup(1)
integrator = MTSIntegrator(4*femtoseconds, [(1,1), (0,2)])
</pre>
For details, see Tuckerman et al., J. Chem. Phys. 97(3) pp. 1990-2001 (1992).
"""
def __init__(self, dt, groups):
"""Create an MTSIntegrator.
Parameters:
- dt (time) The largest (outermost) integration time step to use
- groups (list) A list of tuples defining the force groups. The first element of each
tuple is the force group index, and the second element is the number of times that force
group should be evaluated in one time step.
Parameters
----------
dt : time
The largest (outermost) integration time step to use
groups : list
A list of tuples defining the force groups. The first element of
each tuple is the force group index, and the second element is the
number of times that force group should be evaluated in one time step.
"""
if len(groups) == 0:
raise ValueError("No force groups specified")
......@@ -88,7 +91,7 @@ class MTSIntegrator(CustomIntegrator):
self.addUpdateContextState();
self._createSubsteps(1, groups)
self.addConstrainVelocities();
def _createSubsteps(self, parentSubsteps, groups):
group, substeps = groups[0]
stepsPerParentStep = substeps/parentSubsteps
......
wrappers/python/simtk/unit/unit.py 100644 → 100755
View file @ 1c938ceb
......@@ -54,9 +54,12 @@ class Unit(object):
def __init__(self, base_or_scaled_units):
"""Create a new Unit.
Parameters:
- self (Unit) The newly created Unit.
- base_or_scaled_units (dict) Keys are BaseUnits or ScaledUnits. Values are exponents (numbers).
Parameters
----------
self : Unit
The newly created Unit.
base_or_scaled_units : dict
Keys are BaseUnits or ScaledUnits. Values are exponents (numbers).
"""
# Unit contents are of two types: BaseUnits and ScaledUnits
self._top_base_units = {}
......@@ -389,7 +392,8 @@ class Unit(object):
Strips off any ScaledUnits in the Unit, leaving only BaseUnits.
Parameters
- system: a dictionary of (BaseDimension, BaseUnit) pairs
----------
system : a dictionary of (BaseDimension, BaseUnit) pairs
"""
return system.express_unit(self)
......@@ -583,7 +587,7 @@ class UnitSystem(object):
Parameters
----------
units: ``list``
units : list
List of base units from which to construct the unit system
"""
def __init__(self, units):
......@@ -678,7 +682,7 @@ def is_unit(x):
Returns True if x is a Unit, False otherwise.
Examples
--------
>>> is_unit(16)
False
"""
......
wrappers/python/src/swig_doxygen/OpenMM.i
View file @ 1c938ceb
......@@ -45,7 +45,7 @@ using namespace OpenMM;
%}
%feature("autodoc", "1");
%feature("autodoc", "0");
%nodefaultctor;
%include features.i
......@@ -61,49 +61,3 @@ using namespace OpenMM;
# namespace
__all__ = [k for k in locals().keys() if not (k.endswith('_swigregister') or k.startswith('_'))]
%}
/*
%extend OpenMM::XmlSerializer {
%template(XmlSerializer_serialize_AndersenThermostat) XmlSerializer::serialize<AndersenThermostat>;
%template(XmlSerializer_serialize_RBTorsionForce) XmlSerializer::serialize<RBTorsionForce>;
%template(XmlSerializer_serialize_CMAPTorsionForce) XmlSerializer::serialize<CMAPTorsionForce>;
%template(XmlSerializer_serialize_CMMotionRemover) XmlSerializer::serialize<CMMotionRemover>;
%template(XmlSerializer_serialize_CustomAngleForce) XmlSerializer::serialize<CustomAngleForce>;
%template(XmlSerializer_serialize_CustomBondForce) XmlSerializer::serialize<CustomBondForce>;
%template(XmlSerializer_serialize_CustomExternalForce) XmlSerializer::serialize<CustomExternalForce>;
%template(XmlSerializer_serialize_CustomGBForce) XmlSerializer::serialize<CustomGBForce>;
%template(XmlSerializer_serialize_CustomHbondForce) XmlSerializer::serialize<CustomHbondForce>;
%template(XmlSerializer_serialize_CustomNonbondedForce) XmlSerializer::serialize<CustomNonbondedForce>;
%template(XmlSerializer_serialize_CustomTorsionForce) XmlSerializer::serialize<CustomTorsionForce>;
%template(XmlSerializer_serialize_GBSAOBCForce) XmlSerializer::serialize<GBSAOBCForce>;
%template(XmlSerializer_serialize_GBVIForce) XmlSerializer::serialize<GBVIForce>;
%template(XmlSerializer_serialize_HarmonicAngleForce) XmlSerializer::serialize<HarmonicAngleForce>;
%template(XmlSerializer_serialize_HarmonicBondForce) XmlSerializer::serialize<HarmonicBondForce>;
%template(XmlSerializer_serialize_MonteCarloBarostat) XmlSerializer::serialize<MonteCarloBarostat>;
%template(XmlSerializer_serialize_MonteCarloAnisotropicBarostat) XmlSerializer::serialize<MonteCarloAnisotropicBarostat>;
%template(XmlSerializer_serialize_NonbondedForce) XmlSerializer::serialize<NonbondedForce>;
%template(XmlSerializer_serialize_RBTorsionForce) XmlSerializer::serialize<RBTorsionForce>;
%template(XmlSerializer_serialize_System) XmlSerializer::serialize<System>;
%template(XmlSerializer_deserialize_AndersenThermostat) XmlSerializer::deserialize<AndersenThermostat>;
%template(XmlSerializer_deserialize_RBTorsionForce) XmlSerializer::deserialize<RBTorsionForce>;
%template(XmlSerializer_deserialize_CMAPTorsionForce) XmlSerializer::deserialize<CMAPTorsionForce>;
%template(XmlSerializer_deserialize_CMMotionRemover) XmlSerializer::deserialize<CMMotionRemover>;
%template(XmlSerializer_deserialize_CustomAngleForce) XmlSerializer::deserialize<CustomAngleForce>;
%template(XmlSerializer_deserialize_CustomBondForce) XmlSerializer::deserialize<CustomBondForce>;
%template(XmlSerializer_deserialize_CustomExternalForce) XmlSerializer::deserialize<CustomExternalForce>;
%template(XmlSerializer_deserialize_CustomGBForce) XmlSerializer::deserialize<CustomGBForce>;
%template(XmlSerializer_deserialize_CustomHbondForce) XmlSerializer::deserialize<CustomHbondForce>;
%template(XmlSerializer_deserialize_CustomNonbondedForce) XmlSerializer::deserialize<CustomNonbondedForce>;
%template(XmlSerializer_deserialize_CustomTorsionForce) XmlSerializer::deserialize<CustomTorsionForce>;
%template(XmlSerializer_deserialize_GBSAOBCForce) XmlSerializer::deserialize<GBSAOBCForce>;
%template(XmlSerializer_deserialize_GBVIForce) XmlSerializer::deserialize<GBVIForce>;
%template(XmlSerializer_deserialize_HarmonicAngleForce) XmlSerializer::deserialize<HarmonicAngleForce>;
%template(XmlSerializer_deserialize_HarmonicBondForce) XmlSerializer::deserialize<HarmonicBondForce>;
%template(XmlSerializer_deserialize_MonteCarloBarostat) XmlSerializer::deserialize<MonteCarloBarostat>;
%template(XmlSerializer_deserialize_MonteCarloAnisotropicBarostat) XmlSerializer::deserialize<MonteCarloAnisotropicBarostat>;
%template(XmlSerializer_deserialize_NonbondedForce) XmlSerializer::deserialize<NonbondedForce>;
%template(XmlSerializer_deserialize_RBTorsionForce) XmlSerializer::deserialize<RBTorsionForce>;
%template(XmlSerializer_deserialize_System) XmlSerializer::deserialize<System>;
};
*/
wrappers/python/src/swig_doxygen/swigInputBuilder.py
View file @ 1c938ceb
......@@ -49,8 +49,7 @@ def striphtmltags(s):
s = s.replace('<i>', '_').replace('</i>', '_')
s = s.replace('<b>', '*').replace('</b>', '*')
s = re.sub('\s*(<ul>.*</ul>\s*)', replace_ul_tags, s, flags=re.MULTILINE | re.DOTALL)
s = re.sub('\s*(<ul>.*?</ul>\s*)', replace_ul_tags, s, flags=re.MULTILINE | re.DOTALL)
return s
def trimToSingleSpace(text):
......@@ -146,6 +145,20 @@ def getClassMethodList(classNode, skipMethods):
return methodList
def docstringTypemap(cpptype):
"""Translate a C++ type to Python for inclusion in the Python docstrings.
This doesn't need to be perfectly accurate -- it's not used for generating
the actual swig wrapper code. It's only used for generating the docstrings.
"""
pytype = cpptype
if pytype.startswith('const '):
pytype = pytype[6:]
if pytype.startswith('std::'):
pytype = pytype[5:]
pytype = pytype.strip('&')
return pytype.strip()
class SwigInputBuilder:
def __init__(self,
inputDirname,
......@@ -243,23 +256,43 @@ class SwigInputBuilder:
forceSubclassList.append(shortClassName)
elif baseName == 'OpenMM::Integrator':
integratorSubclassList.append(shortClassName)
self.fOut.write("%factory(OpenMM::Force& OpenMM::System::getForce")
for name in sorted(forceSubclassList):
self.fOut.write(",\n OpenMM::%s" % name)
self.fOut.write(");\n\n")
self.fOut.write("%factory(OpenMM::Force* OpenMM::Force::__copy__")
for name in sorted(forceSubclassList):
self.fOut.write(",\n OpenMM::%s" % name)
self.fOut.write(");\n\n")
self.fOut.write("%factory(OpenMM::Force* OpenMM_XmlSerializer__deserializeForce")
for name in sorted(forceSubclassList):
self.fOut.write(",\n OpenMM::%s" % name)
self.fOut.write(");\n\n")
self.fOut.write("%factory(OpenMM::Integrator* OpenMM::Integrator::__copy__")
for name in sorted(integratorSubclassList):
self.fOut.write(",\n OpenMM::%s" % name)
self.fOut.write(");\n\n")
self.fOut.write("%factory(OpenMM::Integrator* OpenMM_XmlSerializer__deserializeIntegrator")
for name in sorted(integratorSubclassList):
self.fOut.write(",\n OpenMM::%s" % name)
self.fOut.write(");\n\n")
self.fOut.write("%factory(OpenMM::Integrator& OpenMM::Context::getIntegrator")
for name in sorted(integratorSubclassList):
self.fOut.write(",\n OpenMM::%s" % name)
self.fOut.write(");\n\n")
self.fOut.write("%factory(OpenMM::VirtualSite& OpenMM::System::getVirtualSite, OpenMM::TwoParticleAverageSite, OpenMM::ThreeParticleAverageSite, OpenMM::OutOfPlaneSite);\n\n")
self.fOut.write("%factory(OpenMM::Integrator& OpenMM::CompoundIntegrator::getIntegrator")
for name in sorted(integratorSubclassList):
self.fOut.write(",\n OpenMM::%s" % name)
self.fOut.write(");\n\n")
self.fOut.write("%factory(OpenMM::VirtualSite& OpenMM::System::getVirtualSite, OpenMM::TwoParticleAverageSite, OpenMM::ThreeParticleAverageSite, OpenMM::OutOfPlaneSite, OpenMM::LocalCoordinatesSite);\n\n")
self.fOut.write("\n")
def writeGlobalConstants(self):
......@@ -345,7 +378,6 @@ class SwigInputBuilder:
self.fOut.write("\n%s};\n" % INDENT)
if len(enumNodes)>0: self.fOut.write("\n")
def writeMethods(self, classNode):
methodList=getClassMethodList(classNode, self.skipMethods)
......@@ -432,35 +464,44 @@ class SwigInputBuilder:
(shortClassName, memberNode,
shortMethDefinition, methName,
isConstructors, isDestructor, templateType, templateName) = items
paramList=findNodes(memberNode, 'param')
paramList = findNodes(memberNode, 'param')
#write pythonprepend blocks
# write pythonprepend blocks
mArgsstring = getText("argsstring", memberNode)
if self.fOutPythonprepend and \
len(paramList) and \
mArgsstring.find('=0')<0:
key=(shortClassName, methName)
if key in self.configModule.STEAL_OWNERSHIP:
for argNum in self.configModule.STEAL_OWNERSHIP[key]:
if self.SWIG_COMPACT_ARGUMENTS:
argName = 'args[%s]' % argNum
else:
argName = getText('declname', paramList[argNum])
mArgsstring.find('=0') < 0:
text = '''
%pythonprepend OpenMM::{shortClassName}::{methName}{mArgsstring} %{{{{{{0}}
%}}}}'''.format(shortClassName=shortClassName, methName=methName, mArgsstring=mArgsstring)
textInside = ''
key = (shortClassName, methName)
for argNum in self.configModule.STEAL_OWNERSHIP.get(key, []):
if self.SWIG_COMPACT_ARGUMENTS:
argName = 'args[%s]' % argNum
else:
argName = getText('declname', paramList[argNum])
text = '''
%pythonprepend OpenMM::{shortClassName}::{methName}{mArgsstring} %{{
textInside += '''
if not {argName}.thisown:
s = ("the %s object does not own its corresponding OpenMM object"
% self.__class__.__name__)
raise Exception(s)
%}}'''.format(argName=argName, shortClassName=shortClassName, methName=methName, mArgsstring=mArgsstring)
self.fOutPythonprepend.write(text)
raise Exception(s)'''.format(argName=argName)
for argNum in self.configModule.REQUIRE_ORDERED_SET.get(key, []):
if self.SWIG_COMPACT_ARGUMENTS:
argName = 'args[%s]' % argNum
else:
argName = getText('declname', paramList[argNum])
textInside += '''
{argName} = list({argName})'''.format(argName=argName)
if textInside:
self.fOutPythonprepend.write(text.format(textInside))
#write pythonappend blocks
# write pythonappend blocks
if self.fOutPythonappend \
and mArgsstring.find('=0')<0:
key=(shortClassName, methName)
and mArgsstring.find('=0') < 0:
key = (shortClassName, methName)
#print "key %s %s \n" % (shortClassName, methName)
addText=''
returnType = getText("type", memberNode)
......@@ -548,32 +589,58 @@ class SwigInputBuilder:
(shortClassName, memberNode,
shortMethDefinition, methName,
isConstructors, isDestructor, templateType, templateName ) = items
if self.fOutDocstring:
for dNode in findNodes(memberNode, 'detaileddescription'):
dString=""
try:
description=getText('para', dNode)
description.strip()
if description:
dString=description
except IndexError:
pass
params = findNodes(dNode, 'para/parameterlist/parameteritem')
if len(params) > 0:
dString="%s\n Parameters:" % dString
for pNode in params:
argName = getText('parameternamelist/parametername', pNode)
signatureParams = findNodes(memberNode, 'param')
assert len(findNodes(memberNode, 'detaileddescription')) == 1
dNode = findNodes(memberNode, 'detaileddescription')[0]
try:
description=getText('para', dNode)
description.strip()
except IndexError:
description = ''
params = findNodes(dNode, 'para/parameterlist/parameteritem')
paramString = ['Parameters', '----------']
returnString = ['Returns', '-------']
if len(params) > 0:
if len(signatureParams) != len(params):
raise ValueError('docstring in %s.%s does not match the signature' % (shortClassName, methName))
for pNode, pSignatureNode in zip(params, signatureParams):
parameterNameNode = findNodes(pNode, 'parameternamelist/parametername')[0]
argDoc = getText('parameterdescription/para', pNode)
dString="%s\n - %s %s" % (dString, argName, argDoc)
dString.strip()
if dString:
dString=re.sub(r'([^\\])"', r'\g<1>\"', dString)
s = '%%feature("docstring") OpenMM::%s::%s "%s";' \
% (shortClassName, methName, dString)
self.fOutDocstring.write("%s\n" % s)
self.fOutDocstring.write("\n\n")
#print "Done write Docstring info\n"
argName = getNodeText(parameterNameNode)
argType = docstringTypemap(getText('type', pSignatureNode))
isOutput = parameterNameNode.get('direction') == 'out'
if isOutput:
returnString.extend(['%s : %s' % (argName, argType), ' %s' % argDoc])
else:
paramString.extend(['%s : %s' % (argName, argType), ' %s' % argDoc])
returnSection = findNodes(dNode, 'para/simplesect')
if len(returnSection) > 0:
returnNode = returnSection[0]
if returnNode.get('kind') == 'return':
argType = getNodeText(findNodes(memberNode, 'type')[0])
argType = docstringTypemap(argType)
returnString.extend([argType, ' %s' % getNodeText(returnNode).strip()])
dString = '\n'.join(
([description] + [''] if len(description) > 0 else []) +
(paramString + [''] if len(paramString) > 2 else []) +
(returnString if len(returnString) > 2 else [])).strip()
if dString:
dString = re.sub(r'([^\\])"', r'\g<1>\"', dString)
s = '%%feature("docstring") OpenMM::%s::%s "%s";' \
% (shortClassName, methName, dString)
self.fOutDocstring.write("%s\n" % s)
self.fOutDocstring.write("\n\n")
def writeSwigFile(self):
......
wrappers/python/src/swig_doxygen/swigInputConfig.py
View file @ 1c938ceb
......@@ -43,7 +43,6 @@ SKIP_METHODS = [('State',),
('CalcCustomTorsionForceKernel',),
('CalcForcesAndEnergyKernel',),
('CalcGBSAOBCForceKernel',),
('CalcGBVIForceKernel',),
('CalcHarmonicAngleForceKernel',),
('CalcHarmonicBondForceKernel',),
('CalcKineticEnergyKernel',),
......@@ -141,8 +140,15 @@ STEAL_OWNERSHIP = {("Platform", "registerPlatform") : [0],
("CustomHbondForce", "addTabulatedFunction") : [1],
("CustomCompoundBondForce", "addTabulatedFunction") : [1],
("CustomManyParticleForce", "addTabulatedFunction") : [1],
("CompoundIntegrator", "addIntegrator") : [0],
}
REQUIRE_ORDERED_SET = {("CustomNonbondedForce", "addInteractionGroup") : [0, 1],
("CustomNonbondedForce", "setInteractionGroupParameters") : [1, 2],
}
# This is a list of units to attach to return values and method args.
# Indexed by (ClassName, MethodsName)
UNITS = {
......@@ -308,6 +314,7 @@ UNITS = {
("AmoebaWcaDispersionForce", "getShctd") : ( None, ()),
("Context", "getParameter") : (None, ()),
("Context", "getParameters") : (None, ()),
("Context", "getMolecules") : (None, ()),
("CMAPTorsionForce", "getMapParameters") : (None, (None, 'unit.kilojoule_per_mole')),
("CMAPTorsionForce", "getTorsionParameters") : (None, ()),
......@@ -385,14 +392,6 @@ UNITS = {
: (None, ('unit.elementary_charge',
'unit.nanometer', None)),
("GBSAOBCForce", "getSurfaceAreaEnergy") : ('unit.kilojoule_per_mole/unit.nanometer/unit.nanometer', ()),
("GBVIForce", "getBornRadiusScalingMethod") : (None, ()),
("GBVIForce", "getQuinticLowerLimitFactor") : (None, ()),
("GBVIForce", "getQuinticUpperBornRadiusLimit") : ('unit.nanometer', ()),
("GBVIForce", "getBondParameters")
: (None, (None, None, 'unit.nanometer')),
("GBVIForce", "getParticleParameters")
: (None, ('unit.elementary_charge',
'unit.nanometer', 'unit.kilojoule_per_mole')),
("HarmonicAngleForce", "getAngleParameters")
: (None, (None, None, None, 'unit.radian',
'unit.kilojoule_per_mole/(unit.radian*unit.radian)')),
......
wrappers/python/src/swig_doxygen/swig_lib/python/extend.i
View file @ 1c938ceb
%inline %{
typedef int bitmask32t;
%}
%typemap(in) bitmask32t %{
$1 = 0;
#if PY_VERSION_HEX >= 0x03000000
if (PyLong_Check($input)) {
unsigned long u = PyLong_AsUnsignedLongMask($input);
#else
if (PyInt_Check($input)) {
unsigned long u = PyInt_AsUnsignedLongMask($input);
#endif
// 64-bit Windows has 32-bit longs, but other platforms have
// 64-bit longs
$1 = u & 0xffffffff;
} else {
PyErr_SetString(PyExc_ValueError, "in method $symname, argument $argnum could not be converted to type $type");
SWIG_fail;
}
%}
%extend OpenMM::Context {
PyObject *_getStateAsLists(int getPositions,
int getVelocities,
......@@ -5,7 +29,7 @@
int getEnergy,
int getParameters,
int enforcePeriodic,
int groups) {
bitmask32t groups) {
State state;
PyThreadState* _savePythonThreadState = PyEval_SaveThread();
int types = 0;
......@@ -27,54 +51,53 @@
%pythoncode %{
def getState(self,
getPositions=False,
getVelocities=False,
getForces=False,
getEnergy=False,
getParameters=False,
enforcePeriodicBox=False,
groups=-1):
"""
getState(self,
getPositions = False,
getVelocities = False,
getForces = False,
getEnergy = False,
getParameters = False,
enforcePeriodicBox = False,
groups = -1)
-> State
Get a State object recording the current state information stored in this context.
Parameters:
- getPositions (bool=False) whether to store particle positions in the State
- getVelocities (bool=False) whether to store particle velocities in the State
- getForces (bool=False) whether to store the forces acting on particles in the State
- getEnergy (bool=False) whether to store potential and kinetic energy in the State
- getParameter (bool=False) whether to store context parameters in the State
- enforcePeriodicBox (bool=False) if false, the position of each particle will be whatever position is stored in the Context, regardless of periodic boundary conditions. If true, particle positions will be translated so the center of every molecule lies in the same periodic box.
- groups (int=-1) a set of bit flags for which force groups to include when computing forces and energies. Group i will be included if (groups&(1<<i)) != 0. The default value includes all groups.
def getState(self, getPositions=False, getVelocities=False,
getForces=False, getEnergy=False, getParameters=False,
enforcePeriodicBox=False, groups=-1):
"""Get a State object recording the current state information stored in this context.
Parameters
----------
getPositions : bool=False
whether to store particle positions in the State
getVelocities : bool=False
whether to store particle velocities in the State
getForces : bool=False
whether to store the forces acting on particles in the State
getEnergy : bool=False
whether to store potential and kinetic energy in the State
getParameter : bool=False
whether to store context parameters in the State
enforcePeriodicBox : bool=False
if false, the position of each particle will be whatever position
is stored in the Context, regardless of periodic boundary conditions.
If true, particle positions will be translated so the center of
every molecule lies in the same periodic box.
groups : set={0,1,2,...,31}
a set of indices for which force groups to include when computing
forces and energies. The default value includes all groups. groups
can also be passed as an unsigned integer interpreted as a bitmask,
in which case group i will be included if (groups&(1<<i)) != 0.
"""
if getPositions: getP=1
else: getP=0
if getVelocities: getV=1
else: getV=0
if getForces: getF=1
else: getF=0
if getEnergy: getE=1
else: getE=0
if getParameters: getPa=1
else: getPa=0
if enforcePeriodicBox: enforcePeriodic=1
else: enforcePeriodic=0
getP, getV, getF, getE, getPa, enforcePeriodic = map(bool,
(getPositions, getVelocities, getForces, getEnergy, getParameters,
enforcePeriodicBox))
try:
# is the input integer-like?
groups_mask = int(groups)
except TypeError:
if isinstance(groups, set):
# nope, okay, then it should be an set
groups_mask = functools.reduce(operator.or_,
((1<<x) & 0xffffffff for x in groups))
else:
raise TypeError('%s is neither an int nor set' % groups)
(simTime, periodicBoxVectorsList, energy, coordList, velList,
forceList, paramMap) = \
self._getStateAsLists(getP, getV, getF, getE, getPa, enforcePeriodic, groups)
self._getStateAsLists(getP, getV, getF, getE, getPa, enforcePeriodic, groups_mask)
state = State(simTime=simTime,
energy=energy,
coordList=coordList,
......@@ -83,11 +106,11 @@
periodicBoxVectorsList=periodicBoxVectorsList,
paramMap=paramMap)
return state
def setState(self, state):
"""
setState(Context self, State state)
Copy information from a State object into this Context. This restores the Context to
approximately the same state it was in when the State was created. If the State does not include
a piece of information (e.g. positions or velocities), that aspect of the Context is
......@@ -108,7 +131,7 @@
for param in state._paramMap:
self.setParameter(param, state._paramMap[param])
%}
%feature("docstring") createCheckpoint "Create a checkpoint recording the current state of the Context.
This should be treated as an opaque block of binary data. See loadCheckpoint() for more details.
......@@ -185,48 +208,51 @@ Parameters:
getParameters=False,
enforcePeriodicBox=False,
groups=-1):
"""Get a State object recording the current state information about one copy of the system.
Parameters
----------
copy : int
the index of the copy for which to retrieve state information
getPositions : bool=False
whether to store particle positions in the State
getVelocities : bool=False
whether to store particle velocities in the State
getForces : bool=False
whether to store the forces acting on particles in the State
getEnergy : bool=False
whether to store potential and kinetic energy in the State
getParameter : bool=False
whether to store context parameters in the State
enforcePeriodicBox : bool=False
if false, the position of each particle will be whatever position
is stored in the Context, regardless of periodic boundary conditions.
If true, particle positions will be translated so the center of
every molecule lies in the same periodic box.
groups : set={0,1,2,...,31}
a set of indices for which force groups to include when computing
forces and energies. The default value includes all groups. groups
can also be passed as an unsigned integer interpreted as a bitmask,
in which case group i will be included if (groups&(1<<i)) != 0.
"""
getState(self,
copy,
getPositions = False,
getVelocities = False,
getForces = False,
getEnergy = False,
getParameters = False,
enforcePeriodicBox = False,
groups = -1)
-> State
Get a State object recording the current state information about one copy of the system.
Parameters:
- copy (int) the index of the copy for which to retrieve state information
- getPositions (bool=False) whether to store particle positions in the State
- getVelocities (bool=False) whether to store particle velocities in the State
- getForces (bool=False) whether to store the forces acting on particles in the State
- getEnergy (bool=False) whether to store potential and kinetic energy in the State
- getParameter (bool=False) whether to store context parameters in the State
- enforcePeriodicBox (bool=False) if false, the position of each particle will be whatever position is stored in the Context, regardless of periodic boundary conditions. If true, particle positions will be translated so the center of every molecule lies in the same periodic box.
- groups (int=-1) a set of bit flags for which force groups to include when computing forces and energies. Group i will be included if (groups&(1<<i)) != 0. The default value includes all groups.
"""
if getPositions: getP=1
else: getP=0
if getVelocities: getV=1
else: getV=0
if getForces: getF=1
else: getF=0
if getEnergy: getE=1
else: getE=0
if getParameters: getPa=1
else: getPa=0
if enforcePeriodicBox: enforcePeriodic=1
else: enforcePeriodic=0
getP, getV, getF, getE, getPa, enforcePeriodic = map(bool,
(getPositions, getVelocities, getForces, getEnergy, getParameters,
enforcePeriodicBox))
try:
# is the input integer-like?
groups_mask = int(groups)
except TypeError:
if isinstance(groups, set):
groups_mask = functools.reduce(operator.or_,
((1<<x) & 0xffffffff for x in groups))
else:
raise TypeError('%s is neither an int nor set' % groups)
(simTime, periodicBoxVectorsList, energy, coordList, velList,
forceList, paramMap) = \
self._getStateAsLists(copy, getP, getV, getF, getE, getPa, enforcePeriodic, groups)
forceList, paramMap) = \
self._getStateAsLists(copy, getP, getV, getF, getE, getPa, enforcePeriodic, groups_mask)
state = State(simTime=simTime,
energy=energy,
coordList=coordList,
......@@ -299,7 +325,7 @@ Parameters:
ss << inputString;
return OpenMM::XmlSerializer::deserialize<OpenMM::System>(ss);
}
static std::string _serializeForce(const OpenMM::Force* object) {
std::stringstream ss;
OpenMM::XmlSerializer::serialize<OpenMM::Force>(object, "Force", ss);
......@@ -312,7 +338,7 @@ Parameters:
ss << inputString;
return OpenMM::XmlSerializer::deserialize<OpenMM::Force>(ss);
}
static std::string _serializeIntegrator(const OpenMM::Integrator* object) {
std::stringstream ss;
OpenMM::XmlSerializer::serialize<OpenMM::Integrator>(object, "Integrator", ss);
......@@ -327,8 +353,8 @@ Parameters:
}
static std::string _serializeStateAsLists(
const std::vector<Vec3>& pos,
const std::vector<Vec3>& vel,
const std::vector<Vec3>& pos,
const std::vector<Vec3>& vel,
const std::vector<Vec3>& forces,
double kineticEnergy,
double potentialEnergy,
......@@ -341,7 +367,7 @@ Parameters:
OpenMM::XmlSerializer::serialize<OpenMM::State>(&myState, "State", buffer);
return buffer.str();
}
static PyObject* _deserializeStringIntoLists(const std::string &stateAsString) {
std::stringstream ss;
ss << stateAsString;
......@@ -369,7 +395,7 @@ Parameters:
try:
velocities = pythonState.getVelocities().value_in_unit(unit.nanometers/unit.picoseconds)
types |= 2
except:
except:
pass
try:
forces = pythonState.getForces().value_in_unit(unit.kilojoules_per_mole/unit.nanometers)
......@@ -390,14 +416,14 @@ Parameters:
time = pythonState.getTime().value_in_unit(unit.picoseconds)
boxVectors = pythonState.getPeriodicBoxVectors().value_in_unit(unit.nanometers)
string = XmlSerializer._serializeStateAsLists(positions, velocities, forces, kineticEnergy, potentialEnergy, time, boxVectors, params, types)
return string
return string
@staticmethod
def _deserializeState(pythonString):
(simTime, periodicBoxVectorsList, energy, coordList, velList,
forceList, paramMap) = XmlSerializer._deserializeStringIntoLists(pythonString)
state = State(simTime=simTime,
energy=energy,
coordList=coordList,
......@@ -450,6 +476,14 @@ Parameters:
%extend OpenMM::Force {
%pythoncode %{
def __getstate__(self):
serializationString = XmlSerializer.serialize(self)
return serializationString
def __setstate__(self, serializationString):
system = XmlSerializer.deserialize(serializationString)
self.this = system.this
def __deepcopy__(self, memo):
return self.__copy__()
%}
......
wrappers/python/src/swig_doxygen/swig_lib/python/pythoncode.i
View file @ 1c938ceb
......@@ -8,6 +8,8 @@ except ImportError:
import copy
import sys
import math
import functools
import operator
RMIN_PER_SIGMA=math.pow(2, 1/6.0)
RVDW_PER_SIGMA=math.pow(2, 1/6.0)/2.0
if sys.version_info[0] == 2:
......
wrappers/python/src/swig_doxygen/swig_lib/python/typemaps.i
View file @ 1c938ceb
%fragment("Py_StripOpenMMUnits", "header") {
%fragment("Vec3_to_PyVec3", "header") {
/**
* Convert an OpenMM::Vec3 into a Python simtk.openmm.Vec3 object
*
* Returns a new reference.
*/
PyObject* Vec3_to_PyVec3(const OpenMM::Vec3& v) {
static PyObject *__s_mm = NULL;
static PyObject *__s_Vec3 = NULL;
if (__s_mm == NULL) {
__s_mm = PyImport_AddModule("simtk.openmm");
__s_Vec3 = PyObject_GetAttrString(__s_mm, "Vec3");
}
PyObject* tuple = Py_BuildValue("(d,d,d)", v[0], v[1], v[2]);
PyObject* PyVec3 = PyObject_CallObject(__s_Vec3, tuple);
Py_DECREF(tuple);
return PyVec3;
}
}
static PyObject *__s_Quantity = NULL;
static PyObject *__s_md_unit_system_tuple = NULL;
static PyObject *__s_bar_tuple = NULL;
%fragment("Py_StripOpenMMUnits", "header") {
/**
* Strip any OpenMM units of an input PyObject.
*
* This is equivalent to the following Python code
*
* >>> from simtk import unit
* >>> if isinstance(input, unit.Quantity)
* ... if input.is_compatible(unit.bar)
* ... return input.value_in_unit(unit.bar)
* ... return input.value_in_unit_system(unit.md_input_system)
* ... return input
*
* Returns a new reference.
*/
PyObject* Py_StripOpenMMUnits(PyObject *input) {
static PyObject *__s_Quantity = NULL;
static PyObject *__s_md_unit_system_tuple = NULL;
static PyObject *__s_bar_tuple = NULL;
if (__s_Quantity == NULL) {
PyObject* module = NULL;
module = PyImport_ImportModule("simtk.unit");
......@@ -412,19 +446,12 @@ int Py_SequenceToVecVecVecDouble(PyObject* obj, std::vector<std::vector<std::vec
}
%typemap(argout) std::vector<Vec3>& {
int i, n;
PyObject *pyList;
n=(*$1).size();
pyList=PyList_New(n);
PyObject* mm = PyImport_AddModule("simtk.openmm");
PyObject* vec3 = PyObject_GetAttrString(mm, "Vec3");
for (i=0; i<n; i++) {
%typemap(argout, fragment="Vec3_to_PyVec3") std::vector<Vec3>& {
int n = (*$1).size();
PyObject * pyList = PyList_New(n);
for (int i=0; i<n; i++) {
OpenMM::Vec3& v = (*$1).at(i);
PyObject* args = Py_BuildValue("(d,d,d)", v[0], v[1], v[2]);
PyObject* pyVec = PyObject_CallObject(vec3, args);
Py_DECREF(args);
PyObject* pyVec = Vec3_to_PyVec3(v);
PyList_SET_ITEM(pyList, i, pyVec);
}
$result = pyList;
......@@ -436,41 +463,20 @@ int Py_SequenceToVecVecVecDouble(PyObject* obj, std::vector<std::vector<std::vec
}
%typemap(out) Vec3 {
PyObject* mm = PyImport_AddModule("simtk.openmm");
PyObject* vec3 = PyObject_GetAttrString(mm, "Vec3");
PyObject* args = Py_BuildValue("(d,d,d)", ($1)[0], ($1)[1], ($1)[2]);
$result = PyObject_CallObject(vec3, args);
Py_DECREF(args);
%typemap(out, fragment="Vec3_to_PyVec3") Vec3 {
$result = Vec3_to_PyVec3(*$1);
}
%typemap(out) const Vec3& {
PyObject* mm = PyImport_AddModule("simtk.openmm");
PyObject* vec3 = PyObject_GetAttrString(mm, "Vec3");
PyObject* args = Py_BuildValue("(d,d,d)", (*$1)[0], (*$1)[1], (*$1)[2]);
$result = PyObject_CallObject(vec3, args);
Py_DECREF(args);
%typemap(out, fragment="Vec3_to_PyVec3") const Vec3& {
$result = Vec3_to_PyVec3(*$1);
}
/* Convert C++ (Vec3&, Vec3&, Vec3&) object to python tuple or tuples */
%typemap(argout) (Vec3& a, Vec3& b, Vec3& c) {
// %typemap(argout) (Vec3& a, Vec3& b, Vec3& c)
PyObject* mm = PyImport_AddModule("simtk.openmm");
PyObject* vec3 = PyObject_GetAttrString(mm, "Vec3");
PyObject* args1 = Py_BuildValue("(d,d,d)", (*$1)[0], (*$1)[1], (*$1)[2]);
PyObject* args2 = Py_BuildValue("(d,d,d)", (*$2)[0], (*$2)[1], (*$2)[2]);
PyObject* args3 = Py_BuildValue("(d,d,d)", (*$3)[0], (*$3)[1], (*$3)[2]);
PyObject* pyVec1 = PyObject_CallObject(vec3, args1);
PyObject* pyVec2 = PyObject_CallObject(vec3, args2);
PyObject* pyVec3 = PyObject_CallObject(vec3, args3);
Py_DECREF(args1);
Py_DECREF(args2);
Py_DECREF(args3);
Py_DECREF(mm);
Py_DECREF(vec3);
%typemap(argout, fragment="Vec3_to_PyVec3") (Vec3& a, Vec3& b, Vec3& c) {
PyObject* pyVec1 = Vec3_to_PyVec3(*$1);
PyObject* pyVec2 = Vec3_to_PyVec3(*$2);
PyObject* pyVec3 = Vec3_to_PyVec3(*$3);
PyObject *o, *o2, *o3;
o = Py_BuildValue("[N, N, N]", pyVec1, pyVec2, pyVec3);
if ((!$result) || ($result == Py_None)) {
......
wrappers/python/tests/TestCharmmFiles.py
View file @ 1c938ceb
......@@ -4,6 +4,7 @@ from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm.app.element as elem
import warnings
class TestCharmmFiles(unittest.TestCase):
......@@ -90,7 +91,6 @@ class TestCharmmFiles(unittest.TestCase):
def test_NBFIX(self):
"""Tests CHARMM systems with NBFIX Lennard-Jones modifications"""
import warnings
warnings.filterwarnings('ignore', category=CharmmPSFWarning)
psf = CharmmPsfFile('systems/ala3_solv.psf')
crd = CharmmCrdFile('systems/ala3_solv.crd')
......@@ -122,6 +122,34 @@ class TestCharmmFiles(unittest.TestCase):
self.assertEqual(len(list(psf.topology.residues())), 20169)
self.assertEqual(len(list(psf.topology.bonds())), 46634)
def testSystemOptions(self):
""" Test various options in CharmmPsfFile.createSystem """
warnings.filterwarnings('ignore', category=CharmmPSFWarning)
psf = CharmmPsfFile('systems/ala3_solv.psf')
crd = CharmmCrdFile('systems/ala3_solv.crd')
params = CharmmParameterSet('systems/par_all36_prot.prm',
'systems/toppar_water_ions.str')
# Box dimensions (found from bounding box)
psf.setBox(32.7119500*angstroms, 32.9959600*angstroms, 33.0071500*angstroms)
# Check some illegal options
self.assertRaises(ValueError, lambda:
psf.createSystem(params, nonbondedMethod=5))
self.assertRaises(TypeError, lambda:
psf.createSystem(params, nonbondedMethod=PME,
nonbondedCutoff=1*radian)
)
self.assertRaises(TypeError, lambda:
psf.createSystem(params, nonbondedMethod=PME,
switchDistance=1*radian)
)
# Check what should be some legal options
psf.createSystem(params, nonbondedMethod=PME, switchDistance=0.8,
nonbondedCutoff=1.2)
psf.createSystem(params, nonbondedMethod=PME, switchDistance=0.8,
nonbondedCutoff=1.2*nanometer)
def test_ImplicitSolventForces(self):
"""Compute forces for different implicit solvent types, and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
solventType = [HCT, OBC1, OBC2, GBn, GBn2]
......
wrappers/python/tests/TestForceField.py
View file @ 1c938ceb
......@@ -6,6 +6,10 @@ from simtk.unit import *
import simtk.openmm.app.element as elem
import simtk.openmm.app.forcefield as forcefield
import math
if sys.version_info >= (3, 0):
from io import StringIO
else:
from cStringIO import StringIO
class TestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method."""
......@@ -144,7 +148,8 @@ class TestForceField(unittest.TestCase):
context = Context(system, integrator)
context.setPositions(pdb.positions)
state1 = context.getState(getForces=True)
state2 = XmlSerializer.deserialize(open('systems/lysozyme-implicit-forces.xml').read())
with open('systems/lysozyme-implicit-forces.xml') as input:
state2 = XmlSerializer.deserialize(input.read())
numDifferences = 0
for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
diff = norm(f1-f2)
......@@ -197,7 +202,51 @@ class TestForceField(unittest.TestCase):
for i in range(3):
self.assertEqual(vectors[i], self.pdb1.topology.getPeriodicBoxVectors()[i])
self.assertEqual(vectors[i], system.getDefaultPeriodicBoxVectors()[i])
def test_ResidueAttributes(self):
"""Test a ForceField that gets per-particle parameters from residue attributes."""
xml = """
<ForceField>
<AtomTypes>
<Type name="tip3p-O" class="OW" element="O" mass="15.99943"/>
<Type name="tip3p-H" class="HW" element="H" mass="1.007947"/>
</AtomTypes>
<Residues>
<Residue name="HOH">
<Atom name="O" type="tip3p-O" charge="-0.834"/>
<Atom name="H1" type="tip3p-H" charge="0.417"/>
<Atom name="H2" type="tip3p-H" charge="0.417"/>
<Bond from="0" to="1"/>
<Bond from="0" to="2"/>
</Residue>
</Residues>
<NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
<UseAttributeFromResidue name="charge"/>
<Atom type="tip3p-O" sigma="0.315" epsilon="0.635"/>
<Atom type="tip3p-H" sigma="1" epsilon="0"/>
</NonbondedForce>
</ForceField>"""
ff = ForceField(StringIO(xml))
# Build a water box.
modeller = Modeller(Topology(), [])
modeller.addSolvent(ff, boxSize=Vec3(3, 3, 3)*nanometers)
# Create a system and make sure all nonbonded parameters are correct.
system = ff.createSystem(modeller.topology)
nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
atoms = list(modeller.topology.atoms())
for i in range(len(atoms)):
params = nonbonded.getParticleParameters(i)
if atoms[i].element == elem.oxygen:
self.assertEqual(params[0], -0.834*elementary_charge)
self.assertEqual(params[1], 0.315*nanometers)
self.assertEqual(params[2], 0.635*kilojoule_per_mole)
else:
self.assertEqual(params[0], 0.417*elementary_charge)
self.assertEqual(params[1], 1.0*nanometers)
self.assertEqual(params[2], 0.0*kilojoule_per_mole)
class AmoebaTestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method with the AMOEBA forcefield."""
......@@ -270,6 +319,22 @@ class AmoebaTestForceField(unittest.TestCase):
self.assertAlmostEqual(constraints[(0,2)], hoDist)
self.assertAlmostEqual(constraints[(1,2)], hohDist)
def test_Forces(self):
"""Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
pdb = PDBFile('systems/alanine-dipeptide-implicit.pdb')
forcefield = ForceField('amoeba2013.xml', 'amoeba2013_gk.xml')
system = forcefield.createSystem(pdb.topology, polarization='direct')
integrator = VerletIntegrator(0.001)
context = Context(system, integrator)
context.setPositions(pdb.positions)
state1 = context.getState(getForces=True)
with open('systems/alanine-dipeptide-amoeba-forces.xml') as input:
state2 = XmlSerializer.deserialize(input.read())
for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
diff = norm(f1-f2)
self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-3)
if __name__ == '__main__':
unittest.main()
wrappers/python/tests/TestForceGroups.py 0 → 100644
View file @ 1c938ceb
import unittest
import itertools
import simtk.openmm as mm
class TestForceGroups(unittest.TestCase):
def setUp(self):
system = mm.System()
system.addParticle(1.0)
for i in range(32):
force = mm.CustomExternalForce(str(i))
force.addParticle(0, [])
force.setForceGroup(i)
system.addForce(force)
platform = mm.Platform.getPlatformByName('Reference')
context = mm.Context(system, mm.VerletIntegrator(0), platform)
context.setPositions([(0,0,0)])
self.context = context
def test1(self):
n = 32
for (i,j) in itertools.combinations(range(n), 2):
groups = 1<<i | 1<<j
e_0 = self.context.getState(getEnergy=True, groups=groups).getPotentialEnergy()._value
e_1 = self.context.getState(getEnergy=True, groups={i,j}).getPotentialEnergy()._value
e_ref = i+j
self.assertEqual(e_0, e_ref)
self.assertEqual(e_1, e_ref)
def test2(self):
with self.assertRaises(TypeError):
# groups must be an int or set
self.context.getState(getEnergy=True, groups=(1, 2))
def test3(self):
e_0 = self.context.getState(getEnergy=True, groups=-1).getPotentialEnergy()._value
e_ref = sum(range(32))
self.assertEqual(e_0, e_ref)
if __name__ == '__main__':
unittest.main()
wrappers/python/tests/TestModeller.py
View file @ 1c938ceb
......@@ -3,6 +3,10 @@ from validateModeller import *
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
if sys.version_info >= (3, 0):
from io import StringIO
else:
from cStringIO import StringIO
class TestModeller(unittest.TestCase):
""" Test the Modeller class. """
......@@ -397,44 +401,46 @@ class TestModeller(unittest.TestCase):
def test_addSolventNegativeSolvent(self):
""" Test the addSolvent() method; test adding ions to a negatively charged solvent. """
topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers)
# set up modeller with no solvent
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
# add 5 Cl- ions to the original topology
topology_toAdd = Topology()
newChain = topology_toAdd.addChain()
for i in range(5):
topology_toAdd.addResidue('CL', newChain)
residues = [residue for residue in topology_toAdd.residues()]
for i in range(5):
topology_toAdd.addAtom('Cl',Element.getBySymbol('Cl'), residues[i])
positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
modeller.add(topology_toAdd, positions_toAdd)
modeller.addSolvent(self.forcefield, ionicStrength=1.0*molar)
topology_after = modeller.getTopology()
for neutralize in (True, False):
# set up modeller with no solvent
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
water_count = 0
sodium_count = 0
chlorine_count = 0
for residue in topology_after.residues():
if residue.name=='HOH':
water_count += 1
elif residue.name=='NA':
sodium_count += 1
elif residue.name=='CL':
chlorine_count += 1
# add 5 Cl- ions to the original topology
topology_toAdd = Topology()
newChain = topology_toAdd.addChain()
for i in range(5):
topology_toAdd.addResidue('CL', newChain)
residues = [residue for residue in topology_toAdd.residues()]
for i in range(5):
topology_toAdd.addAtom('Cl',Element.getBySymbol('Cl'), residues[i])
positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
modeller.add(topology_toAdd, positions_toAdd)
modeller.addSolvent(self.forcefield, ionicStrength=1.0*molar, neutralize=neutralize)
topology_after = modeller.getTopology()
total_water_ions = water_count+sodium_count+chlorine_count
expected_ion_fraction = 1.0*molar/(55.4*molar)
expected_ions = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5
self.assertEqual(sodium_count, expected_ions)
self.assertEqual(chlorine_count, expected_ions)
water_count = 0
sodium_count = 0
chlorine_count = 0
for residue in topology_after.residues():
if residue.name=='HOH':
water_count += 1
elif residue.name=='NA':
sodium_count += 1
elif residue.name=='CL':
chlorine_count += 1
total_water_ions = water_count+sodium_count+chlorine_count
expected_ion_fraction = 1.0*molar/(55.4*molar)
expected_chlorine = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5
expected_sodium = expected_chlorine if neutralize else expected_chlorine-5
self.assertEqual(sodium_count, expected_sodium)
self.assertEqual(chlorine_count, expected_chlorine)
def test_addSolventPositiveSolvent(self):
""" Test the addSolvent() method; test adding ions to a positively charged solvent. """
......@@ -442,42 +448,44 @@ class TestModeller(unittest.TestCase):
topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 3.5, 3.5)*nanometers)
# set up modeller with no solvent
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
for neutralize in (True, False):
# set up modeller with no solvent
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
# add 5 Na+ ions to the original topology
topology_toAdd = Topology()
newChain = topology_toAdd.addChain()
for i in range(5):
topology_toAdd.addResidue('NA', newChain)
residues = [residue for residue in topology_toAdd.residues()]
for i in range(5):
topology_toAdd.addAtom('Na',Element.getBySymbol('Na'), residues[i])
positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
# positions_toAdd doesn't need to change
modeller.add(topology_toAdd, positions_toAdd)
modeller.addSolvent(self.forcefield, ionicStrength=1.0*molar)
topology_after = modeller.getTopology()
# add 5 Na+ ions to the original topology
topology_toAdd = Topology()
newChain = topology_toAdd.addChain()
for i in range(5):
topology_toAdd.addResidue('NA', newChain)
residues = [residue for residue in topology_toAdd.residues()]
for i in range(5):
topology_toAdd.addAtom('Na',Element.getBySymbol('Na'), residues[i])
positions_toAdd = [Vec3(1.0,1.2,1.5), Vec3(1.7,1.0,1.4), Vec3(1.5,2.0,1.0),
Vec3(2.0,2.0,2.0), Vec3(2.0,1.5,1.0)]*nanometers
# positions_toAdd doesn't need to change
modeller.add(topology_toAdd, positions_toAdd)
modeller.addSolvent(self.forcefield, ionicStrength=1.0*molar, neutralize=neutralize)
topology_after = modeller.getTopology()
water_count = 0
sodium_count = 0
chlorine_count = 0
for residue in topology_after.residues():
if residue.name=='HOH':
water_count += 1
elif residue.name=='NA':
sodium_count += 1
elif residue.name=='CL':
chlorine_count += 1
water_count = 0
sodium_count = 0
chlorine_count = 0
for residue in topology_after.residues():
if residue.name=='HOH':
water_count += 1
elif residue.name=='NA':
sodium_count += 1
elif residue.name=='CL':
chlorine_count += 1
total_water_ions = water_count+sodium_count+chlorine_count
expected_ion_fraction = 1.0*molar/(55.4*molar)
expected_ions = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5
self.assertEqual(sodium_count, expected_ions)
self.assertEqual(chlorine_count, expected_ions)
total_water_ions = water_count+sodium_count+chlorine_count
expected_ion_fraction = 1.0*molar/(55.4*molar)
expected_sodium = math.floor((total_water_ions-10)*expected_ion_fraction/2+0.5)+5
expected_chlorine = expected_sodium if neutralize else expected_sodium-5
self.assertEqual(sodium_count, expected_sodium)
self.assertEqual(chlorine_count, expected_chlorine)
def test_addSolventIons(self):
""" Test the addSolvent() method with all possible choices for positive and negative ions. """
......@@ -915,6 +923,135 @@ class TestModeller(unittest.TestCase):
self.assertEqual(1, len(ep))
def test_addVirtualSites(self):
"""Test adding extra particles defined by virtual sites."""
xml = """
<ForceField>
<AtomTypes>
<Type name="C" class="C" element="C" mass="10"/>
<Type name="N" class="N" element="N" mass="10"/>
<Type name="O" class="O" element="O" mass="10"/>
<Type name="V" class="V" mass="0.0"/>
</AtomTypes>
<Residues>
<Residue name="Test">
<Atom name="C" type="C"/>
<Atom name="N" type="N"/>
<Atom name="O" type="O"/>
<Atom name="V1" type="V"/>
<Atom name="V2" type="V"/>
<Atom name="V3" type="V"/>
<Atom name="V4" type="V"/>
<VirtualSite type="average2" index="3" atom1="0" atom2="1" weight1="0.7" weight2="0.3"/>
<VirtualSite type="average3" index="4" atom1="0" atom2="1" atom3="2" weight1="0.2" weight2="0.3" weight3="0.5"/>
<VirtualSite type="outOfPlane" index="5" atom1="0" atom2="1" atom3="2" weight12="0.1" weight13="-0.2" weightCross="0.8"/>
<VirtualSite type="localCoords" index="6" atom1="0" atom2="1" atom3="2" wo1="0.1" wo2="0.5" wo3="0.4" wx1="1" wx2="-0.6" wx3="-0.4" wy1="0.1" wy2="0.9" wy3="-1" p1="-0.5" p2="0.4" p3="1.1"/>
</Residue>
</Residues>
</ForceField>"""
ff = ForceField(StringIO(xml))
# Create the three real atoms.
topology = Topology()
chain = topology.addChain()
residue = topology.addResidue('Test', chain)
topology.addAtom('C', element.carbon, residue)
topology.addAtom('N', element.nitrogen, residue)
topology.addAtom('V', element.oxygen, residue)
# Add the virtual sites.
modeller = Modeller(topology, [Vec3(0.1, 0.2, 0.3), Vec3(1.0, 0.9, 0.8), Vec3(1.5, 1.1, 0.7)]*nanometers)
modeller.addExtraParticles(ff)
top = modeller.topology
pos = modeller.positions
# Check that the correct particles were added.
self.assertEqual(len(pos), 7)
for atom, elem in zip(top.atoms(), [element.carbon, element.nitrogen, element.oxygen, None, None, None, None]):
self.assertEqual(elem, atom.element)
# Check that the positions were calculated correctly.
system = ff.createSystem(top)
integ = VerletIntegrator(1.0)
context = Context(system, integ)
context.setPositions(pos)
context.computeVirtualSites()
state = context.getState(getPositions=True)
for p1, p2 in zip (pos, state.getPositions()):
self.assertVecAlmostEqual(p1.value_in_unit(nanometers), p2.value_in_unit(nanometers), 1e-6)
def test_multiSiteIon(self):
"""Test adding extra particles whose positions are determined based on bonds."""
xml = """
<ForceField>
<AtomTypes>
<Type name="Zn" class="Zn" element="Zn" mass="53.380"/>
<Type name="DA" class="DA" mass="3.0"/>
</AtomTypes>
<Residues>
<Residue name="ZN">
<Atom name="ZN" type="Zn"/>
<Atom name="D1" type="DA"/>
<Atom name="D2" type="DA"/>
<Atom name="D3" type="DA"/>
<Atom name="D4" type="DA"/>
<Bond from="0" to="2"/>
<Bond from="0" to="1"/>
<Bond from="0" to="3"/>
<Bond from="0" to="4"/>
<Bond from="1" to="2"/>
<Bond from="1" to="3"/>
<Bond from="1" to="4"/>
<Bond from="2" to="4"/>
<Bond from="2" to="3"/>
<Bond from="3" to="4"/>
</Residue>
</Residues>
<HarmonicBondForce>
<Bond class1="DA" class2="Zn" length="0.09" k="535552.0"/>
<Bond class1="DA" class2="DA" length="0.147" k="535552.0"/>
</HarmonicBondForce>
</ForceField>"""
ff = ForceField(StringIO(xml))
# Create two zinc atoms.
topology = Topology()
chain = topology.addChain()
residue = topology.addResidue('ZN', chain)
topology.addAtom('ZN', element.zinc, residue)
residue = topology.addResidue('ZN', chain)
topology.addAtom('ZN', element.zinc, residue)
# Add the extra particles.
modeller = Modeller(topology, [Vec3(0.5, 1.0, 1.5), Vec3(2.0, 2.0, 0.0)]*nanometers)
modeller.addExtraParticles(ff)
top = modeller.topology
pos = modeller.positions
# Check that the correct particles were added.
self.assertEqual(len(pos), 10)
for i, atom in enumerate(top.atoms()):
self.assertEqual(element.zinc if i in (0,5) else None, atom.element)
# Check that the positions in the first residue are reasonable.
center = Vec3(0.5, 1.0, 1.5)*nanometers
self.assertEqual(center, modeller.positions[0])
for i in range(1, 5):
for j in range(i):
dist = norm(pos[i]-pos[j])
expectedDist = 0.09 if j == 0 else 0.147
self.assertTrue(dist > (expectedDist-0.01)*nanometers and dist < (expectedDist+0.01)*nanometers)
def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
scale = max(1.0, norm(p1),)
for i in range(3):
......
wrappers/python/tests/TestPickle.py
View file @ 1c938ceb
......@@ -3,13 +3,14 @@ from validateConstraints import *
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm
import simtk.openmm.app.element as elem
import simtk.openmm.app.forcefield as forcefield
import copy
import pickle
class TestPickle(unittest.TestCase):
"""Pickling / deepcopy of OpenMM state and integrator objects."""
"""Pickling / deepcopy of OpenMM objects."""
def setUp(self):
"""Set up the tests by loading the input pdb files and force field
......@@ -26,28 +27,46 @@ class TestPickle(unittest.TestCase):
self.pdb2 = PDBFile('systems/alanine-dipeptide-implicit.pdb')
self.forcefield2 = ForceField('amber99sb.xml', 'amber99_obc.xml')
def check_copy(self, object, object_copy):
"""Check that an object's copy is an accurate replica."""
# Check class name is same.
self.assertEqual(object.__class__.__name__, object_copy.__class__.__name__)
# Check serialized contents are the same.
self.assertEqual(XmlSerializer.serialize(object), XmlSerializer.serialize(object_copy))
def test_deepcopy(self):
"""Test that serialization/deserialization works (via deepcopy)."""
# Create system, integrator, and state.
system = self.forcefield1.createSystem(self.pdb1.topology)
integrator = VerletIntegrator(2*femtosecond)
context = Context(system, integrator)
context.setPositions(self.pdb1.positions)
state = context.getState(getPositions=True, getForces=True, getEnergy=True)
system2 = copy.deepcopy(system)
integrator2 = copy.deepcopy(integrator)
state2 = copy.deepcopy(state)
str_state = pickle.dumps(state)
str_integrator = pickle.dumps(integrator)
state3 = pickle.loads(str_state)
context.setState(state3)
#
# Test deepcopy
#
self.check_copy(system, copy.deepcopy(system))
self.check_copy(integrator, copy.deepcopy(integrator))
self.check_copy(state, copy.deepcopy(state))
for force_index in range(system.getNumForces()):
force = system.getForce(force_index)
force_copy = copy.deepcopy(force)
self.check_copy(force, force_copy)
del context, integrator
#
# Test pickle
#
self.check_copy(system, pickle.loads(pickle.dumps(system)))
self.check_copy(integrator, pickle.loads(pickle.dumps(integrator)))
self.check_copy(state, pickle.loads(pickle.dumps(state)))
for force_index in range(system.getNumForces()):
force = system.getForce(force_index)
force_copy = pickle.loads(pickle.dumps(force))
self.check_copy(force, force_copy)
if __name__ == '__main__':
unittest.main()
......
wrappers/python/tests/TestSwigWrappers.py 0 → 100644
View file @ 1c938ceb
import unittest
import simtk.openmm as mm
class TestSwigWrappers(unittest.TestCase):
def test_1(self):
# This tests for a refcounting bug in the swig wrappers
# that was previously problematic.
# See https://github.com/pandegroup/openmm/issues/1214
for cycle in range(10):
system = mm.System()
system.getDefaultPeriodicBoxVectors()
if __name__ == '__main__':
unittest.main()
wrappers/python/tests/TestTopology.py 0 → 100644
View file @ 1c938ceb
import sys
import unittest
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm.app.element as elem
if sys.version_info >= (3, 0):
from io import StringIO
else:
from cStringIO import StringIO
class TestTopology(unittest.TestCase):
"""Test the Topology object"""
def check_pdbfile(self, pdbfilename, natoms, nres, nchains):
"""Check that a PDB file has the specified number of atoms, residues, and chains."""
pdb = PDBFile(pdbfilename)
top = pdb.topology
self.assertEqual(pdb.topology.getNumAtoms(), natoms)
self.assertEqual(pdb.topology.getNumResidues(), nres)
self.assertEqual(pdb.topology.getNumChains(), nchains)
def test_getters(self):
"""Test getters for number of atoms, residues, chains."""
self.check_pdbfile('systems/1T2Y.pdb', 271, 25, 1)
if __name__ == '__main__':
unittest.main()
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