Amber masthead
Filler image AmberTools21 Amber20 Manuals Tutorials Force Fields Contacts History
Filler image

Useful links:

Amber Home
Download Amber
Installation
Amber Citations
GPU Support
Updates
Mailing Lists
For Educators
File Formats
8 Chemical Reactions and Equilibria
 

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Vinícius Cruzeiro. 2018


Section 2: Minimization, heating, and equilibration

In this section we will describe how to perform minimization, heating, and equilibration on the system built in Section 1.

These are three general setup phases that come before the main production C(pH,E)MD simulations (common also with many other types of simulations):

  1. Minimize system to relax bad contacts (jump here)
  2. Heat the system slowly to the target temperature (jump here)
  3. Equilibrate the system at the target temperature (jump here)

Minimizing the system

In this first stage, we will minimize the system applying restraint on backbone atoms and on some heme atoms. The initial charge distribution on the titratable residues corresponds to the default states on leaprc.constph and leaprc.conste in tleap. Thus, we begin with the heme group in the oxidized state, and both propionates and the glutamate in the deprotonated state.

The input files for both implicit and explicit solvent simulations are shown below. The ntmin flag represents the minimization method. When ntmin is 1, the steepest descent method is going to be used during the initial ncyc cycles followed by the conjugate gradient method. The maxcyc flag represents the total number of cycles.

min.implicit.mdin min.explicit.mdin
Energy Minimization Stage in Implicit Solvent
&cntrl
  imin=1,            ! Perform an energy minimization
  ntb=0,             ! Boundaries (=0, non-periodic)
  cut=1000.0,        ! Cutoff
  ntmin=1,           ! Method of minimization
  ncyc=100,          ! Number of steepest descent cycles
  maxcyc=4000,       ! Maximum number of minimization cycles
  igb=2,             ! GB model
  saltcon=0.1,       ! Salt concentration
  ntr=1,             ! Turn on positional restraints
  restraintmask='@CA,C,O,N|(:HEH,PRN)@FE,NA,NB,NC,ND,C3D,C2A,C3B,C2C,CA,CB',
  restraint_wt=10.0, ! 10 kcal/mol.A**2 restraint force constant
/
Energy Minimization Stage in Explicit Solvent
&cntrl
  imin=1,            ! Perform an energy minimization
  ntb=1,             ! Boundaries (=1, constant volume)
  cut=8.0,           ! Cutoff
  ntmin=1,           ! Method of minimization
  ncyc=1000,         ! Number of steepest descent cycles
  maxcyc=5000,       ! Maximum number of minimization cycles
  igb=0,             ! GB model (=0, explicit solvent)
  saltcon=0.1,       ! Salt concentration
  ntr=1,             ! Turn on positional restraints
  restraintmask='@CA,C,O,N|(:HEH,PRN)@FE,NA,NB,NC,ND,C3D,C2A,C3B,C2C,CA,CB',
  restraint_wt=10.0, ! 10 kcal/mol.A**2 restraint force constant
/

For implicit solvent simulation, the command

pmemd -O -i min.implicit.mdin -p mp8_is.prmtop -c mp8_is.rst7 -inf mp8_is.min.mdinfo -o mp8_is.min.mdout -r mp8_is.min.rst7 -ref mp8_is.rst7

will generate the following files: mp8_is.min.mdout and mp8_is.min.rst7.

     

For explicit solvent simulation, the command

pmemd -O -i min.explicit.mdin -p mp8_es.new.prmtop -c mp8_es.rst7 -inf mp8_es.min.mdinfo -o mp8_es.min.mdout -r mp8_es.min.rst7 -ref mp8_es.rst7

will generate the following files: mp8_es.min.mdout and mp8_es.min.rst7.

The minimization may also be performed using sander instead of pmemd. It is a good idea to inspect the mdout files to evaluate how the energies behave during minimization. This can be done with the aid of the mdout_analyzer.py AmberTool. For instructions on how to do this, please check here. We should see the enegies dropping as the structure relaxes.

Heating the system

We will now perform a heating simulation at constant volume on the minimized structures. Here, we will heat the system slowly from 10 to 300 K during the initial 0.6 ns, followed by 2.4 ns of simulation with a target temperature of 300 K. We keep the restraints on the backbone atoms, but with a weaker force constant in comparison to the one used during minimization. The input files for implicit and explicit solvent simulations are shown below:

heat.implicit.mdin heat.explicit.mdin
Heating Stage in Implicit Solvent
&cntrl
  ntt=3,             ! Temperature scaling (=3, Langevin dynamics)
  gamma_ln=5.0,      ! Collision frequency of the Langevin dynamics in ps-1
  ntc=2,             ! SHAKE constraints (=2, hydrogen bond lengths constrained)
  ntf=2,             ! Force evaluation (=2, hydrogen bond interactions omitted)
  ntb=0,             ! Boundaries (=0, non-periodic)
  cut=1000.0,        ! Cutoff
  dt=0.002,          ! The time step in picoseconds
  nstlim=1500000,    ! Number of MD steps to be performed
  ig=-1,             ! Random seed (=-1, get a number from current date and time)
  ntwr=10000,        ! Restart file written every ntwr steps
  ntwx=10000,        ! Trajectory file written every ntwx steps
  ntpr=10000,        ! The mdout and mdinfo files written every ntpr steps
  ioutfm=1,          ! Trajectory file format (=1, Binary NetCDF)
  igb=2,             ! GB model
  saltcon=0.1,       ! Salt concentration
  ntr=1,             ! Turn on positional restraints
  restraintmask='@CA,C,O,N|(:HEH,PRN)@FE,NA,NB,NC,ND,C3D,C2A,C3B,C2C,CA,CB',
  restraint_wt=1.0,  ! 1 kcal/mol.A**2 restraint force constant
  nmropt=1,          ! To vary target temperature below
/
&wt
  TYPE='TEMP0', ISTEP1=0, ISTEP2=300000
  VALUE1=10.00, VALUE2=300.0,
/
&wt TYPE='END'
/
Heating Stage in Explicit Solvent
&cntrl
  ntt=3,             ! Temperature scaling (=3, Langevin dynamics)
  gamma_ln=5.0,      ! Collision frequency of the Langevin dynamics in ps-1
  ntc=2,             ! SHAKE constraints (=2, hydrogen bond lengths constrained)
  ntf=2,             ! Force evaluation (=2, hydrogen bond interactions omitted)
  ntb=1,             ! Boundaries (=1, constant volume)
  cut=8.0,           ! Cutoff
  dt=0.002,          ! The time step in picoseconds
  nstlim=1500000,    ! Number of MD steps to be performed
  ig=-1,             ! Random seed (=-1, get a number from current date and time)
  ntwr=10000,        ! Restart file written every ntwr steps
  ntwx=10000,        ! Trajectory file written every ntwx steps
  ntpr=10000,        ! The mdout and mdinfo files written every ntpr steps
  ioutfm=1,          ! Trajectory file format (=1, Binary NetCDF)
  iwrap=1,           ! Translate water molecules into the original simulation box
  igb=0,             ! GB model (=0, explicit solvent)
  saltcon=0.1,       ! Salt concentration
  ntr=1,             ! Turn on positional restraints
  restraintmask='@CA,C,O,N|(:HEH,PRN)@FE,NA,NB,NC,ND,C3D,C2A,C3B,C2C,CA,CB',
  restraint_wt=1.0,  ! 1 kcal/mol.A**2 restraint force constant
  nmropt=1,          ! To vary target temperature below
/
&wt
  TYPE='TEMP0', ISTEP1=0, ISTEP2=300000
  VALUE1=10.00, VALUE2=300.0,
/
&wt TYPE='END'
/

For implicit solvent simulation, the command

pmemd.cuda -O -i heat.implicit.mdin -p mp8_is.prmtop -c mp8_is.min.rst7 -x mp8_is.heat.nc -inf mp8_is.heat.mdinfo -o mp8_is.heat.mdout -r mp8_is.heat.rst7 -ref mp8_is.min.rst7

will generate the files: mp8_is.heat.mdout, mp8_is.heat.nc and mp8_is.heat.rst7.

     

For explicit solvent simulation, the command

pmemd.cuda -O -i heat.explicit.mdin -p mp8_es.new.prmtop -c mp8_es.min.rst7 -x mp8_es.heat.nc -inf mp8_es.heat.mdinfo -o mp8_es.heat.mdout -r mp8_es.heat.rst7 -ref mp8_es.min.rst7

will generate the files: mp8_es.heat.mdout, mp8_es.heat.nc and mp8_es.heat.rst7.

These simulations may also be performed using pmemd, pmemd.MPI, sander or sander.MPI instead of using the GPU-accelerated code with pmemd.cuda, however, the calculations may take a significantly longer time to finish.

Equilibrating the system

In the equilibration stage, we will let the system equilibrate more at the target temperature of 300 K. Here, the backbone restraints are still in place, however, with a weaker force constant in comparison to the heating stage. The explicit solvent simulation is performed at constant pressure to allow the system density to be stabilized. The simulations are performed for 10 ns in implicit solvent and 8 ns in explicit solvent. Follows the input files for the simulations using both solvent models:

equi.implicit.mdin equi.explicit.mdin
Equilibration Stage in Implicit Solvent
&cntrl
  ntt=3,             ! Temperature scaling (=3, Langevin dynamics)
  gamma_ln=10.0,     ! Collision frequency of the Langevin dynamics in ps-1
  ntc=2,             ! SHAKE constraints (=2, hydrogen bond lengths constrained)
  ntf=2,             ! Force evaluation (=2, hydrogen bond interactions omitted)
  ntb=0,             ! Boundaries (=0, non-periodic)
  cut=1000.0,        ! Cutoff
  dt=0.002,          ! The time step in picoseconds
  nstlim=5000000,    ! Number of MD steps to be performed
  ig=-1,             ! Random seed (=-1, get a number from current date and time)
  ntwr=10000,        ! Restart file written every ntwr steps
  ntwx=10000,        ! Trajectory file written every ntwx steps
  ntpr=10000,        ! The mdout and mdinfo files written every ntpr steps
  ioutfm=1,          ! Trajectory file format (=1, Binary NetCDF)
  igb=2,             ! GB model
  saltcon=0.1,       ! Salt concentration
  irest=1,           ! Flag to restart the simulation
  ntx=5,             ! Initial condition (=5, coord. and veloc. read from the inpcrd file)
  ntr=1,             ! Turn on positional restraints
  restraintmask='@CA,C,O,N|(:HEH,PRN)@FE,NA,NB,NC,ND,C3D,C2A,C3B,C2C,CA,CB',
  restraint_wt=0.1,  ! 0.1 kcal/mol.A**2 restraint force constant
/
Equilibration Stage in Explicit Solvent
&cntrl
  ntt=3,             ! Temperature scaling (=3, Langevin dynamics)
  gamma_ln=2.0,      ! Collision frequency of the Langevin dynamics in ps-1
  ntc=2,             ! SHAKE constraints (=2, hydrogen bond lengths constrained)
  ntf=2,             ! Force evaluation (=2, hydrogen bond interactions omitted)
  ntb=2,             ! Boundaries (=2, constant pressure)
  ntp=1,             ! Constant pressure dynamics (=1, isotropic position scaling)
  barostat=2,        ! Barostat to control the pressure (=2, Monte Carlo)
  pres0=1.0,         ! Target pressure (1 bar)
  taup=1.0,          ! Pressure relaxation time in ps
  cut=8.0,           ! Cutoff
  dt=0.002,          ! The time step in picoseconds
  nstlim=4000000,    ! Number of MD steps to be performed
  ig=-1,             ! Random seed (=-1, get a number from current date and time)
  ntwr=10000,        ! Restart file written every ntwr steps
  ntwx=10000,        ! Trajectory file written every ntwx steps
  ntpr=10000,        ! The mdout and mdinfo files written every ntpr steps
  ioutfm=1,          ! Trajectory file format (=1, Binary NetCDF)
  iwrap=1,           ! Translate water molecules into the original simulation box
  igb=0,             ! GB model (=0, explicit solvent)
  saltcon=0.1,       ! Salt concentration
  irest=1,           ! Flag to restart the simulation
  ntx=5,             ! Initial condition (=5, coord. and veloc. read from the inpcrd file)
  ntr=1,             ! Turn on positional restraints
  restraintmask='@CA,C,O,N|(:HEH,PRN)@FE,NA,NB,NC,ND,C3D,C2A,C3B,C2C,CA,CB',
  restraint_wt=0.1,  ! 0.1 kcal/mol.A**2 restraint force constant
/

For implicit solvent simulation, the command

pmemd.cuda -O -i equi.implicit.mdin -p mp8_is.prmtop -c mp8_is.heat.rst7 -x mp8_is.equi.nc -inf mp8_is.equi.mdinfo -o mp8_is.equi.mdout -r mp8_is.equi.rst7 -ref mp8_is.heat.rst7

will generate the files: mp8_is.equi.mdout, mp8_is.equi.nc and mp8_is.equi.rst7.

     

For explicit solvent simulation, the command

pmemd.cuda -O -i equi.explicit.mdin -p mp8_es.new.prmtop -c mp8_es.heat.rst7 -x mp8_es.equi.nc -inf mp8_es.equi.mdinfo -o mp8_es.equi.mdout -r mp8_es.equi.rst7 -ref mp8_es.heat.rst7

will generate the files: mp8_es.equi.mdout, mp8_es.equi.nc and mp8_es.equi.rst7.

These simulations may also be performed using pmemd, pmemd.MPI, sander or sander.MPI.

In explicit solvent we perform one additional equilibration at constant volume to let the system further relax at the new density. This simulation is executed for 20 ns and at 300 K. The input file for this simulation is:

equi2.explicit.mdin
Second Equilibration Stage in Explicit Solvent
&cntrl
  ntt=3,             ! Temperature scaling (=3, Langevin dynamics)
  gamma_ln=2.0,      ! Collision frequency of the Langevin dynamics in ps-1
  ntc=2,             ! SHAKE constraints (=2, hydrogen bond lengths constrained)
  ntf=2,             ! Force evaluation (=2, hydrogen bond interactions omitted)
  ntb=1,             ! Boundaries (=1, constant volume)
  cut=8.0,           ! Cutoff
  dt=0.002,          ! The time step in picoseconds
  nstlim=10000000,   ! Number of MD steps to be performed
  ig=-1,             ! Random seed (=-1, get a number from current date and time)
  ntwr=10000,        ! Restart file written every ntwr steps
  ntwx=10000,        ! Trajectory file written every ntwx steps
  ntpr=10000,        ! The mdout and mdinfo files written every ntpr steps
  ioutfm=1,          ! Trajectory file format (=1, Binary NetCDF)
  iwrap=1,           ! Translate water molecules into the original simulation box
  igb=0,             ! GB model (=0, explicit solvent)
  saltcon=0.1,       ! Salt concentration
  irest=1,           ! Flag to restart the simulation
  ntx=5,             ! Initial condition (=5, coord. and veloc. read from the inpcrd file)
  ntr=1,             ! Turn on positional restraints
  restraintmask='@CA,C,O,N|(:HEH,PRN)@FE,NA,NB,NC,ND,C3D,C2A,C3B,C2C,CA,CB',
  restraint_wt=0.1,  ! 0.1 kcal/mol.A**2 restraint force constant
/

The command

pmemd.cuda -O -i equi2.explicit.mdin -p mp8_es.new.prmtop -c mp8_es.equi.rst7 -x mp8_es.equi2.nc -inf mp8_es.equi2.mdinfo -o mp8_es.equi2.mdout -r mp8_es.equi2.rst7 -ref mp8_es.equi.rst7

will generate the files mp8_es.equi2.mdout, mp8_es.equi2.nc and mp8_es.equi2.rst7. This simulation may also be performed with pmemd, pmemd.MPI, sander or sander.MPI instead of using the GPU-accelerated code with pmemd.cuda.

Now we have the equilibrated restart files mp8_is.equi.rst7 and mp8_es.equi2.rst7 that will be used as input in the production simulations in the next section.


Click here to go to Section 3


Click here to go back to Section 1


(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Vinícius Cruzeiro. 2018