Merge pull request #2531 from peastman/params

Updated simulation parameters in documentation and examples

docs-source/usersguide/application.rst

@@ -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()
@@ -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)
@@ -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:
@@ -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()
@@ -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()
@@ -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
@@ -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
 ===============
@@ -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
@@ -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
@@ -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
@@ -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::

docs-source/usersguide/library.rst

@@ -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)
 
@@ -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);

examples/HelloSodiumChloride.cpp

@@ -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)
 
@@ -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);

examples/HelloSodiumChlorideInC.c

@@ -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)  */
 
@@ -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);

examples/HelloSodiumChlorideInFortran.f90

@@ -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)
     
@@ -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.
@@ -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)

examples/simulateAmber.py

@@ -6,7 +6,7 @@ from sys import stdout
 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:

examples/simulateCharmm.py

@@ -20,9 +20,8 @@ params = CharmmParameterSet('charmm22.rtf', 'charmm22.par')
 # 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()

examples/simulateGromacs.py

@@ -6,7 +6,7 @@ from sys import stdout
 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()

examples/simulatePdb.py

@@ -6,7 +6,7 @@ from sys import stdout
 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()