ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/owl/trunk/modelit/scripts/refine.sh
Revision: 1299
Committed: Tue Jan 11 17:38:33 2011 UTC (8 years, 8 months ago) by hstehr
File size: 9235 byte(s)
Log Message:
adding Model-It server files to repository
Line File contents
1 #!/bin/sh
2 # Script for automated simulated annealing structure refinement using Gromacs
3
4 # CONSTANTS
5 BOXTYPE="dodecahedron"
6 BOXSIDES="0.5"
7
8 MDPDIR=/project/StruPPi/Software/gromacs-3.3.3/mdp_files
9 # energy minimization parameters (used for ions too)
10 MDP_EM=$MDPDIR/em.mdp
11 # position restrained equilibration parameters
12 MDP_PR=$MDPDIR/pr.mdp
13 # MD simulation parameters
14 MDP_MD=$MDPDIR/md.mdp
15
16 SOLVENTGROUP=12
17
18 DEFAULT_NP=1
19
20 # END CONSTANTS
21
22 # command line parameters
23 inPdb=""
24 outDir=""
25 np=$DEFAULT_NP
26 begT="-1"
27 endT="-1"
28
29 while getopts i:o:p:b:e: opt
30 do
31 case "$opt" in
32 i) inPdb="$OPTARG";;
33 o) outDir="$OPTARG";;
34 p) np="$OPTARG";;
35 b) begT="$OPTARG";;
36 e) endT="$OPTARG";;
37 esac
38 done
39
40 if [ -z "$inPdb" ] || [ -z "$outDir" ]
41 then
42 echo ""
43 echo "Usage: $0 "
44 echo " -i <in pdb>"
45 echo " -o <out dir> "
46 echo " [-p] <number of processes>"
47 echo " [-b] <begin time (ns)>"
48 echo " [-e] <end time (ns)>"
49 echo "Use begin and end time (ns) so that the md simulation will output to different files every "
50 echo "1ns, instead of using 1 big file for all output."
51 echo "Use absolute paths to run on cluster!!"
52 echo ""
53 exit 1
54 fi
55
56 # gromacs dir
57 arch=`uname -m`
58 gmxDir=/project/StruPPi/Software/gromacs-3.3.3/i686/bin
59 if [ "$arch" == "x86_64" ]
60 then
61 gmxDir=/project/StruPPi/Software/gromacs-3.3.3/x86_64/bin
62 fi
63
64 # mdrun command (use mpi or not)
65 mdrunCmd="$gmxDir/mdrun"
66 if [ "$np" -gt "1" ]
67 then
68 mdrunCmd="mpirun -np $np $gmxDir/mdrun_mpi"
69 fi
70
71 # mdp files
72
73 mdpEM=$MDP_EM
74 mdpPR=$MDP_PR
75 mdpMD=$MDP_MD
76 if [ -e "$outDir/em.mdp" ]
77 then
78 mdpEM="em.mdp"
79 echo "Using $mdpEM in $outDir as parameter file for EM step"
80 fi
81 if [ -e "$outDir/pr.mdp" ]
82 then
83 mdpPR="pr.mdp"
84 echo "Using $mdpPR in $outDir as parameter file for PR step"
85 fi
86 if [ -e "$outDir/md.mdp" ]
87 then
88 mdpMD="md.mdp"
89 echo "Using $mdpMD in $outDir as parameter file for MD step"
90 fi
91
92
93 outLog=out.log
94 cmdLog=cmd.log
95
96 basename=`basename $inPdb .pdb`
97
98 # gro files
99 groFile=$basename.gro
100 boxFile=$basename.box.gro
101 watFile=$basename.wat.gro
102 ionFile=$basename.ion.gro
103 emFile=$basename.em.gro
104 prFile=$basename.pr.gro
105
106 # final pdb file
107 mdPdbFile=$basename.md.pdb
108
109 # top files
110 topFile=$basename.top
111
112
113 # FUNCTIONS
114
115 # checkExitStatus, takes 2 parameters:
116 # exit value
117 # program name
118 checkExitStatus () {
119 exitVal="$1"
120 prog="$2"
121 if [ "$exitVal" -ne "0" ]
122 then
123 echo "$prog failed. Revise log file $outLog. Exiting"
124 exit 1
125 fi
126 }
127
128 # runSimulation, takes 5 parameters:
129 # mdp file
130 # in gro file
131 # in top file
132 # suffix of the input gro file (without the .gro), e.g. for T0123.pr.gro the suffix is "pr"
133 # suffix of the output files, e.g. if input is T0123.pr.gro and we specify "md" as suffix output, we'd get output files like T0123.md.gro or T0123.md.xtc
134 runSimulation () {
135 if [ -z "$5" ]
136 then
137 echo "Not enough parameters were passed to runSimulation. This is a bug in this script... exiting"
138 exit 1
139 fi
140
141 # globals: outLog, cmdLog, gmxDir, mdrunCmd, np
142
143 # parameters
144 mdp="$1"
145 gro="$2"
146 top="$3"
147 suffixIn="$4"
148 suffixOut="$5"
149
150 # locals
151 bn1=`basename $gro .gro`
152 bn2=`basename $bn1 .$suffixIn`
153 tpr=$bn2.$suffixOut.tpr
154 trr=$bn2.$suffixOut.trr
155 outGro=$bn2.$suffixOut.gro
156 edr=$bn2.$suffixOut.edr
157 xtc=$bn2.$suffixOut.xtc
158 log=$bn2.$suffixOut.log
159
160 # grompp creates a mdout.mdp file by default, which is (I think) useless
161 # If there is already one in the outdir it will back it up giving the backed up file a version number.
162 # But gromacs won't backup up beyond version 128: because of that we remove the file here before running grompp
163 rm -f mdout.mdp
164
165 echo -e "##########\n# preparing run with input $bn1.gro and output $outGro (grompp)\n##########" >> $outLog
166 echo "$gmxDir/grompp -f $mdp -c $gro -p $top -o $tpr -np $np" >> $cmdLog
167 $gmxDir/grompp -f $mdp -c $gro -p $top -o $tpr -np $np 1>> $outLog 2>&1
168
169 echo -e "##########\n# mdrun with input $bn1.gro and output $outGro (mdrun)\n##########" >> $outLog
170 echo "$mdrunCmd -s $tpr -o $trr -c $outGro -e $edr -x $xtc -g $log -np $np" >> $cmdLog
171 $mdrunCmd -s $tpr -o $trr -c $outGro -e $edr -x $xtc -g $log -np $np 1>> $outLog 2>&1
172
173 checkExitStatus $? mdrun
174 }
175
176 # addIons, takes 7 parameters:
177 # positive charge
178 # negative charge
179 # mdp file
180 # in gro file
181 # in top file
182 # suffix of the input gro file (without the .gro), e.g. for T0123.pr.gro the suffix is "pr"
183 # suffix of the output files, e.g. if input is T0123.pr.gro and we specify "md" as suffix output, we'd get output files like T0123.md.gro or T0123.md.xtc
184 addIons () {
185 if [ -z "$7" ]
186 then
187 echo "Not enough parameters were passed to addIons. This is a bug in this script... exiting"
188 exit 1
189 fi
190 # globals: SOLVENTGROUP, mdpEM, outLog, cmdLog, gmxDir
191
192 # parameters
193 poscharge="$1"
194 negcharge="$2"
195 mdp="$3"
196 gro="$4"
197 top="$5"
198 suffixIn="$6"
199 suffixOut="$7"
200
201 #locals
202 bn1=`basename $gro .gro`
203 bn2=`basename $bn1 .$suffixIn`
204 tpr=$bn2.$suffixOut.tpr
205 outGro=$bn2.$suffixOut.gro
206 log=$bn2.$suffixOut.log
207 ionOpt="-nname CL- -nn $negcharge -pname NA+ -np $poscharge"
208
209 echo "Adding $negcharge CL- and $poscharge NA+"
210 echo -e "##########\n# preparing to add ions (grompp)\n##########" >> $outLog
211 echo "$gmxDir/grompp -f $mdp -c $gro -p $top -o $tpr" >> $cmdLog
212 $gmxDir/grompp -f $mdp -c $gro -p $top -o $tpr 1>> $outLog 2>&1
213 echo -e "##########\n# adding ions (genion)\n##########" >> $outLog
214 echo "echo $SOLVENTGROUP | $gmxDir/genion -s $tpr -o $outGro -p $top $ionOpt -g $log" >> $cmdLog
215 echo $SOLVENTGROUP | $gmxDir/genion -s $tpr -o $outGro -p $top $ionOpt -g $log 1>> $outLog 2>&1
216
217 checkExitStatus $? genion
218 }
219
220
221 # END FUNCTIONS
222
223 # gromacs assumes that its output directory is the same as the current, a few things don't work well if that's not the case
224 # thus we change to the outDir before starting to run any commands,
225 # that means we use directly the file names without paths for all input/output files, except for $inPdb that can be in another directory
226 firstChar=`echo $inPdb | cut -c 1` # first character of inPdb
227 if [ "$firstChar" != "/" ] # i.e. inPdb is not an absoute path
228 then
229 inPdb=`pwd`/$inPdb
230 fi
231
232
233 cd $outDir
234 echo "" > $outLog
235 echo "" > $cmdLog
236
237
238 # 1. convert pdb to gro and topology files
239 echo "Converting pdb to gro"
240 echo -e "##########\n# pdb2gmx\n##########" >> $outLog
241 echo "$gmxDir/pdb2gmx -f $inPdb -o $groFile -p $topFile -ff G43a1 -ignh" >> $cmdLog
242 $gmxDir/pdb2gmx -f $inPdb -o $groFile -p $topFile -ff G43a1 -ignh 1>> $outLog 2>&1
243 checkExitStatus $? pdb2gmx
244 charge=`cat $outLog | grep "^Total charge " | perl -pe 's/^Total charge (-?\d+)\.\d+ e/$1/'`
245
246 # 2. center the molecule within a box (with specified type and space on the sides)
247 echo "Creating box"
248 echo -e "##########\n# creating box (editconf)\n##########" >> $outLog
249 echo "$gmxDir/editconf -bt $BOXTYPE -f $groFile -o $boxFile -c -d $BOXSIDES" >> $cmdLog
250 $gmxDir/editconf -bt $BOXTYPE -f $groFile -o $boxFile -c -d $BOXSIDES 1>> $outLog 2>&1
251 checkExitStatus $? editconf
252
253 # 3. add solvent water to the box (-cp is solute, -cs solvent)
254 echo "Adding solvent to box"
255 echo -e "##########\n# adding solvent (genbox)\n##########" >> $outLog
256 echo "$gmxDir/genbox -cp $boxFile -cs spc216.gro -o $watFile -p $topFile" >> $cmdLog
257 $gmxDir/genbox -cp $boxFile -cs spc216.gro -o $watFile -p $topFile 1>> $outLog 2>&1
258 checkExitStatus $? genbox
259
260 # 4. adding ions if needed
261 poscharge=0
262 negcharge=0
263 if [ "$charge" -lt "0" ]
264 then
265 # no easy way to get abs value, we use bc instead of expr
266 poscharge=`echo "$charge/1*-1 + 1"|bc`
267 negcharge=1
268 else
269 if [ "$charge" -gt "0" ]
270 then
271 poscharge=1
272 negcharge=`expr $charge + 1`
273 else
274 if [ "$charge" -eq "0" ]
275 then
276 poscharge=1
277 negcharge=1
278 fi
279 fi
280 fi
281 if [ "$poscharge" -eq "0" ]
282 then
283 # if poscharge is still 0, then none of the if conditions above executed --> $charge was not set
284 echo "Charge value couldn't be parsed. Exiting"
285 exit 1
286 fi
287 echo "Total charge $charge"
288 addIons $poscharge $negcharge $mdpEM $watFile $topFile wat ion
289
290 # 5. minimization
291 echo "Energy minimization"
292 runSimulation $mdpEM $ionFile $topFile ion em
293
294 # 6. equilibration
295 echo "Position restrained equilibration"
296 runSimulation $mdpPR $emFile $topFile em pr
297
298 # 7. MD simulation
299 echo "Molecular dynamics simulation"
300 lastSuffix=""
301 if [ "$begT" -eq "-1" ] || [ "$endT" -eq "-1" ]
302 then
303 runSimulation $mdpMD $prFile $topFile pr md
304 lastSuffix="md"
305 else
306 echo "Splitting output per ns"
307 tStr=`printf %04.0f $begT`ns
308 echo "Running simulation up to $tStr"
309 runSimulation $mdpMD $prFile $topFile pr $tStr.md
310 endTminus1=`expr $endT - 1`
311 for bt in `seq $begT $endTminus1`
312 do
313 btStr=`printf %04.0f $bt`ns
314 et=`expr $bt + 1`
315 etStr=`printf %04.0f $et`ns
316 echo "Running simulation up to $etStr"
317 runSimulation $mdpMD $basename.$btStr.md.gro $topFile $btStr.md $etStr.md
318 lastSuffix="$etStr.md"
319 done
320
321 fi
322
323
324 # 8. converting the final md gro file to pdb
325 mdFile=$basename.$lastSuffix.gro
326 tprMD=$basename.$lastSuffix.tpr
327
328 echo "Converting $mdFile to $mdPdbFile"
329 echo -e "##########\n# converting gro to pdb (trjconv)\n##########" >> $outLog
330 echo "echo 1 | $gmxDir/trjconv -f $mdFile -o $mdPdbFile -s $tprMD" >> $cmdLog
331 echo 1 | $gmxDir/trjconv -f $mdFile -o $mdPdbFile -s $tprMD 1>> $outLog 2>&1
332 checkExitStatus $? trjconv

Properties

Name Value
svn:executable *