#!/bin/bash

# This script will create several benzene molecules and describe their interactions,
# using the Newtonian approach of the ClassicalDynamics program.
# Apart from bash (the main language) this script also uses bc and python for calculations.

# It consists of several functions containing the code for a part of the output file.
# Creation of this file is controlled by the "main" function, where the desired number of molecules can be submitted.
# A test calculation of your desired value, without creating the output, can be done by entering '1234' for N, then the number in question.
# Of course, the program can also be run without the interaction. Uncomment "N", ""Generate_cbs" (bottom line) and comment "main".
# The function of the important variables is explained where they're used.

################################ Startup variables: ####################################################
OutputFile="Water_Benzene.cbs" # The file to which the output should be directed

mN="100" # The maximum number of benzene and water molucules that may be calculated, 
mNO="100" # to prevent long cpu usage.

dBenzene="8.979212737211E-04" #This is the density factor in atomic units, equal to 0.8765 kg/l in SI.
dWater="4.943301883245E-03" #Density factor in a.u.m./Bohr^3, equal to 0.9982 kg/l.
dMixture="4.241232703061E-3" #Density factor of a 2,5M mixture in a.u.m./Bohr^3, equal to 0.9806 kg/l.
########################################################################################################

function Main() {	
# This is the main function, controlling the input process and starting the output process if required.

if [ $1 = "" ] ; then 
  Input_vars # Provide a prompt to enter the desired values for the required variables.
  Check_vars # Checks whether the values are in the allowed range.
elif [ $1 = "-*"] ; then # Route taken when an option is specified.
  case $1 in
    -h*)
      Help_Info
      ;;
    -d*)
      echo "Running in Direct mode."
      Read_vars
      Check_vars
      ;;
    -t*]
      echo "Running in Test mode"
      Read_vars
      Test_Calculation
      ;;
    -s*)
      Silent_Mode
      ;;
    -*)
      Help_Info      
  esac
else
  Help_Info
fi

}

function Input_vars {
# This prompts the user to enter the desired values for "N" and "NO", stored in "aN" and "aNO".

clear
echo "This script will generate a ClassicalDynamics input file with approximatly N benzene molecules and NO water molecules per benzene."
echo ""
echo "Fill in N:"
declare -i aN
read aN
echo "Fill in NO:"
declare -i aNO
read aNO

}


function Read_vars {
# Secures the variables from the command arguments for further use.
	
aN="$2"
aNO="$3"

}

function Check_vars {
# Checks whether the desired variables are in the approved range, then alter the check value.
# The outcome of the check value determines if the outputfile is created.

check="0"
		
if [ $aN -ge 1 -a $aN -le $mN ] ; then
  check=$[$check+0]
elif [ $aN -gt $mN ] ; then
  check=$[$check+1]
else
  check=$[$check+4]
fi

if [ $aNO -ge 0 -a $aNO -le $mN ] ; then
  check=$[$check+0]
elif [ $aNO -gt $mN ] ; then
  check=$[$check+2]
else
  check=$[$check+5]
fi

case $check in
  0)
    Generate_cbs # Check has been succesfull, the outputfile will be created.
    ;;
  1)
    echo "The number of benzene molecules may not exceed $mN," # Message displayed when you exceed the limitation value for "aN".
    echo "since creating $(echo "print $aN*($aN-1)*18" | python) benzene-benzene and"
    echo "$(echo "print pow($aN,2)*12*$aNO" | python) benzene-water interactions require lengthy calculations."
	;;
  2)
    echo "The number of water molecules per benzene may not exceed $mNO," # Message displayed when you exceed the limitation value for "aNO".
    echo "since creating $(echo "print $aN*($aN-1)*18" | python) benzene-benzene and"
    echo "$(echo "print pow($aN,2)*12*$aNO" | python) benzene-water interactions require lengthy calculations."
	;;  
  3)
    echo "The number of benzene molecules may not exceed $mN, and the number of water molecules per benzene may not exceed $mNO," 
    echo "since creating $(echo "print $aN*($aN-1)*18" | python) benzene-benzene and" # Message displayed when you exceed the limitation value for "aN" & "aNO".
    echo "$(echo "print pow($aN,2)*12*$aNO" | python) benzene-water interactions require lengthy calculations."
	;; 
  4)
	echo "Use a valid number for N." # Message displayed when no nonzero integer is filled in for "aN".
    ;;
  5)  
    echo "Use a valid number for NO." # Message displayed when no nonzero integer is filled in for "aNO".
    ;;
  6)
    echo "The values for \"N\" and \"NO\" are out of range or invalid."
    ;;
  9)
    echo "The values for \"N\" and \"NO\" are invalid."
    ;;
  *)
    echo "Error checking the variable values."
    ;;
esac

}

function Help_Info {
echo "Syntax is:"
echo "Water_Benzene_Generator [option] [argument1] [argument2]"
echo "[option] can be: -h or -help or help for info, -d or -direct for direct mode and -t or -test to review the calculations."
echo "In direct mode, [argument1] = desired number of benzene molecules and [argument2] = the number of water molecules per benzene."
echo "The same applies to the test mode.
echo "When no or no valid option is given, the script will run in standard mode."
}

function Test_Calculation() {
#Function to test the calculations done with "aN" and "aNO".
#It asks for a new value to test and prints a list of the variables calculated. 
  
  Perform_Calculations
  
  echo "$Nd, $nX,$nY,$nZ"
  echo "$dXYZ, $BoxSize, $N"
  echo "$Od, $oX,$oY,$oZ"
  echo "$odXYZ, $OBoxSize, $O"
  echo "$tdXYZ, $TBoxSize, $LBoxSize"

}

function Generate_cbs() {
#This function writes the outputfile.
#It displays a dot or x when a function is completed, so you can follow the progress.

printf "\nCreating $OutputFile\n"\
"\n"
complete="0"

Perform_Calculations
sleep 1
complete=$[$complete+5]
printf " $complete percent completed.\n"

Generate_Header > $OutputFile
sleep 1
complete=$[$complete+5]
printf "$complete\n"

Generate_Atoms >> $OutputFile
sleep 1
complete=$[$complete+30]
printf "$complete\n"

Generate_Molecules >> $OutputFile
sleep 1
complete=$[$complete+10]
printf "$complete\n"

Benzene_Molecular_Interactions >> $OutputFile
sleep 1
complete=$[$complete+10]
printf "$complete\n"

Water_Molecular_Interactions >> $OutputFile
sleep 1
complete=$[$complete+10]
printf "$complete\n"

Benzene_Water_Interactions >> $OutputFile
sleep 1
complete=$[$complete+25]
printf "$complete\n"

Generate_Footer >> $OutputFile
sleep 1
complete=$[$complete+5]
printf "$complete percent\n"\
"\nFile $OutputFile is built!\n"

}

function Perform_Calculations() {	
# This function performs the dimension calculations.
# "Nd" is the cubic-root of "aN", which is the desired number of molecules.
# "nX", "nY" and "nZ" are the rounded, integer values of "Nd".
# These are used to determine the actual number of molecules to be created.
# "dXYZ" is the fraction of space each molecule gets in the x,y and z direction.
# "BoxSize" is "N" divided by the density.
# The cubic root is used to determine the length of the sides.
# "N" is the total number of molecules, thus the product of the number in x,y and z direction.

Nd=$(echo "print pow($aN,1.0/3.0)" | python)
nX=$(echo "print int($Nd)" | python)	
nY=$(echo "print int($Nd)" | python)	 
nZ=$(echo "print int($Nd)" | python)
dXYZ=$(echo "print (1/$Nd)" | python)	

N=$(echo "print ($nX*$nY*$nZ)" | python)
aO=$[$aNO*$N]

Od=$(echo "print pow($aO,1.0/3.0)" | python)
oX=$(echo "print int($Od)" | python)	
oY=$(echo "print int($Od)" | python)	 
oZ=$(echo "print int($Od)" | python)
odXYZ=$(echo "print (1/$Od)" | python)

O=$(echo "print ($oX*$oY*$oZ)" | python)

tdXYZ=$(echo "print (1/($Nd+$Od))" | python)
BoxSize=$(echo "print pow((($N)/$dBenzene),1.0/3.0)" | python)
OBoxSize=$(echo "print pow((($O)/$dWater),1.0/3.0)" | python)
TBoxSize=$(echo "print pow((($N+$O)/$dMixture),1.0/3.0)" | python)	
LBoxSize=$(echo "print ($TBoxSize+10)" | python)	

}

function Generate_Header() {
# This script creates the header of the output file, containing the description.
# "LBoxSize" is the lenght of 1 axis of the cubic box, defined as "BoxSize" + 5.

printf "# This input file for ClassicalDynamics simualates the harmonical bend,\n"\
"# harmonical stretch and torsional interactions in benzene molecules,\n"\
"# the 12-6 Lennard-Jones interactions between the carbon atoms of the benzene molecules,\n"\
"# the harmonical bend, harmonical stretch and coulomb interactions in and between water molecules,\n"\
"# and the 12-6 Lennard-Jones interactions between the carbon atoms of the benzene molecules and the oxygen atoms of the water molecules.\n"\
"\n"\
"# Define some colors:\n"\
"	color id=0 rgb=[0.0,1.0,1.0] #!CB name=\"Aqua\"\n"\
"	color id=1 rgb=[1.0,0.0,0.0] #!CB name=\"Red\"\n"\
"	color id=2 rgb=[1.0,1.0,1.0] #!CB name=\"White\"\n"\
"\n"\
"# Create the box:\n"\
"	box periodic=[$TBoxSize,$TBoxSize,$TBoxSize]\n"\
"\n"

}

function Generate_Atoms() {
# This is the script that calculates the position for each molecule.
# Lowercase "n" is the number of the benzene molecule the script is currently describing. It ranges from 0 to "N" - 1.
# "i", "j" and "k" are local variables used within loops, to separate different molecules and their individual atoms.
# "nX", "nY" and "nZ" are the total number of benzene molecules on each axis. Their product equals "N".
# "X", "Y" and "Z" are the positions in their respective axes where the molecule is placed. 
# The formula with "BoxSize" calculates positions so that each molecule has the same amount of space.
# "BoxSize" is the scale factor, "dXYZ" the density of benzene.
# The center of mass is placed in the middle of the designated space.

printf "# Construct the particles and define the intramolecular interactions:\n"\
"# Generating a model of $N benzene molecules with density equal to 0.8765 kg/l ... \n"

n="0"
while [ $n -lt $N ]; do
   i="0"
   while [ $i -lt $nX ]; do
      X=$(echo "print $BoxSize*(($i+0.5)*$dXYZ-0.5)" | python)
      j="0"
      while [ $j -lt $nY ]; do	
         Y=$(echo "print $BoxSize*(($j+0.5)*$dXYZ-0.5)" | python)	
         k="0"
         while [ $k -lt $nZ ]; do
            Z=$(echo "print $BoxSize*(($k+0.5)*$dXYZ-0.5)" | python)
            Benzene_Particle_Code
            #printf "n=$n, i=$i,j=$j,k=$k, X=$X,Y=$Y,Z=$Z\n";	#T his line is for testing pourposes. 
            n=$[$n+1]	# Uncomment it and comment the line above to see the values of "N","n","i","j"&"k".
            k=$[$k+1];
         done
      j=$[$j+1];
      done
   i=$[$i+1];
   done;
done

printf "# Adding $NO water molecules per benzene, $O in total.\n"\
"# They are placed in a model with a density equal to 0.9982 kg/l.\n"

n="0"
no="0"
o="0"
while [ $o -lt $O ]; do
   ii="0"
   while [ $ii -lt $oX ]; do
      X=$(echo "print $OBoxSize*(($ii+0.5)*$odXYZ-0.5)" | python)
      jj="0"
      while [ $jj -lt $oY ]; do	
         Y=$(echo "print $OBoxSize*(($jj+0.5)*$odXYZ-0.5)" | python)	
         kk="0"
         while [ $kk -lt $oZ ]; do
            Z=$(echo "print $OBoxSize*(($kk+0.5)*$odXYZ-0.5)" | python)
            Water_Particle_Code
            #printf "n=$n, no=$no, o=$o, ii=$ii,jj=$jj,kk=$kk, X=$X,Y=$Y,Z=$Z\n";
            if [ $no -ge $[$NO-1] ] ; then
              no="0"
              n=$[$n+1]
            else        
              no=$[$no+1]
            fi
            kk=$[$kk+1]
            o=$[$o+1];
         done
     jj=$[$jj+1];
      done
   ii=$[$ii+1];
   done;
done

}

function Benzene_Particle_Code() {
# This script contains the actual code for each molecule and is called on by "Generate_Molecules".
# The variabeles "Ang","Lx", etc. ensure the random orientation of the benzene molecules.
# The "C" variables calculate the necessary id numbers for the atoms in 1 molecule.
# The script uses them to specify the intramolecular interactions.
# Molecule number 1 has id's 0..11, number 2 id's range from 100..111 and so on.

Ang=$(echo "($RANDOM % 360)" | bc)	
Lx=$(echo "c=(1000-($RANDOM%2000)); scale=4; (c/1000)" | bc)
Ly=$(echo "c=(1000-($RANDOM%2000)); scale=4; (c/1000)" | bc)
Lz=$(echo "c=(1000-($RANDOM%2000)); scale=4; (c/1000)" | bc)

C1=$[100*$n+0]
C2=$[100*$n+1]	
C3=$[100*$n+2]	
C4=$[100*$n+3]
C5=$[100*$n+4]
C6=$[100*$n+5]
H1=$[100*$n+6]
H2=$[100*$n+7]
H3=$[100*$n+8]
H4=$[100*$n+9]
H5=$[100*$n+10]
H6=$[100*$n+11]

# Here the center of mass is shifted to the new position.
# And the random orientation parameters are translated here.

printf "# Molecule $n\n"\
"	translate x=$X y=$Y z=$Z\n"\
"	rotate deg=$Ang axis=[$Lx,$Ly,$Lz]\n"\
"	# Carbon\n"\
"	particle id=$C1 c=0 m=21894.71 q=-0.115 r=1.46 x=[4.5,0.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n C1\"\n"\
"	particle id=$C2 x=[2.25,-4.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n C2\"\n"\
"	particle id=$C3 x=[-2.25,-4.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n C3\"\n"\
"	particle id=$C4 x=[-4.5,0.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n C4\"\n"\
"	particle id=$C5 x=[-2.25,4.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n C5\"\n"\
"	particle id=$C6 x=[2.25,4.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n C6\"\n"\
"	# Hydrogen\n"\
"	particle id=$H1 c=2 m=1827.81 q=0.115 r=0.701 x=[9.0,0.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n H1\"\n"\
"	particle id=$H2 x=[4.5,-8.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n H2\"\n"\
"	particle id=$H3 x=[-4.5,-8.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n H3\"\n"\
"	particle id=$H4 x=[-9.0,0.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n H4\"\n"\
"	particle id=$H5 x=[-4.5,8.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n H5\"\n"\
"	particle id=$H6 x=[4.5,8.0,0.0] p=[0.0,0.0,0.0] #!CB name=\"$n H6\"\n"

}

function Water_Particle_Code() {

Ang=$(echo "($RANDOM % 360)" | bc)	
Lx=$(echo "c=(1000-($RANDOM%2000)); scale=4; (c/1000)" | bc)
Ly=$(echo "c=(1000-($RANDOM%2000)); scale=4; (c/1000)" | bc)
Lz=$(echo "c=(1000-($RANDOM%2000)); scale=4; (c/1000)" | bc)

oO=$[$[100*$n]+$[3*$no]+15]
oH1=$[$[100*$n]+$[3*$no]+16]	
oH2=$[$[100*$n]+$[3*$no]+17]

printf "# Molecule $n $no ($o)\n"\
"	translate x=$X y=$Y z=$Z\n"\
"	rotate deg=$Ang axis=[$Lx,$Ly,$Lz]\n"\
"	particle id=$oO c=1 m=29165.1 q=-0.8476 r=1.38 x=[0.0,1.5,1.0] p=[0.0,0.0,0.0] #!CB name=\"$n $no O1 \"\n"\
"	particle id=$oH1 c=2 m=1827.81 q=0.4238 r=0.701 x=[0.0,0.75,3.75] p=[0.0,0.0,0.0] #!CB name=\"$n $no H1\"\n"\
"	particle id=$oH2 x=[0.0,1.5,3.75] p=[0.0,0.0,0.0] #!CB name=\"$n $no H2\"\n"

}

function Generate_Molecules() {

printf "\n# Creating the molecules by defining the intramolecular interactions.\n"\

n="0"
while [ $n -lt $N ]; do
   Benzene_Internal_Interactions 
   n=$[$n+1];
done

printf "\n"

n="0"
o="0"
while [ $o -lt $O ]; do
   no="0"
   while [ $no -lt $NO ]; do
      Water_Internal_Interactions
      CI="$CI $oO $oH1 $oH2"
      LJO="$LJO $oO"
      no=$[$no+1]
      o=$[$o+1];
   done
   n=$[$n+1];
done

}

function Benzene_Internal_Interactions() {

C1=$[100*$n+0]
C2=$[100*$n+1]	
C3=$[100*$n+2]	
C4=$[100*$n+3]
C5=$[100*$n+4]
C6=$[100*$n+5]
H1=$[100*$n+6]
H2=$[100*$n+7]
H3=$[100*$n+8]
H4=$[100*$n+9]
H5=$[100*$n+10]
H6=$[100*$n+11]

printf "\n# Benzene $n Internal interactions:\n"\
"	interaction harmonic f=0.3505 r0=2.06 { $C1 $H1 } #!CB name=\"HS $n C1H1\"\n"\
"	interaction harmonic f=0.3505 r0=2.06 { $C2 $H2 } #!CB name=\"HS $n C2H2\"\n"\
"	interaction harmonic f=0.3505 r0=2.06 { $C3 $H3 } #!CB name=\"HS $n C3H3\"\n"\
"	interaction harmonic f=0.3505 r0=2.06 { $C4 $H4 } #!CB name=\"HS $n C4H4\"\n"\
"	interaction harmonic f=0.3505 r0=2.06 { $C5 $H5 } #!CB name=\"HS $n C5H5\"\n"\
"	interaction harmonic f=0.3505 r0=2.06 { $C6 $H6 } #!CB name=\"HS $n C6H6\"\n"\
"	interaction harmonic f=0.4792 r0=2.642 { $C1 $C2 } #!CB name=\"HS $n C1C2\"\n"\
"	interaction harmonic f=0.4792 r0=2.642 { $C2 $C3 } #!CB name=\"HS $n C1C2\"\n"\
"	interaction harmonic f=0.4792 r0=2.642 { $C3 $C4 } #!CB name=\"HS $n C1C2\"\n"\
"	interaction harmonic f=0.4792 r0=2.642 { $C4 $C5 } #!CB name=\"HS $n C1C2\"\n"\
"	interaction harmonic f=0.4792 r0=2.642 { $C5 $C6 } #!CB name=\"HS $n C1C2\"\n"\
"	interaction harmonic f=0.4792 r0=2.642 { $C6 $C1 } #!CB name=\"HS $n C1C2\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C2 $C1 $H1 } #!CB name=\"HB $n C2CH1\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C3 $C2 $H2 } #!CB name=\"HB $n C3CH2\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C4 $C3 $H3 } #!CB name=\"HB $n C4CH3\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C5 $C4 $H4 } #!CB name=\"HB $n C5CH4\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C6 $C5 $H5 } #!CB name=\"HB $n C6CH5\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C1 $C6 $H6 } #!CB name=\"HB $n C1CH6\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C6 $C1 $H1 } #!CB name=\"HB $n C6CH1\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C5 $C6 $H6 } #!CB name=\"HB $n C5CH6\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C4 $C5 $H5 } #!CB name=\"HB $n C4CH5\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C3 $C4 $H4 } #!CB name=\"HB $n C3CH4\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C2 $C3 $H3 } #!CB name=\"HB $n C2CH3\"\n"\
"	interaction bending f=0.285 rad=2.094395102 { $C1 $C2 $H2 } #!CB name=\"HB $n C1CH2\"\n"\
"	interaction bending f=1.173 rad=2.094395102 { $C1 $C2 $C3 } #!CB name=\"HB $n C1-C3\"\n"\
"	interaction bending f=1.173 rad=2.094395102 { $C2 $C3 $C4 } #!CB name=\"HB $n C2-C4\"\n"\
"	interaction bending f=1.173 rad=2.094395102 { $C3 $C4 $C5 } #!CB name=\"HB $n C3-C5\"\n"\
"	interaction bending f=1.173 rad=2.094395102 { $C4 $C5 $C6 } #!CB name=\"HB $n C4-C6\"\n"\
"	interaction bending f=1.173 rad=2.094395102 { $C5 $C6 $C1 } #!CB name=\"HB $n C5-C1\"\n"\
"	interaction bending f=1.173 rad=2.094395102 { $C6 $C1 $C2 } #!CB name=\"HB $n C6-C2\"\n"\
"	interaction torsional n=3.0 f=0.0004464 deg=180.0 { $C1 $C2 $C3 $H3 } #!CB name=\"PT $n C123H3\"\n"\
"	interaction torsional n=3.0 f=0.0004464 deg=180.0 { $C2 $C3 $C4 $H4 } #!CB name=\"PT $n C234H4\"\n"\
"	interaction torsional n=3.0 f=0.0004464 deg=180.0 { $C3 $C4 $C5 $H5 } #!CB name=\"PT $n C345H5\"\n"\
"	interaction torsional n=3.0 f=0.0004464 deg=180.0 { $C4 $C5 $C6 $H6 } #!CB name=\"PT $n C456H6\"\n"\
"	interaction torsional n=3.0 f=0.0004464 deg=180.0 { $C5 $C6 $C1 $H1 } #!CB name=\"PT $n C561H1\"\n"\
"	interaction torsional n=3.0 f=0.0004464 deg=180.0 { $C6 $C1 $C2 $H2 } #!CB name=\"PT $n C612H2\"\n"\

}

function Water_Internal_Interactions() {

oO=$[$[100*$n]+$[3*$no]+15]
oH1=$[$[100*$n]+$[3*$no]+16]	
oH2=$[$[100*$n]+$[3*$no]+17]

printf "# Water molecule $n $no ($o) internal interactions:\n"\
"	interaction harmonic f=1.74616112472 r0=1.88972598858 { $oO $oH1 } #!CB name=\"HS $o O1H1\"\n"\
"	interaction harmonic f=1.74616112472 r0=1.88972598858 { $oO $oH2 } #!CB name=\"HS $o O1H2\"\n"\
"	interaction bending f=0.159147835842 deg=109.47 { $oH1 $oO $oH2 } #!CB name=\"HB $o H1O1H2\"\n"

}

function Benzene_Molecular_Interactions() {
# This funtion generates the intermolecular interactions.
# At this moment only the 12-6 Lennard-Jones is created between atom "Cni" and "Cjk".
# The outer loop runs until the total number of molecules is created.
# "i" goes from 0 to 5, for the 6 carbom atoms in each molecule.
# "j" takes the values 1.."N", since the interactions are with all the other molecules.
# "k" also varies from 0 to 5, accounting for 6 atoms.

printf "\n# Benzene intermolecular interactions\n"\

n="0"
while [ $n -lt $N ]; do	
   i="0"
   while [ $i -lt 6 ]; do
   Cni=$[100*$n+$i]
      j=$[$n+1]
      while [ $j -lt $N ]; do
         k="0"
         while [ $k -lt 6 ]; do
            Cjk=$[100*$j+$k]
            printf "	interaction lennardjones f=0.00106 r0=6.3697 { $Cni $Cjk } #!CB name=\"LJ C$Cni C$Cjk\"\n"
            k=$[$k+1];
         done
         j=$[$j+1];
      done
      i=$[$i+1];
   done
   n=$[$n+1];
done

}

function Water_Molecular_Interactions() {
# Creates the coulomb attraction between all water atoms,
# and the LJ interaction between all the oxygen atoms in water.
# The strings containing all the necessary id's, "CI" and "LJO", are created during molecule generation.

printf "\n# Intermolecular interactions\n"\
"	interaction coulomb f=1.0 {$CI } #!CB name=\"Coulomb All Water Molecules\"\n"\
"	interaction lennardjones f=9.887E-4 r0=5.96 {$LJO } #!CB name=\"LJ All Oxygen\"\n"\
"\n"

}

function Benzene_Water_Interactions() {
# This function generates the LJ interactions between the C and H of benzene and the O in water.

n="0"  
while [ $n -lt $N ]; do
   i="0"
   while [ $i -lt 6 ]; do
      nn="0"
      while [ $nn -lt $N ]; do 
         Cni=$[100*$n+$i]
         Hni=$[100*$n+$i+6]    
         no="0"
         while [ $no -lt $NO ]; do
            Ono=$[100*$nn+3*$no+15] 
            printf "	interaction lennardjones f=6.644E-4 r0=6.334 { $Cni $Ono } #!CB name=\"LJ C$Cni O$Ono\"\n	interaction lennardjones f=4.349E-4 r0=5.266 { $Hni $Ono } #!CB name=\"LJ H$Hni O$Ono\"\n"
            no=$[$o+1];
         done
         nn=$[$no+1];
      done
      i=$[$i+1];
   done
   n=$[$n+1];
done

}

function Generate_Footer() {
# This function creates the footer, containing the parameters for the dynamics simulation.

printf "\n# Constraints:\n"\
"\n"\
"# Conformation:\n"\
"	conformation n=10000 error=1.0E-4 maxstep=2.5\n"\
"\n"\
"# Temperature:\n"\
"	temperature k=3.166829E-6 constant=293.15\n"\
"\n"\
"# Dynamics:\n"\
"	dynamics dt=4.13411 tend=4.1341105461E7 t=0.0 error=3.67493E-6 snapshots=10\n"\
"\n"\
"# End input file"\
"\n"\
;
}

Main # Here the main function is initialized, to take control of the file creation.

function Silent_Mode {

#N="2"	# Capital "N" is the total number of benzene molecule the script is going to create.
#NO="27"
#O=$[$NO*$N]
		
}