(5)在带有周期性界限的显式的溶剂中运行最小化和MD：
a.在poly(A)－poly(T)下平衡并运行MD
注：理想的截断值在8－10埃，对于有周期边界的截断值不应小于8埃。
最小化 1——溶质分子的固定（仅仅最小化水和离子）
有两种方法："BELLY"和"positional restraints"，文章中应用的是后一种。
建立输入文件：polyAT_wat_min1.in
polyA-polyT 10-mer: initial minimization solvent + ions
&cntrl
imin = 1, //IMIN = 1: Minimization is turned on (no MD)
maxcyc = 1000, //MAXCYC = 1000: Conduct a total of 1,000 steps of minimization.
ncyc = 500, //NCYC = 500: Initially do 500 steps of steepest descent minimization followed by 500 steps (MAXCYC - NCYC) steps of conjugate gradient minimization.
ntb = 1, //NTB = 1: Use constant volume periodic boundaries (PME is always "on" when NTB>0).
ntr = 1, //NTR = 1: Use position restraints based on the GROUP input given in the input file.
cut = 10.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.)
/
Hold the DNA fixed
500.0 //(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.)
RES 1 20
END
END
运行：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
Input files: polyAT_wat_min1.in, polyAT_wat.prmtop, polyAT_wat.inpcrd
Output files: polyAT_wat_min1.out, polyAT_wat_min1.rst
最小化2－－整个系统的最小化: （这步要用到上一步的最小化的文件：polyAT_wat_min1.rst ）
在这里没有时间的限制运行2500步。输入文件：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 -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
Input files: polyAT_wat_min2.in, polyAT_wat.prmtop, polyAT_wat_min1.rst
Output files: polyAT_wat_min2.out, polyAT_wat_min2.rst
产生两个PDB文件:一个只有水，另一个是最小化的文件:
ambpdb -p polyAT_wat.prmtop < polyAT_wat.inpcrd > polyAT_wat.pdb
ambpdb -p polyAT_wat.prmtop < polyAT_wat_min2.rst > polyAT_wat_min2.pdb
用VMD观看.pdb
vmd
File -> New Molecule (load the first pdb polyAT_wat.pdb)
We can make the DNA stand out from the water better:
Graphics -> Representations
Hit "Create Rep" This will create a second representation of our molecule
In the "Selected Atoms" box type: all not water
Then under "Drawing Method" select "Ribbons"
Next select the representation that still has the style "Lines"
And change the drawing method for this to "Points"
拟合两个.pdb
vmd polyAT_wat.pdb
File -> New Molecule Load the polyAT_wat_min2.pdb file
Turn off the water to make it easier to compare the DNA's:
Select Graphics->Representations
In the Selected Atoms box enter: all not water
Click "Apply"
Now under the Selected Molecule drop down box at the top of the Graphical Representations window select molecule 0: polyAT_wat.pdb and repeat the process above.
在溶质上加热的MD：
equilibration protocol is allow our system to heat up from 0 K to 300 K
NTT=1 维持温度不变；NTT=3 which is significantly better at maintaining and equalizing the system temperature
在阶段一位置限制被用作输入文件：polyAT_wat_md1.in
polyA-polyT 10-mer: 20ps MD with res on DNA
&cntrl
imin = 0, //IMIN = 0: Minimization is turned off (run molecular dynamics)
irest = 0, ntx = 1, //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, //NTB = 1: Use constant volume periodic boundaries (PME is always "on" when NTB>0).
cut = 10.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, //NTC = 2, NTF = 2: SHAKE should be turned on and used to constrain bonds involving hydrogen.
ntf = 2,
tempi = 0.0, //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.
temp0 = 300.0,
ntt = 3, //NTT = 3, GAMMA_LN = 1.0: The langevin dynamics should be used to control the temperature using a collision frequency of 1.0 ps
gamma_ln = 1.0,
nstlim = 10000, dt = 0.002 //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 //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.
/
Keep DNA fixed with weak restraints
10.0
RES 1 20
END
END
运行：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
Input files: polyAT_wat_md1.in, polyAT_wat.prmtop, polyAT_wat_min2.rst
Output files: polyAT_wat_md1.out,polyAT_wat_md1.rst,polyAT_wat_md1.mdcrd.gz 平衡整个系统MD：
polyAT_wat_md2.in
polyA-polyT 10-mer: 100ps MD
&cntrl
imin = 0, //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 2ps should be used (TAUP=2.0).
cut = 10.0, ntr = 0, //CUT = 10.0: Use a cut off of 10 angstroms.NTR = 0: We are no longer using positional restraints.
ntc = 2, ntf = 2, //NTC = 2, NTF = 2: SHAKE should be turned on and used to constrain bonds involving hydrogen.
tempi = 300.0, temp0 = 300.0, //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, //NTT = 3, GAMMA_LN = 1.0: The Langevin dynamics should be used to control the temperature using a collision frequency of 1.0 ps
nstlim = 50000, dt = 0.002, //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 //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.
/
用上一阶段的限制文件：polyAT_wat_md1.rst
开始：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
Input files: polyAT_wat_md2.in, polyAT_wat.prmtop, polyAT_wat_md1.rst Output files: polyAT_wat_md2.out,polyAT_wat_md2.rst,(polyAT_wat_md2.mdcrd
(6)练习：ADNA
a.建立结构：
a-dna.nab写入：
molecule m;
m = fd_helix( "adna", "aaaaaaaaaa", "dna" );
putpdb( "a-dna.nab.pdb", m, "-wwpdb");
Run nab to create our pdb file:
$AMBEROME/exe/nab a-dna.nab
./a.out
产生文件：a-dna.nab.pdb
b.创建sander的输入文件：
运行：xleap
>source leaprc.ffSB99
>adna = loadpdb a-dna.nab.pdb
中性化系统：
>addions adna Na+ 0
保存中性化的文件：
>saveamberparm adna a-dna_charges_only.prmtop a-dna_charges_only.inpcrd
输出文件： a-dna_charges_only.prmtop, a-dna_charges_only.inpcrd
给系统加8埃buffer的溶剂：
>solvateoct adna TIP3PBOX 8.0
>edit adna
保存加水的文件：
>saveamberparm adna a-dna_wat.prmtop a-dna_wat.inpcrd
c.最小化系统
Stage 1: 500 steps of steepest descent followed by 500 steps of conjugate gradient minimization, with a 500 kcal/mol restraint force on the DNA molecule (Res 1-20)
a-dna_min1.in (固定DNA)
A-DNA 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
sander -O -i a-dna_min1.in -o a-dna_min1.out -p a-dna_wat.prmtop -c a-dna_wat.inpcrd -r a-dna_min1.rst -ref a-dna_wat.inpcrd
Input files: a-dna_min1.in, a-dna_wat.prmtop, a-dna_wat.inpcrd
Output files: a-dna_min1.out,a-dna_min1.rst
Stage 2: 1,000 steps of steepest descent followed by 1,500 steps of conjugate gradient minimization with no restraints
a-dna_min2.in （用到上一步的 a-dna_min1.rst ）
A-DNA 10-mer: initial minimization solvent + ions
&cntrl
imin = 1,
maxcyc = 2500,
ncyc = 1000,
ntb = 1,
ntr = 0,
cut = 10.0
/
sander -O -i a-dna_min2.in -o a-dna_min2.out -p a-dna_wat.prmtop -c a-dna_min1.rst -r a-dna_min2.rst
Input files: a-dna_min2.in, a-dna_wat.prmtop, a-dna_min1.rst
Output files: a-dna_min2.out,a-dna_min2.rst
d.初始平衡
The next stage is to start running molecular dynamics. We will initially run 20 ps of heating/equilibration at constant volume in the same fashion as we did for the B-DNA explicit water simulation. This will be done using the same protocols:
i.e. 20 ps at constant volume periodic boundaries with weak (10 kcal/mol) restraints on the DNA. An initial temperature of 0 K, final temperature of 300 K, Langevin temperature controls, SHAKE constraints on hydrogen atoms and a 2 fs time step.
To save disk space we will only write to our output and mdcrd files every 250 steps (0.5 ps). We will also compress our mdcrd file upon completion using gzip. This will give us a compression ratio of about 60 %. This is important because the next step
will be to run 1.8 ns of MD which can create some very big mdcrd files:
a-dna_md1.in
A-DNA 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 = 250, ntwx = 250, ntwr = 10000
/
Keep DNA fixed with weak restraints
10.0
RES 1 20
END
END
运行：sander -O -i a-dna_md1.in -o a-dna_md1.out -p a-dna_wat.prmtop -c a-dna_min2.rst -r a-dna_md1.rst -ref a-dna_min2.rst -x a-dna_md1.mdcrd
Input files: a-dna_md1.in, a-dna_wat.prmtop, a-dna_min2.rst
Output files: a-dna_md1.out,a-dna_md1.rst, a-dna_md1.mdcrd.gz
e.进一步平衡：
In the next stage of our simulation we are going to run a further 1.8 ns of MD. Rather than run a separate equilibration and production phase we will, for convenience, combine the two since our production conditions
will be identical to our final equilibration conditions. Here we will use the same conditions as we did for the B-DNA explicit solvent equilibration. Namely we will run without restraints on the DNA. We shall maintain
the temperature at 300 K using the Langevin thermostat, run with constant pressure (1 atm) periodic boundaries, use SHAKE constraints on the hydrogen atoms and a 2 fs time step.
Since this simulation will take a long time to run it is probably not a good idea to run all 1.8 ns in one shot. Thus we will actually run 9 x 200 ps simulations with each successive simulation continuing on from the previous one.
In this way if a job crashes we only lose the progress of that job and can restart from where we got to. The input file will be the same in each case:
a-dna_md_1800ps.in
A-DNA 10-mer: 200ps of 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 = 100000, dt = 0.002,
ntpr = 250, ntwx = 250, ntwr = 10000
/
We will run the 9 simulations one after the other using the restrt file from the previous run as the input for the next run. While you could run each of these jobs manually using a command line as below
sander -O -i a-dna_md_1800ps.in -o a-dna_md2.out -p a-dna_wat.prmtop -c a-dna_md1.rst -r a-dna_md2.rst -x a-dna_md2.mdcrd
因为要运行循环的运行9次,所以最好写一个脚本：