The input file..

# Monte Carlo Simulations

Any line that begins with a hashtag (#) will be assumed to be a comment, by towhee.

inputformat

Here, we specify the input format for our configuration file. One option is ‘Towhee’, while the other is ‘Lammps’

ensemble

Enter the ensemble you are working in, for general property calculation (say, density) 'npt' is an obvious choice (density is measured at constant P/T ). Also different ensembles require different input parameters, check the manual before digressing

temperature

Since we chose an NPT ensemble, we will have to specify the constant Temperature, Pressure and the number of molecules
Let’s stick with the room temperature (in Kelvin): 298.15

pressure

Again, normal atmospheric pressure will suffice (in kPa): 101.325

nmolty

nmolty is the number of different types of molecule in the system.
As an example, in case our system has “water+ionic liquid+co2”, then nmolty will be 3

nmolectyp

nmolectyp refers to the number of molecules of each type in the system.
If our system has three different molecule, the we would specify three values (separated by spaces): 500 200 100

numboxes

Number of box. A box can be thought of as an ensemble or a phase.
1 is sufficient for normal property calculations, while you may need two boxes if you have two different phases say vapor and liquid.

stepstyle

Stepstyle can be 'moves' , where each move is 1 MC step, or it can be ‘cycle’. In 1 cycle there are N moves, where N is the number of olecules in the system.

nstep

nstep is the number of steps the MC simulation will run for. A normal simulation might consist of an equilibration stage of 25000 steps and another 20000 steps for production run.

printfreq

Frequency for printing the result in the output file. A printfreq of 200 will print the Pressure, density, etc. values every 200 steps.

blocksize

Towhee doesn’t calculate the density after each move. Rather it collects a number of steps as a block and calculates the average value of properties in that block. 200 is sufficient.

moviefreq

As far as I remember, I used the analyse_movie function once, to visualize two phase equilibria in VLE. Never used it again. However, it won’t be a good idea to ignore it all together, assign: 200

backupfreq

Highly useful quantity, if your simulation stops due to power failure, this is your only hope. PRO TIP : Keep all the frequency values same, this will avoid any confusion later on. Recommended value: 200

runoutput

An option that let’s you select the kind of information you want in your output file. 'full', as we don’t wan’t to miss out on any kind of information

pdb_output_freq

Towhee generates a pdb file of the system after every certain moves. 200 it is. (Increase this value in case you need more data for analysis)

loutlammps

LAMMPS is a molecular dynamics software. This option lets you generate file which will serve as LAMMPS input file '.true.' is helpful, and doesn't cost much.

pressurefreq

Virial pressure is the system of pressure calculated using statistical mechanics. These P values are helpful in further analysis of data. Pressure calculation takes time, so if you don’t really need them you can increase it. 200 isn’t a bad idea for eqb run.

trmaxdispfreq

If your basics of MC simulation are clear, you will understand what it is. The max displacement value you provide has to be changed to achieve the required acceptance rate. 20 means that it will refresh the trmax value after every 20 step.

volmaxdispfreq

Same logic, extended for Volume move. Lower value is helpful in equilibrating the system more elegantly. Assign: 20

potentialstyle

There are two options here, 'internal' or 'external'. 'External' is used when you have an external code to compute the energy of the system. Internal is used otherwise. We stick to: internal

ffnumber

The number of different force field files you are using. Though, you should try to use a single force field in your simulation, but in case it is Tip3P and OPlS-AA, well you need to specify: 2.

ff_filename

Next, you need to specify the directory where the forcefield files are located. A better idea is to copy the force field file from the towhee directory to your current directory, this will save you the trouble of copying the exact location of the file. Specify the full path, something like:
/towheebase/ForceFields/towhee_ff_Charmm22

classical_potential

Again, there are a number of options to choose the potential energy function. For our purposes, we stick to the ubiquitous:
Lennardjones

classical_mixrule

Mixing rules help to determine the epsilon and sigma parameters for intramolecular interactions. The choice of the mixrule should be consistent with your force field. For TraPPE force fields used here, we can stick to (arithmetic mean of σ, geometric mean of ε):
'Lorentz-Berthelot'

lshift

This option truncates the potential energy function at the cutoff and makes it equal to zero. We won't be doing that. Instead of arbitrarily bringing the potential energy function to zero we will apply a tail correction: .false.

ltailc

This option helps to apply tail corrections that ensure a sort of continuity while applying the cut offs: .true.

rmin

Don't confuse it with the cutoff value. Instead, this one is the inner cutoff or the closest distance two atoms can come together. (Lj goes to infinity for too small distances). This value ensures that mathematically we don't encounter infinity, instead we reject the move or avoid it somehow. The manual suggests a value of (in A): 1.0

rcut

The potential energy calculation cutoff value. This value needs to be handled carefully. Too short and your simulations won't give good results, too big and the computation time takes a jump. Ensure that this value is smaller than half of box length.
12

rcutin

This value is required for the configurational bias moves. Although the rcut value can be used, however, a smaller cutoff for the cbmc move allows a faster algorithm. What value to set? Manual suggests a value of 5 for non-columbic system and 10 otherwise. 10

electrostatic_form

Next, you need to specify whether you wish to calculate the columbic interactions. Depending on the system being studied, you can turn this one on or off. Our molecule is polar and charges do play a significant role. Hence, we will stick to coulombic interactions.
coulomb

coulombstyle

There are a couple of options. There is a third one too, but that is not recommended by the users themselves. We will stick with: ewald_fixed_kmax

kalp

Value used in the Ewald sum to compute alpha: 5.6

kmax

Maximum number of inverse space vectors to use in any dimension for the Ewald sum: 5

dielect

The dielectric constant used when computing coulombic interactions: 1

nfield

Use this option if you wish to apply any external field, for example, Wall, umbrella, etc. We aren't getting into those stuffs here.
0

linit

This option is crucial. Not for running the simulation from scratch, but for continuing a simulation after it crashes, or the power cuts off, or just dies. Remember, linit = true implies that the simulation is being started from scratch (which is the case here); incase you wish to continue a simulation that has dies. Look here. .true.

initboxtype

Towhee offers different ways to generate the configuration. We will stick to the option: dimensions

initstyle

This one is required if you set the initboxtype to dimensions. Here we specify that the initial coordinates should be generated using configurational bias method. You can specify other methods too. Won't hurt. full cbmc

initlattice

All these options can be read in detail in manual, they are just trying to specify how the initial set of molecules will be generated. No worries if you don't understand. simple cubic

initmol

Here we specify again the number of molecules to be filled in each box. 500

inix iniy iniz

Rule of thumb: take the cube root of the number of molecules. Round it up to the nearest bigger integer (i.e. 8.1 becomes 9). This is the value you need to specify here, thrice. For example, cube root of 500 is 7.sth , so we specify: 8 8 8

hmatrix

Here we define the box size or length. Enter the x, y, z values in a matriz notation. For example:
22.6 0.0 0.0
0.0 22.6 0.0
0.0 0.0 22.6

# A side note

Next we specify the different monte carlo moves to be taken into account and also their respective probablities of occurence. All these pm* values denote probability.

pmvol

Probability of performing a volume move. For our NPT ensemble, we need to specify the probability for the box: 0.01

pmcb

Probablity for performing a configurational bias move. This move is the reason why towhee came into existence. Ensure that you write the cumulative value, that is, the value you specify here should be probability of cbmc move plus probability of volume move. We specify a value of: 0.33

pmtracm

Probablity for a translation move. Note that we specify the cumulative probabilty, this means that the value we specify here should be the sum of pvol, pmcb and finally the pmtracm. pmtracm should be kept a bit high, around 33% should be good: 0.67

pmrotate

Again, we add the cumulative number. Rotational move probability is usually kept same as the translational value. Note, this is the last monte carlo move we are specifying so we need to ensure that the sum of all probability is 1. 1.0

# other probabilties

There are a plethora of monte carlo moves that can be specified for monte carlo simulations. In this tutorial, we are keeping ourselves close to the elementary or the basic moves.

input_style

Next we need to tell towhee about the molecules that we are using in our system, how the atoms in the molecules are connected and what are the charges that they carry. Remember this section is to be repeated for each molecule in the system. These are called: 'basic connectivity map'

nunit

Specify the number of atoms in the first molecule. Ethanol has 4 (CH3, CH2, O and H): 4

nmaxcbmc

Set this to nunit. Simple. 4

lpdbnames

None of our concern. .false.

forcefield

Specify the force field we are going to use for our simulatoin. Ensure you don't mistype any letter here: TraPPE-UA

charge_assignment

Manual if your molecule has charge, none if it doesn't, else read manual. 'manual'

unit ntype qqatom

Next you need to tell about each of the atoms in the molecule. There are three sections for each atom. The first one is unit ntype qqatom. unit refers to atom number, it can be anything however, you need to be consistent. This will help towhee to remember which atom is connected to which one. ntype refers to atom type from the force field file. You can have a look here to get some help in deciding the atom types. qqatom refers to charge, if any, in the atom.
1 'CT3' -0.27d0

vibration

Here you specify the connectivity of the atom to other atoms. First you will tell how many atoms a single atom is connected to , then you will specify the atom numbers (nunit) of all those atoms.
4
2 5 6 7

improper torsion

This one refers to torsion terms if any. You can ignore this parameter if your force field file asks you to. 0

# the end

Repeat the above section for each atom in the molecule, and we are good to go.