applications/utilities/waterBoxer.in

#!@PERLINTERP@ -w

# program that builds water boxes

# author    = "Chris Fennell
# version   = "$Revision: 1.4 $"
# date      = "$Date: 2006-10-02 23:54:30 $"
# 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;

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

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

# 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 "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;
} else {
    $fileName = 'freshWater.md';
}
if ($opt_m){
  die "$fileName cannot be \"water.md\"\n\tPlease choose a different name\n" if $fileName eq 'water.md';
  $waterFileHandle = 'WATERMD';
} else {
  $waterFileHandle = 'OUTFILE';
}
if ($opt_r){
  $doRandomize = $opt_r;
}

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

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

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

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

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

writeOutFile();
print "Water box \"$fileName\" generated.\n";

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