Per-residue energy decompositions in NAMD
** 原文地址 之前一篇文章中，通过计算范德华和静电相互作用，发现amber和namd的计算结果一致。在本文中，将通过namd来对轨迹进行分析，将范德华和静电相互作用分解到受体的每一个残基上，方法相当于将受体分解为单个的残基pdb，分别计算之后再加和。当然，这完全可以由amber的相应模块实现。此处贴出本文章，主要是认为这种实现方式很有趣，如果让我来做的话，我该是要扑街了。 Calculate Decompositions Using Namd **
The following script uses vmd to create pdb files. Each atom is labeled with a 0 (excluded), 1(species #1) or 2(species #2) in beta column. namd will compute the interaction energy between species #1 and #2. A separate pdb file is created for each residue (species #2). The ligand is always species #1. We are using amber inputs in this example. The namd input file loops over each pdb file to compute the interaction energy for every timestep. run.005.vmd.namd.energy.decomp_sw.csh: 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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135#! /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/user05/lingling/zzz.programs/vmd-1.9.1/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 # Default
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"
/nfs/user03/fochtman/zzz.programs/NAMD_2.9_Source/Linux-x86_64-g++/namd2 07_to_08pr.${count01}.in > 07_to_08pr.${count01}.out
EOF
qsub fp.namd.${count01}.csh
@ count01 = $count01 + 1
###########################################
###########################################
end
exit
** Energy Decomposition Analysis **
This script using awk to calculate the mean and standard deviation. run.006.vmd.namd.energy.decomp_sw_analysis.csh: 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
31
32
33
34
35
36
37#! /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