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

#!/usr/bin/env perl

# program that builds water boxes

# author    = "Chris Fennell
# version   = "$Revision: 1.6 $"
# date      = "$Date: 2009-10-22 19:12:14 $"
# 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:x:y:z:');

# 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";
    print "  -x real    : dimension of the box along the x direction\n";
    print "  -y real    : dimension of the box along the y direction\n";
    print "  -z real    : dimension of the box along the z direction\n\n";
    print "Note: you can only use values of x, y, or z that are smaller\n";
    print "      than the derived box length for a given density and\n";
    print "      number of molecules.\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_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_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 ($waterName eq 'DPD') {
  # DPD waters are stand-ins for 4 water molecules
  $density = $density * 0.25;
}
if ($waterName eq 'CG2') {
  # CG2 waters are stand-ins for 2 water molecules
  $density = $density * 0.5;
}

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";
  }
}
if (defined($opt_x)){
  $nothingSelected = 0;
  if ($opt_x =~ /^[0-9]/) {
    $boxx = $opt_x;
  } else {
    die "Error: the value for '-x' ($opt_x) is not a valid number\n       Please choose a positive real # value\n";
  }
}
if (defined($opt_y)){
  $nothingSelected = 0;
  if ($opt_y =~ /^[0-9]/) {
    $boxy = $opt_y;
  } else {
    die "Error: the value for '-y' ($opt_y) is not a valid number\n       Please choose a positive real # value\n";
  }
}
if (defined($opt_z)){
  $nothingSelected = 0;
  if ($opt_z =~ /^[0-9]/) {
    $boxz = $opt_z;
  } else {
    die "Error: the value for '-z' ($opt_z) is not a valid number\n       Please choose a positive real # value\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 + $tolerance);
  $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;

if ($boxx != 0) {
  if ($boxLength < $boxx) {
    print "Computed box length is smaller than requested x axis.  Use more\n";
    die "molecules.";    
  }
} else {
  $boxx = $boxLength;
}
if ($boxy != 0) {
  if ($boxLength < $boxy) {
    print "Computed box length is smaller than requested y axis.  Use more\n";
    die "molecules.";    
  }
} else {
  $boxy = $boxLength;
}
if ($boxz != 0) {
  if ($boxLength < $boxz) {
    print "Computed box length is smaller than requested z axis.  Use more\n";
    die "molecules.";    
  }
} else {
  $boxz = $boxLength;
}

$nx = int($boxx / $cellLength);
$ny = int($boxy / $cellLength);
$nz = int($boxz / $cellLength);

if ($lattice == 0) {
  $nMol = 4 * $nx * $ny * $nz;
} else {
  $nMol = $nx * $ny * $nz;
}

$newDensity = $nMol * $densityConvert / ($boxx*$boxy*$boxz);

if (abs($newDensity-$density) > $tolerance) {
  print "Resetting density to $newDensity to make chosen box sides work out\n";
}
$cellLengthX = $boxx/$nx;
$cellLengthY = $boxy/$ny;
$cellLengthZ = $boxz/$nz;

$cell2X = $cellLengthX*0.5;
$cell2Y = $cellLengthY*0.5;
$cell2Z = $cellLengthZ*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] = $cell2Y;
  $zCorr[1] = $cell2Z;
# molecule 2
  $xCorr[2] = $cell2X;
  $yCorr[2] = $cell2Y;
  $zCorr[2] = 0.0;
# molecule 3
  $xCorr[3] = $cell2X;
  $yCorr[3] = 0.0;
  $zCorr[3] = $cell2Z;
# assemble the lattice
  $counter = 0;
  for ($z = 0; $z < $nz; $z++) {
    for ($y = 0; $y < $ny; $y++) {
      for ($x = 0; $x < $nx; $x++) {
        for ($uc = 0; $uc < 4; $uc++) {
          $xCorr[$uc+$counter] = $xCorr[$uc] + $cellLengthX*$x;
          $yCorr[$uc+$counter] = $yCorr[$uc] + $cellLengthY*$y;
          $zCorr[$uc+$counter] = $zCorr[$uc] + $cellLengthZ*$z;
        }
        $counter = $counter + 4;
      }
    }
  }
  
} elsif ($lattice == 1) {
  # build the unit cell
  # molecule 0
  $xCorr[0] = $cell2X;
  $yCorr[0] = $cell2Y;
  $zCorr[0] = $cell2Z;
#assemble the lattice
  $counter = 0;
  for ($z = 0; $z < $nz; $z++) {
    for ($y = 0; $y < $ny; $y++) {
      for ($x = 0; $x < $nx; $x++) {
        $xCorr[$counter] = $xCorr[0] + $cellLengthX*$x;
        $yCorr[$counter] = $yCorr[0] + $cellLengthY*$y;
        $zCorr[$counter] = $zCorr[0] + $cellLengthZ*$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*$boxx;
    $yCorr[$i] -= 0.5*$boxy;
    $zCorr[$i] -= 0.5*$boxz;
    
    $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 {
  if ($boxy < $boxx) {
    $bm = $boxy;
  } else {
    $bm = $boxx;
  } 
  if ($boxz < $bm) {
    $bm = $boxz;
  }
  $boxLength2 = 0.5*$bm;
  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: {{ $boxx, 0, 0 }, { 0, $boxy, 0 }, { 0, 0, $boxz }}
    </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();
  printCG2();
  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);
  }

}"
}

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);
  }

}"
}

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);
  }

}"
}

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);
  }

}"
}

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);
  }

}"
}

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);
  }

}"
}

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);
  }

}"
}

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);
  }

}"
}

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

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