applications/utilities/waterBoxer.in

#!@PERLINTERP@ -w

# program that builds water boxes

# author    = "Chris Fennell
# version   = "$Revision: 1.2 $"
# date      = "$Date: 2006-09-01 22:58:33 $"
# 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;

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

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

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

# set some variables to be used in the code
$fileName = $ARGV[0];
if (defined($fileName)){
} else {
    $fileName = 'waterBox';
}
if ($opt_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 (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();

if ($opt_m || $invalidWater){
  printWaterMD();
}

# 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\";
dampingAlpha = $alpha;
cutoffRadius = $cutoff;

targetTemp = 300;
targetPressure = 1.0;

tauThermostat = 1e3;
tauBarostat = 1e4;

dt = 2.0;
runTime = 1e3;

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

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

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

sub printFrameData{
print OUTFILE 
"  <Snapshot>
    <FrameData>
        Time: 0
        Hmat: {{ $boxLength, 0, 0 }, { 0, $boxLength, 0 }, { 0, 0, $boxLength }}
    </FrameData>\n";
}

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

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