MD simulation of BPTI

/6PTI.pdb (if want to include cystal water, use: "-cyrstalwater", if want ssbond, use: "-ssbond") 6pti.pdb>6pti_min.pdb
on sierra:
edit hosts.txt 6 /home/zwei/work/MD/BPTI/6PTI/REX/VA
dont use

and start.txt

native structure as starting structures -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 -inx 1:1000 -ref 6pti_min.pdb -dir varex

load into vmd
vmd -e ~/toolset/loadrex.tcl -args varex

create RMSD plot varex


6pti.pdb was copied from /BPTI/6PTI/6pti.pdb -out one 6pti.pdb> 32 xxx( 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: -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

6pti_va_rmsd_grid 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)
eigplot_rmsd.txt and eigresult_rmsd.txt will be created.

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
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: -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 -inx 1:100 -ref 6pti.pdb -dir gbrex_va_resd
vmd -e ~/toolset/loadrex.tcl -ags gbrex_va_resd drawcell

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

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*:


for ((i=5;i<=100;i+=5)); do for ((j=-120;j<=120;j+=20)); do -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


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 gbrex_va_resd_dihe


when plot 2d:

awk '{if ($3<9999) {print $0} }' whamoutput-2d

to delete probability=99999


wordom_0.22-rc2.x86-64 -iA ~/toolset/corr.inp -imol final.pdb -itrj traj-aa1.dcd
"corrtest" will be created under current directory 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

make sure 6pti.pdb and 6pti.psf( -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
then click play, until the end of trajectory, click stop
then look file amino_stat.txt under current directory.
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


wordom -iA ~/toolset/psn.inp -imol 6pti.pdb -itrj 6pti_equi.dcd avgpsn-example 57 avgHUB