/6PTI.pdb
solvate.sh (if want to include cystal water, use: "-cyrstalwater", if want ssbond, use: "-ssbond")
minCHARMM.pl 6pti.pdb>6pti_min.pdb
on sierra:
edit hosts.txt
172.16.0.1 6 /home/zwei/work/MD/BPTI/6PTI/REX/VA
....
dont use 172.16.0.9
....
and start.txt
6pti_min.txt
6pti_min.txt
....
native structure as starting structures
aarex.pl -n 1000 -hosts hosts.txt -par archive,natpdb=6pti.pdb -temp 96:300:850 -log rex.log -charmmlog charmm.log -elog energy.log -dir varex -f start.txt&
create trajectories
trajarchive2dcd.pl -inx 1:1000 -ref 6pti_min.pdb -dir varex
load into vmd
vmd -e ~/toolset/loadrex.tcl -args varex
create RMSD plot
rmsd_grid.sh varex
6pti.pdb was copied from /BPTI/6PTI/6pti.pdb
genseq.pl -out one 6pti.pdb>6pti.one
genini.sh 6pti.one 32 xxx(genini.sh blocked the termini)
if any one of them has energy greater than 0, replace them with backup ones, 32,31,30..
in one word make sure we have good initial structures from 1 to 30.
Replica Exchange simulation was performed on medusa:
aarex.pl -n 1000 -par archive,natpdb=6pti.pdb -temp 30:300:850 -log rex.log -charmmlog charmm.log -elog energy.log -mdpar cmap,blocked -dir gbrex_va 6pti.{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30}.pdb
RMSD_ENER.pl 1:1000 30
rmsd_ener.txt will be created under current directory:
(XXX_aa1_1.pdb) (rmsd) (potential energy)
(XXX_aa1_2.pdb) (rmsd) (potential energy)
.....
(XXX_aa30_1000.pdb) (rmsd) (potential energy)
eig.pl
eigplot_rmsd.txt and eigresult_rmsd.txt will be created.
gnuplot
gnuplot>set dgrid3d 100,100
gnuplot>set pm3d
gnuplot>set hidden3d
gnuplot>set contour both
gnuplot>set title "6pti_temp_rex"
gnuplot>set xlabel "PC1"
gnuplot>set ylabel "PC2"
gnuplot>set zlabel "Potential Energy(kcal/mol)" rotate by 90
gnuplot>set cntrparam levels incr -800,50,0
gnuplot>splot "eigplot_rmsd.txt" u 2:3:4 w lines notitle
set terminal postscript eps enhanced color size 3.5,2.62
set output "bpti_landscape.eps"
1D Umbrella replica exchange simulation
cond_create.sh
will create ./cond_resd
bias resd
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=57,atom2=CA
300 force=1.0,target=15,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=57,atom2=CA
......
300 force=1.0,target=205,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=57,atom2=CA
300 force=1.0,target=215,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=57,atom2=CA
then:
aarex.pl -n 100 -condfile cond_resd -par archive,natpdb=6pti.pdb -log rex.log -charmmlog charmm.log -elog energy.log -mdpar cmap,blocked -dir gbrex_va_resd 6pti.pdb
trajarchive2dcd.pl -inx 1:100 -ref 6pti.pdb -dir gbrex_va_resd
vmd -e ~/toolset/loadrex.tcl -ags gbrex_va_resd drawcell
NOTE:
in temperature replica exchange simulation, the format of energy.log under ./aa$ is
prod/0/1:322.340000: 0 0.00 328.55 -394.98 731.56 -1126.54 -9.48 -1586.33 0.00 0.00 0.00 0.000000 0.000000
but in conditional replica exchange simulation the format of energy.log turn into:
prod/0/1:300.000000:15: 0 0.00 301.14 -978.93 670.54 -1649.47 -182.40 -1904.98 0.00 0.00 0.00 0.000000 0.000000
RMSD_ENER.pl 1:100 25
resd_rmsd_ener.txt will be created.
wham1d analysis:
~/softwares/wham_2.0.9/wham/wham 5 215 400 0.0001 300 0 whaminput1d whamoutput1d
2D Umbrella replica exchange simulation
This time I add a dihedral angle of 25.CA 26.CA 27.CA 28.CA as second reaction coordinate.
cond_create.sh
bias resd
bias dihe
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=57,atom2=CA,force=0.05,target=-150,sel1=25.CA,sel2=26.CA,sel3=27.CA,sel4=28.CA
300 force=1.0,target=5,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=57,atom2=CA,force=0.05,target=-100,sel1=25.CA,sel2=26.CA,sel3=27.CA,sel4=28.CA
..........
300 force=1.0,target=245,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=57,atom2=CA,force=0.05,target=100,sel1=25.CA,sel2=26.CA,sel3=27.CA,sel4=28.CA
300 force=1.0,target=245,seg1=PRO0,res1=1,atom1=CA,seg2=PRO0,res2=57,atom2=CA,force=0.05,target=150,sel1=25.CA,sel2=26.CA,sel3=27.CA,sel4=28.CA
output data.r*d*:
cond=0
for ((i=5;i<=100;i+=5)); do for ((j=-120;j<=120;j+=20)); do rexinfo.pl -inx 1:100 -condsel $cond run,val1,val2 -dir gbrex_va_resd_dihe_new >data.r${i}d$j; cond=`expr $cond + 1|bc`; done; done
create wham2dinput
for ((i=5;i<=100;i+=5)); do for ((j=-120;j<=120;j+=20)); do echo "data.r${i}d$j $i $j 1.0 0.0015">>wham2dinput; done; done
NOTE:
in 2D conditional replica exchange simulation the format of energy.log turn into:
prod/0/1:300.000000:15:: 0 0.00 301.14 -978.93 670.54 -1649.47 -182.40 -1904.98 0.00 0.00 0.00 0.000000 0.000000
rmsd_grid.sh gbrex_va_resd_dihe
when plot 2d:
awk '{if ($3<9999) {print $0} }' whamoutput-2d
to delete probability=99999
CORR ANALYSIS:
wordom_0.22-rc2.x86-64 -iA ~/toolset/corr.inp -imol final.pdb -itrj traj-aa1.dcd
"corrtest" will be created under current directory
matrix_trans.pl corrtest
"corrmodi" will be created under current directory
gnuplot>set pm3d map
gnuplot>set xrange [0:56]
gnuplot>set yrange[0:56]
gnuplot>set cbrange[-1:1]
gnuplot>plot "corrmodi" matrix with image
SMD:
make sure 6pti.pdb and 6pti.psf(genPSF.pl -xplor ) exist.
vmd -dispdev text -e solvate_ionize.tcl -args 6pti
6pti_wb_ionized.pdb 6pi_wb_ionized.psf 6pti_wb_ionized.ref 6pti.namd will be created.
namd2 +p num 6pti.namd >equi.log
6pti_equi.coor will be created.
edit pull_v.namd
change coordinates,structure,cellBasis,cellOrigin,fixedAtomsFile
edit pull_v.tcl
set fix num
set smd num
set numatoms
array set fix_ini
array set smd_ini
then namd2 +p num pull_v.namd>pull_v.log
BPTI has total 4 turns:1CA--25CA,36CA--17CA,36CA--46CA, 43CA--57CA
pull_v.namd is the same, but tcl are different: pull_v_endtoend.tcl,pull_v_1_5.tcl, pull_v_17_36.tcl,...
Water Bridge:
vmd 6pti_wb_ionized.psf 5000.dcd
source ~/toolset/tmp.tcl
start_protein_bridging
then click play, until the end of trajectory, click stop
stop_protein_bridging
then look file amino_stat.txt under current directory.
gnuplot
gnuplot>plot "amino_stat.txt" matrix w image
if I want to look the value greater than some criteria:
awk '{for(i=1;i<=57;i++) if ($i > 2000) {print FNR ":" $0}}' amino_stat.txt
load trajecoties of each cluster(t.1,t.2,....)
vmd -e ~/toolset/loadcluster.tcl -args t.1
draw color yellow
draw text {25 -10 0} "t.2" size 3
HUB:
wordom -iA ~/toolset/psn.inp -imol 6pti.pdb -itrj 6pti_equi.dcd
avg_trans.pl avgpsn-example 57
avg_hub_gnu.sh avgHUB