waterBoxer

#!@Python3_EXECUTABLE@
 
__author__ = "Chris Fennell and Dan Gezelter"
__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
__license__ = "OpenMD"
 
import sys
import os
import math
import random
import argparse
from argparse import RawDescriptionHelpFormatter
 
def usage():
    print(__doc__)
 
# function to calculate the arccosine
def acos(rad):
    return math.atan2(math.sqrt(1 - rad * rad), rad)
 
# Function to open the output file
def open_output_file(fileName):
    f = open(fileName, 'w')
    if os.path.exists(fileName):
        return f
    else:
        sys.exit("Error: can't open file {}\n".format(fileName))
 
# Function to find the cutoff
def find_cutoff(cutoff, boxx, boxy, boxz):
    bm = min(boxx, boxy, boxz)
    boxLength2 = 0.5 * bm
    if boxLength2 <= cutoff:
        cutoff = int(boxLength2)
 
# Functions to print the water model types
def print_cl(outfile):
    outfile.write("""
molecule{
  name = "Cl-";
  atom[0]{
    type = "Cl-";
    position(0.0, 0.0, 0.0);
  }
}""")
 
def print_na(outfile):
    outfile.write("""
molecule{
  name = "Na+";
  atom[0]{
    type = "Na+";
    position(0.0, 0.0, 0.0);
  }
}""")
 
def print_ssd_e(outfile):
    outfile.write("""
molecule{
  name = "SSD_E";
  atom[0]{
    type = "SSD_E";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}""")
 
def print_ssd_rf(outfile):
    outfile.write("""
molecule{
  name = "SSD_RF";
  atom[0]{
    type = "SSD_RF";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}""")
 
def print_ssd(outfile):
    outfile.write("""
molecule{
  name = "SSD";
  atom[0]{
    type = "SSD";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}""")
 
def print_ssdq(outfile):
    outfile.write("""
molecule{
  name = "SSDQ";
  atom[0]{
    type = "SSDQ";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}""")
 
def print_ssdqo(outfile):
    outfile.write("""
molecule{
  name = "SSDQO";
  atom[0]{
    type = "SSDQO";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}""")
 
def print_ssd1(outfile):
    outfile.write("""
molecule{
  name = "SSD1";
  atom[0]{
    type = "SSD1";
    position( 0.0, 0.0, 0.0 );
    orientation( 0.0, 0.0, 0.0 );
  }
}""")
 
def print_tip3p(outfile):
    outfile.write("""
molecule{
  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);
  }
}""")
 
def print_tip4p(outfile):
    outfile.write("""
molecule{
  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);
  }
}""")
 
def print_tip4p_ice(outfile):
    outfile.write("""
molecule{
  name = "TIP4P-Ice";
  atom[0]{
    type = "O_TIP4P-Ice";
    position( 0.0, 0.0, -0.06556 );
  }
  atom[1]{
    type = "H_TIP4P-Ice";
    position( 0.0, 0.75695, 0.52032 );
  }
  atom[2]{
    type = "H_TIP4P-Ice";
    position( 0.0, -0.75695, 0.52032 );
  }
  atom[3]{
    type = "EP_TIP4P-Ice";
    position( 0.0, 0.0, 0.09214 );
  }
  rigidBody[0]{
    members(0, 1, 2, 3);
  }
}""")
 
def print_tip4p_2005(outfile):
    outfile.write("""
molecule{
  name = "TIP4P-2005";
  atom[0]{
    type = "O_TIP4P-2005";
    position( 0.0, 0.0, -0.06556 );
  }
  atom[1]{
    type = "H_TIP4P-2005";
    position( 0.0, 0.75695, 0.52032 );
  }
  atom[2]{
    type = "H_TIP4P-2005";
    position( 0.0, -0.75695, 0.52032 );
  }
  atom[3]{
    type = "EP_TIP4P-2005";
    position( 0.0, 0.0, 0.08904 );
  }
  rigidBody[0]{
    members(0, 1, 2, 3);
  }
}""")
 
def print_tip4pew(outfile):
    outfile.write("""
molecule{
  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);
  }
}""")
 
def print_tip5p(outfile):
    outfile.write("""
molecule{
  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);
  }
}""")
 
def print_tip5pe(outfile):
    outfile.write("""
molecule{
  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);
  }
}""")
 
def print_spce(outfile):
    outfile.write("""
molecule{
  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);
  }
}""")
 
def print_spc(outfile):
    outfile.write("""
molecule{
  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);
  }
}""")
 
def print_spc_hw(outfile):
    outfile.write("""
molecule{
  name = "SPC-HW";
  atom[0]{
    type = "O_SPC-HW";
    position( 0.0, 0.0, -0.06461 );
  }
  atom[1]{
    type = "D_SPC-HW";
    position( 0.0, 0.81649, 0.51275 );
  }
  atom[2]{
    type = "D_SPC-HW";
    position( 0.0, -0.81649, 0.51275 );
  }
  rigidBody[0]{
    members(0, 1, 2);
  }
}""")
 
def print_ne6(outfile):
    outfile.write("""
molecule{
  name = "NE6";
  atom[0]{
    type = "O_NE6";
    position( 0.0, 0.0, 0.0 );
  }
  atom[1]{
    type = "H_NE6";
    position( 0.0, 0.576029, 0.79283665 );
  }
  atom[2]{
    type = "H_NE6";
    position( 0.0, -0.576029, 0.79283665 );
  }
  atom[3]{
    type = "EP_NE6";
    position( 0.0, 0.23, 0.0 );
  }
  atom[4]{
    type = "LP_NE6";
    position( 0.732813007, -0.50364843, 0.0  );
  }
  atom[5]{
    type = "LP_NE6";
    position( -0.732813007, -0.50364843, 0.0  );
  }
  rigidBody[0]{
    members(0, 1, 2, 3, 4, 5);
  }
}""")
 
def print_tip3p_fb(outfile):
    outfile.write("""
molecule{
  name = "TIP3P-FB";
  atom[0]{
    type = "O_TIP3P-FB";
    position( 0.0, 0.0, -0.066424 );
  }
  atom[1]{
    type = "H_TIP3P-FB";
    position( 0.0, 0.819341, 0.527225 );
  }
  atom[2]{
    type = "H_TIP3P-FB";
    position( 0.0, -0.819341, 0.527225 );
  }
  rigidBody[0]{
    members(0, 1, 2);
  }
}""")
 
def print_tip4p_fb(outfile):
    outfile.write("""
molecule{
  name = "TIP4P-FB";
  atom[0]{
    type = "O_TIP4P-FB";
    position( 0.0, 0.0, -0.0655549 );
  }
  atom[1]{
    type = "H_TIP4P-FB";
    position( 0.0, 0.75695, 0.520327 );
  }
  atom[2]{
    type = "H_TIP4P-FB";
    position( 0.0, -0.75695, 0.520327 );
  }
  atom[3]{
    type = "EP_TIP4P-FB";
    position( 0.0, 0.0, 0.0397151 );
  }
  rigidBody[0]{
    members(0, 1, 2, 3);
  }
}""")
 
def print_opc(outfile):
    outfile.write("""
molecule{
  name = "OPC";
  atom[0]{
    type = "O_OPC";
    position( 0.0, 0.0, -0.0603651 );
  }
  atom[1]{
    type = "H_OPC";
    position( 0.0, 0.685582, 0.479134 );
  }
  atom[2]{
    type = "H_OPC";
    position( 0.0, -0.685582, 0.479134 );
  }
  atom[3]{
    type = "EP_OPC";
    position( 0.0, 0.0, 0.0990349 );
  }
  rigidBody[0]{
    members(0, 1, 2, 3);
  }
}""")
 
def print_opc3(outfile):
    outfile.write("""
molecule{
  name = "OPC3";
  atom[0]{
    type = "O_OPC3";
    position( 0.0, 0.0, -0.0632382 );
  }
  atom[1]{
    type = "H_OPC3";
    position( 0.0, 0.799262, 0.501939 );
  }
  atom[2]{
    type = "H_OPC3";
    position( 0.0, -0.799262, 0.501939 );
  }
  rigidBody[0]{
    members(0, 1, 2);
  }
}""")
    
def print_dpd(outfile):
    outfile.write("""
molecule{
  name = "DPD";
  atom[0]{
    type = "DPD";
    position( 0.0, 0.0, 0.0 );
  }
}""")
 
def print_cg2(outfile):
    outfile.write("""
molecule{
  name = "CG2";
  atom[0]{
    type = "CG2";
    position( 0.0, 0.0, 0.0 );
  }
}""")
 
# Function to print fake water model
def print_fake_water(waterName, outfile):
    outfile.write("\n\nmolecule{{\n  name = \"{}\";\n  atom[0]{{\n    type = \"{}\";\n    position(0.0, 0.0, 0.0);\n  }}\n}}".format(waterName, waterName))
 
# Function to print the water model
def print_water_model(waterCase, outfile):
    if waterCase == 0: print_cl(outfile)
    elif waterCase == 1: print_na(outfile)
    elif waterCase == 2: print_ssd_e(outfile)
    elif waterCase == 3: print_ssd_rf(outfile)
    elif waterCase == 4: print_ssd(outfile)
    elif waterCase == 5: print_ssd1(outfile)
    elif waterCase == 6: print_tip3p(outfile)
    elif waterCase == 7: print_tip4p(outfile)
    elif waterCase == 8: print_tip4pew(outfile)
    elif waterCase == 9: print_tip5p(outfile)
    elif waterCase == 10: print_tip5pe(outfile)
    elif waterCase == 11: print_spce(outfile)
    elif waterCase == 12: print_spc(outfile)
    elif waterCase == 13: print_dpd(outfile)
    elif waterCase == 14: print_cg2(outfile)
    elif waterCase == 15: print_ssdq(outfile)
    elif waterCase == 16: print_ssdqo(outfile)
    elif waterCase == 17: print_tip4p_ice(outfile)
    elif waterCase == 18: print_tip4p_2005(outfile)
    elif waterCase == 19: print_spc_hw(outfile)
    elif waterCase == 20: print_ne6(outfile)
    elif waterCase == 21: print_tip3p_fb(outfile)
    elif waterCase == 22: print_tip4p_fb(outfile)
    elif waterCase == 23: print_opc(outfile)
    elif waterCase == 24: print_opc3(outfile)
 
# Function to validate the water model
def validate_water(waterName):
    water_models = {
        'Cl-': 0, 'Na+': 1, 'SSD_E': 2, 'SSD_RF': 3, 'SSD': 4, 'SSD1': 5,
        'TIP3P': 6, 'TIP4P': 7, 'TIP4P-Ew': 8, 'TIP5P': 9, 'TIP5P-E': 10,
        'SPCE': 11, 'SPC': 12, 'DPD': 13, 'CG2': 14, 'SSDQ': 15, 'SSDQO': 16,
        'TIP4P-Ice': 17, 'TIP4P-2005': 18, 'SPC-HW': 19, 'NE6': 20,
        'TIP3P-FB': 21, 'TIP4P-FB': 22, 'OPC': 23, 'OPC3': 24
    }
    waterCase = water_models.get(waterName, -1)
    invalidWater = False
    if waterCase == -1:
        invalidWater = True
    return (waterCase, invalidWater)
 
# Function to print the meta data
def print_meta_data(outfile, doWaterInc, waterName, waterCase, invalidWater,
                        nMol, cutoff, alpha):
    outfile.write("  <MetaData>\n")
 
    if doWaterInc:
        outfile.write('#include "water.inc"')
    else:
        print_water_model(waterCase, outfile)
    
    if invalidWater:
        print_fake_water(waterName, outfile)
    
    outfile.write("\n\n")
    outfile.write("component{\n")
    outfile.write("  type = \"{}\";\n".format(waterName))
    outfile.write("  nMol = {};\n".format(nMol))
    outfile.write("}\n")
    outfile.write("ensemble = \"NVT\";\n")
    outfile.write("forceField = \"Water\";\n")
    outfile.write("electrostaticSummationMethod = \"shifted_force\";\n")
    outfile.write("electrostaticScreeningMethod = \"damped\";\n")
    outfile.write("cutoffRadius = {};\n".format(cutoff))
    outfile.write("dampingAlpha = {};\n".format(alpha))
    outfile.write("\n")
    outfile.write("targetTemp = 300;\n")
    outfile.write("targetPressure = 1.0;\n")
    outfile.write("\n")
    outfile.write("tauThermostat = 1e3;\n")
    outfile.write("tauBarostat = 1e4;\n")
    outfile.write("\n")
    outfile.write("dt = 2.0;\n")
    outfile.write("runTime = 1e3;\n")
    outfile.write("\n")
    outfile.write("tempSet = \"true\";\n")
    outfile.write("thermalTime = 10;\n")
    outfile.write("sampleTime = 100;\n")
    outfile.write("statusTime = 2;\n")
    outfile.write("  </MetaData>\n")
 
# Function to print the frame data
def print_frame_data(outfile, boxx, boxy, boxz):
    outfile.write("  <Snapshot>\n")
    outfile.write("    <FrameData>\n")
    outfile.write("        Time: 0\n")
    outfile.write("        Hmat: {{{{ {}, 0, 0 }}, {{ 0, {}, 0 }}, {{ 0, 0, {} }}}}\n".format(boxx, boxy, boxz))
    outfile.write("    </FrameData>\n")
        
# Function to write the output file
def write_output_file(outfile, nMol, xCorr, yCorr, zCorr, cutoff, alpha,
                          boxx, boxy, boxz, waterName, waterCase, invalidWater,
                          doRandomize, doWaterInc):
    # Write out the header
    outfile.write("<OpenMD version=2>\n")
    find_cutoff(cutoff, boxx, boxy, boxz)
    print_meta_data(outfile, doWaterInc, waterName, waterCase, invalidWater,
                        nMol, cutoff, alpha)
    print_frame_data(outfile, boxx, boxy, boxz)
    outfile.write("    <StuntDoubles>\n")
 
    # Shift the box center to the origin and write out the coordinates
    for i in range(nMol):
        xCorr[i] -= 0.5 * boxx
        yCorr[i] -= 0.5 * boxy
        zCorr[i] -= 0.5 * boxz
 
        q0 = 1.0
        q1 = q2 = q3 = 0.0
 
        if doRandomize:
            cosTheta = 2.0 * random.random() - 1.0
            theta = acos(cosTheta)
            phi = 2.0 * math.pi * random.random()
            psi = 2.0 * math.pi * random.random()
 
            q0 = math.cos(0.5 * theta) * math.cos(0.5 * (phi + psi))
            q1 = math.sin(0.5 * theta) * math.cos(0.5 * (phi - psi))
            q2 = math.sin(0.5 * theta) * math.sin(0.5 * (phi - psi))
            q3 = math.cos(0.5 * theta) * math.sin(0.5 * (phi + psi))
 
        outfile.write("{}\tpq\t{} {} {} {} {} {} {}\n".format(i, xCorr[i], yCorr[i], zCorr[i], q0, q1, q2, q3))
 
    outfile.write("    </StuntDoubles>\n")
    outfile.write("  </Snapshot>\n")
    outfile.write("</OpenMD>\n")
 
def main(argv):
 
    parser = argparse.ArgumentParser(
        description='builds water boxes',
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=
'''
Note: you can only use values of x, y, or z that are smaller
      than the derived box length for a given density and
      number of molecules.
 
Example:
       waterBoxer -d 0.997 -x 20 -y 20 -z 20 -w TIP4P-Ice -o tip4pWater.omd
''')
    parser.add_argument("-d", "--density=", action="store", dest='density',
                        type=float, help="density in g/cm^3",
                        required=False, default=1.0)
    parser.add_argument("-l", "--lattice=", choices=[0,1],
                        action="store", dest='lattice',
                        type=int, default=0, required=False,
                        help="0 - face centered cubic, 1 - simple cubic")
    parser.add_argument("-c", "--rcut=", action="store", dest="rcut",
                        type=float, default=3.2, required=False,
                        help="default overlap cutoff in angstroms")
    parser.add_argument("-o", "--output=", action="store", dest='fileName', 
                        help="output file name", required=False,
                        default="freshWater.omd")
    parser.add_argument("-w", "--water=", action="store", dest="waterName", 
                        help="name of the water StuntDouble",
                        required=False, default="SPCE")
    parser.add_argument("-x", "--boxx=", action="store", required=True,
                        dest='boxx', type=float,
                        help="dimension of the box along the x direction")
    parser.add_argument("-y", "--boxy=", action="store", required=True,
                        dest='boxy', type=float,
                        help="dimension of the box along the y direction")
    parser.add_argument("-z", "--boxz=", action="store", required=True,
                        dest='boxz', type=float,
                        help="dimension of the box along the z direction")
    parser.add_argument("-r", "--randomize",
                        action=argparse.BooleanOptionalAction, default=False,
                        type=bool, dest='doRandomize',
                        help="randomize orientations")
    parser.add_argument("-m", "--printInclude",
                        action=argparse.BooleanOptionalAction, default=False,
                        type=bool, dest='doWaterInc',
                        help="print out a water.inc file (file with all water models)")
 
    if len(sys.argv) == 1:
        parser.print_help()
        sys.exit(2)
    args = parser.parse_args()
    
    density = args.density
    lattice = args.lattice
    rcut = args.rcut
    fileName = args.fileName
    waterName = args.waterName
    boxx = args.boxx
    boxy = args.boxy
    boxz = args.boxz
    doRandomize = args.doRandomize
    doWaterInc = args.doWaterInc
 
    tolerance = 1.0E-8
    mass = 2.99151E-23  # mass of H2O in grams
    cm3ToAng3 = 1E24   # convert cm^3 to angstroms^3
    densityConvert = mass * cm3ToAng3
    cutoff = 9
    alpha = 0.18
    invalidWater = False
    waterCase = -1
 
    if doWaterInc:
        if fileName == 'water.inc':
            sys.exit("Error: {} cannot be \"water.inc\"\n       Please choose a different name\n".format(fileName))
    if density <= 0:
        usage()
        sys.exit("Error: the value for '-d' ({}) is not a valid number\n       Please choose a positive real # value\n".format(density))
    if rcut <= 0:
        usage()
        sys.exit("Error: the value for '-c' ({}) is not a valid number\n       Please choose a positive real # value\n".format(arg))
    if lattice != 0 and lattice != 1:
        usage()
        sys.exit("Error: the '-l' value ({}) is not a valid number\n       Please choose 0 or 1\n".format(arg))
    if boxx <= 0:
        usage()
        sys.exit("Error: the value for '-x' ({}) is not a valid number\n       Please choose a positive real # value\n".format(boxx))
    if boxy <= 0:
        usage()
        sys.exit("Error: the value for '-y' ({}) is not a valid number\n       Please choose a positive real # value\n".format(boxy))
    if boxz <= 0:
        usage()
        sys.exit("Error: the value for '-z' ({}) is not a valid number\n       Please choose a positive real # value\n".format(boxz))
 
    outfile = open_output_file(fileName)
 
    # Set magic lattice numbers
    if lattice == 0:
        a = (4 * 18.01 / (0.602 * density)) ** (1.0 / 3.0)
        acut = rcut * math.sqrt(2.0)
        if acut > a:
            a = acut
            newDensity = 4.0 * 18.01 / (0.602 * a ** 3)
            print("using density of {} to match cutoff value ({})\n".format(newDensity, rcut))
    elif lattice == 1:
        a = (18.01 / (0.602 * density)) ** (1.0 / 3.0)
        acut = rcut
        if acut > a:
            a = acut
            newDensity = 18.01 / (0.602 * a ** 3)
            print("using density of {} to match cutoff value ({})\n".format(newDensity, rcut))
 
    nxFloat = boxx / a
    nx = int(nxFloat + nxFloat / abs(nxFloat * 2))
    nyFloat = boxy / a
    ny = int(nyFloat + nyFloat / abs(nyFloat * 2))
    nzFloat = boxz / a
    nz = int(nzFloat + nzFloat / abs(nzFloat * 2))
 
    anew = min(boxx / nx, boxy / ny, boxz / nz)
 
    if anew < acut:
        anew = acut
 
    nxFloat = boxx / anew
    nx = int(nxFloat + nxFloat / abs(nxFloat * 2))
    nyFloat = boxy / anew
    ny = int(nyFloat + nyFloat / abs(nyFloat * 2))
    nzFloat = boxz / anew
    nz = int(nzFloat + nzFloat / abs(nzFloat * 2))
 
    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 {} to make chosen box sides work out\n".format(newDensity))
 
    cellLengthX = boxx / nx
    cellLengthY = boxy / ny
    cellLengthZ = boxz / nz
 
    cell2X = cellLengthX * 0.5
    cell2Y = cellLengthY * 0.5
    cell2Z = cellLengthZ * 0.5
 
    xCorr = []
    yCorr = []
    zCorr = []
 
    if lattice == 0:
        xref = [0.0, 0.0, cell2X, cell2X]
        yref = [0.0, cell2Y, cell2Y, 0.0]
        zref = [0.0, cell2Z, 0.0, cell2Z]
    
        counter = 0
        for z in range(nz):
            for y in range(ny):
                for x in range(nx):
                    for uc in range(4):
                        xCorr.append(xref[uc] + cellLengthX * x)
                        yCorr.append(yref[uc] + cellLengthY * y)
                        zCorr.append(zref[uc] + cellLengthZ * z)
                    counter += 4
    elif lattice == 1:
        xref = [cell2X]
        yref = [cell2Y]
        zref = [cell2Z]
 
        counter = 0
        for z in range(nz):
            for y in range(ny):
                for x in range(nx):
                    xCorr.append(xref[uc] + cellLengthX * x)
                    yCorr.append(yref[uc] + cellLengthY * y)
                    zCorr.append(zref[uc] + cellLengthZ * z)
                    counter += 1
 
    (waterCase, invalidWater) = validate_water(waterName)
    if invalidWater:
        print("Warning: '{}' is not a recognized water model name.\n".format(waterName))
        print("         Use the '-m' option to generate a 'water.inc' with the\n")
        print("         recognized water model geometries.\n\n")
    if waterName == 'DPD':
        density *= 0.25
    elif waterName == 'CG2':
        density *= 0.5
 
    write_output_file(outfile, nMol, xCorr, yCorr, zCorr, cutoff, alpha,
                          boxx, boxy, boxz, waterName, waterCase, invalidWater,
                          doRandomize, doWaterInc)
    outfile.close()
 
    print("The water box \"{}\" was generated.\n".format(fileName))
 
    if doWaterInc:
        try:
            outfile = open("./water.inc", "w")
        except:
            sys.exit("Error: can't open file {}\n".format("water.inc"))
 
        outfile.write("#ifndef _WATER_INC_\n#define _WATER_INC_\n")
        print_cl(outfile)
        print_na(outfile)
        print_ssd_e(outfile)
        print_ssd_rf(outfile)
        print_ssd(outfile)
        print_ssd1(outfile)
        print_tip3p(outfile)
        print_tip4p(outfile)
        print_tip4pew(outfile)
        print_tip5p(outfile)
        print_tip5pe(outfile)
        print_spce(outfile)
        print_spc(outfile)
        print_dpd(outfile)
        print_cg2(outfile)
        print_ssdq(outfile)
        print_ssdqo(outfile)
        print_tip4p_ice(outfile)
        print_tip4p_2005(outfile)
        print_spc_hw(outfile)
        print_ne6(outfile)
        print_tip3p_fb(outfile)
        print_tip4p_fb(outfile)
        print_opc(outfile)
        print_opc3(outfile)
        outfile.write("\n\n#endif")
 
        print("The file \"water.inc\" was generated for inclusion in \"{}\"\n".format(fileName))
 
        outfile.close()
 
if __name__ == "__main__":
    main(sys.argv[1:])