[ViewVC] Annotation of: OpenMD/trunk/src/applications/utilities/waterBoxer

#!/usr/bin/env perl

# program that builds water boxes

# author    = "Chris Fennell
# version   = "$Revision: 1.1 $"
# date      = "$Date: 2006-10-10 14:08:10 $"
# 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;
$invalidWater = 0;
$waterCase = -1;
$nothingSelected = 1;

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

# just to test opt_h
$opt_h = "true" if $#ARGV >= 0;

# our option output 
if ($opt_h){
    print "waterBoxer: builds water boxes\n\n";
    print "usage: waterBoxer [-hmrv] [-d density] [-l lattice] [-n # waters]\n";
    print "\t[-o file name] [-w water name] \n\n";
    print "  -h : show this message\n";
    print "  -m : print out a water.md file (file with all water models)\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 "  -o char    : output file name\n";
    print "                 (default: freshWater.md)\n";
    print "  -w char    : name of the water stunt double\n";
    print "                 (default: SPCE)\n\n";
    print "Example:\n";
    die   "   waterBoxer -d 0.997 -n 864 -w SSD_RF -o ssdrfWater.md\n";
}

# set some variables to be used in the code
if (defined($opt_o)){
  $fileName = $opt_o;
  $nothingSelected = 0;
} else {
    $fileName = 'freshWater.md';
}
if ($opt_m){
  die "Error: $fileName cannot be \"water.md\"\n       Please choose a different name\n" if $fileName eq 'water.md';
  $waterFileHandle = 'WATERMD';
  $nothingSelected = 0;
} else {
  $waterFileHandle = 'OUTFILE';
}
if ($opt_r){
  $doRandomize = $opt_r;
  $nothingSelected = 0;
}

if (defined($opt_w)){
  $waterName = $opt_w;
  $nothingSelected = 0;
} else {
  $waterName = 'SPCE';
}
validateWater();
if ($invalidWater == 1){
  print "Warning: \'$waterName\' is not a recognized water model name.\n";
  print "         Use the \'-m\' option to generate a \'water.md\' with the\n";
  print "         recognized water model geometries.\n\n"; 
}
if (defined($opt_d)){
  $nothingSelected = 0;
  if ($opt_d =~ /^[0-9]/) {
    $density = $opt_d;
  } else {
    die "Error: the value for '-d' ($opt_d) is not a valid number\n       Please choose a positive real # value\n";
  }
}
if (defined($opt_l)){
  $nothingSelected = 0;
  if ($opt_l =~ /^[0-9]/) {
    $lattice = $opt_l;
    if ($lattice != 0 && $lattice != 1){
      die "Error: the '-l' value ($opt_l) is not a valid number\n       Please choose 0 or 1\n";
    }
  } else {
    die "Error: the '-l' value ($opt_l) is not a valid number\n       Please choose 0 or 1\n";
  }
}
if (defined($opt_n)){
  $nothingSelected = 0;
  if ($opt_n =~ /^[0-9]/) {
    $nMol = $opt_n;
  } else { 
    die "Error: the '-n' value ($opt_n) is not a valid number\n       Please choose a non-negative integer\n";
  }
}

# open the file writer
open(OUTFILE, ">./$fileName") || die "Error: 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 fcc lattice.\n";
    print "         The number of molecules has been increased to the next magic number ($newMol).\n\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 simple cubic lattice.\n";
    print "         The number of molecules has been increased to the next magic number ($newMol).\n\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();
print "The water box \"$fileName\" was generated.\n";

if ($opt_m){
  printWaterMD();
  print "The file \"water.md\" was generated for inclusion in \"$fileName\"\n";
}

if ($nothingSelected == 1) {
    print "(For help, use the \'-h\' option.)\n";
}


# 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>\n";

# print the water model or includes
  if ($opt_m){
    print OUTFILE "#include \"water.md\"";
  } else {
    printWaterModel();
  }
  printFakeWater() if $invalidWater == 1;

# now back to the metaData output
  print OUTFILE "\n\ncomponent{
  type = \"$waterName\";
  nMol = $nMol;
}

ensemble = NVE;
forceField = \"DUFF\";
electrostaticSummationMethod = \"shifted_force\";
electrostaticScreeningMethod = \"damped\";
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";
}

sub printWaterMD {
  open(WATERMD, ">./water.md") || die "Error: can't open file water.md\n";
  $waterFileHandle = 'WATERMD';

  print WATERMD "#ifndef _WATER_MD_\n#define _WATER_MD_\n";
  printCl();
  printNa();
  printSSD_E();
  printSSD_RF();
  printSSD();
  printSSD1();
  printTRED();
  printTIP3P();
  printTIP4P();
  printTIP4PEw();
  printTIP5P();
  printTIP5PE();
  printSPCE();
  printSPC();
  printDPD();
  print WATERMD "\n\n#endif";
}

sub printCl {
  print $waterFileHandle "\n\nmolecule{
  name = \"Cl-\";
  
  atom[0]{
    type = \"Cl-\";
    position(0.0, 0.0, 0.0);
  }
}"
}

sub printNa {
  print $waterFileHandle "\n\nmolecule{
  name = \"Na+\";
  
  atom[0]{
    type = \"Na+\";
    position(0.0, 0.0, 0.0);
  }
}"
}

sub printSSD_E {
  print $waterFileHandle "\n\nmolecule{
  name = \"SSD_E\";
  
  atom[0]{
    type = \"SSD_E\";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}"
}

sub printSSD_RF {
  print $waterFileHandle "\n\nmolecule{
  name = \"SSD_RF\";
  
  atom[0]{
    type = \"SSD_RF\";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}"
}

sub printSSD {
  print $waterFileHandle "\n\nmolecule{
  name = \"SSD\";
  
  atom[0]{
    type = \"SSD\";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}"
}

sub printSSD1 {
  print $waterFileHandle "\n\nmolecule{
  name = \"SSD1\";
  
  atom[0]{
    type = \"SSD1\";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}"
}

sub printTRED {
  print $waterFileHandle "\n\nmolecule{
  name = \"TRED\";
  
  atom[0]{
    type = \"TRED\";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
  atom[1]{
    type = \"EP_TRED\";
    position( 0.0, 0.0, 0.5 );
  }

  rigidBody[0]{
    members(0, 1);
  }

  cutoffGroup{
    members(0, 1);
  }
}"
}

sub printTIP3P {
  print $waterFileHandle "\n\nmolecule{
  name = \"TIP3P\";
  
  atom[0]{
    type = \"O_TIP3P\";
    position( 0.0, 0.0, -0.06556 );
  }
  atom[1]{
    type = \"H_TIP3P\";
    position( 0.0, 0.75695, 0.52032 );
  }
  atom[2]{
    type = \"H_TIP3P\";
    position( 0.0, -0.75695, 0.52032 );
  }

  rigidBody[0]{
    members(0, 1, 2);
  }

  cutoffGroup{
    members(0, 1, 2);
  }
}"
}

sub printTIP4P {
  print $waterFileHandle "\n\nmolecule{
  name = \"TIP4P\";
  
  atom[0]{
    type = \"O_TIP4P\";
    position( 0.0, 0.0, -0.06556 );
  }
  atom[1]{
    type = \"H_TIP4P\";
    position( 0.0, 0.75695, 0.52032 );
  }
  atom[2]{
    type = \"H_TIP4P\";
    position( 0.0, -0.75695, 0.52032 );
  }
  atom[3]{
    type = \"EP_TIP4P\";
    position( 0.0, 0.0, 0.08444 );
  }

  rigidBody[0]{
    members(0, 1, 2, 3);
  }

  cutoffGroup{
    members(0, 1, 2, 3);
  }
}"
}

sub printTIP4PEw {
  print $waterFileHandle "\n\nmolecule{
  name = \"TIP4P-Ew\";
  
  atom[0]{
    type = \"O_TIP4P-Ew\";
    position( 0.0, 0.0, -0.06556 );
  }
  atom[1]{
    type = \"H_TIP4P-Ew\";
    position( 0.0, 0.75695, 0.52032 );
  }
  atom[2]{
    type = \"H_TIP4P-Ew\";
    position( 0.0, -0.75695, 0.52032 );
  }
  atom[3]{
    type = \"EP_TIP4P-Ew\";
    position( 0.0, 0.0, 0.05944 );
  }
  
  rigidBody[0]{
    members(0, 1, 2, 3);
  }

  cutoffGroup{
    members(0, 1, 2, 3);
  }
}"
}

sub printTIP5P {
  print $waterFileHandle "\n\nmolecule{
  name = \"TIP5P\";
  
  atom[0]{
    type = \"O_TIP5P\";
    position( 0.0, 0.0, -0.06556 );
  }
  atom[1]{
    type = \"H_TIP5P\";
    position( 0.0, 0.75695, 0.52032 );
  }
  atom[2]{
    type = \"H_TIP5P\";
    position( 0.0, -0.75695, 0.52032 );
  }
  atom[3]{
    type = \"EP_TIP5P\";
    position( 0.57154, 0.0, -0.46971 );
  }
  atom[4]{
    type = \"EP_TIP5P\";
    position( -0.57154, 0.0, -0.46971 );
  }
  
  rigidBody[0]{
    members(0, 1, 2, 3, 4);
  }

  cutoffGroup{
    members(0, 1, 2, 3, 4);
  }
}"
}

sub printTIP5PE {
  print $waterFileHandle "\n\nmolecule{
  name = \"TIP5P-E\";
  
  atom[0]{
    type = \"O_TIP5P-E\";
    position( 0.0, 0.0, -0.06556 );
  }
  atom[1]{
    type = \"H_TIP5P\";
    position( 0.0, 0.75695, 0.52032 );
  }
  atom[2]{
    type = \"H_TIP5P\";
    position( 0.0, -0.75695, 0.52032 );
  }
  atom[3]{
    type = \"EP_TIP5P\";
    position( 0.57154, 0.0, -0.46971 );
  }
  atom[4]{
    type = \"EP_TIP5P\";
    position( -0.57154, 0.0, -0.46971 );
  }
  
  rigidBody[0]{
    members(0, 1, 2, 3, 4);
  }

  cutoffGroup{
    members(0, 1, 2, 3, 4);
  }
}"
}

sub printSPCE {
print $waterFileHandle "\n\nmolecule{
  name = \"SPCE\";
  
  atom[0]{
    type = \"O_SPCE\";
    position( 0.0, 0.0, -0.06461 );
  }
  atom[1]{
    type = \"H_SPCE\";
    position( 0.0, 0.81649, 0.51275 );
  }
  atom[2]{
    type = \"H_SPCE\";
    position( 0.0, -0.81649, 0.51275 );
  }
  
  rigidBody[0]{
    members(0, 1, 2);
  }

  cutoffGroup{
    members(0, 1, 2);
  }
}"
}

sub printSPC {
print $waterFileHandle "\n\nmolecule{
  name = \"SPC\";
  
  atom[0]{
    type = \"O_SPC\";
    position( 0.0, 0.0, -0.06461 );
  }
  atom[1]{
    type = \"H_SPC\";
    position( 0.0, 0.81649, 0.51275 );
  }
  atom[2]{
    type = \"H_SPC\";
    position( 0.0, -0.81649, 0.51275 );
  }
  
  rigidBody[0]{
    members(0, 1, 2);
  }

  cutoffGroup{
    members(0, 1, 2);
  }
}"
}

sub printDPD {
  print $waterFileHandle "\n\nmolecule{
  name = \"DPD\";
  
  atom[0]{
    type = \"DPD\";
    position(0.0, 0.0, 0.0);
  }
}"
}


sub printFakeWater {
  print $waterFileHandle "\n\nmolecule{
  name = \"$waterName\";
  
  atom[0]{
    type = \"$waterName\";
    position(0.0, 0.0, 0.0);
  }
}"
}


sub validateWater {
  if    ($waterName eq 'Cl-')      { $waterCase = 0;    }
  elsif ($waterName eq 'Na+')      { $waterCase = 1;    }
  elsif ($waterName eq 'SSD_E')    { $waterCase = 2;    }
  elsif ($waterName eq 'SSD_RF')   { $waterCase = 3;    }
  elsif ($waterName eq 'SSD')      { $waterCase = 4;    }
  elsif ($waterName eq 'SSD1')     { $waterCase = 5;    }
  elsif ($waterName eq 'TIP3P')    { $waterCase = 6;    }
  elsif ($waterName eq 'TIP4P')    { $waterCase = 7;    }
  elsif ($waterName eq 'TIP4P-Ew') { $waterCase = 8;    }
  elsif ($waterName eq 'TIP5P')    { $waterCase = 9;    }
  elsif ($waterName eq 'TIP5P-E')  { $waterCase = 10;   }
  elsif ($waterName eq 'SPCE')     { $waterCase = 11;   }
  elsif ($waterName eq 'SPC')      { $waterCase = 12;   }
  elsif ($waterName eq 'DPD')      { $waterCase = 13;   }
  else                             { $invalidWater = 1; }
}

sub printWaterModel {
  if    ($waterCase == 0)  { printCl();      }
  elsif ($waterCase == 1)  { printNa();      }
  elsif ($waterCase == 2)  { printSSD_E();   }
  elsif ($waterCase == 3)  { printSSD_RF();  }
  elsif ($waterCase == 4)  { printSSD();     }
  elsif ($waterCase == 5)  { printSSD1();    }
  elsif ($waterCase == 6)  { printTIP3P();   }
  elsif ($waterCase == 7)  { printTIP4P();   }
  elsif ($waterCase == 8)  { printTIP4PEw(); }
  elsif ($waterCase == 9)  { printTIP5P();   }
  elsif ($waterCase == 10) { printTIP5PE();  }
  elsif ($waterCase == 11) { printSPCE();    }
  elsif ($waterCase == 12) { printSPC();     }
  elsif ($waterCase == 13) { printDPD();     }
}
Revision:	1063
Committed:	Tue Oct 10 14:08:10 2006 UTC (19 years ago) by chrisfen
File size:	15684 byte(s)
Log Message:	who needs .in files for perl? using env
#	User	Rev	Content
1	chrisfen	1063	#!/usr/bin/env perl
2
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 $"
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			$invalidWater = 0;
26			$waterCase = -1;
27			$nothingSelected = 1;
28
29			# get our options
30			getopts('hmrvd:l:n:o:w:');
31
32			# just to test opt_h
33			$opt_h = "true" if $#ARGV >= 0;
34
35			# our option output
36			if ($opt_h){
37			print "waterBoxer: builds water boxes\n\n";
38			print "usage: waterBoxer [-hmrv] [-d density] [-l lattice] [-n # waters]\n";
39			print "\t[-o file name] [-w water name] \n\n";
40			print " -h : show this message\n";
41			print " -m : print out a water.md file (file with all water models)\n";
42			print " -r : randomize orientations\n";
43			print " -v : verbose output\n\n";
44			print " -d real : density in g/cm^3\n";
45			print " (default: 1)\n";
46			print " -l integer : 0 - face centered cubic, 1 - simple cubic\n";
47			print " (default: 0)\n";
48			print " -n integer : # of water molecules\n";
49			print " (default: 500)\n";
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";
54			print "Example:\n";
55			die " waterBoxer -d 0.997 -n 864 -w SSD_RF -o ssdrfWater.md\n";
56			}
57
58			# set some variables to be used in the code
59			if (defined($opt_o)){
60			$fileName = $opt_o;
61			$nothingSelected = 0;
62			} else {
63			$fileName = 'freshWater.md';
64			}
65			if ($opt_m){
66			die "Error: $fileName cannot be \"water.md\"\n Please choose a different name\n" if $fileName eq 'water.md';
67			$waterFileHandle = 'WATERMD';
68			$nothingSelected = 0;
69			} else {
70			$waterFileHandle = 'OUTFILE';
71			}
72			if ($opt_r){
73			$doRandomize = $opt_r;
74			$nothingSelected = 0;
75			}
76
77			if (defined($opt_w)){
78			$waterName = $opt_w;
79			$nothingSelected = 0;
80			} else {
81			$waterName = 'SPCE';
82			}
83			validateWater();
84			if ($invalidWater == 1){
85			print "Warning: \'$waterName\' is not a recognized water model name.\n";
86			print " Use the \'-m\' option to generate a \'water.md\' with the\n";
87			print " recognized water model geometries.\n\n";
88			}
89			if (defined($opt_d)){
90			$nothingSelected = 0;
91			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			}
96			}
97			if (defined($opt_l)){
98			$nothingSelected = 0;
99			if ($opt_l =~ /^[0-9]/) {
100			$lattice = $opt_l;
101			if ($lattice != 0 && $lattice != 1){
102			die "Error: the '-l' value ($opt_l) is not a valid number\n Please choose 0 or 1\n";
103			}
104			} else {
105			die "Error: the '-l' value ($opt_l) is not a valid number\n Please choose 0 or 1\n";
106			}
107			}
108			if (defined($opt_n)){
109			$nothingSelected = 0;
110			if ($opt_n =~ /^[0-9]/) {
111			$nMol = $opt_n;
112			} else {
113			die "Error: the '-n' value ($opt_n) is not a valid number\n Please choose a non-negative integer\n";
114			}
115			}
116
117			# open the file writer
118			open(OUTFILE, ">./$fileName") \|\| die "Error: can't open file $fileName\n";
119
120			# check to set magic lattice numbers
121			if ($lattice == 0){
122			$crystalNumReal = ($nMol/4.0)**(1.0/3.0);
123			$crystalNum = int($crystalNumReal + $tolerance);
124			$remainder = $crystalNumReal - $crystalNum;
125
126			# if crystalNumReal wasn't an integer, we bump the crystal to the next
127			# magic number
128			if ($remainder > $tolerance){
129			$crystalNum = $crystalNum + 1;
130			$newMol = 4 * $crystalNum**3;
131			print "Warning: The number chosen ($nMol) failed to build a clean fcc lattice.\n";
132			print " The number of molecules has been increased to the next magic number ($newMol).\n\n";
133			$nMol = $newMol;
134			}
135			} elsif ($lattice == 1){
136			$crystalNumReal = ($nMol/1.0)**(1.0/3.0);
137			$crystalNum = int($crystalNumReal);
138			$remainder = $crystalNumReal - $crystalNum;
139
140			# again, if crystalNumReal wasn't an integer, we bump the crystal to the next
141			# magic number
142			if ($remainder > $tolerance){
143			$crystalNum = $crystalNum + 1;
144			$newMol = $crystalNum**3;
145			print "Warning: The number chosen ($nMol) failed to build a clean simple cubic lattice.\n";
146			print " The number of molecules has been increased to the next magic number ($newMol).\n\n";
147			$nMol = $newMol;
148			}
149			}
150
151			# now we can start building the crystals
152			$boxLength = ($nMol$densityConvert/$density)*(1.0/3.0);
153			$cellLength = $boxLength / $crystalNum;
154			$cell2 = $cellLength*0.5;
155
156			if ($lattice == 0) {
157			# build the unit cell
158			# molecule 0
159			$xCorr[0] = 0.0;
160			$yCorr[0] = 0.0;
161			$zCorr[0] = 0.0;
162			# molecule 1
163			$xCorr[1] = 0.0;
164			$yCorr[1] = $cell2;
165			$zCorr[1] = $cell2;
166			# molecule 2
167			$xCorr[2] = $cell2;
168			$yCorr[2] = $cell2;
169			$zCorr[2] = 0.0;
170			# molecule 3
171			$xCorr[3] = $cell2;
172			$yCorr[3] = 0.0;
173			$zCorr[3] = $cell2;
174			# assemble the lattice
175			$counter = 0;
176			for ($z = 0; $z < $crystalNum; $z++) {
177			for ($y = 0; $y < $crystalNum; $y++) {
178			for ($x = 0; $x < $crystalNum; $x++) {
179			for ($uc = 0; $uc < 4; $uc++) {
180			$xCorr[$uc+$counter] = $xCorr[$uc] + $cellLength*$x;
181			$yCorr[$uc+$counter] = $yCorr[$uc] + $cellLength*$y;
182			$zCorr[$uc+$counter] = $zCorr[$uc] + $cellLength*$z;
183			}
184			$counter = $counter + 4;
185			}
186			}
187			}
188
189			} elsif ($lattice == 1) {
190			# build the unit cell
191			# molecule 0
192			$xCorr[0] = $cell2;
193			$yCorr[0] = $cell2;
194			$zCorr[0] = $cell2;
195			#assemble the lattice
196			$counter = 0;
197			for ($z = 0; $z < $crystalNum; $z++) {
198			for ($y = 0; $y < $crystalNum; $y++) {
199			for ($x = 0; $x < $crystalNum; $x++) {
200			$xCorr[$counter] = $xCorr[0] + $cellLength*$x;
201			$yCorr[$counter] = $yCorr[0] + $cellLength*$y;
202			$zCorr[$counter] = $zCorr[0] + $cellLength*$z;
203
204			$counter++;
205			}
206			}
207			}
208			}
209
210			writeOutFile();
211			print "The water box \"$fileName\" was generated.\n";
212
213			if ($opt_m){
214			printWaterMD();
215			print "The file \"water.md\" was generated for inclusion in \"$fileName\"\n";
216			}
217
218			if ($nothingSelected == 1) {
219			print "(For help, use the \'-h\' option.)\n";
220			}
221
222
223			# this marks the end of the main program, below is subroutines
224
225			sub acos {
226			my ($rad) = @_;
227			my $ret = atan2(sqrt(1 - $rad*$rad), $rad);
228			return $ret;
229			}
230
231			sub writeOutFile {
232			# write out the header
233			print OUTFILE "<OOPSE version=4>\n";
234			findCutoff();
235			findAlpha();
236			printMetaData();
237			printFrameData();
238			print OUTFILE " <StuntDoubles>\n";
239
240			# shift the box center to the origin and write out the coordinates
241			for ($i = 0; $i < $nMol; $i++) {
242			$xCorr[$i] -= 0.5*$boxLength;
243			$yCorr[$i] -= 0.5*$boxLength;
244			$zCorr[$i] -= 0.5*$boxLength;
245
246			$q0 = 1.0;
247			$q1 = 0.0;
248			$q2 = 0.0;
249			$q3 = 0.0;
250
251			if ($doRandomize == 1){
252			$cosTheta = 2.0*rand() - 1.0;
253			$theta = acos($cosTheta);
254			$phi = 2.03.14159265359rand();
255			$psi = 2.03.14159265359rand();
256
257			$q0 = cos(0.5$theta)cos(0.5*($phi + $psi));
258			$q1 = sin(0.5$theta)cos(0.5*($phi - $psi));
259			$q2 = sin(0.5$theta)sin(0.5*($phi - $psi));
260			$q3 = cos(0.5$theta)sin(0.5*($phi + $psi));
261			}
262
263			print OUTFILE "$i\tpq\t$xCorr[$i] $yCorr[$i] $zCorr[$i] ";
264			print OUTFILE "$q0 $q1 $q2 $q3\n";
265			}
266
267			print OUTFILE " </StuntDoubles>\n </Snapshot>\n</OOPSE>\n";
268			}
269
270			sub printMetaData {
271			print OUTFILE " <MetaData>\n";
272
273			# print the water model or includes
274			if ($opt_m){
275			print OUTFILE "#include \"water.md\"";
276			} else {
277			printWaterModel();
278			}
279			printFakeWater() if $invalidWater == 1;
280
281			# now back to the metaData output
282			print OUTFILE "\n\ncomponent{
283			type = \"$waterName\";
284			nMol = $nMol;
285			}
286
287			ensemble = NVE;
288			forceField = \"DUFF\";
289			electrostaticSummationMethod = \"shifted_force\";
290			electrostaticScreeningMethod = \"damped\";
291			cutoffRadius = $cutoff;
292
293			targetTemp = 300;
294			targetPressure = 1.0;
295
296			tauThermostat = 1e3;
297			tauBarostat = 1e4;
298
299			dt = 2.0;
300			runTime = 1e3;
301
302			tempSet = \"true\";
303			thermalTime = 10;
304			sampleTime = 100;
305			statusTime = 2;
306			</MetaData>\n";
307			}
308
309			sub findCutoff {
310			$boxLength2 = 0.5*$boxLength;
311			if ($boxLength2 > $cutoff){
312			# the default is good
313			} else {
314			$cutoff = int($boxLength2);
315			}
316			}
317
318			sub findAlpha {
319			$alpha = $alphaInt - $cutoff*$alphaSlope;
320			}
321
322			sub printFrameData {
323			print OUTFILE
324			" <Snapshot>
325			<FrameData>
326			Time: 0
327			Hmat: {{ $boxLength, 0, 0 }, { 0, $boxLength, 0 }, { 0, 0, $boxLength }}
328			</FrameData>\n";
329			}
330
331			sub printWaterMD {
332			open(WATERMD, ">./water.md") \|\| die "Error: can't open file water.md\n";
333			$waterFileHandle = 'WATERMD';
334
335			print WATERMD "#ifndef _WATER_MD_\n#define _WATER_MD_\n";
336			printCl();
337			printNa();
338			printSSD_E();
339			printSSD_RF();
340			printSSD();
341			printSSD1();
342			printTRED();
343			printTIP3P();
344			printTIP4P();
345			printTIP4PEw();
346			printTIP5P();
347			printTIP5PE();
348			printSPCE();
349			printSPC();
350			printDPD();
351			print WATERMD "\n\n#endif";
352			}
353
354			sub printCl {
355			print $waterFileHandle "\n\nmolecule{
356			name = \"Cl-\";
357
358			atom[0]{
359			type = \"Cl-\";
360			position(0.0, 0.0, 0.0);
361			}
362			}"
363			}
364
365			sub printNa {
366			print $waterFileHandle "\n\nmolecule{
367			name = \"Na+\";
368
369			atom[0]{
370			type = \"Na+\";
371			position(0.0, 0.0, 0.0);
372			}
373			}"
374			}
375
376			sub printSSD_E {
377			print $waterFileHandle "\n\nmolecule{
378			name = \"SSD_E\";
379
380			atom[0]{
381			type = \"SSD_E\";
382			position( 0.0, 0.0, 0.0 );
383			orientation( 0.0, 0.0, 0.0 );
384			}
385			}"
386			}
387
388			sub printSSD_RF {
389			print $waterFileHandle "\n\nmolecule{
390			name = \"SSD_RF\";
391
392			atom[0]{
393			type = \"SSD_RF\";
394			position( 0.0, 0.0, 0.0 );
395			orientation( 0.0, 0.0, 0.0 );
396			}
397			}"
398			}
399
400			sub printSSD {
401			print $waterFileHandle "\n\nmolecule{
402			name = \"SSD\";
403
404			atom[0]{
405			type = \"SSD\";
406			position( 0.0, 0.0, 0.0 );
407			orientation( 0.0, 0.0, 0.0 );
408			}
409			}"
410			}
411
412			sub printSSD1 {
413			print $waterFileHandle "\n\nmolecule{
414			name = \"SSD1\";
415
416			atom[0]{
417			type = \"SSD1\";
418			position( 0.0, 0.0, 0.0 );
419			orientation( 0.0, 0.0, 0.0 );
420			}
421			}"
422			}
423
424			sub printTRED {
425			print $waterFileHandle "\n\nmolecule{
426			name = \"TRED\";
427
428			atom[0]{
429			type = \"TRED\";
430			position( 0.0, 0.0, 0.0 );
431			orientation( 0.0, 0.0, 0.0 );
432			}
433			atom[1]{
434			type = \"EP_TRED\";
435			position( 0.0, 0.0, 0.5 );
436			}
437
438			rigidBody[0]{
439			members(0, 1);
440			}
441
442			cutoffGroup{
443			members(0, 1);
444			}
445			}"
446			}
447
448			sub printTIP3P {
449			print $waterFileHandle "\n\nmolecule{
450			name = \"TIP3P\";
451
452			atom[0]{
453			type = \"O_TIP3P\";
454			position( 0.0, 0.0, -0.06556 );
455			}
456			atom[1]{
457			type = \"H_TIP3P\";
458			position( 0.0, 0.75695, 0.52032 );
459			}
460			atom[2]{
461			type = \"H_TIP3P\";
462			position( 0.0, -0.75695, 0.52032 );
463			}
464
465			rigidBody[0]{
466			members(0, 1, 2);
467			}
468
469			cutoffGroup{
470			members(0, 1, 2);
471			}
472			}"
473			}
474
475			sub printTIP4P {
476			print $waterFileHandle "\n\nmolecule{
477			name = \"TIP4P\";
478
479			atom[0]{
480			type = \"O_TIP4P\";
481			position( 0.0, 0.0, -0.06556 );
482			}
483			atom[1]{
484			type = \"H_TIP4P\";
485			position( 0.0, 0.75695, 0.52032 );
486			}
487			atom[2]{
488			type = \"H_TIP4P\";
489			position( 0.0, -0.75695, 0.52032 );
490			}
491			atom[3]{
492			type = \"EP_TIP4P\";
493			position( 0.0, 0.0, 0.08444 );
494			}
495
496			rigidBody[0]{
497			members(0, 1, 2, 3);
498			}
499
500			cutoffGroup{
501			members(0, 1, 2, 3);
502			}
503			}"
504			}
505
506			sub printTIP4PEw {
507			print $waterFileHandle "\n\nmolecule{
508			name = \"TIP4P-Ew\";
509
510			atom[0]{
511			type = \"O_TIP4P-Ew\";
512			position( 0.0, 0.0, -0.06556 );
513			}
514			atom[1]{
515			type = \"H_TIP4P-Ew\";
516			position( 0.0, 0.75695, 0.52032 );
517			}
518			atom[2]{
519			type = \"H_TIP4P-Ew\";
520			position( 0.0, -0.75695, 0.52032 );
521			}
522			atom[3]{
523			type = \"EP_TIP4P-Ew\";
524			position( 0.0, 0.0, 0.05944 );
525			}
526
527			rigidBody[0]{
528			members(0, 1, 2, 3);
529			}
530
531			cutoffGroup{
532			members(0, 1, 2, 3);
533			}
534			}"
535			}
536
537			sub printTIP5P {
538			print $waterFileHandle "\n\nmolecule{
539			name = \"TIP5P\";
540
541			atom[0]{
542			type = \"O_TIP5P\";
543			position( 0.0, 0.0, -0.06556 );
544			}
545			atom[1]{
546			type = \"H_TIP5P\";
547			position( 0.0, 0.75695, 0.52032 );
548			}
549			atom[2]{
550			type = \"H_TIP5P\";
551			position( 0.0, -0.75695, 0.52032 );
552			}
553			atom[3]{
554			type = \"EP_TIP5P\";
555			position( 0.57154, 0.0, -0.46971 );
556			}
557			atom[4]{
558			type = \"EP_TIP5P\";
559			position( -0.57154, 0.0, -0.46971 );
560			}
561
562			rigidBody[0]{
563			members(0, 1, 2, 3, 4);
564			}
565
566			cutoffGroup{
567			members(0, 1, 2, 3, 4);
568			}
569			}"
570			}
571
572			sub printTIP5PE {
573			print $waterFileHandle "\n\nmolecule{
574			name = \"TIP5P-E\";
575
576			atom[0]{
577			type = \"O_TIP5P-E\";
578			position( 0.0, 0.0, -0.06556 );
579			}
580			atom[1]{
581			type = \"H_TIP5P\";
582			position( 0.0, 0.75695, 0.52032 );
583			}
584			atom[2]{
585			type = \"H_TIP5P\";
586			position( 0.0, -0.75695, 0.52032 );
587			}
588			atom[3]{
589			type = \"EP_TIP5P\";
590			position( 0.57154, 0.0, -0.46971 );
591			}
592			atom[4]{
593			type = \"EP_TIP5P\";
594			position( -0.57154, 0.0, -0.46971 );
595			}
596
597			rigidBody[0]{
598			members(0, 1, 2, 3, 4);
599			}
600
601			cutoffGroup{
602			members(0, 1, 2, 3, 4);
603			}
604			}"
605			}
606
607			sub printSPCE {
608			print $waterFileHandle "\n\nmolecule{
609			name = \"SPCE\";
610
611			atom[0]{
612			type = \"O_SPCE\";
613			position( 0.0, 0.0, -0.06461 );
614			}
615			atom[1]{
616			type = \"H_SPCE\";
617			position( 0.0, 0.81649, 0.51275 );
618			}
619			atom[2]{
620			type = \"H_SPCE\";
621			position( 0.0, -0.81649, 0.51275 );
622			}
623
624			rigidBody[0]{
625			members(0, 1, 2);
626			}
627
628			cutoffGroup{
629			members(0, 1, 2);
630			}
631			}"
632			}
633
634			sub printSPC {
635			print $waterFileHandle "\n\nmolecule{
636			name = \"SPC\";
637
638			atom[0]{
639			type = \"O_SPC\";
640			position( 0.0, 0.0, -0.06461 );
641			}
642			atom[1]{
643			type = \"H_SPC\";
644			position( 0.0, 0.81649, 0.51275 );
645			}
646			atom[2]{
647			type = \"H_SPC\";
648			position( 0.0, -0.81649, 0.51275 );
649			}
650
651			rigidBody[0]{
652			members(0, 1, 2);
653			}
654
655			cutoffGroup{
656			members(0, 1, 2);
657			}
658			}"
659			}
660
661			sub printDPD {
662			print $waterFileHandle "\n\nmolecule{
663			name = \"DPD\";
664
665			atom[0]{
666			type = \"DPD\";
667			position(0.0, 0.0, 0.0);
668			}
669			}"
670			}
671
672
673			sub printFakeWater {
674			print $waterFileHandle "\n\nmolecule{
675			name = \"$waterName\";
676
677			atom[0]{
678			type = \"$waterName\";
679			position(0.0, 0.0, 0.0);
680			}
681			}"
682			}
683
684
685			sub validateWater {
686			if ($waterName eq 'Cl-') { $waterCase = 0; }
687			elsif ($waterName eq 'Na+') { $waterCase = 1; }
688			elsif ($waterName eq 'SSD_E') { $waterCase = 2; }
689			elsif ($waterName eq 'SSD_RF') { $waterCase = 3; }
690			elsif ($waterName eq 'SSD') { $waterCase = 4; }
691			elsif ($waterName eq 'SSD1') { $waterCase = 5; }
692			elsif ($waterName eq 'TIP3P') { $waterCase = 6; }
693			elsif ($waterName eq 'TIP4P') { $waterCase = 7; }
694			elsif ($waterName eq 'TIP4P-Ew') { $waterCase = 8; }
695			elsif ($waterName eq 'TIP5P') { $waterCase = 9; }
696			elsif ($waterName eq 'TIP5P-E') { $waterCase = 10; }
697			elsif ($waterName eq 'SPCE') { $waterCase = 11; }
698			elsif ($waterName eq 'SPC') { $waterCase = 12; }
699			elsif ($waterName eq 'DPD') { $waterCase = 13; }
700			else { $invalidWater = 1; }
701			}
702
703			sub printWaterModel {
704			if ($waterCase == 0) { printCl(); }
705			elsif ($waterCase == 1) { printNa(); }
706			elsif ($waterCase == 2) { printSSD_E(); }
707			elsif ($waterCase == 3) { printSSD_RF(); }
708			elsif ($waterCase == 4) { printSSD(); }
709			elsif ($waterCase == 5) { printSSD1(); }
710			elsif ($waterCase == 6) { printTIP3P(); }
711			elsif ($waterCase == 7) { printTIP4P(); }
712			elsif ($waterCase == 8) { printTIP4PEw(); }
713			elsif ($waterCase == 9) { printTIP5P(); }
714			elsif ($waterCase == 10) { printTIP5PE(); }
715			elsif ($waterCase == 11) { printSPCE(); }
716			elsif ($waterCase == 12) { printSPC(); }
717			elsif ($waterCase == 13) { printDPD(); }
718			}