AmberTools23 Amber22 Manuals Tutorials Force Fields Contacts History

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

## 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.