Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updated simulation parameters in documentation and examples #2531

Merged
merged 1 commit into from
Jan 20, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
35 changes: 19 additions & 16 deletions docs-source/usersguide/application.rst
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ steps.
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
Expand Down Expand Up @@ -210,14 +210,14 @@ convenient and less error-prone. We could have equivalently specified
The units system will be described in more detail later, in Section :ref:`units-and-dimensional-analysis`.
::

integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)

This line creates the integrator to use for advancing the equations of motion.
It specifies a :class:`BAOABLangevinIntegrator`, which performs Langevin dynamics,
and assigns it to a variable called :code:`integrator`\ . It also specifies
the values of three parameters that are specific to Langevin dynamics: the
simulation temperature (300 K), the friction coefficient (1 ps\ :sup:`-1`\ ), and
the step size (0.002 ps).
the step size (0.004 ps).
::

simulation = Simulation(pdb.topology, system, integrator)
Expand Down Expand Up @@ -295,7 +295,7 @@ found in OpenMM’s :file:`examples` folder with the name :file:`simulateAmber.p
inpcrd = AmberInpcrdFile('input.inpcrd')
system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
constraints=HBonds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(prmtop.topology, system, integrator)
simulation.context.setPositions(inpcrd.positions)
if inpcrd.boxVectors is not None:
Expand Down Expand Up @@ -389,7 +389,7 @@ with the name :file:`simulateGromacs.py`.
includeDir='/usr/local/gromacs/share/gromacs/top')
system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
constraints=HBonds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(top.topology, system, integrator)
simulation.context.setPositions(gro.positions)
simulation.minimizeEnergy()
Expand Down Expand Up @@ -453,7 +453,7 @@ on the :class:`CharmmPsfFile`.
params = CharmmParameterSet('charmm22.rtf', 'charmm22.prm')
system = psf.createSystem(params, nonbondedMethod=NoCutoff,
nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(psf.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
Expand Down Expand Up @@ -981,9 +981,10 @@ Value Meaning

The main reason to use constraints is that it allows one to use a larger
integration time step. With no constraints, one is typically limited to a time
step of about 1 fs for typical biomolecular force fields like AMBER or CHARMM. With :code:`HBonds` constraints, this can be increased
to about 2 fs. With :code:`HAngles`\ , it can be further increased to 3.5 or
4 fs.
step of about 1 fs for typical biomolecular force fields like AMBER or CHARMM.
With :code:`HBonds` constraints, this can be increased to about 2 fs for Verlet
dynamics, or about 4 fs for Langevin dynamics. With :code:`HAngles`\ , it can
sometimes be increased even further.

Regardless of the value of this parameter, OpenMM makes water molecules
completely rigid, constraining both their bond lengths and angles. You can
Expand All @@ -997,7 +998,9 @@ step size, typically to about 0.5 fs.

.. note::

The AMOEBA forcefield is intended to be used without constraints.
The AMOEBA forcefield is designed to be used without constraints, so by
default OpenMM makes AMOEBA water flexible. You can still force it to be
rigid by specifying :code:`rigidWater=True`.

Heavy Hydrogens
===============
Expand All @@ -1012,7 +1015,7 @@ optionally tell OpenMM to increase the mass of hydrogen atoms. For example,
This applies only to hydrogens that are bonded to heavy atoms, and any mass
added to the hydrogen is subtracted from the heavy atom. This keeps their total
mass constant while slowing down the fast motions of hydrogens. When combined
with constraints (typically :code:`constraints=AllBonds`\ ), this allows a
with constraints (typically :code:`constraints=AllBonds`\ ), this often allows a
further increase in integration step size.

Integrators
Expand All @@ -1028,13 +1031,13 @@ BAOAB Langevin Integrator
In the examples of the previous sections, we used Langevin integration:
::

integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)

The three parameter values in this line are the simulation temperature (300 K),
the friction coefficient (1 ps\ :sup:`-1`\ ), and the step size (0.002 ps). You
the friction coefficient (1 ps\ :sup:`-1`\ ), and the step size (0.004 ps). You
are free to change these to whatever values you want. Be sure to specify units
on all values. For example, the step size could be written either as
:code:`0.002*picoseconds` or :code:`2*femtoseconds`\ . They are exactly
:code:`0.004*picoseconds` or :code:`4*femtoseconds`\ . They are exactly
equivalent.

Langevin Integrator
Expand Down Expand Up @@ -1155,7 +1158,7 @@ previous section:
system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
constraints=HBonds)
system.addForce(MonteCarloBarostat(1*bar, 300*kelvin))
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
...

The parameters of the Monte Carlo barostat are the pressure (1 bar) and
Expand Down Expand Up @@ -1747,7 +1750,7 @@ coordinates. Here is the code to do it:
system.addForce(force)
for i in range(system.getNumParticles()):
force.addParticle(i, [])
integrator = BAOABLangevinIntegrator(300*kelvin, 91/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 91/picosecond, 0.004*picoseconds)
...

.. caption::
Expand Down
4 changes: 2 additions & 2 deletions docs-source/usersguide/library.rst
Original file line number Diff line number Diff line change
Expand Up @@ -1297,7 +1297,7 @@ existing MD program.
static const double SolventDielectric = 80.; // typical for water
static const double SoluteDielectric = 2.; // typical for protein

static const double StepSizeInFs = 2; // integration step size (fs)
static const double StepSizeInFs = 4; // integration step size (fs)
static const double ReportIntervalInFs = 50; // how often to issue PDB frame (fs)
static const double SimulationTimeInPs = 100; // total simulation time (ps)

Expand Down Expand Up @@ -1631,7 +1631,7 @@ along with the handle :code:`omm`\ , back to the calling function.
// best available Platform. Initialize the configuration from the default
// positions we collected above. Initial velocities will be zero but could
// have been set here.
omm->integrator = new OpenMM::LangevinIntegrator(temperature,
omm->integrator = new OpenMM::BAOABLangevinIntegrator(temperature,
frictionInPs,
stepSizeInFs * OpenMM::PsPerFs);
omm->context = new OpenMM::Context(*omm->system, *omm->integrator);
Expand Down
4 changes: 2 additions & 2 deletions examples/HelloSodiumChloride.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ static const double FrictionInPerPs = 91.; // collisions per picosecond
static const double SolventDielectric = 80.; // typical for water
static const double SoluteDielectric = 2.; // typical for protein

static const double StepSizeInFs = 2; // integration step size (fs)
static const double StepSizeInFs = 4; // integration step size (fs)
static const double ReportIntervalInFs = 50; // how often to issue PDB frame (fs)
static const double SimulationTimeInPs = 100; // total simulation time (ps)

Expand Down Expand Up @@ -249,7 +249,7 @@ myInitializeOpenMM( const MyAtomInfo atoms[],
// best available Platform. Initialize the configuration from the default
// positions we collected above. Initial velocities will be zero but could
// have been set here.
omm->integrator = new OpenMM::LangevinIntegrator(temperature, frictionInPs,
omm->integrator = new OpenMM::BAOABLangevinIntegrator(temperature, frictionInPs,
stepSizeInFs * OpenMM::PsPerFs);
omm->context = new OpenMM::Context(*omm->system, *omm->integrator);
omm->context->setPositions(initialPosInNm);
Expand Down
4 changes: 2 additions & 2 deletions examples/HelloSodiumChlorideInC.c
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ static const double FrictionInPerPs = 91.; /*collisions per ps*/
static const double SolventDielectric = 80.; /*typical for water */
static const double SoluteDielectric = 2.; /*typical for protein */

static const double StepSizeInFs = 2; /*integration step size (fs) */
static const double StepSizeInFs = 4; /*integration step size (fs) */
static const double ReportIntervalInFs = 50; /*how often for PDB frame (fs)*/
static const double SimulationTimeInPs = 100; /*total simulation time (ps) */

Expand Down Expand Up @@ -252,7 +252,7 @@ myInitializeOpenMM( const MyAtomInfo atoms[],
* best available Platform. Initialize the configuration from the default
* positions we collected above. Initial velocities will be zero but could
* have been set here. */
omm->integrator = (OpenMM_Integrator*)OpenMM_LangevinIntegrator_create(
omm->integrator = (OpenMM_Integrator*)OpenMM_BAOABLangevinIntegrator_create(
temperature, frictionInPerPs,
stepSizeInFs * OpenMM_PsPerFs);
omm->context = OpenMM_Context_create(omm->system, omm->integrator);
Expand Down
8 changes: 4 additions & 4 deletions examples/HelloSodiumChlorideInFortran.f90
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ MODULE MyAtomInfo
parameter(SoluteDielectric = 2) !typical for protein

real*8 StepSizeInFs, ReportIntervalInFs, SimulationTimeInPs
parameter(StepSizeInFs = 2) !integration step size (fs)
parameter(StepSizeInFs = 4) !integration step size (fs)
parameter(ReportIntervalInFs = 50) !how often for PDB frame (fs)
parameter(SimulationTimeInPs = 100) !total simulation time (ps)

Expand Down Expand Up @@ -171,7 +171,7 @@ SUBROUTINE myInitializeOpenMM(ommHandle, platformName)
! These are the objects we'll create here thare are stored in the
! Context for later access. Don't forget to delete them at the end.
type (OpenMM_System) system
type (OpenMM_LangevinIntegrator) langevin
type (OpenMM_BAOABLangevinIntegrator) langevin
type (OpenMM_Context) context

! These are temporary OpenMM objects used and discarded here.
Expand Down Expand Up @@ -236,11 +236,11 @@ SUBROUTINE myInitializeOpenMM(ommHandle, platformName)
! best available Platform. Initialize the configuration from the default
! positions we collected above. Initial velocities will be zero but could
! have been set here.
call OpenMM_LangevinIntegrator_create(langevin, &
call OpenMM_BAOABLangevinIntegrator_create(langevin, &
Temperature, FrictionInPerPs, &
StepSizeInFs * OpenMM_PsPerFs)

! Convert LangevinIntegrator to generic Integrator type for this call.
! Convert BAOABLangevinIntegrator to generic Integrator type for this call.
call OpenMM_Context_create(context, system, &
transfer(langevin, OpenMM_Integrator(0)))
call OpenMM_Context_setPositions(context, initialPosInNm)
Expand Down
2 changes: 1 addition & 1 deletion examples/simulateAmber.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
prmtop = AmberPrmtopFile('input.prmtop')
inpcrd = AmberInpcrdFile('input.inpcrd')
system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(prmtop.topology, system, integrator)
simulation.context.setPositions(inpcrd.positions)
if inpcrd.boxVectors is not None:
Expand Down
5 changes: 2 additions & 3 deletions examples/simulateCharmm.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,8 @@
# http://mackerell.umaryland.edu/CHARMM_ff_params.html

# Instantiate the system
system = psf.createSystem(params, nonbondedMethod=NoCutoff,
nonbondedCutoff=None)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
system = psf.createSystem(params, nonbondedMethod=NoCutoff, constraints=HBonds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(psf.topology, system, integrator)
simulation.context.setPositions(pdb.getPositions())
simulation.minimizeEnergy()
Expand Down
2 changes: 1 addition & 1 deletion examples/simulateGromacs.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
gro = GromacsGroFile('input.gro')
top = GromacsTopFile('input.top', periodicBoxVectors=gro.getPeriodicBoxVectors())
system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(top.topology, system, integrator)
simulation.context.setPositions(gro.positions)
simulation.minimizeEnergy()
Expand Down
2 changes: 1 addition & 1 deletion examples/simulatePdb.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
pdb = PDBFile('input.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
integrator = BAOABLangevinIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
Expand Down