1 |
#!/usr/bin/env perl |
2 |
|
3 |
# program that reads in a water box and solute (md), merges the two systems, |
4 |
# and deletes overlapping molecules |
5 |
|
6 |
# author = "Dan Gezelter" |
7 |
# version = "$Revision: 1.1 $" |
8 |
# date = "$Date: 2007-01-11 21:46:14 $" |
9 |
# copyright = "Copyright (c) 2007 by the University of Notre Dame" |
10 |
# license = "OOPSE" |
11 |
|
12 |
|
13 |
use Getopt::Std; |
14 |
|
15 |
$d2tolerance = 7.5625; #distance to start cutting |
16 |
$fileName = 'solvatedSystem.md'; |
17 |
$startSnap = 0; |
18 |
$startFrame = 0; |
19 |
$startStunts = 0; |
20 |
$startMeta = 0; |
21 |
$nSolvent = 0; |
22 |
|
23 |
# get our options |
24 |
getopts('hd:v:u:o:'); |
25 |
|
26 |
# if we don't have a filename, drop to -h |
27 |
$opt_h = 'true' if $#ARGV != -1; |
28 |
|
29 |
# our option output |
30 |
if ($opt_h){ |
31 |
print "solvator: carves a solute void in an oopse water box\n\n"; |
32 |
print "usage: solvator [-h] [-d distance tolerance] [-v file name (solvent)]"; |
33 |
print " [-u file name (solute)] [-o file name (output)]\n\n"; |
34 |
print " -h : show this message\n\n"; |
35 |
print " -d real : overlap removal distance\n"; |
36 |
print " (default: 2.75 ang)\n"; |
37 |
print " -u char : solute input file (oopse .md file)\n"; |
38 |
print " -v char : solvent input file (oopse .md file)\n"; |
39 |
print " -o char : carved solvent output file (oopse .md format)\n"; |
40 |
print " (default: 'solvatedSystem.md')\n"; |
41 |
print "Example:\n"; |
42 |
die " md-solvator -v freshWater.md -u mySolute.md -o mySystem.md\n"; |
43 |
} |
44 |
|
45 |
# set some variables to be used in the code |
46 |
if (defined($opt_v)){ |
47 |
$solventName = $opt_v; |
48 |
} else { |
49 |
die "Error: No solvent box specified\n Please select a solvent box using the -v flag\n"; |
50 |
} |
51 |
|
52 |
$soluteFileName = $opt_u if defined($opt_u); |
53 |
$fileName = $opt_o if defined($opt_o); |
54 |
|
55 |
if (!defined($opt_u)){ |
56 |
die "Error: No solute file specifed\n Please select a solute file with the -u flag\n"; |
57 |
} |
58 |
|
59 |
if (defined($opt_d)){ |
60 |
if ($opt_d =~ /^[0-9]/) { |
61 |
$dval = $opt_d; |
62 |
$d2tolerance = $dval*$dval; |
63 |
} else { |
64 |
die "Error: The '-d' value ($opt_d) is not a valid number\n Please choose a positive real # value\n"; |
65 |
} |
66 |
} |
67 |
|
68 |
# open the file writer |
69 |
open(OUTFILE, ">./$fileName") || die "Error: can't open file $fileName\n"; |
70 |
|
71 |
# open and read the pdb or xyz file and get solute coordinates |
72 |
open(SOLUTEFILE, "$soluteFileName") || die "Error: Can't open file $soluteFileName\n"; |
73 |
|
74 |
|
75 |
while (<SOLUTEFILE>){ |
76 |
$startSnap = 0 if /\/Snapshot/; |
77 |
$startFrame = 0 if /\/FrameData/; |
78 |
$startStunts = 0 if /\/StuntDoubles/; |
79 |
$startMeta = 0 if /\/MetaData/; |
80 |
|
81 |
# save the MetaData lines |
82 |
push(@soluteMetaLines, $_) if $startMeta == 1; |
83 |
|
84 |
if ($startSnap == 1){ |
85 |
# save the frame data |
86 |
push(@soluteFrameData, $_) if $startFrame == 1; |
87 |
if (/Hmat/){ |
88 |
@line = split; |
89 |
$soluteHxx = $line[2]; |
90 |
chop($soluteHxx); |
91 |
$soluteHyy = $line[8]; |
92 |
chop($soluteHyy); |
93 |
$soluteHzz = $line[14]; |
94 |
} |
95 |
if ($startStunts == 1){ |
96 |
@line = split; |
97 |
$x_val = $line[2] - ($soluteHxx*round($line[2]/$soluteHxx)); |
98 |
$y_val = $line[3] - ($soluteHyy*round($line[3]/$soluteHyy)); |
99 |
$z_val = $line[4] - ($soluteHzz*round($line[4]/$soluteHzz)); |
100 |
|
101 |
push(@solute_x, $x_val); |
102 |
push(@solute_y, $y_val); |
103 |
push(@solute_z, $z_val); |
104 |
$saveLine = ""; |
105 |
for ($i = 0; $i <= $#line; $i++){ |
106 |
$saveLine = join(' ', $saveLine, $line[$i]); |
107 |
} |
108 |
push(@soluteStuntDoubles, $saveLine); |
109 |
$soluteCount++; |
110 |
} |
111 |
} |
112 |
$startSnap = 1 if /Snapshot/; |
113 |
$startFrame = 1 if /FrameData/; |
114 |
$startStunts = 1 if /StuntDoubles/; |
115 |
$startMeta = 1 if /MetaData/; |
116 |
# check again since "/value" contains "value" |
117 |
$startSnap = 0 if /\/Snapshot/; |
118 |
$startFrame = 0 if /\/FrameData/; |
119 |
$startStunts = 0 if /\/StuntDoubles/; |
120 |
$startMeta = 0 if /\/MetaData/; |
121 |
} |
122 |
|
123 |
$solventCount = $soluteCount; |
124 |
# okay, let's open the solvent box and write a carved file |
125 |
open(SOLVENT, "$solventName") || die "Error: Can't open file $solventName\n"; |
126 |
|
127 |
|
128 |
while (<SOLVENT>){ |
129 |
$startSnap = 0 if /\/Snapshot/; |
130 |
$startFrame = 0 if /\/FrameData/; |
131 |
$startStunts = 0 if /\/StuntDoubles/; |
132 |
$startMeta = 0 if /\/MetaData/; |
133 |
|
134 |
# save the MetaData lines |
135 |
push(@metaLines, $_) if $startMeta == 1; |
136 |
|
137 |
if ($startSnap == 1){ |
138 |
# save the frame data |
139 |
push(@frameData, $_) if $startFrame == 1; |
140 |
if (/Hmat/){ |
141 |
@line = split; |
142 |
$hxx = $line[2]; |
143 |
chop($hxx); |
144 |
$hyy = $line[8]; |
145 |
chop($hyy); |
146 |
$hzz = $line[14]; |
147 |
} |
148 |
if ($startStunts == 1){ |
149 |
@line = split; |
150 |
|
151 |
# wrap positions back into the box |
152 |
$x_val = $line[2] - ($hxx*round($line[2]/$hxx)); |
153 |
$y_val = $line[3] - ($hyy*round($line[3]/$hyy)); |
154 |
$z_val = $line[4] - ($hzz*round($line[4]/$hzz)); |
155 |
|
156 |
$saveFlag = 1; |
157 |
for($i=0; $i<=$#solute_x; $i++){ |
158 |
$diff_x = $x_val - $solute_x[$i]; |
159 |
$diff_y = $y_val - $solute_y[$i]; |
160 |
$diff_z = $z_val - $solute_z[$i]; |
161 |
$dist2 = $diff_x*$diff_x + $diff_y*$diff_y + $diff_z*$diff_z; |
162 |
if ($dist2 < $d2tolerance) { |
163 |
$saveFlag = 0; |
164 |
last; |
165 |
} |
166 |
} |
167 |
if ($saveFlag == 1){ |
168 |
$saveLine = "$solventCount\t$line[1]\t$line[2]"; |
169 |
for ($i = 3; $i <= $#line; $i++){ |
170 |
$saveLine = join(' ', $saveLine, $line[$i]); |
171 |
} |
172 |
push(@goodSolventMolecules, $saveLine); |
173 |
$solventCount++; |
174 |
} |
175 |
} |
176 |
} |
177 |
$startSnap = 1 if /Snapshot/; |
178 |
$startFrame = 1 if /FrameData/; |
179 |
$startStunts = 1 if /StuntDoubles/; |
180 |
$startMeta = 1 if /MetaData/; |
181 |
# check again since "/value" contains "value" |
182 |
$startSnap = 0 if /\/Snapshot/; |
183 |
$startFrame = 0 if /\/FrameData/; |
184 |
$startStunts = 0 if /\/StuntDoubles/; |
185 |
$startMeta = 0 if /\/MetaData/; |
186 |
} |
187 |
|
188 |
$nSolvent = $#goodSolventMolecules + 1; |
189 |
|
190 |
# write out the final file |
191 |
writeOutFile(); |
192 |
print "The solvated system \"$fileName\" was generated.\n"; |
193 |
|
194 |
|
195 |
sub writeOutFile { |
196 |
# write out the header |
197 |
print OUTFILE "<OOPSE version=4>\n"; |
198 |
printMetaData(); |
199 |
printFrameData(); |
200 |
print OUTFILE " <StuntDoubles>\n"; |
201 |
# now print out the atoms |
202 |
for ($i=0; $i<=$#solute_x; $i++){ |
203 |
print OUTFILE "$soluteStuntDoubles[$i]\n"; |
204 |
} |
205 |
for ($i=0; $i<=$#goodSolventMolecules; $i++) { |
206 |
print OUTFILE "$goodSolventMolecules[$i]\n"; |
207 |
} |
208 |
print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OOPSE>\n"; |
209 |
} |
210 |
|
211 |
sub printMetaData() { |
212 |
print OUTFILE " <MetaData>\n"; |
213 |
|
214 |
writeSoluteDescription(); |
215 |
|
216 |
for ($i = 0; $i<=$#metaLines; $i++) { |
217 |
# reset the number of solvent molecules |
218 |
if ($metaLines[$i] =~ /nMol/) { |
219 |
$metaLines[$i] = " nMol = $nSolvent;\n"; |
220 |
} |
221 |
|
222 |
print OUTFILE "$metaLines[$i]"; |
223 |
} |
224 |
|
225 |
print OUTFILE " </MetaData>\n"; |
226 |
} |
227 |
|
228 |
sub writeSoluteDescription { |
229 |
# include a solute model description in the meta data region |
230 |
|
231 |
for ($i = 0; $i<=$#soluteMetaLines; $i++) { |
232 |
print OUTFILE "$soluteMetaLines[$i]"; |
233 |
} |
234 |
print "The solute model definition was included in '$fileName'.\n"; |
235 |
} |
236 |
|
237 |
sub printFrameData { |
238 |
print OUTFILE " <Snapshot>\n <FrameData>\n"; |
239 |
|
240 |
for($i = 0; $i<=$#frameData; $i++) { |
241 |
print OUTFILE "$frameData[$i]"; |
242 |
} |
243 |
|
244 |
print OUTFILE " </FrameData>\n"; |
245 |
} |
246 |
|
247 |
sub round { |
248 |
return int( $_[0] + 0.5 * ($_[0] <=> 0) ); |
249 |
} |