3# program that builds water boxes
5# author = "Chris Fennell
6# copyright = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
10use List::Util qw[min max];
14$mass = 2.99151E-23; # mass of H2O in grams
15$cm3ToAng3 = 1E24; # convert cm^3 to angstroms^3
16$densityConvert = $mass*$cm3ToAng3;
28getopts('hmrvd:l:c:o:w:x:y:z:');
31$opt_h = "true" if $#ARGV >= 0;
35 print "waterBoxer: builds water boxes\n\n";
36 print "usage: waterBoxer [-hmrv] [-d density] [-l lattice] [-n # waters]\n";
37 print "\t[-o file name] [-w water name] \n\n";
38 print " -h : show this message\n";
39 print " -m : print out a water.inc file (file with all water models)\n";
40 print " -r : randomize orientations\n";
41 print " -v : verbose output\n\n";
42 print " -d real : density in g/cm^3\n";
43 print " (default: 1)\n";
44 print " -l integer : 0 - face centered cubic, 1 - simple cubic\n";
45 print " (default: 0)\n";
48 print " -c real : default overlap cutoff in angstroms\n";
49 print " (default: 3.3)\n";
50 print " -o char : output file name\n";
51 print " (default: freshWater.omd)\n";
52 print " -w char : name of the water stunt double\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";
61 die " waterBoxer -d 0.997 -w SSD_RF -o ssdrfWater.omd\n";
64# set some variables to be used in the code
69 $fileName = 'freshWater.omd';
72 die "Error: $fileName cannot be \"water.inc\"\n Please choose a different name\n" if $fileName eq 'water.inc';
73 $waterFileHandle = 'WATERMD';
76 $waterFileHandle = 'OUTFILE';
79 $doRandomize = $opt_r;
85 if ($opt_d =~ /^[0-9]/) {
88 die "Error: the value for '-d' ($opt_d) is not a valid number\n Please choose a positive real # value\n";
93 if ($opt_c =~ /^[0-9]/) {
96 die "Error: the value for '-c' ($opt_c) is not a valid number\n Please choose a positive real # value\n";
102 $nothingSelected = 0;
107if ($invalidWater == 1){
108 print "Warning: \'$waterName\' is not a recognized water model name.\n";
109 print " Use the \'-m\' option to generate a \'water.inc\' with the\n";
110 print " recognized water model geometries.\n\n";
112if ($waterName eq 'DPD') {
113 # DPD waters are stand-ins for 4 water molecules
114 $density = $density * 0.25;
116if ($waterName eq 'CG2') {
117 # CG2 waters are stand-ins for 2 water molecules
118 $density = $density * 0.5;
122 $nothingSelected = 0;
123 if ($opt_l =~ /^[0-9]/) {
125 if ($lattice != 0 && $lattice != 1){
126 die "Error: the '-l' value ($opt_l) is not a valid number\n Please choose 0 or 1\n";
129 die "Error: the '-l' value ($opt_l) is not a valid number\n Please choose 0 or 1\n";
133 $nothingSelected = 0;
134 if ($opt_x =~ /^[0-9]/) {
137 die "Error: the value for '-x' ($opt_x) is not a valid number\n Please choose a positive real # value\n";
140 die "Error: the size of the box '-x' must be set\n";
143 $nothingSelected = 0;
144 if ($opt_y =~ /^[0-9]/) {
147 die "Error: the value for '-y' ($opt_y) is not a valid number\n Please choose a positive real # value\n";
150 die "Error: the size of the box '-y' must be set\n";
153 $nothingSelected = 0;
154 if ($opt_z =~ /^[0-9]/) {
157 die "Error: the value for '-z' ($opt_z) is not a valid number\n Please choose a positive real # value\n";
160 die "Error: the size of the box '-z' must be set\n";
163if ($nothingSelected == 1) {
164 die "(For help, use the \'-h\' option.)\n";
167# open the file writer
168open(OUTFILE, ">./$fileName") || die "Error: can't open file $fileName\n";
170# check to set magic lattice numbers
172 $a = (4*18.01/(0.602*$density))**(1.0/3.0);
173 $acut = $rcut * sqrt(2.0);
176 $newDensity = 4.0*18.01/(0.602*$a**3);
177 print "using density of $newDensity to match cutoff value ($rcut)\n";
179} elsif ($lattice == 1){
180 $a = (18.01/(0.602*$density))**(1.0/3.0);
184 $newDensity = 18.01/(0.602*$a**3);
185 print "using density of $newDensity to match cutoff value ($rcut)\n";
189$nx = int($nxFloat + $nxFloat/abs($nxFloat*2));
191$ny = int($nyFloat + $nyFloat/abs($nyFloat*2));
193$nz = int($nzFloat + $nzFloat/abs($nzFloat*2));
195$anew = min($boxx/$nx, $boxy/$ny, $boxz/$nz);
201$nxFloat = $boxx/$anew;
202$nx = int($nxFloat + $nxFloat/abs($nxFloat*2));
203$nyFloat = $boxy/$anew;
204$ny = int($nyFloat + $nyFloat/abs($nyFloat*2));
205$nzFloat = $boxz/$anew;
206$nz = int($nzFloat + $nzFloat/abs($nzFloat*2));
209 $nMol = 4 * $nx * $ny * $nz;
211 $nMol = $nx * $ny * $nz;
214$newDensity = $nMol * $densityConvert / ($boxx*$boxy*$boxz);
216if (abs($newDensity-$density) > $tolerance) {
217 print "Resetting density to $newDensity to make chosen box sides work out\n";
219$cellLengthX = $boxx/$nx;
220$cellLengthY = $boxy/$ny;
221$cellLengthZ = $boxz/$nz;
223$cell2X = $cellLengthX*0.5;
224$cell2Y = $cellLengthY*0.5;
225$cell2Z = $cellLengthZ*0.5;
245# assemble the lattice
247 for ($z = 0; $z < $nz; $z++) {
248 for ($y = 0; $y < $ny; $y++) {
249 for ($x = 0; $x < $nx; $x++) {
250 for ($uc = 0; $uc < 4; $uc++) {
251 $xCorr[$uc+$counter] = $xCorr[$uc] + $cellLengthX*$x;
252 $yCorr[$uc+$counter] = $yCorr[$uc] + $cellLengthY*$y;
253 $zCorr[$uc+$counter] = $zCorr[$uc] + $cellLengthZ*$z;
255 $counter = $counter + 4;
260} elsif ($lattice == 1) {
261 # build the unit cell
268 for ($z = 0; $z < $nz; $z++) {
269 for ($y = 0; $y < $ny; $y++) {
270 for ($x = 0; $x < $nx; $x++) {
271 $xCorr[$counter] = $xCorr[0] + $cellLengthX*$x;
272 $yCorr[$counter] = $yCorr[0] + $cellLengthY*$y;
273 $zCorr[$counter] = $zCorr[0] + $cellLengthZ*$z;
282print "The water box \"$fileName\" was generated.\n";
286 print "The file \"water.inc\" was generated for inclusion in \"$fileName\"\n";
291# this marks the end of the main program, below is subroutines
295 my $ret = atan2(sqrt(1 - $rad*$rad), $rad);
300 # write out the header
301 print OUTFILE "<OpenMD version=1>\n";
305 print OUTFILE " <StuntDoubles>\n";
307 # shift the box center to the origin and write out the coordinates
308 for ($i = 0; $i < $nMol; $i++) {
309 $xCorr[$i] -= 0.5*$boxx;
310 $yCorr[$i] -= 0.5*$boxy;
311 $zCorr[$i] -= 0.5*$boxz;
318 if ($doRandomize == 1){
319 $cosTheta = 2.0*rand() - 1.0;
320 $theta = acos($cosTheta);
321 $phi = 2.0*3.14159265359*rand();
322 $psi = 2.0*3.14159265359*rand();
324 $q0 = cos(0.5*$theta)*cos(0.5*($phi + $psi));
325 $q1 = sin(0.5*$theta)*cos(0.5*($phi - $psi));
326 $q2 = sin(0.5*$theta)*sin(0.5*($phi - $psi));
327 $q3 = cos(0.5*$theta)*sin(0.5*($phi + $psi));
330 print OUTFILE "$i\tpq\t$xCorr[$i] $yCorr[$i] $zCorr[$i] ";
331 print OUTFILE "$q0 $q1 $q2 $q3\n";
334 print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OpenMD>\n";
338 print OUTFILE " <MetaData>\n";
340# print the water model or includes
342 print OUTFILE "#include \"water.inc\"";
346 printFakeWater() if $invalidWater == 1;
348# now back to the metaData output
349 print OUTFILE "\n\ncomponent{
350 type = \"$waterName\";
355forceField = \"Water\";
356electrostaticSummationMethod = \"shifted_force\";
357electrostaticScreeningMethod = \"damped\";
358cutoffRadius = $cutoff;
359switchingRadius = $cutoff;
360dampingAlpha = $alpha;
387 $boxLength2 = 0.5*$bm;
388 if ($boxLength2 > $cutoff){
389 # the default is good
391 $cutoff = int($boxLength2);
400 Hmat: {{ $boxx, 0, 0 }, { 0, $boxy, 0 }, { 0, 0, $boxz }}
405 open(WATERMD, ">./water.inc") || die "Error: can't open file water.inc\n";
406 $waterFileHandle = 'WATERMD';
408 print WATERMD "#ifndef _WATER_INC_\n#define _WATER_INC_\n";
431 print WATERMD "\n\n#endif";
435 print $waterFileHandle "\n\nmolecule{
440 position(0.0, 0.0, 0.0);
446 print $waterFileHandle "\n\nmolecule{
451 position(0.0, 0.0, 0.0);
457 print $waterFileHandle "\n\nmolecule{
462 position( 0.0, 0.0, 0.0 );
463 orientation( 0.0, 0.0, 0.0 );
469 print $waterFileHandle "\n\nmolecule{
474 position( 0.0, 0.0, 0.0 );
475 orientation( 0.0, 0.0, 0.0 );
481 print $waterFileHandle "\n\nmolecule{
486 position( 0.0, 0.0, 0.0 );
487 orientation( 0.0, 0.0, 0.0 );
493 print $waterFileHandle "\n\nmolecule{
498 position( 0.0, 0.0, 0.0 );
499 orientation( 0.0, 0.0, 0.0 );
505 print $waterFileHandle "\n\nmolecule{
510 position( 0.0, 0.0, 0.0 );
511 orientation( 0.0, 0.0, 0.0 );
517 print $waterFileHandle "\n\nmolecule{
522 position( 0.0, 0.0, 0.0 );
523 orientation( 0.0, 0.0, 0.0 );
529 print $waterFileHandle "\n\nmolecule{
534 position( 0.0, 0.0, 0.0 );
535 orientation( 0.0, 0.0, 0.0 );
539 position( 0.0, 0.0, 0.5 );
550 print $waterFileHandle "\n\nmolecule{
555 position( 0.0, 0.0, -0.06556 );
559 position( 0.0, 0.75695, 0.52032 );
563 position( 0.0, -0.75695, 0.52032 );
574 print $waterFileHandle "\n\nmolecule{
579 position( 0.0, 0.0, -0.06556 );
583 position( 0.0, 0.75695, 0.52032 );
587 position( 0.0, -0.75695, 0.52032 );
591 position( 0.0, 0.0, 0.08444 );
602 print $waterFileHandle "\n\nmolecule{
603 name = \"TIP4P-Ice\";
606 type = \"O_TIP4P-Ice\";
607 position( 0.0, 0.0, -0.06556 );
610 type = \"H_TIP4P-Ice\";
611 position( 0.0, 0.75695, 0.52032 );
614 type = \"H_TIP4P-Ice\";
615 position( 0.0, -0.75695, 0.52032 );
618 type = \"EP_TIP4P-Ice\";
619 position( 0.0, 0.0, 0.09214 );
630 print $waterFileHandle "\n\nmolecule{
634 type = \"O_TIP4P-Ew\";
635 position( 0.0, 0.0, -0.06556 );
638 type = \"H_TIP4P-Ew\";
639 position( 0.0, 0.75695, 0.52032 );
642 type = \"H_TIP4P-Ew\";
643 position( 0.0, -0.75695, 0.52032 );
646 type = \"EP_TIP4P-Ew\";
647 position( 0.0, 0.0, 0.05944 );
658 print $waterFileHandle "\n\nmolecule{
659 name = \"TIP4P-2005\";
662 type = \"O_TIP4P-2005\";
663 position( 0.0, 0.0, -0.06556 );
666 type = \"H_TIP4P-2005\";
667 position( 0.0, 0.75695, 0.52032 );
670 type = \"H_TIP4P-2005\";
671 position( 0.0, -0.75695, 0.52032 );
674 type = \"EP_TIP4P-2005\";
675 position( 0.0, 0.0, 0.08904 );
686 print $waterFileHandle "\n\nmolecule{
691 position( 0.0, 0.0, -0.06556 );
695 position( 0.0, 0.75695, 0.52032 );
699 position( 0.0, -0.75695, 0.52032 );
703 position( 0.57154, 0.0, -0.46971 );
707 position( -0.57154, 0.0, -0.46971 );
711 members(0, 1, 2, 3, 4);
718 print $waterFileHandle "\n\nmolecule{
722 type = \"O_TIP5P-E\";
723 position( 0.0, 0.0, -0.06556 );
727 position( 0.0, 0.75695, 0.52032 );
731 position( 0.0, -0.75695, 0.52032 );
735 position( 0.57154, 0.0, -0.46971 );
739 position( -0.57154, 0.0, -0.46971 );
743 members(0, 1, 2, 3, 4);
750 print $waterFileHandle "\n\nmolecule{
755 position( 0.0, 0.0, 0.0 );
759 position( 0.0, 0.576029, 0.79283665 );
763 position( 0.0, -0.576029, 0.79283665 );
767 position( 0.0, 0.23, 0.0 );
771 position( 0.732813007, -0.50364843, 0.0 );
775 position( -0.732813007, -0.50364843, 0.0 );
779 members(0, 1, 2, 3, 4, 5);
785print $waterFileHandle "\n\nmolecule{
790 position( 0.0, 0.0, -0.06461 );
794 position( 0.0, 0.81649, 0.51275 );
798 position( 0.0, -0.81649, 0.51275 );
809print $waterFileHandle "\n\nmolecule{
814 position( 0.0, 0.0, -0.06461 );
818 position( 0.0, 0.81649, 0.51275 );
822 position( 0.0, -0.81649, 0.51275 );
834print $waterFileHandle "\n\nmolecule{
839 position( 0.0, 0.0, -0.06461 );
843 position( 0.0, 0.81649, 0.51275 );
847 position( 0.0, -0.81649, 0.51275 );
858 print $waterFileHandle "\n\nmolecule{
863 position(0.0, 0.0, 0.0);
869 print $waterFileHandle "\n\nmolecule{
874 position(0.0, 0.0, 0.0);
880 print $waterFileHandle "\n\nmolecule{
881 name = \"$waterName\";
884 type = \"$waterName\";
885 position(0.0, 0.0, 0.0);
892 if ($waterName eq 'Cl-') { $waterCase = 0; }
893 elsif ($waterName eq 'Na+') { $waterCase = 1; }
894 elsif ($waterName eq 'SSD_E') { $waterCase = 2; }
895 elsif ($waterName eq 'SSD_RF') { $waterCase = 3; }
896 elsif ($waterName eq 'SSD') { $waterCase = 4; }
897 elsif ($waterName eq 'SSD1') { $waterCase = 5; }
898 elsif ($waterName eq 'TIP3P') { $waterCase = 6; }
899 elsif ($waterName eq 'TIP4P') { $waterCase = 7; }
900 elsif ($waterName eq 'TIP4P-Ew') { $waterCase = 8; }
901 elsif ($waterName eq 'TIP5P') { $waterCase = 9; }
902 elsif ($waterName eq 'TIP5P-E') { $waterCase = 10; }
903 elsif ($waterName eq 'SPCE') { $waterCase = 11; }
904 elsif ($waterName eq 'SPC') { $waterCase = 12; }
905 elsif ($waterName eq 'DPD') { $waterCase = 13; }
906 elsif ($waterName eq 'CG2') { $waterCase = 14; }
907 elsif ($waterName eq 'SSDQ') { $waterCase = 15; }
908 elsif ($waterName eq 'SSDQO') { $waterCase = 16; }
909 elsif ($waterName eq 'TIP4P-Ice'){ $waterCase = 17; }
910 elsif ($waterName eq 'TIP4P-2005'){ $switchCase = 18; }
911 elsif ($waterName eq 'SPC-HW') { $waterCase == 19; }
912 elsif ($waterName eq 'NE6') { $waterCase == 20; }
913 else { $invalidWater = 1; }
917 if ($waterCase == 0) { printCl(); }
918 elsif ($waterCase == 1) { printNa(); }
919 elsif ($waterCase == 2) { printSSD_E(); }
920 elsif ($waterCase == 3) { printSSD_RF(); }
921 elsif ($waterCase == 4) { printSSD(); }
922 elsif ($waterCase == 5) { printSSD1(); }
923 elsif ($waterCase == 6) { printTIP3P(); }
924 elsif ($waterCase == 7) { printTIP4P(); }
925 elsif ($waterCase == 8) { printTIP4PEw(); }
926 elsif ($waterCase == 9) { printTIP5P(); }
927 elsif ($waterCase == 10) { printTIP5PE(); }
928 elsif ($waterCase == 11) { printSPCE(); }
929 elsif ($waterCase == 12) { printSPC(); }
930 elsif ($waterCase == 13) { printDPD(); }
931 elsif ($waterCase == 14) { printCG2(); }
932 elsif ($waterCase == 15) { printSSDQ(); }
933 elsif ($waterCase == 16) { printSSDQO(); }
934 elsif ($waterCase == 17) { printTIP4PIce();}
935 elsif ($waterCase == 18) { printTIP4P2005();}
936 elsif ($waterCase == 19) { printSPCHW(); }
937 elsif ($waterCase == 20) { printNE6(); }