understanding openMM minimization, equilibration and production

[Gudrun Gygli]

I have recently switched to openMM coming from GROMACS for simulations of proteins and ligands in explicit solvent
.
In GROMACS, energy minimization, equilibration and production parameters are stored and read in ".mdp"-files, for example:
Energy minimization:

Code: Select all

; LINES STARTING WITH ';' ARE COMMENTS from tutorial gromacs
title		  	 	 = Minimization	; Title of run

; Parameters describing what to do, when to stop and what to save
integrator	  	 	 = steep		; Algorithm (steep = steepest descent minimization)
emtol			   	 = 1000.0  	; Stop minimization when the maximum force < 10.0 kJ/mol
emstep      	   	 = 0.01      ; Energy step size
nsteps		   		 = 50000	  	; Maximum number of (minimization) steps to perform
nstxout-compressed	 = 500		; save coordinates every 1 ps (to xtc file)

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist		    = 1		        ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet
ns_type		    = grid		    ; Method to determine neighbor list (simple, grid)
rlist		    = 1.2		    ; Cut-off for making neighbor list (short range forces)
coulombtype	    = PME		    ; Treatment of long range electrostatic interactions
rcoulomb	    = 1.2		    ; long range electrostatic cut-off
vdwtype         = cutoff
vdw-modifier    = force-switch
rvdw-switch     = 1.0
rvdw		    = 1.2		    ; long range Van der Waals cut-off
pbc             = xyz 		    ; Periodic Boundary Conditions
DispCorr        = no

All this is now just one line of code in openMM:

Code: Select all

simulation.minimizeEnergy()

And according to the docs (http://docs.openmm.org/7.1.0/api-python ... ation.html), minimizeEnergy() only accepts 2 parameters: tolerance and maxIterations.

SO if I understand correctly, I can define only one integrator for my entire simulation using for example :

Code: Select all

integrator = LangevinIntegrator(
                        temperature*unit.kelvin,       # Temperature of heat bath
                        1.0/unit.picoseconds,  # Friction coefficient
                        2.0*unit.femtoseconds, # Time step
)

which I then call when defining my simulation, including other parameters that I define:

Code: Select all

simulation = Simulation(complex.topology, system, integrator, platform)

I am having trouble wrapping my head around this because in GROMACS I have previously used different integrators for minimization and equilibration/production (steepest descent and leap-frog/leap-frog, respectively), and I am now wondering if this is possible/recommended to do in openMM.

Any thoughts/input on this are welcome.

[Peter Eastman]

That's correct. The fact that GROMACS represents energy minimization as "integration" is kind of an implementation quirk.

OpenMM does provide ways of switching integrators part way through a simulation (for example with CompoundIntegrator), but that's a more obscure feature only used in less common simulation protocols. In most cases, you only use a single integrator.

[Gudrun Gygli]

So in the example I posted, openMM would use Langevin Dynamics and minimize until the default values for minimizeEnergy() are reached?

Or does it ignore the Langevin integrator and use a default energy minimization algorithm?

I might "just" have a problem with the terminology...

I cannot find online how to set a leapfrog integrator in openMM, like I do in the GROMACS mdp posted in the end.

I found how to set PME in openMM (https://gpantel.github.io/assets/PDF/Op ... torial.pdf):

Code: Select all

# use PME for long-range electrostatics, cutoff for short-range interactions 
# constrain H-bonds with RATTLE, constrain water with SETTLE 
system = topology.createSystem(parameters, nonbondedMethod=app.PME,  
    nonbondedCutoff=1.2*unit.nanometers, constraints=app.HBonds, rigidWater=True,  
    ewaldErrorTolerance=0.0005)

This tutorial (https://gpantel.github.io/assets/PDF/Op ... torial.pdf) also says:

Code: Select all

...
# There is only a Langevin and Andersen thermostat in OpenMM in the pre-compiled distribution 
# However, you can define custom integrators (more on this later) 
print('Constructing integrator') 
integrator = mm.LangevinIntegrator(325*unit.kelvin, 1.0/unit.picosecond, 2.0*unit.femtosecond)

This was in 2019, now it is 2021, is there a berendsen thermostat implemented now?

Where can I find a list of available integrators, thermostats and cutoff schemes for electrostatic interactions? (is there such a list?)

In case it matters: I am using openmm 7.5.1

Example of a GROMACS mdp file for production:

Code: Select all

;title       = Protein  NPT run

; Run parameters
integrator  = md        ; leap-frog integrator
nsteps      = 250000000  ; 0.002 * 250000000 = 500000 ps
dt          = 0.002     ; 2 fs
; Output control
nstcomm	    = 1
nstxout     = 500000       ; save coordinates every 1000 ps
nstvout     = 50000       ; save velocities every 100 ps
nstenergy   = 50000       ; save energies every 100 ps
nstlog      = 100000      ; update log file every 200 ps
nstxout-compressed = 2500 ; save coordinates every 5 ps in .xtc format
energygrps  = protein non-protein
; Bond parameters
continuation    = yes           ; dynamics run
constraint_algorithm = lincs    ; holonomic constraints 
constraints     = h-bonds       ; h-bonds constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
ns_type     = grid      ; search neighboring grid cells
nstlist     = 10        ; 20 fs
rlist       = 1.4       ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.4       ; short-range electrostatic cutoff (in nm)
vdwtype     = cut-off
rvdw        = 1.4       ; short-range van der Waals cutoff (in nm)
cutoff-scheme   = Verlet

; Electrostatics
coulombtype     = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order       = 4         ; cubic interpolation
fourierspacing  = 0.16      ; grid spacing for FFT
ewald_rtol      = 1e-5
optimize_fft    = yes
; Berendsen temperature coupling is on
tcoupl      = V-rescale                 ; thermostat
tc-grps     = Protein non-protein   	; coupling group
tau_t       = 0.1  0.1                  ; time constant, in ps
ref_t       = 298  298                  ; reference temperature in K
; Pressure coupling is on
pcoupl      = Berendsen         ; pressure coupling is on for NPT
pcoupltype  = isotropic         ; uniform scaling of box vectors
tau_p       = 0.6               ; time constant, in ps
ref_p       = 1.0               ; reference pressure, in bar
compressibility = 4.5e-5        ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc         = xyz       ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = no        ; assign velocities from Maxwell distribution
gen_temp    = 298       ; temperature for Maxwell distribution
gen_seed    = -1        ; generate a random seed

[Peter Eastman]

minimizeEnergy() uses an L-BFGS optimizer. It does that regardless of what integrator you're using.

I cannot find online how to set a leapfrog integrator in openMM

You specify an integrator object, and that determines the algorithm. Most simulations use either LangevinMiddleIntegrator (for constant temperature) or VerletIntegrator (for constant energy). Both of those use leapfrog algorithms. See http://docs.openmm.org/latest/userguide ... ntegrators for brief descriptions of all the standard integrators, or http://docs.openmm.org/latest/userguide ... ators.html for more detail.