Difference between revisions of "Per-residue energy decompositions"

From Rizzo_Lab
Jump to: navigation, search
(→‎Per-residue energy decompositions)
(→‎Using Namd)
Line 1: Line 1:
  
−
== Using Namd ==
+
The Molecular Mechanics energy function is a pair-wise additive function which can be broken into per-residue components.
 +
 
 +
We call these plots molecular interaction footprints.
 +
 
 +
== Calculate Decompositions Using Namd ==
 +
 
 +
 
 +
The following script uses vmd to create pdb files label every atom with the
 +
 
 +
 
 +
run.005.vmd.namd.energy.decomp_sw.csh:
 +
 
 +
#! /bin/csh -fe
 +
 +
set mountdir        = "/nfs/user03/chembio/AMBER_Project"
 +
 +
set start          = 1
 +
set stop            = 35
 +
set offset          = 1
 +
 +
set workdir        = "${mountdir}/006.FP.VMD.NAMD"
 +
 +
cd        ${workdir}
 +
 +
###########################################
 +
 +
set res_start = ${start}
 +
set res_stop  = ${stop}
 +
@ loop_stop = ${res_stop} + 1
 +
 +
cat << EOF > ptraj.in
 +
trajin ${mountdir}/004.ptraj/1qp6.rec.trj.stripfit
 +
trajout 1qp6.rec.stripfit.dcd charmm
 +
EOF
 +
 +
ptraj ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.prm7 ptraj.in > ptraj.out
 +
 +
 +
###########################################
 +
cat <<EOFK >pair.maker.tcl
 +
#opens the file that has the atoms that are to be set-up for pair-interaction calculations
 +
mol load parm7 ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.prm7 rst7 ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.rst7
 +
 +
for {set resid_num $res_start} {\$resid_num < $loop_stop} {incr resid_num} {
 +
set all [atomselect top all]
 +
\$all set beta 0
 +
set case1 [atomselect top "protein and resid 36 to 70"]
 +
\$case1 set beta 1
 +
set case2 [atomselect top "protein and resid \$resid_num"]
 +
\$case2 set beta 2
 +
\$all writepdb 1qp6.rec.gas.pair.\$resid_num.pdb
 +
}
 +
quit
 +
EOFK
 +
###########################################
 +
# copy DCD files required in the calculations
 +
# to the current directory
 +
 +
# make a special PDB files defining what the two species will be for pair interaciton calcs
 +
  /nfs/user03/sudipto/vmd-1.9/bin/vmd -dispdev text -e pair.maker.tcl        >> zzz.pair.maker.out
 +
 +
echo "MAKE INPUT FILES"
 +
date
 +
 +
set res_start = ${start}
 +
set res_stop  = ${stop}
 +
@ loop_stop = ${res_stop} + 1
 +
 +
set start_residue = $res_start
 +
set loop_stop  = $loop_stop
 +
set count01 = $start_residue
 +
 +
###########################################
 +
###########################################
 +
while ( $count01 < $loop_stop )
 +
###########################################
 +
cat <<EOFL >07_to_08pr.$count01.in
 +
# initial config
 +
temperature    0
 +
 +
# output params
 +
outputname      1qp6.rec.gas.$count01
 +
binaryoutput    no
 +
 +
# integrator params
 +
timestep        1.0
 +
 +
# force field params
 +
amber          on
 +
parmfile        ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.prm7
 +
ambercoor      ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.rst7
 +
readexclusions  yes  # exclusions are taken from parmfile
 +
exclude        scaled1-4
 +
1-4scaling      0.833333  # =1/1.2, default is 1.0
 +
scnb            2.0  # Defalt
 +
switching      off
 +
cutoff          999.0
 +
 +
# Atoms in group 1 have a 1 in the B column; group 2 has a 2
 +
pairInteraction        on
 +
pairInteractionFile    1qp6.rec.gas.pair.$count01.pdb
 +
pairInteractionCol      B
 +
pairInteractionGroup1  1
 +
pairInteractionGroup2  2
 +
 +
# First frame saved was frame 1000.
 +
set ts 1000
 +
 +
coorfile open dcd 1qp6.rec.stripfit.dcd
 +
 +
# Read all frames until nonzero is returned.
 +
while { ![coorfile read] } {
 +
  # Set firstTimestep so our energy output has the correct TS.
 +
  firstTimestep \$ts
 +
  # Compute energies and forces, but don't try to move the atoms.
 +
  run 0
 +
  incr ts 1000
 +
}
 +
coorfile close
 +
EOFL
 +
###########################################
 +
 +
 +
cat << EOF > fp.namd.${count01}.csh
 +
#! /bin/csh
 +
#PBS -l walltime=24:00:00
 +
#PBS -V
 +
 +
cd ${workdir}
 +
 +
hostname
 +
date
 +
echo "START CALCS"
 +
 +
namd2 07_to_08pr.${count01}.in > 07_to_08pr.${count01}.out
 +
 +
EOF
 +
 +
qsub fp.namd.${count01}.csh
 +
 +
 +
@ count01 = $count01 + 1
 +
###########################################
 +
###########################################
 +
end
 +
 +
exit
 +
 
 +
 
 +
Analysis script using ack which calculates the mean and standard deviation
 +
 
 +
run.006.vmd.namd.energy.decomp_sw_analysis.csh:
 +
#! /bin/csh -fe
 +
 +
set mountdir        = "/nfs/user03/chembio/AMBER_Project"
 +
 +
set start          = 1
 +
set stop            = 35
 +
set offset          = 1
 +
 +
set workdir        = "${mountdir}/006.FP.VMD.NAMD"
 +
 +
cd        ${workdir}
 +
 +
set res_start = ${start}
 +
set res_stop  = ${stop}
 +
@ loop_stop = ${res_stop} + 1
 +
 +
set start_residue = $res_start
 +
set loop_stop  = $loop_stop
 +
set count01 = $start_residue
 +
 +
echo "resid,mean,std" > vdw.fp.output.csv
 +
echo "resid,mean,std" > es.fp.output.csv
 +
 +
###########################################
 +
while ( $count01 < $loop_stop )
 +
 +
## Gets VDW
 +
awk ' BEGIN{count=0;sum=0;sum2=0;} /^ENERGY:/{count=count+1;sum=sum+$8;sum2 = sum2 + $8^2} END{mean=sum/count; std= ((sum2/count)- (mean^2))^(1/2); print "'${count01}'," mean "," std } ' 07_to_08pr.${count01}.out >> vdw.fp.output.csv
 +
 +
##Get ES
 +
awk ' BEGIN{count=0;sum=0;sum2=0;} /^ENERGY:/{count=count+1;sum=sum+$7;sum2 = sum2 + $7^2} END{mean=sum/count; std= ((sum2/count)- (mean^2))^(1/2); print "'${count01}'," mean "," std } ' 07_to_08pr.${count01}.out >> es.fp.output.csv
 +
 +
@ count01 = $count01 + 1
 +
 +
end # while
 +
 +
exit

Revision as of 17:03, 5 May 2011

The Molecular Mechanics energy function is a pair-wise additive function which can be broken into per-residue components.

We call these plots molecular interaction footprints.

Calculate Decompositions Using Namd

The following script uses vmd to create pdb files label every atom with the


run.005.vmd.namd.energy.decomp_sw.csh: 

#! /bin/csh -fe

set mountdir        = "/nfs/user03/chembio/AMBER_Project"

set start           = 1
set stop            = 35
set offset          = 1

set workdir         = "${mountdir}/006.FP.VMD.NAMD"

cd        ${workdir}

###########################################

set res_start = ${start}
set res_stop  = ${stop}
@ loop_stop = ${res_stop} + 1

cat << EOF > ptraj.in
trajin ${mountdir}/004.ptraj/1qp6.rec.trj.stripfit
trajout 1qp6.rec.stripfit.dcd charmm
EOF

ptraj ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.prm7 ptraj.in > ptraj.out


###########################################
cat <<EOFK >pair.maker.tcl
#opens the file that has the atoms that are to be set-up for pair-interaction calculations
mol load parm7 ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.prm7 rst7 ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.rst7

for {set resid_num $res_start} {\$resid_num < $loop_stop} {incr resid_num} {
set all [atomselect top all]
\$all set beta 0
set case1 [atomselect top "protein and resid 36 to 70"]
\$case1 set beta 1
set case2 [atomselect top "protein and resid \$resid_num"]
\$case2 set beta 2
\$all writepdb 1qp6.rec.gas.pair.\$resid_num.pdb
}
quit
EOFK
###########################################
# copy DCD files required in the calculations
# to the current directory

# make a special PDB files defining what the two species will be for pair interaciton calcs
  /nfs/user03/sudipto/vmd-1.9/bin/vmd -dispdev text -e pair.maker.tcl         >> zzz.pair.maker.out

echo "MAKE INPUT FILES"
date

set res_start = ${start}
set res_stop  = ${stop}
@ loop_stop = ${res_stop} + 1

set start_residue = $res_start
set loop_stop  = $loop_stop
set count01 = $start_residue

###########################################
###########################################
while ( $count01 < $loop_stop )
###########################################
cat <<EOFL >07_to_08pr.$count01.in
# initial config
temperature     0

# output params
outputname      1qp6.rec.gas.$count01
binaryoutput    no

# integrator params
timestep        1.0

# force field params
amber           on
parmfile        ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.prm7
ambercoor       ${mountdir}/002.TLEAP/1qp6.rec.gas.leap.rst7
readexclusions  yes  # exclusions are taken from parmfile
exclude         scaled1-4
1-4scaling      0.833333  # =1/1.2, default is 1.0
scnb            2.0  # Defalt
switching       off
cutoff          999.0

# Atoms in group 1 have a 1 in the B column; group 2 has a 2
pairInteraction         on
pairInteractionFile     1qp6.rec.gas.pair.$count01.pdb
pairInteractionCol      B
pairInteractionGroup1   1
pairInteractionGroup2   2

# First frame saved was frame 1000.
set ts 1000

coorfile open dcd 1qp6.rec.stripfit.dcd

# Read all frames until nonzero is returned.
while { ![coorfile read] } {
  # Set firstTimestep so our energy output has the correct TS.
  firstTimestep \$ts
  # Compute energies and forces, but don't try to move the atoms.
  run 0
  incr ts 1000
}
coorfile close
EOFL
###########################################


cat << EOF > fp.namd.${count01}.csh
#! /bin/csh
#PBS -l walltime=24:00:00
#PBS -V

cd ${workdir}

hostname
date
echo "START CALCS"

namd2 07_to_08pr.${count01}.in > 07_to_08pr.${count01}.out

EOF

qsub fp.namd.${count01}.csh


@ count01 = $count01 + 1
###########################################
###########################################
end

exit


Analysis script using ack which calculates the mean and standard deviation

run.006.vmd.namd.energy.decomp_sw_analysis.csh: 

#! /bin/csh -fe

set mountdir        = "/nfs/user03/chembio/AMBER_Project"

set start           = 1
set stop            = 35
set offset          = 1

set workdir         = "${mountdir}/006.FP.VMD.NAMD"

cd        ${workdir}

set res_start = ${start}
set res_stop  = ${stop}
@ loop_stop = ${res_stop} + 1

set start_residue = $res_start
set loop_stop  = $loop_stop
set count01 = $start_residue

echo "resid,mean,std" > vdw.fp.output.csv
echo "resid,mean,std" > es.fp.output.csv

###########################################
while ( $count01 < $loop_stop )

## Gets VDW
awk ' BEGIN{count=0;sum=0;sum2=0;} /^ENERGY:/{count=count+1;sum=sum+$8;sum2 = sum2 + $8^2} END{mean=sum/count; std= ((sum2/count)- (mean^2))^(1/2); print "'${count01}'," mean "," std } ' 07_to_08pr.${count01}.out >> vdw.fp.output.csv

##Get ES
awk ' BEGIN{count=0;sum=0;sum2=0;} /^ENERGY:/{count=count+1;sum=sum+$7;sum2 = sum2 + $7^2} END{mean=sum/count; std= ((sum2/count)- (mean^2))^(1/2); print "'${count01}'," mean "," std } ' 07_to_08pr.${count01}.out >> es.fp.output.csv

@ count01 = $count01 + 1

end # while

exit
Retrieved from "https://ringo.ams.stonybrook.edu/index.php?title=Per-residue_energy_decompositions&oldid=6425"