applications/utilities/waterBoxer.in

#!@PERLINTERP@ -w

# program that builds water boxes

# author    = "Chris Fennell
# version   = "$Revision: 1.1 $"
# date      = "$Date: 2006-09-01 21:08:03 $"
# copyright = "Copyright (c) 2006 by the University of Notre Dame"
# license   = "OOPSE"

use Getopt::Std;

$tolerance = 1.0E-8;
$mass = 2.99151E-23; # mass of H2O in grams
$cm3ToAng3 = 1E24;   # convert cm^3 to angstroms^3
$densityConvert = $mass*$cm3ToAng3;
$lattice = 0;
$nMol = 500;
$density = 1.0;
$doRandomize = 0;
$cutoff = 12;
$alpha = 0.2125;
$alphaInt = 0.5125;
$alphaSlope = 0.025;

# get our options
getopts('hrvd:l:n:w:');

# if we don't have a filename, drop to -h
$opt_h = 'true' if $#ARGV != 0;

# our option output 
if ($opt_h){
    print "waterBoxer: builds water boxes\n\n";
    print "usage: waterBoxer [-hv] [-d density] [-l lattice] [-n # waters]\n";
    print "\t[-w water name] [file name]\n\n";
    print "  -h : show this message\n";
    print "  -r : randomize orientations\n";
    print "  -v : verbose output\n\n";
    print "  -d real    : density in g/cm^3\n";
    print "                 (default: 1)\n";
    print "  -l integer : 0 - face centered cubic, 1 - simple cubic\n";
    print "                 (default: 0)\n";
    print "  -n integer : # of water molecules\n";
    print "                 (default: 500)\n";
    print "  -w char    : name of the water stunt double\n";
    print "                 (default: SPCE)\n";
    print "Example:\n";
    die   "   waterBoxer -d 0.997 -n 864 -w SSD_RF ssdrfWater.md\n";
}

# set some variables to be used in the code
$fileName = $ARGV[0];
if (defined($fileName)){
} else {
    $fileName = 'waterBox';
}
if ($opt_r){
  $doRandomize = $opt_r;
}
if (defined($opt_w)){
  $waterName = $opt_w;
} else {
  $waterName = 'SPCE';
}

if (defined($opt_d)){
  if ($opt_d =~ /^[0-9]/) {
    $density = $opt_d;
  } else {
    die "\t-d value ($opt_d) is not a valid number\n\tPlease choose a positive real # value\n";
  }
}
if (defined($opt_l)){
  if ($opt_l =~ /^[0-9]/) {
    $lattice = $opt_l;
    if ($lattice != 0 && $lattice != 1){
      die "\t-l value ($opt_l) is not a valid number\n\tPlease choose 0 or 1\n";
    }
  } else {
    die "\t-l value ($opt_l) is not a valid number\n\tPlease choose 0 or 1\n";
  }
}
if (defined($opt_n)){
  if ($opt_n =~ /^[0-9]/) {
    $nMol = $opt_n;
  } else { 
    die "\t-n value ($opt_n) is not a valid number\n\tPlease choose a non-negative integer\n";
  }
}

# open the file writer
open(OUTFILE, ">./$fileName") || die "\tError: can't open file $fileName\n";

# check to set magic lattice numbers
if ($lattice == 0){
  $crystalNumReal = ($nMol/4.0)**(1.0/3.0);
  $crystalNum = int($crystalNumReal + $tolerance);
  $remainder = $crystalNumReal - $crystalNum;
  
  # if crystalNumReal wasn't an integer, we bump the crystal to the next
  # magic number
  if ($remainder > $tolerance){
    $crystalNum = $crystalNum + 1;
    $newMol = 4 * $crystalNum**3;
    print "WARNING: The number chosen ($nMol) failed to build a clean\n";
    print "fcc lattice. The number of molecules has been increased to\n";
    print "the next magic number ($newMol).\n";
    $nMol = $newMol;
  }
} elsif ($lattice == 1){
  $crystalNumReal = ($nMol/1.0)**(1.0/3.0);
  $crystalNum = int($crystalNumReal);
  $remainder = $crystalNumReal - $crystalNum;
  
  # again, if crystalNumReal wasn't an integer, we bump the crystal to the next
  # magic number
  if ($remainder > $tolerance){
    $crystalNum = $crystalNum + 1;
    $newMol = $crystalNum**3;
    print "WARNING: The number chosen ($nMol) failed to build a clean\n";
    print "simple cubic lattice. The number of molecules has been\n";
    print "increased to the next magic number ($newMol).\n";
    $nMol = $newMol;
  }
}

# now we can start building the crystals
$boxLength = ($nMol*$densityConvert/$density)**(1.0/3.0);
$cellLength = $boxLength / $crystalNum;
$cell2 = $cellLength*0.5;

if ($lattice == 0) {
# build the unit cell
# molecule 0
  $xCorr[0] = 0.0;
  $yCorr[0] = 0.0;
  $zCorr[0] = 0.0;
# molecule 1
  $xCorr[1] = 0.0;
  $yCorr[1] = $cell2;
  $zCorr[1] = $cell2;
# molecule 2
  $xCorr[2] = $cell2;
  $yCorr[2] = $cell2;
  $zCorr[2] = 0.0;
# molecule 3
  $xCorr[3] = $cell2;
  $yCorr[3] = 0.0;
  $zCorr[3] = $cell2;
# assemble the lattice
  $counter = 0;
  for ($z = 0; $z < $crystalNum; $z++) {
    for ($y = 0; $y < $crystalNum; $y++) {
      for ($x = 0; $x < $crystalNum; $x++) {
        for ($uc = 0; $uc < 4; $uc++) {
          $xCorr[$uc+$counter] = $xCorr[$uc] + $cellLength*$x;
          $yCorr[$uc+$counter] = $yCorr[$uc] + $cellLength*$y;
          $zCorr[$uc+$counter] = $zCorr[$uc] + $cellLength*$z;
        }
        $counter = $counter + 4;
      }
    }
  }
  
} elsif ($lattice == 1) {
# build the unit cell
# molecule 0
  $xCorr[0] = $cell2;
  $yCorr[0] = $cell2;
  $zCorr[0] = $cell2;
#assemble the lattice
  $counter = 0;
  for ($z = 0; $z < $crystalNum; $z++) {
    for ($y = 0; $y < $crystalNum; $y++) {
      for ($x = 0; $x < $crystalNum; $x++) {
        $xCorr[$counter] = $xCorr[0] + $cellLength*$x;
        $yCorr[$counter] = $yCorr[0] + $cellLength*$y;
        $zCorr[$counter] = $zCorr[0] + $cellLength*$z;
        
        $counter++;
      }
    }
  }
}

writeOutFile();

# this marks the end of the main program, below is subroutines

sub acos {
  my ($rad) = @_;
  my $ret = atan2(sqrt(1 - $rad*$rad), $rad);
  return $ret;
}

sub writeOutFile{
  # write out the header
  print OUTFILE "<OOPSE version=4>\n";
  findCutoff();
  findAlpha();
  printMetaData();
  printFrameData();
  print OUTFILE "    <StuntDoubles>\n";
  
  # shift the box center to the origin and write out the coordinates
  for ($i = 0; $i < $nMol; $i++) {
    $xCorr[$i] -= 0.5*$boxLength;
    $yCorr[$i] -= 0.5*$boxLength;
    $zCorr[$i] -= 0.5*$boxLength;
    
    $q0 = 1.0;
    $q1 = 0.0;
    $q2 = 0.0;
    $q3 = 0.0;
    
    if ($doRandomize == 1){
      $cosTheta = 2.0*rand() - 1.0;
      $theta = acos($cosTheta);
      $phi = 2.0*3.14159265359*rand();
      $psi = 2.0*3.14159265359*rand();
      
      $q0 = cos(0.5*$theta)*cos(0.5*($phi + $psi));
      $q1 = sin(0.5*$theta)*cos(0.5*($phi - $psi));
      $q2 = sin(0.5*$theta)*sin(0.5*($phi - $psi));
      $q3 = cos(0.5*$theta)*sin(0.5*($phi + $psi));
    }
    
    print OUTFILE "$i\tpq\t$xCorr[$i] $yCorr[$i] $zCorr[$i] ";
    print OUTFILE "$q0 $q1 $q2 $q3\n";
  }

  print OUTFILE "    </StuntDoubles>\n  </Snapshot>\n</OOPSE>\n";
}

sub printMetaData {
print OUTFILE"  <MetaData>
#include \"water.md\"

component{
  type = \"$waterName\";
  nMol = $nMol;
}


ensemble = NVE;
forceField = \"DUFF\";
electrostaticSummationMethod = \"shifted_force\";
electrostaticScreeningMethod = \"damped\";
dampingAlpha = $alpha;
cutoffRadius = $cutoff;

targetTemp = 300;
targetPressure = 1.0;

tauThermostat = 1e3;
tauBarostat = 1e4;

dt = 2.0;
runTime = 1e3;

tempSet = \"true\";
thermalTime = 10;
sampleTime = 100;
statusTime = 2;
  </MetaData>\n";
}

sub findCutoff{
  $boxLength2 = 0.5*$boxLength;
  if ($boxLength2 > $cutoff){
    # the default is good
  } else {
    $cutoff = int($boxLength2);
  }
}

sub findAlpha{
  $alpha = $alphaInt - $cutoff*$alphaSlope;
}

sub printFrameData{
print OUTFILE 
"  <Snapshot>
    <FrameData>
        Time: 0
        Hmat: {{ $boxLength, 0, 0 }, { 0, $boxLength, 0 }, { 0, 0, $boxLength }}
    </FrameData>\n";
}
Revision:	1034
Committed:	Fri Sep 1 21:08:03 2006 UTC (18 years, 8 months ago) by chrisfen
File size:	7261 byte(s)
Log Message:	added the waterBoxer script
#	User	Rev	Content
1	chrisfen	1034	#!@PERLINTERP@ -w
2
3			# program that builds water boxes
4
5			# author = "Chris Fennell
6			# version = "$Revision: 1.1 $"
7			# date = "$Date: 2006-09-01 21:08:03 $"
8			# copyright = "Copyright (c) 2006 by the University of Notre Dame"
9			# license = "OOPSE"
10
11			use Getopt::Std;
12
13			$tolerance = 1.0E-8;
14			$mass = 2.99151E-23; # mass of H2O in grams
15			$cm3ToAng3 = 1E24; # convert cm^3 to angstroms^3
16			$densityConvert = $mass*$cm3ToAng3;
17			$lattice = 0;
18			$nMol = 500;
19			$density = 1.0;
20			$doRandomize = 0;
21			$cutoff = 12;
22			$alpha = 0.2125;
23			$alphaInt = 0.5125;
24			$alphaSlope = 0.025;
25
26			# get our options
27			getopts('hrvd:l:n:w:');
28
29			# if we don't have a filename, drop to -h
30			$opt_h = 'true' if $#ARGV != 0;
31
32			# our option output
33			if ($opt_h){
34			print "waterBoxer: builds water boxes\n\n";
35			print "usage: waterBoxer [-hv] [-d density] [-l lattice] [-n # waters]\n";
36			print "\t[-w water name] [file name]\n\n";
37			print " -h : show this message\n";
38			print " -r : randomize orientations\n";
39			print " -v : verbose output\n\n";
40			print " -d real : density in g/cm^3\n";
41			print " (default: 1)\n";
42			print " -l integer : 0 - face centered cubic, 1 - simple cubic\n";
43			print " (default: 0)\n";
44			print " -n integer : # of water molecules\n";
45			print " (default: 500)\n";
46			print " -w char : name of the water stunt double\n";
47			print " (default: SPCE)\n";
48			print "Example:\n";
49			die " waterBoxer -d 0.997 -n 864 -w SSD_RF ssdrfWater.md\n";
50			}
51
52			# set some variables to be used in the code
53			$fileName = $ARGV[0];
54			if (defined($fileName)){
55			} else {
56			$fileName = 'waterBox';
57			}
58			if ($opt_r){
59			$doRandomize = $opt_r;
60			}
61			if (defined($opt_w)){
62			$waterName = $opt_w;
63			} else {
64			$waterName = 'SPCE';
65			}
66
67			if (defined($opt_d)){
68			if ($opt_d =~ /^[0-9]/) {
69			$density = $opt_d;
70			} else {
71			die "\t-d value ($opt_d) is not a valid number\n\tPlease choose a positive real # value\n";
72			}
73			}
74			if (defined($opt_l)){
75			if ($opt_l =~ /^[0-9]/) {
76			$lattice = $opt_l;
77			if ($lattice != 0 && $lattice != 1){
78			die "\t-l value ($opt_l) is not a valid number\n\tPlease choose 0 or 1\n";
79			}
80			} else {
81			die "\t-l value ($opt_l) is not a valid number\n\tPlease choose 0 or 1\n";
82			}
83			}
84			if (defined($opt_n)){
85			if ($opt_n =~ /^[0-9]/) {
86			$nMol = $opt_n;
87			} else {
88			die "\t-n value ($opt_n) is not a valid number\n\tPlease choose a non-negative integer\n";
89			}
90			}
91
92			# open the file writer
93			open(OUTFILE, ">./$fileName") \|\| die "\tError: can't open file $fileName\n";
94
95			# check to set magic lattice numbers
96			if ($lattice == 0){
97			$crystalNumReal = ($nMol/4.0)**(1.0/3.0);
98			$crystalNum = int($crystalNumReal + $tolerance);
99			$remainder = $crystalNumReal - $crystalNum;
100
101			# if crystalNumReal wasn't an integer, we bump the crystal to the next
102			# magic number
103			if ($remainder > $tolerance){
104			$crystalNum = $crystalNum + 1;
105			$newMol = 4 * $crystalNum**3;
106			print "WARNING: The number chosen ($nMol) failed to build a clean\n";
107			print "fcc lattice. The number of molecules has been increased to\n";
108			print "the next magic number ($newMol).\n";
109			$nMol = $newMol;
110			}
111			} elsif ($lattice == 1){
112			$crystalNumReal = ($nMol/1.0)**(1.0/3.0);
113			$crystalNum = int($crystalNumReal);
114			$remainder = $crystalNumReal - $crystalNum;
115
116			# again, if crystalNumReal wasn't an integer, we bump the crystal to the next
117			# magic number
118			if ($remainder > $tolerance){
119			$crystalNum = $crystalNum + 1;
120			$newMol = $crystalNum**3;
121			print "WARNING: The number chosen ($nMol) failed to build a clean\n";
122			print "simple cubic lattice. The number of molecules has been\n";
123			print "increased to the next magic number ($newMol).\n";
124			$nMol = $newMol;
125			}
126			}
127
128			# now we can start building the crystals
129			$boxLength = ($nMol$densityConvert/$density)*(1.0/3.0);
130			$cellLength = $boxLength / $crystalNum;
131			$cell2 = $cellLength*0.5;
132
133			if ($lattice == 0) {
134			# build the unit cell
135			# molecule 0
136			$xCorr[0] = 0.0;
137			$yCorr[0] = 0.0;
138			$zCorr[0] = 0.0;
139			# molecule 1
140			$xCorr[1] = 0.0;
141			$yCorr[1] = $cell2;
142			$zCorr[1] = $cell2;
143			# molecule 2
144			$xCorr[2] = $cell2;
145			$yCorr[2] = $cell2;
146			$zCorr[2] = 0.0;
147			# molecule 3
148			$xCorr[3] = $cell2;
149			$yCorr[3] = 0.0;
150			$zCorr[3] = $cell2;
151			# assemble the lattice
152			$counter = 0;
153			for ($z = 0; $z < $crystalNum; $z++) {
154			for ($y = 0; $y < $crystalNum; $y++) {
155			for ($x = 0; $x < $crystalNum; $x++) {
156			for ($uc = 0; $uc < 4; $uc++) {
157			$xCorr[$uc+$counter] = $xCorr[$uc] + $cellLength*$x;
158			$yCorr[$uc+$counter] = $yCorr[$uc] + $cellLength*$y;
159			$zCorr[$uc+$counter] = $zCorr[$uc] + $cellLength*$z;
160			}
161			$counter = $counter + 4;
162			}
163			}
164			}
165
166			} elsif ($lattice == 1) {
167			# build the unit cell
168			# molecule 0
169			$xCorr[0] = $cell2;
170			$yCorr[0] = $cell2;
171			$zCorr[0] = $cell2;
172			#assemble the lattice
173			$counter = 0;
174			for ($z = 0; $z < $crystalNum; $z++) {
175			for ($y = 0; $y < $crystalNum; $y++) {
176			for ($x = 0; $x < $crystalNum; $x++) {
177			$xCorr[$counter] = $xCorr[0] + $cellLength*$x;
178			$yCorr[$counter] = $yCorr[0] + $cellLength*$y;
179			$zCorr[$counter] = $zCorr[0] + $cellLength*$z;
180
181			$counter++;
182			}
183			}
184			}
185			}
186
187			writeOutFile();
188
189			# this marks the end of the main program, below is subroutines
190
191			sub acos {
192			my ($rad) = @_;
193			my $ret = atan2(sqrt(1 - $rad*$rad), $rad);
194			return $ret;
195			}
196
197			sub writeOutFile{
198			# write out the header
199			print OUTFILE "<OOPSE version=4>\n";
200			findCutoff();
201			findAlpha();
202			printMetaData();
203			printFrameData();
204			print OUTFILE " <StuntDoubles>\n";
205
206			# shift the box center to the origin and write out the coordinates
207			for ($i = 0; $i < $nMol; $i++) {
208			$xCorr[$i] -= 0.5*$boxLength;
209			$yCorr[$i] -= 0.5*$boxLength;
210			$zCorr[$i] -= 0.5*$boxLength;
211
212			$q0 = 1.0;
213			$q1 = 0.0;
214			$q2 = 0.0;
215			$q3 = 0.0;
216
217			if ($doRandomize == 1){
218			$cosTheta = 2.0*rand() - 1.0;
219			$theta = acos($cosTheta);
220			$phi = 2.03.14159265359rand();
221			$psi = 2.03.14159265359rand();
222
223			$q0 = cos(0.5$theta)cos(0.5*($phi + $psi));
224			$q1 = sin(0.5$theta)cos(0.5*($phi - $psi));
225			$q2 = sin(0.5$theta)sin(0.5*($phi - $psi));
226			$q3 = cos(0.5$theta)sin(0.5*($phi + $psi));
227			}
228
229			print OUTFILE "$i\tpq\t$xCorr[$i] $yCorr[$i] $zCorr[$i] ";
230			print OUTFILE "$q0 $q1 $q2 $q3\n";
231			}
232
233			print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OOPSE>\n";
234			}
235
236			sub printMetaData {
237			print OUTFILE" <MetaData>
238			#include \"water.md\"
239
240			component{
241			type = \"$waterName\";
242			nMol = $nMol;
243			}
244
245
246			ensemble = NVE;
247			forceField = \"DUFF\";
248			electrostaticSummationMethod = \"shifted_force\";
249			electrostaticScreeningMethod = \"damped\";
250			dampingAlpha = $alpha;
251			cutoffRadius = $cutoff;
252
253			targetTemp = 300;
254			targetPressure = 1.0;
255
256			tauThermostat = 1e3;
257			tauBarostat = 1e4;
258
259			dt = 2.0;
260			runTime = 1e3;
261
262			tempSet = \"true\";
263			thermalTime = 10;
264			sampleTime = 100;
265			statusTime = 2;
266			</MetaData>\n";
267			}
268
269			sub findCutoff{
270			$boxLength2 = 0.5*$boxLength;
271			if ($boxLength2 > $cutoff){
272			# the default is good
273			} else {
274			$cutoff = int($boxLength2);
275			}
276			}
277
278			sub findAlpha{
279			$alpha = $alphaInt - $cutoff*$alphaSlope;
280			}
281
282			sub printFrameData{
283			print OUTFILE
284			" <Snapshot>
285			<FrameData>
286			Time: 0
287			Hmat: {{ $boxLength, 0, 0 }, { 0, $boxLength, 0 }, { 0, 0, $boxLength }}
288			</FrameData>\n";
289			}