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, 9 months ago) by hstehr
File size: 9235 byte(s)
Log Message:
adding Model-It server files to repository
Line User Rev File contents
1 hstehr 1299 #!/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 *