Amber TUTORIAL B1: Simulating a DNA polyA-polyT Decamer

Section 1: Introduction

The input files required (using their default file names):

• prmtop - a file containing a description of the molecular topology and the necessary force field parameters.
• inpcrd (or a restrt from a previous run) - a file containing a description of the atom coordinates and optionally velocities and current periodic box dimensions.
• mdin - the sander input file consisting of a series of namelists and control variables that determine the options and type of simulation to be run.

The approximate order of this tutorial will be as follows:

1. Create the prmtop and inpcrd files: This is a description of how to generate the initial structure and set up the molecular topology/parameter and coordinate files necessary for performing minimization or dynamics with sander.
2. An introduction to minimization and molecular dynamics. Run short MD simulations in-vacuo. Perform basic analysis such as calculating root-mean-squared deviations (RMSd) and plotting various energy terms as a function of time. Visualizing results with VMD.
3. Minimization and molecular dynamics in implicit solvent: Setting up and running equilibration and production minimization and molecular dynamics simulations for our DNA model using the Born implicit solvent model.
4. Minimization and molecular dynamics in explicit solvent: Setting up and running equilibration and production simulations for our DNA model using TIP3P explicit water.

Section 2: Setting up duplex DNA: polyA-polyT

2.1) Generating the coordinates of the model structure

2.1.1) Creating our DNA duplex using NAB

構建NAB輸入文件nuc.nab

$$$ nuc.nab

molecule m; m = fd_helix( "abdna", "aaaaaaaaaa", "dna" ); putpdb( "nuc.pdb", m, "-wwpdb");

Running NAB

nab nuc.nab ./a.out (Note: type "./a.exe" for  Cygwin/Windows)

2.1.2) Loading the structure into Leap

如下命令能夠分步再xleap中或tleap中運行，也可把下列命令放到名爲 tleap.in 的文件中，用命令 tleap -f tleap.in 執行，下面命令會產生三種不一樣的模型，詳情看註釋內容

source leaprc.DNA.bsc1 source leaprc.water.tip3p dna1= loadpdb "nuc.pdb"
 # An in vacuo model of the poly(A)-poly(T) structure (named polyAT_vac) saveamberparm dna1 polyAT_vac.prmtop polyAT_vac.inpcrd
 # An in vacuo model of the poly(A)-poly(T) structure with explicit counterions (named polyAT_cio) addions dna1 Na+ 0 saveamberparm dna1 polyAT_cio.prmtop polyAT_cio.inpcrd
 # A TIP3P (water) solvated model of the poly(A)-poly(T) structure in a periodic box (named polyAT_wat) dna2 = copy dna1 solvatebox dna1 TIP3PBOX 8.0 solvateoct dna2 TIP3PBOX 8.0 saveamberparm dna2 polyAT_wat.prmtop polyAT_wat.inpcrd
 quit

其中，能夠運行如下命令增長對leap的理解：

edit dna1　　　　　　　　　　　　# 可在xleap狀態下查看目標結構 help loadpdb　　　　　　　　　　# loadpdb幫助文檔 help saveamberparm help addions help solvatebox help solvateoct list　　　　　　　　　　　　　　# 列舉當前力場下載入的全部原子，分子，離子及水模型

Section 3: Running Minimization and Molecular Dynamics (in vacuo)

This section of the tutorial will consist of 3 stages:

1. Relaxing the system prior to MD
2. Molecular dynamics at constant temperature
3. Analyzing the results

3.1) Relaxing the System Prior to MD

The basic usage for sander is as follows:

sander [-O] -i mdin -o mdout -p prmtop -c inpcrd -r restrt [-ref refc] [-x mdcrd] [-v mdvel] [-e mden] [-inf mdinfo]

• Arguments in []'s are optional
• If an argument is not specified, the default name will be used.
• -O    overwrite all output files (the default behavior is to quit if any output files already exist)
• -i      the name of the input file (which describes the simulation options), mdin by default.
• -o     the name of the output file, mdout by default.
• -p     the parameter/topology file, prmtop by default.
• -c     the set of initial coordinates for this run, inpcrd by default.
• -r     the final set of coordinates from this MD or minimization run, restrt by default.
• -ref  reference coordinates for positional restraints, if this option is specified in the input file, refc by default.
• -x    the molecular dynamics trajectory file (if running MD), mdcrd by default.
• -v    the molecular dynamics velocities file (if running MD), mdvel by default.
• -e    a summary file of the energies (if running MD), mden by default.
• -inf  a summary file written every time energy information is printed in the output file for the current step of the minimization of MD, useful for checking on the progress of a simulation, mdinfo by default.

3.1.1) Building the mdin input file

由於先前步驟產生的默認結構並非能量最優的，內部可能存在原子間的conflict和overlap，這可能使MD模擬變得不穩定，因此最好再MD以前先對初始結構進行優化(minimization)

$$$   polyAT_vac_init_min.in

polyA-polyT 10-mer: initial minimization prior to MD &cntrl imin = 1, maxcyc = 500, ncyc = 250, ntb = 0, igb = 0, cut = 12
 /

上面爲指定sander輸入的namelist 變量文件，除最開始的說明信息外，此文件一般以&cntrl或&ewald開頭，以「/」結尾，下面一段英文是對polyAT_vac_init_min.in文件的論述；

　　To turn on minimization, we specify IMIN = 1. We want a fairly short minimization since we don't actually need to reach the minima, but just want to move away from any local maxima, so we select 500 steps of minimization by specifying MAXCYC = 500. Sandersupports a number of different minimization algorithms, the most commonly used being steepest descent and conjugate gradient. The steepest descent algorithm is good for quickly removing the largest strains in the system but it also converges slowly when close to a minima. Here the conjugate gradient method is more efficient. The use of these two algorithms can be controlled using the NCYC flag. If NCYC < MAXCYC, sander will use the steepest descent algorithm for the first NCYC steps before switching to the conjugate gradient algorithm for the remaining (MAXCYC - NCYC) steps. In this case we will run an equal number of steps with each algorithm so we set NCYC = 250. Since sander assumes that the system is periodic by default we need to explicitly turn this off (NTB = 0). In this simulation we will be using a constant dielectric and not an implicit (or explicit) solvent model so we set IGB = 0 (no generalized born solvation model). This is the default so we  don't strictly need to specify this, but I will include it here so that we can see what differences we have in the input file when we switch on implicit solvent later. We also need to choose a value for the non-bonded cut off. A larger cut off introduces less error in the non-bonded force evaluation but increases the computational complexity and thus calculation time. 12 angstroms would seem like a reasonable tradeoff for gas phase so that is what we will initially use (CUT = 12).

3.1.2) Running sander for the first time

指定好namelist文件後，執行sander：

sander -O -i polyAT_vac_init_min.in -o polyAT_vac_init_min.out -c polyAT_vac.inpcrd -p polyAT_vac.prmtop -r polyAT_vac_init_min.rst

Input files: polyAT_vac_init_min.in, polyAT_vac.inpcrd, polyAT_vac.prmtop
Output files: polyAT_vac_init_min.out, polyAT_vac_init_min.rst

3.1.3) Creating PDB files from the AMBER coordinate files

使用ambpdb命令能夠從rst文件中提取pdb信息：

ambpdb -p polyAT_vac.prmtop -c polyAT_vac_init_min.rst > polyAT_vac_init_min.pdb

3.2) Running MD in-vacuo (100ps)

經過比較兩種不一樣的cut值，咱們能夠更加明確cut的意義：

$$$ polyAT_vac_md1_12Acut.in

10-mer DNA MD in-vacuo, 12 angstrom cut off &cntrl imin = 0, ntb = 0, igb = 0, ntpr = 100, ntwx = 100, ntt = 3, gamma_ln = 1.0, tempi = 300.0, temp0 = 300.0 nstlim = 100000, dt = 0.001, cut = 12.0
 /

 

$$$ polyAT_vac_md1_nocut.in

10-mer DNA MD in-vacuo, infinite cut off &cntrl imin = 0, ntb = 0, igb = 0, ntpr = 100, ntwx = 100, ntt = 3, gamma_ln = 1.0, tempi = 300.0, temp0 = 300.0 nstlim = 100000, dt = 0.001, cut = 999
 /

 To run a molecular dynamics simulation with sander, we need to turn off minimization (IMIN=0). Since we are running in-vacuo we also need to disable periodicity (NTB=0) and set IGB=0 since we are not using implicit solvent. For these two examples we will write information to the output file and trajectory coordinates file every 100 steps (NTPR=100, NTWX=100). The CUT variable specifies the cut off range for the long range non-bonded interactions. Here, two different values for the cutoff will be used, one run with a cut off of 12 angstroms (CUT = 12.0) and one run without a cutoff. To run without a cutoff we simply set CUT to be larger than the extent of the system (e.g. CUT = 999). For temperature regulation we will use the Langevin thermostat (NTT=3) to maintain the temperature of our system at 300 K. This temperature control method uses Langevin dynamics with a collision frequency given by GAMMA_LN and is significantly more efficient at equilibrating the system temperature than the Berendsen temperature coupling scheme (NTT=1). For simplicity we will run this entire simulation with NTT=3 and GAMMA_LN=1. We shall also set the initial and final temperatures to 300 K (TEMPI=300.0, TEMP0=300.0) which will mean our system's temperature should remain around 300 K. In these two examples we will run a total of 100,000 steps each (NSTLIM=100000) with a 1 fs time step (DT=0.001) giving simulation lengths of 100 ps (100,000 x 1fs).

運行真空狀態下MD（也可用sander.MP或pmemd或pmemd.MPI）：

CUT=12.0:

sander -O -i polyAT_vac_md1_12Acut.in -o polyAT_vac_md1_12Acut.out -c polyAT_vac_init_min.rst -p polyAT_vac.prmtop -r polyAT_vac_md1_12Acut.rst -x polyAT_vac_md1_12Acut.mdcrd

CUT=999:

sander -O -i polyAT_vac_md1_nocut.in -o polyAT_vac_md1_nocut.out -c polyAT_vac_init_min.rst -p polyAT_vac.prmtop -r polyAT_vac_md1_nocut.rst -x polyAT_vac_md1_nocut.mdcrd

實時觀測程序輸出：

tail -f polyAT_vac_md1_12Acut.out

3.3) Analyzing the results

3.3.1) Extracting the energies, etc. from the mdout file

使用process_mdout.perl腳本，從mdout文件中提取能量勢能等信息：

mkdir polyAT_vac_md1_12Acut cd polyAT_vac_md1_12Acut
process_mdout.perl ../polyAT_vac_md1_12Acut.out

 

mkdir polyAT_vac_md1_nocut cd polyAT_vac_md1_nocut process_mdout.perl ../polyAT_vac_md1_nocut.out

用grace對各個能量進行繪圖：

xmgrace ./polyAT_vac_md1_12Acut/summary.EPTOT ./polyAT_vac_md1_nocut/summary.EPTOT

 圖表以下：

 

3.3.2) Calculating the RMSd vs. time

計算12Acut RMSD變化狀況：

trajin polyAT_vac_md1_12Acut.mdcrd rms first mass out polyAT_vac_md1_12Acut.rms time 0.1

計算nocut RMSD變化狀況：

trajin polyAT_vac_md1_nocut.mdcrd rms first mass out polyAT_vac_md1_nocut.rms time 0.1

where trajin specifies the name of the trajectory file to process, rms first mass tells cpptraj that we want it to calculate a MASS weighted, RMS fit to the FIRST structure. out specifies the name of the file to write the results to and time 0.1 tells cpptraj that each frame of the coordinate file represents 0.1 ps. Since we have two simulations we have two input files, with different trajectory and output files in each.

使用cpptraj程序：

$AMBERHOME/bin/cpptraj -p polyAT_vac.prmtop -i polyAT_vac_md1_12Acut.calc_rms $AMBERHOME/bin/cpptraj -p polyAT_vac.prmtop -i polyAT_vac_md1_nocut.calc_rms

使用grace做圖：

>xmgrace polyAT_vac_md1_12Acut.rms polyAT_vac_md1_nocut.rms

 圖表以下：

 

Section 4: Running Minimization and MD (in implicit solvent)

與第三節相同的是，在正式作MD以前須要對體系進行優化，不一樣的是，這節的運行環境再也不是真空(vacuo)，而是在隱式溶劑模型中(implicit solvent model)，這種模型相對於真空更加與真實環境相近，同時運算量也要大不少，效果比真空環境好，可是不如準確溶劑模型(explicit solvent model)，可是explicit solvent計算資源耗費更大。

Use of the Born solvation model is controlled by the IGB flag. The default value for this flag is 0 which corresponds to no generalized Born term being used. Alternative values are IGB=1 which corresponds to the "standard" pairwise generalized Born model using the default radii set up by LEaP (the method we will be using here), IGB=2 a modified GB model developed by A. Onufriev, D. Bashford and D.A.Case and IGB=5 which is essentially the same as IGB=2 but with alternative parameters. For further details see the AMBER manual.

4.1) Relaxing the System Prior to MD

generalized Born method (IGB=1), CUT=12.0

$$$ polyAT_gb_init_min.in

polyA-polyT 10-mer: initial minimization prior to MD GB model &cntrl imin = 1, maxcyc = 500, ncyc = 250, ntb = 0, igb = 1, cut = 12
 /

運行sander(也可以使用pmemd):

sander -O -i polyAT_gb_init_min.in -o polyAT_gb_init_min.out -c polyAT_vac.inpcrd -p polyAT_vac.prmtop -r polyAT_gb_init_min.rst

分析優化後pdb結構：

ambpdb -p polyAT_vac.prmtop -c polyAT_gb_init_min.rst > polyAT_gb_init_min.pdb

4.2) Running MD with generalized Born solvation

$$$   polyAT_gb_md1_12Acut.in　　　　　　　　　　　　　　# 12.0 angstrom long range cutoff, Generalized Born

10-mer DNA MD Generalized Born, 12 angstrom cut off &cntrl imin = 0, ntb = 0, igb = 1, ntpr = 100, ntwx = 100, ntt = 3, gamma_ln = 1.0, tempi = 300.0, temp0 = 300.0 nstlim = 100000, dt = 0.001, cut = 12.0
 /

 

$$$ polyAT_gb_md1_nocut.in　　　　　　　　　　　　　　# no long range cutoff, Generalized Born

10-mer DNA MD Generalized Born, infinite cut off &cntrl imin = 0, ntb = 0, igb = 1, ntpr = 100, ntwx = 100, ntt = 3, gamma_ln = 1.0, tempi = 300.0, temp0 = 300.0 nstlim = 100000, dt = 0.001, cut = 999
 /

We are using the same settings as before, namely we turn off minimization (IMIN=0). We disable periodicity (NTB=0) but this time set IGB=1 since we want to use the Born implicit solvent model. We will again write information to the output file and trajectory coordinates file every 100 steps (NTPR=100, NTWX=100). Two different values for the cutoff will be used, one run will be with a cutoff of 12 angstroms (CUT = 12.0) and one run will be without a cutoff (CUT = 999). For temperature regulation we will use the Langevin thermostat (NTT=3, GAMMA_LN=1) to maintain the temperature of our system at 300 K.

We shall also set the initial and final temperatures to 300 K (TEMPI=300.0, TEMP0=300.0) which will mean our system's temperature should remain around 300 K. Again we will run the two examples for a total of 100,000 steps each (NSTLIM=100000) with a 1 fs time step (DT=0.001) giving simulation lengths of 100 ps (100,000 x 1fs).

As mentioned in the previous section when running production simulations with ntt=2 or 3 you should always change the random seed (ig) between restarts. You can achieve this automatically by setting ig=-1 in the ctrl namelist. We do not do this here in the tutorial for reproducibility but it is generally recommended that you do this in your calculations. 

 

運行100ns MD（也可用sander.MP或pmemd或pmemd.MPI）：

sander -O -i polyAT_gb_md1_12Acut.in -o polyAT_gb_md1_12Acut.out -c polyAT_gb_init_min.rst -p polyAT_vac.prmtop -r polyAT_gb_md1_12Acut.rst -x polyAT_gb_md1_12Acut.mdcrd

 

sander -O -i polyAT_gb_md1_nocut.in -o polyAT_gb_md1_nocut.out -c polyAT_gb_init_min.rst -p polyAT_vac.prmtop -r polyAT_gb_md1_nocut.rst -x polyAT_gb_md1_nocut.mdcrd

實時觀測程序輸出：

tail -f polyAT_gb_md1_12Acut.out

4.3) Analyzing the results

使用process_mdout.perl腳本，從mdout文件中提取能量勢能等信息：

mkdir polyAT_gb_md1_12Acut mkdir polyAT_gb_md1_nocut cd polyAT_gb_md1_12Acut process_mdout.perl ../polyAT_gb_md1_12Acut.out cd ../polyAT_gb_md1_nocut process_mdout.perl ../polyAT_gb_md1_nocut.out cd .. xmgrace ./polyAT_gb_md1_12Acut/summary.EPTOT ./polyAT_gb_md1_nocut/summary.EPTOT

結果如圖:

計算RMSD變化狀況：

$$$ polyAT_gb_md1_12Acut.calc_rms trajin polyAT_gb_md1_12Acut.mdcrd rms first mass out polyAT_gb_md1_12Acut.rms time 0.1

 

$$$ polyAT_gb_md1_nocut.calc_rms trajin polyAT_gb_md1_nocut.mdcrd rms first mass out polyAT_gb_md1_nocut.rms time 0.1

 使用cpptraj分析RMSD變化：

$AMBERHOME/bin/cpptraj -p polyAT_vac.prmtop -i polyAT_gb_md1_12Acut.calc_rms $AMBERHOME/bin/cpptraj -p polyAT_vac.prmtop -i polyAT_gb_md1_nocut.calc_rms
xmgrace polyAT_gb_md1_12Acut.rms polyAT_gb_md1_nocut.rms

 grace做圖：

 

Section 5: Running Minimization and MD (in explicit solvent)

5.1) Equilibrating and running MD with solvated poly(A)-poly(T)

由於準確溶劑化模型的緣由，這次對模型的優化分爲兩個階段，第一階段保持溶質固定，優化溶劑及離子；第二階段優化整個體系：

5.1.1) Minimization Stage 1 - Holding the solute fixed

$$$ polyAT_wat_min1.in polyA-polyT 10-mer: initial minimization solvent + ions &cntrl imin = 1, maxcyc = 1000, ncyc = 500, ntb = 1, ntr = 1, cut = 10.0
 / Hold the DNA fixed
500.0 RES 1 20 END END

The meaning of each of the terms are as follows:

• IMIN = 1: Minimization is turned on (no MD)
• MAXCYC = 1000: Conduct a total of 1,000 steps of minimization.
• NCYC = 500: Initially do 500 steps of steepest descent minimization followed by 500 steps (MAXCYC - NCYC) steps of conjugate gradient minimization.
• NTB = 1: Use constant volume periodic boundaries (PME is always "on" when NTB > 0).
• CUT = 10.0: Use a cutoff of 10 angstroms. (This is probably overkill here, we could get away with 8.0 or 9.0 angstroms if we wanted but using a larger cutoff does not harm the accuracy of the results. It just makes the calculations slightly more costly.)
• NTR = 1: Use position restraints based on the GROUP input given in the input file. GROUP input is described in detail in the appendix of the user's manual. In this example, we use a force constant of 500 kcal mol^-1 angstrom^-2 and restrain residues 1 through 20 (the DNA). This means that the water and counterions are free to move. The GROUP input is the last 5 lines of the input file:

• Hold the DNA fixed
  500.0 RES 1 20 END END

  Note that whenever you run using the GROUP option in the input you should carefully check the top of the output file to make sure you've selected as many atoms as you thought you did. (Note also that 500 kcal/mol-A**2 is a very large force constant, much larger than is needed. You can use this for minimization, but for dynamics, stick to much smaller values, say around 10.)

運行sander(pmemd)進行第一階段優化：

sander -O -i polyAT_wat_min1.in -o polyAT_wat_min1.out -p polyAT_wat.prmtop -c polyAT_wat.inpcrd -r polyAT_wat_min1.rst -ref polyAT_wat.inpcrd

5.1.2) Minimization Stage 2 - Minimizing the entire system

$$$ polyAT_wat_min2.in polyA-polyT 10-mer: initial minimization whole system &cntrl imin = 1, maxcyc = 2500, ncyc = 1000, ntb = 1, ntr = 0, cut = 10.0
 /

運行sander(pmemd)進行第二階段優化：

sander -O -i polyAT_wat_min2.in -o polyAT_wat_min2.out -p polyAT_wat.prmtop -c polyAT_wat_min1.rst -r polyAT_wat_min2.rst

提出pdb，查看結構差異：

$AMBERHOME/bin/ambpdb -p polyAT_wat.prmtop < polyAT_wat.inpcrd > polyAT_wat.pdb $AMBERHOME/bin/ambpdb -p polyAT_wat.prmtop < polyAT_wat_min2.rst > polyAT_wat_min2.pdb

注：ambpdb出現Error: Could not read restart atoms/time.時，命令變動爲 ambpdb -p polyAT_wat.prmtop -c polyAT_wat_min2.rst > polyAT_wat_min2.pdb .

5.1.3) Molecular Dynamics (heating) with restraints on the solute

$$$ polyAT_wat_md1.in polyA-polyT 10-mer: 20ps MD with res on DNA &cntrl imin = 0, irest = 0, ntx = 1, ntb = 1, cut = 10.0, ntr = 1, ntc = 2, ntf = 2, tempi = 0.0, temp0 = 300.0, ntt = 3, gamma_ln = 1.0, nstlim = 10000, dt = 0.002 ntpr = 100, ntwx = 100, ntwr = 1000
 / Keep DNA fixed with weak restraints 10.0 RES 1 20 END END

As mentioned in Section 3, when running production simulations with ntt=2/3 you should always change the random seed (ig) between restarts. If you are using AMBER 10 (bugfix.26 or later) or AMBER 11 or later then you can achieve this automatically by setting ig=-1 in the ctrl namelist. We do not do this here in the tutorial for reproducibility but it is generally recommended that you do this in your calculations.

The meaning of each of the terms are as follows:

• IMIN = 0: Minimization is turned off (run molecular dynamics)
• IREST = 0, NTX = 1: We are generating random initial velocities from a Boltzmann distribution and only read in the coordinates from the inpcrd. In other words, this is the first stage of our molecular dynamics. Later we will change these values to indicate that we want to restart a molecular dynamics run from where we left off.
• NTB = 1: Use constant volume periodic boundaries (PME is always "on" when NTB>0).
• CUT = 10.0: Use a cutoff of 10 angstroms.
• NTR = 1: Use position restraints based on the GROUP input given in the input file. In this case we will restrain the DNA with a force constant of 10 angstroms.
• NTC = 2, NTF = 2: SHAKE should be turned on and used to constrain bonds involving hydrogen.
• TEMPI = 0.0, TEMP0 = 300.0: We will start our simulation with a temperature, derived from the kinetic energy, of 0 K and we will allow it to heat up to 300 K. The system should be maintained, by adjusting the kinetic energy, as 300 K.
• NTT = 3, GAMMA_LN = 1.0: The Langevin dynamics should be used to control the temperature using a collision frequency of 1.0 ps^-1.
• NSTLIM = 10000, DT = 0.002: We are going to run a total of 10,000 molecular dynamics steps with a time step of 2 fs per step, possible since we are now using SHAKE, to give a total simulation time of 20 ps.
• NTPR = 100, NTWX = 100, NTWR = 1000: Write to the output file (NTPR) every 100 steps (200 fs), to the trajectory file (NTWX) every 100 steps and write a restart file (NTWR), in case our job crashes and we want to restart it, every 1,000 steps.

 運行sander(pmemd, pmemd.MPI etc.)進行限制性pre-MD (heating):

sander -O -i polyAT_wat_md1.in -o polyAT_wat_md1.out -p polyAT_wat.prmtop -c polyAT_wat_min2.rst -r polyAT_wat_md1.rst -x polyAT_wat_md1.mdcrd -ref polyAT_wat_min2.rst

 

5.1.6) Running MD Equilibration on the whole system

$$$ polyAT_wat_md2.in polyA-polyT 10-mer: 100ps MD &cntrl imin = 0, irest = 1, ntx = 7, ntb = 2, pres0 = 1.0, ntp = 1, taup = 2.0, cut = 10.0, ntr = 0, ntc = 2, ntf = 2, tempi = 300.0, temp0 = 300.0, ntt = 3, gamma_ln = 1.0, nstlim = 50000, dt = 0.002, ntpr = 100, ntwx = 100, ntwr = 1000
 /

The meaning of each of the terms are as follows:

• IMIN = 0: Minimization is turned off (run molecular dynamics)
• IREST = 1, NTX = 7: We want to restart our MD simulation where we left off after the 20 ps of constant volume simulation. IREST tells sander that we want to restart a simulation, so the time is not reset to zero but will start at 20 ps. Previously we have had NTX set at the default of 1 which meant only the coordinates were read from the restrt file. This time, however, we want to continue from where we finished so we set NTX = 7 which means the coordinates, velocities and box information will be read from a formatted (ASCII) restrt file.
• NTB = 2, PRES0 = 1.0, NTP = 1, TAUP = 2.0: Use constant pressure periodic boundary with an average pressure of 1 atm (PRES0). Isotropic position scaling should be used to maintain the pressure (NTP=1) and a relaxation time of 2 ps should be used (TAUP=2.0).
• CUT = 10.0: Use a cut off of 10 angstroms.
• NTR = 0: We are no longer using positional restraints.
• NTC = 2, NTF = 2: SHAKE should be turned on and used to constrain bonds involving hydrogen.
• TEMPI = 300.0, TEMP0 = 300.0: Our system was already heated to 300 K during the first stage of MD so here it will start at 300 K and should be maintained at 300 K.
• NTT = 3, GAMMA_LN = 1.0: The Langevin dynamics should be used to control the temperature using a collision frequency of 1.0 ps^-1.
• NSTLIM = 50000, DT = 0.002: We are going to run a total of 50,000 molecular dynamics steps with a time step of 2 fs per step, possible since we are now using SHAKE, to give a total simulation time of 100 ps.
• NTPR = 100, NTWX = 100, NTWR = 1000: Write to the output file (NTPR) every 100 steps (200 fs), to the trajectory file (NTWX) every 100 steps and write a restart file (NTWR), in case our job crashes and we want to restart it, every 1,000 steps.

 對整個體系進行MD running (sander, sander.MPI, pmemd, pmemd.MPI)

sander -O -i polyAT_wat_md2.in -o polyAT_wat_md2.out -p polyAT_wat.prmtop -c polyAT_wat_md1.rst -r polyAT_wat_md2.rst -x polyAT_wat_md2.mdcrd

 

5.2) Analyzing the results to test the equilibration

5.2.1) Analyzing the output files

mkdir analysis cd analysis process_mdout.perl ../polyAT_wat_md1.out ../polyAT_wat_md2.out xmgrace summary.EPTOT summary.EKTOT summary.ETOT
xmgrace summary.TEMP
xmgrace summary.PRES

 

vi summary.VOLUME d100 (removes the current (first) line plus the next 100 lines) <esc>:wq summary.VOLUME.modified xmgrace summary.VOLUME.modified

5.2.2) Analyzing the trajectory

trajin polyAT_wat_md1.mdcrd trajin polyAT_wat_md2.mdcrd rms first out polyAT_wat_backbone.rms @P,O3',O5',C3',C4',C5' time 0.2
 $AMBERHOME/bin/cpptraj -p polyAT_wat.prmtop -i polyAT_wat_calc_backbone_rms.in xmgrace polyAT_wat_backbone.rms