(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):
- Minimize system to relax bad contacts (jump here)
- Heat the system slowly to the target temperature (jump here)
- 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 |
For explicit solvent simulation, the command |
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 |
For explicit solvent simulation, the command |
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 |
For explicit solvent simulation, the command |
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 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