| 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 |
|
|
