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

#!@PERL_EXECUTABLE@

# program that builds water boxes

# author    = "Chris Fennell
# version   = "$Revision$"
# date      = "$Date$"
# copyright = "Copyright (c) 2006 by the University of Notre Dame"
# license   = "OpenMD"

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 = 9;
$alpha = 0.18;
$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 "<OpenMD version=1>\n";
  findCutoff();
  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</OpenMD>\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 = \"Amber\";
electrostaticSummationMethod = \"shifted_force\";
electrostaticScreeningMethod = \"damped\";
cutoffRadius = $cutoff;
switchingRadius = $cutoff;
dampingAlpha = $alpha;

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