3 |
|
# program that builds water boxes |
4 |
|
|
5 |
|
# author = "Chris Fennell |
6 |
< |
# version = "$Revision: 1.1 $" |
7 |
< |
# date = "$Date: 2006-10-10 14:08:10 $" |
6 |
> |
# version = "$Revision: 1.2 $" |
7 |
> |
# date = "$Date: 2007-01-09 03:42:14 $" |
8 |
|
# copyright = "Copyright (c) 2006 by the University of Notre Dame" |
9 |
|
# license = "OOPSE" |
10 |
|
|
27 |
|
$nothingSelected = 1; |
28 |
|
|
29 |
|
# get our options |
30 |
< |
getopts('hmrvd:l:n:o:w:'); |
30 |
> |
getopts('hmrvd:l:n:o:w:x:y:z:'); |
31 |
|
|
32 |
|
# just to test opt_h |
33 |
|
$opt_h = "true" if $#ARGV >= 0; |
50 |
|
print " -o char : output file name\n"; |
51 |
|
print " (default: freshWater.md)\n"; |
52 |
|
print " -w char : name of the water stunt double\n"; |
53 |
< |
print " (default: SPCE)\n\n"; |
53 |
> |
print " (default: SPCE)\n"; |
54 |
> |
print " -x real : dimension of the box along the x direction\n"; |
55 |
> |
print " -y real : dimension of the box along the y direction\n"; |
56 |
> |
print " -z real : dimension of the box along the z direction\n\n"; |
57 |
> |
print "Note: you can only use values of x, y, or z that are smaller\n"; |
58 |
> |
print " than the derived box length for a given density and\n"; |
59 |
> |
print " number of molecules.\n\n"; |
60 |
|
print "Example:\n"; |
61 |
|
die " waterBoxer -d 0.997 -n 864 -w SSD_RF -o ssdrfWater.md\n"; |
62 |
|
} |
80 |
|
$nothingSelected = 0; |
81 |
|
} |
82 |
|
|
83 |
+ |
if (defined($opt_d)){ |
84 |
+ |
$nothingSelected = 0; |
85 |
+ |
if ($opt_d =~ /^[0-9]/) { |
86 |
+ |
$density = $opt_d; |
87 |
+ |
} else { |
88 |
+ |
die "Error: the value for '-d' ($opt_d) is not a valid number\n Please choose a positive real # value\n"; |
89 |
+ |
} |
90 |
+ |
} |
91 |
|
if (defined($opt_w)){ |
92 |
|
$waterName = $opt_w; |
93 |
|
$nothingSelected = 0; |
100 |
|
print " Use the \'-m\' option to generate a \'water.md\' with the\n"; |
101 |
|
print " recognized water model geometries.\n\n"; |
102 |
|
} |
103 |
< |
if (defined($opt_d)){ |
104 |
< |
$nothingSelected = 0; |
105 |
< |
if ($opt_d =~ /^[0-9]/) { |
92 |
< |
$density = $opt_d; |
93 |
< |
} else { |
94 |
< |
die "Error: the value for '-d' ($opt_d) is not a valid number\n Please choose a positive real # value\n"; |
95 |
< |
} |
103 |
> |
if ($waterName eq 'DPD') { |
104 |
> |
# DPD waters are stand-ins for 4 water molecules |
105 |
> |
$density = $density * 0.25; |
106 |
|
} |
107 |
+ |
|
108 |
|
if (defined($opt_l)){ |
109 |
|
$nothingSelected = 0; |
110 |
|
if ($opt_l =~ /^[0-9]/) { |
124 |
|
die "Error: the '-n' value ($opt_n) is not a valid number\n Please choose a non-negative integer\n"; |
125 |
|
} |
126 |
|
} |
127 |
+ |
if (defined($opt_x)){ |
128 |
+ |
$nothingSelected = 0; |
129 |
+ |
if ($opt_x =~ /^[0-9]/) { |
130 |
+ |
$boxx = $opt_x; |
131 |
+ |
} else { |
132 |
+ |
die "Error: the value for '-x' ($opt_x) is not a valid number\n Please choose a positive real # value\n"; |
133 |
+ |
} |
134 |
+ |
} |
135 |
+ |
if (defined($opt_y)){ |
136 |
+ |
$nothingSelected = 0; |
137 |
+ |
if ($opt_y =~ /^[0-9]/) { |
138 |
+ |
$boxy = $opt_y; |
139 |
+ |
} else { |
140 |
+ |
die "Error: the value for '-y' ($opt_y) is not a valid number\n Please choose a positive real # value\n"; |
141 |
+ |
} |
142 |
+ |
} |
143 |
+ |
if (defined($opt_z)){ |
144 |
+ |
$nothingSelected = 0; |
145 |
+ |
if ($opt_z =~ /^[0-9]/) { |
146 |
+ |
$boxz = $opt_z; |
147 |
+ |
} else { |
148 |
+ |
die "Error: the value for '-z' ($opt_z) is not a valid number\n Please choose a positive real # value\n"; |
149 |
+ |
} |
150 |
+ |
} |
151 |
|
|
152 |
+ |
|
153 |
|
# open the file writer |
154 |
|
open(OUTFILE, ">./$fileName") || die "Error: can't open file $fileName\n"; |
155 |
|
|
187 |
|
# now we can start building the crystals |
188 |
|
$boxLength = ($nMol*$densityConvert/$density)**(1.0/3.0); |
189 |
|
$cellLength = $boxLength / $crystalNum; |
154 |
– |
$cell2 = $cellLength*0.5; |
190 |
|
|
191 |
+ |
if ($boxx != 0) { |
192 |
+ |
if ($boxLength < $boxx) { |
193 |
+ |
print "Computed box length is smaller than requested x axis. Use more\n"; |
194 |
+ |
die "molecules."; |
195 |
+ |
} |
196 |
+ |
} else { |
197 |
+ |
$boxx = $boxLength; |
198 |
+ |
} |
199 |
+ |
if ($boxy != 0) { |
200 |
+ |
if ($boxLength < $boxy) { |
201 |
+ |
print "Computed box length is smaller than requested y axis. Use more\n"; |
202 |
+ |
die "molecules."; |
203 |
+ |
} |
204 |
+ |
} else { |
205 |
+ |
$boxy = $boxLength; |
206 |
+ |
} |
207 |
+ |
if ($boxz != 0) { |
208 |
+ |
if ($boxLength < $boxz) { |
209 |
+ |
print "Computed box length is smaller than requested z axis. Use more\n"; |
210 |
+ |
die "molecules."; |
211 |
+ |
} |
212 |
+ |
} else { |
213 |
+ |
$boxz = $boxLength; |
214 |
+ |
} |
215 |
+ |
|
216 |
+ |
$nx = int($boxx / $cellLength); |
217 |
+ |
$ny = int($boxy / $cellLength); |
218 |
+ |
$nz = int($boxz / $cellLength); |
219 |
+ |
|
220 |
|
if ($lattice == 0) { |
221 |
+ |
$nMol = 4 * $nx * $ny * $nz; |
222 |
+ |
} else { |
223 |
+ |
$nMol = $nx * $ny * $nz; |
224 |
+ |
} |
225 |
+ |
|
226 |
+ |
$newDensity = $nMol * $densityConvert / ($boxx*$boxy*$boxz); |
227 |
+ |
|
228 |
+ |
if (abs($newDensity-$density) > $tolerance) { |
229 |
+ |
print "Resetting density to $newDensity to make chosen box sides work out\n"; |
230 |
+ |
} |
231 |
+ |
$cellLengthX = $boxx/$nx; |
232 |
+ |
$cellLengthY = $boxy/$ny; |
233 |
+ |
$cellLengthZ = $boxz/$nz; |
234 |
+ |
|
235 |
+ |
$cell2X = $cellLengthX*0.5; |
236 |
+ |
$cell2Y = $cellLengthY*0.5; |
237 |
+ |
$cell2Z = $cellLengthZ*0.5; |
238 |
+ |
|
239 |
+ |
if ($lattice == 0) { |
240 |
|
# build the unit cell |
241 |
|
# molecule 0 |
242 |
|
$xCorr[0] = 0.0; |
244 |
|
$zCorr[0] = 0.0; |
245 |
|
# molecule 1 |
246 |
|
$xCorr[1] = 0.0; |
247 |
< |
$yCorr[1] = $cell2; |
248 |
< |
$zCorr[1] = $cell2; |
247 |
> |
$yCorr[1] = $cell2Y; |
248 |
> |
$zCorr[1] = $cell2Z; |
249 |
|
# molecule 2 |
250 |
< |
$xCorr[2] = $cell2; |
251 |
< |
$yCorr[2] = $cell2; |
250 |
> |
$xCorr[2] = $cell2X; |
251 |
> |
$yCorr[2] = $cell2Y; |
252 |
|
$zCorr[2] = 0.0; |
253 |
|
# molecule 3 |
254 |
< |
$xCorr[3] = $cell2; |
254 |
> |
$xCorr[3] = $cell2X; |
255 |
|
$yCorr[3] = 0.0; |
256 |
< |
$zCorr[3] = $cell2; |
256 |
> |
$zCorr[3] = $cell2Z; |
257 |
|
# assemble the lattice |
258 |
|
$counter = 0; |
259 |
< |
for ($z = 0; $z < $crystalNum; $z++) { |
260 |
< |
for ($y = 0; $y < $crystalNum; $y++) { |
261 |
< |
for ($x = 0; $x < $crystalNum; $x++) { |
259 |
> |
for ($z = 0; $z < $nx; $z++) { |
260 |
> |
for ($y = 0; $y < $ny; $y++) { |
261 |
> |
for ($x = 0; $x < $nz; $x++) { |
262 |
|
for ($uc = 0; $uc < 4; $uc++) { |
263 |
< |
$xCorr[$uc+$counter] = $xCorr[$uc] + $cellLength*$x; |
264 |
< |
$yCorr[$uc+$counter] = $yCorr[$uc] + $cellLength*$y; |
265 |
< |
$zCorr[$uc+$counter] = $zCorr[$uc] + $cellLength*$z; |
263 |
> |
$xCorr[$uc+$counter] = $xCorr[$uc] + $cellLengthX*$x; |
264 |
> |
$yCorr[$uc+$counter] = $yCorr[$uc] + $cellLengthY*$y; |
265 |
> |
$zCorr[$uc+$counter] = $zCorr[$uc] + $cellLengthZ*$z; |
266 |
|
} |
267 |
|
$counter = $counter + 4; |
268 |
|
} |
270 |
|
} |
271 |
|
|
272 |
|
} elsif ($lattice == 1) { |
273 |
< |
# build the unit cell |
274 |
< |
# molecule 0 |
275 |
< |
$xCorr[0] = $cell2; |
276 |
< |
$yCorr[0] = $cell2; |
277 |
< |
$zCorr[0] = $cell2; |
273 |
> |
# build the unit cell |
274 |
> |
# molecule 0 |
275 |
> |
$xCorr[0] = $cell2X; |
276 |
> |
$yCorr[0] = $cell2Y; |
277 |
> |
$zCorr[0] = $cell2Z; |
278 |
|
#assemble the lattice |
279 |
|
$counter = 0; |
280 |
< |
for ($z = 0; $z < $crystalNum; $z++) { |
281 |
< |
for ($y = 0; $y < $crystalNum; $y++) { |
282 |
< |
for ($x = 0; $x < $crystalNum; $x++) { |
283 |
< |
$xCorr[$counter] = $xCorr[0] + $cellLength*$x; |
284 |
< |
$yCorr[$counter] = $yCorr[0] + $cellLength*$y; |
285 |
< |
$zCorr[$counter] = $zCorr[0] + $cellLength*$z; |
280 |
> |
for ($z = 0; $z < $nx; $z++) { |
281 |
> |
for ($y = 0; $y < $ny; $y++) { |
282 |
> |
for ($x = 0; $x < $nz; $x++) { |
283 |
> |
$xCorr[$counter] = $xCorr[0] + $cellLengthX*$x; |
284 |
> |
$yCorr[$counter] = $yCorr[0] + $cellLengthY*$y; |
285 |
> |
$zCorr[$counter] = $zCorr[0] + $cellLengthZ*$z; |
286 |
|
|
287 |
|
$counter++; |
288 |
|
} |
322 |
|
|
323 |
|
# shift the box center to the origin and write out the coordinates |
324 |
|
for ($i = 0; $i < $nMol; $i++) { |
325 |
< |
$xCorr[$i] -= 0.5*$boxLength; |
326 |
< |
$yCorr[$i] -= 0.5*$boxLength; |
327 |
< |
$zCorr[$i] -= 0.5*$boxLength; |
325 |
> |
$xCorr[$i] -= 0.5*$boxx; |
326 |
> |
$yCorr[$i] -= 0.5*$boxy; |
327 |
> |
$zCorr[$i] -= 0.5*$boxz; |
328 |
|
|
329 |
|
$q0 = 1.0; |
330 |
|
$q1 = 0.0; |
390 |
|
} |
391 |
|
|
392 |
|
sub findCutoff { |
393 |
< |
$boxLength2 = 0.5*$boxLength; |
393 |
> |
if ($boxy < $boxx) { |
394 |
> |
$bm = $boxy; |
395 |
> |
} else { |
396 |
> |
$bm = $boxx; |
397 |
> |
} |
398 |
> |
if ($boxz < $bm) { |
399 |
> |
$bm = $boxz; |
400 |
> |
} |
401 |
> |
$boxLength2 = 0.5*$bm; |
402 |
|
if ($boxLength2 > $cutoff){ |
403 |
|
# the default is good |
404 |
|
} else { |
415 |
|
" <Snapshot> |
416 |
|
<FrameData> |
417 |
|
Time: 0 |
418 |
< |
Hmat: {{ $boxLength, 0, 0 }, { 0, $boxLength, 0 }, { 0, 0, $boxLength }} |
418 |
> |
Hmat: {{ $boxx, 0, 0 }, { 0, $boxy, 0 }, { 0, 0, $boxz }} |
419 |
|
</FrameData>\n"; |
420 |
|
} |
421 |
|
|