IBIsCO
It is Boltzmann Inversion software for Coarse Graining Simulations

> Home
> Input files
        > control
        > coordinate
        > interaction
> IBIsCO algorithm
> IBIsCO units
> Boltzmann Inversion
> News
> Authors
> Contact Info.



Users Manual

Introduction

Input files

  1. control

This is the file that contains the general parameters for the simulation (time step, ensemble, thermodynamic conditions ).

Example control file

Ensemble

NVT

(NVE, NVT, NPT)

Temperature

300

(Temperature (K))

Pressure

101.3

(Pressure (kPa))

Natoms

4800

(Number of atoms)

Nsteps

1000

(Number of time steps)

DT

10

(Time step in fs)

TAUT

0.15D0

(Temperature relaxation time in ps)

TAUP

5.0D0

(Pressure relaxation time in ps)

BETA

1.0D-6

(Isothermal compressibility (1/kPa))

Cutoff

1.60

(The cutoff distance for non-bonded interaction in nm)

Neighbor_list_cutoff

1.60

(The distance in nm up to which pairs of particles are included in the neighbor list)

Update_neighbor_list

5

(Set intervals of neighbor list updates)

Nsampling

100

(Interval of sampling of quantities)

Ntrajectory

100

(Number of time steps between storing configuration)

Halt_Drift

1000

(Interval at which the net drift of the system is reset to zero)

Naverage

100

(Number of time steps between storing average data and restart file)

Non-Bonded

4

(Use non-bonded potential on 1..4 OR 1..5? 4=1..4, 5=1..5)

Interaction

GAUSSIAN

(Use Gaussian function for bond and bend interactions or tables? (GAUSSIAN, TABLE))

Restart_velocities

NO

(Do you want to restart the initial velocities (use Boltzmann distribution)? YES or NO)

# Total number of processors that you need to run the program will be NPX*NPY*NPZ+1

NPX

2

(Number of processors in X direction, THIS NUMBER SHOULD NOT BE LESS THAN TWO)

NPY

2

(Number of processors in Y direction, THIS NUMBER SHOULD NOT BE LESS THAN TWO)

NPZ

2

(Number of processors in Z direction, THIS NUMBER SHOULD NOT BE LESS THAN TWO)

Boundary

1.80

(Width of boundary region to receive the neighbor coordinates, Boundary >= Cutoff)

 

    1. Ensemble: specify which ensemble you want to use for simulation. Current option include three types micro canonical (NVE), canonical (NVT) and isobaric-isothermal (NPT) ensemble.

 

    1. Temperature: initial temperature in Kelvin (attention this temperature is used for thermostat and also for calculating the forces and potentials by Bolzmann inversion).

 

    1. Pressure: initial pressure in kPa (target pressure in NPT ensemble).

 

    1. Natoms: total number of atoms included in the system.

 

    1. Nsteps: the total number of time steps of the simulation.

 

    1. DT: the simulation time step in ps.

 

    1. TAUT: thermostat coupling time in ps.

 

    1. TAUP: pressure relaxation time in ps.

 

    1. BETA: isothermal compressibility (1/kPa).

 

    1. Cutoff: the cutoff distance for non-bonded interaction in nm.

 

    1. Neighbor_list_cutoff: the distance in nm up to which pairs of particles are included in the neighbor list.

 

    1. Update_neighbor_list: number of time steps between updating the neighbor list.

 

    1. Nsampling: number of time steps between printing the instantaneous values.

 

    1. Ntrajectory: number of time steps between writing the trajectory file.

 

    1. Halt_Drift: number of time steps between putting the net momentum equal to zero.

 

    1. Naverage: number of time steps between calculating the averages and printing the restart file.

 

    1. Non-Bonded: in this program there are two possible definitions of non-bonded pairs within the same molecules: atoms separated by more than three bonds ('1..5' interactions and above) or atoms separated by more than two bonds ('1..4' interactions and above).

 

    1. Interaction: in this program there are two options to define the bond and bend interactions. You can prepare Gaussian file and specify the Gaussian functions that are fitting to the bond length and angles distributions, or you can prepare the potential tables for bond and bend interactions.

 

    1. Reset_velocities: in this option you clarify that program must use the velocities are given in the initial coordinate file (NO) or initialize the particle velocities according to a Maxwell-Boltzmann distribution (YES).

 

    1. NPX: number of processors in X direction. This number must not be less that two.

 

    1. NPY: number of processors in Y direction. This number must not be less that two.

 

    1. NPZ: number of processors in Z direction. This number should not be less that two. In this program, one of the processors is assigned to write the output files. Then total number of processors that you need to run the program will be NPX+NPY+NPZ+1.

 

    1. Boundary: This program is using a domain decomposition strategy to use a parallel machine. To calculate the non-bonded interaction each processor needs to receive the boundary region from neighbor processors. In this case the width of boundary region is our cutoff, but some time after mapping from atomistic to coarse graining supper atoms are so far, and to calculate the bonded part of forces and potentials (angle and torsion) you need to use a bigger number for boundary width to receive all information from boundary processors.

 

  1. coordinate

This file contains the coordinates, the velocities and the connectivity table.

Example coordinate file

***** Hexane *****

Time     0      (ps)

***** Box Length(X, Y, Z in nanometer) *****

   6.0000       6.0000       6.0000

**************************************************************************************

********** Record for each atom is in the form:-                             **********

********** Index Atom_Type No._of_bonds X Y Z (coords.in m)       **********

********** Vx Vy Vz (in m/s) Indices_of_bonded_atoms                 **********

**************************************************************************************

      800       Total_No._of_molecules

      6       Atoms_in_Molecule_No.       1

   1

1

1

0.1110827430000

-3.0459108170000

-2.3974193558000

 

   0.0000000000000000E+000

0.0000000000000000E+000

0.0000000000000000E+000

2

 

   2

2

2

0.252678727000

-3.062108242000

-2.453072034500

 

   0.0000000000000000E+000

0.0000000000000000E+000

0.0000000000000000E+000

1

3

 

   3

2

2

0.2573729000000

-3.1743108940000

-2.5569831960000

 

   0.0000000000000000E+000

0.0000000000000000E+000

0.0000000000000000E+000

2

4

 

   4

2

2

0.40422115400000

-3.2018947220000

-2.5899031410000

 

   0.0000000000000000E+000

0.0000000000000000E+000

0.0000000000000000E+000

3

5

 

   5

2

2

0.41202051600000

-3.3455466990000

-2.6419826890000

 

   0.0000000000000000E+000

0.0000000000000000E+000

0.0000000000000000E+000

4

6

 

   6

1

1

0.34237269500000

-3.3405625230000

-2.7781198950000

 

   0.0000000000000000E+000

0.0000000000000000E+000

0.0000000000000000E+000

5

      6       Atoms_in_Molecule_No.       2

   7

1

1

-0.6092115250000

0.47064512800000

0.14301834600000

 

   0.0000000000000000E+000

0.0000000000000000E+000

0.0000000000000000E+000

8

   ...

 

    1. Comment

 

    1. Time: previous simulated time

 

    1. Comment

 

    1. Box length for X, Y, and Z directions in nanometer

 

    1. Comment (5 lines)

 

    1. Total number of molecules.

 

    1. Total number of atom in the molecule + Comment (Atoms_in_Molecule_No.) + index of the molecule.

 

    1. For each atom in the molecule, specify the index of the atom, type of the atom (correspond to the atom type in the interaction file), number of other atoms connected to this one, coordinate (Sx, Sy, Sz), velocity (Vx, Vy, Vz) and indices of the connected atoms. (The connectors (bonds) are specified twice: one time for each of the two atoms forming the connection.).

 

  1. interaction

This is the file that contains types of atoms, mass and different type of bonds, angles, torsions and non-bonded interactions. In this file you have to specify the name of the potential files for bonds and angles, if you want to use the table for these interactions, and also the name of the tabulated potential file for the different type of non-bonded interactions.

Example interaction file

******* CG trig2 bulk numerical potential 1=bead 1st gen. 2=bead 2nd gen.

      2    No._of_different_atom_types

**** Type

Label

Mass/(g/mole)

     1

CH3

15.023475

     2

CH2

14.01565

      2    No._of_different_types_of_bonds              (Distance : nanometer, Potential: kJ/mol)

****Type_1

Type_2

pot_file

     1

2

a.pot

     2

2

a.pot

      2   No._of_different_types_of_bond_angles       (Angle : Degree, Potential : kJ/mol)

****Type_1

Type_2

Type_3

pot_file

     1

2

2

a-b.pot

     2

2

2

a-b.pot

      2   No._of_different_types_of_torsions              (Angle : Degree, Potential: kJ/mol)

****Type_1

Type_2

Type_3

Type_4

pot_file

     1

2

2

2

a-b-b.pot

     2

2

2

2

a-b-b.pot

      3    No._of_non-bonded_interactions              (Distance : nanometer, Potential: kJ/mol)

****Type_1

Type_2

nb_pot_file

     1

1

aa.pot

     2

2

bb.pot

     1

2

ab.pot

 

    1. Arbitrary title

 

    1. The number of different atom type (attention the expression after this number is a keyword and must not be changed).

 

    1. Comment

 

    1. For each atom type, specify the arbitrary type number, atom label and mass (in g/mole).

 

    1. The number of different of types of bonds

 

    1. Comment

 

    1. Specify the atom types for each type of bond. If you want to use the potential table for bond interaction, you should put the name of files after the atom types. If you don't want to use tabulated potential, don't put any things after the bond types.

 

    1. The number of different types of bond angles.

 

    1. Comment

 

    1. For each type of angles, specify the atom types that form the angle. If you want to use the potential table for bend interactions, put the names of corresponding files after the atom types. If you don't want to use tabulated potential, don't put any things after the angle types.

 

    1. The number of different types of torsions.

 

    1. Comment

 

    1. For each type of torsions, specify the atom types and name of the potential table.

 

    1. The number of non-bonded interactions.

 

    1. Comment

 

    1. The atom types that are belong to each kind of non-bonded interactions, and the name of potential table.

 

IBIsCO algorithm

The position of atoms on the continuous chain is saving on the S(t) arrays. The procedure in IBIsCO is as follows:

  1. Calculate the scaled positions R(t) of the atoms of the chain within the primary cell, i.e. all the components of the R(t) lie between Box-Length/2, using the following sequence of operations:

R(t) = S(t) 2/BOX

R'(t) = S(t) - 2 INT(S(t)/2)

R(t) = (R'(t) - 2 INT(R'(t))) BOX/2

where by INT we refer to the FORTRAN function which returns the integer part of the expression in the bracket.

  1. Make the neighbour list (if necessary).

If the size of the box is bigger than three times the neighbour list cut off, we use a link cell algorithm to make the neighbour list. The simulation box is divided into cubic cells with the neighbour list cut off length, and instead of the looking for the neighbour atoms in the whole box, searching is restricted to the neighbour cells. If the box size is smaller, the neighbour list update is performed by searching all possible atom pairs. All bonded atoms are excluded from the neighbour list.

  1. Calculate the bonded and non-bonded interactions.

The non-bonded part of the force is calculating by looping over the neighbour list. The contribution of the force fij(t), between atoms i and j to the virial tensor is calculating by:

V(t)Pnb(t) =



i



j>i

Rij(t) fij(t)

where Rij = Ri - Rj is the relative poison vector of the atoms i and j, and the double sum is over the all interaction pairs. At this point we are storing the total non-bonded force on each atom in a second vector fnb(t). To calculate the bonded part of the force, we are making a loop over all the atoms, and for each atom, over all bonds, angles and torsions to which the atom belongs. The bonded part of the virial tensor is calculated using the position of atoms on the continuous chain:

V(t)Pb(t) =

N

i=1

Si(t) [ fi(t) - fnb i(t)]

where fi(t) is now the total force on a atom. For the virial, the forces due to the non-bonded part are subtracted as their contribution has already been calculated.

  1. Update the positions and velocity of the atoms.

If constant temperature conditions are used (NVT, NPT):
First, from the temperature at
t - ∆t, T(t - ∆t), the temperature scale factor λ is obtained.

λ = [1 + ∆t/τT(T0/T(t - ∆t) - 1)]1/2

where τT is the relaxation time which controls the rate at which the system tends towards the desired temperature, T0. λ is chosen equal to 1 in NVE simulations.
Second, calculate the new positions using the leapfrog algorithm:

vi(t + 1/2∆t) = [vi(t - 1/2∆t) + fi(t)/mi ∆t] λ
si(t + ∆t) = si(t) + vi(t + 1/2∆t)∆t

vi(t) is calculated from the vi(t - 1/2∆t) and the vi(t + 1/2∆t):

vi(t) = (vi(t + 1/2∆t) + vi(t - 1/2∆t))/2

The on-step temperature is calculating by means of vi(t):

T(t) = [

N

i=1

mi v2i(t)]/(kBNf)

where Nf is the total number of degree of freedom in the system.
The on-step kinetic contribution to the pressure tensor is calculating by
vi(t):

V(t)Pkin(t) =

N

i=1

mi(t) vi(t) vi(t)

Now, we can calculate the total pressure tensor:

P(t) = Pkin(t) + Pb(t) + Pnb(t)

  1. Reset the net momentum to zero.
  2. Calculate the new volume and scale the positions (in NPT ensemble).


In IBIsCO, box lengths and continuous positions of the atoms are scaled with same factor. The pressure scaling factor is:

μ = [1 + ∆t/τP β(P(t) - P0)]1/3

where τP is the pressure relaxation time, β is the isothermal compressibility and P0 is the target pressure.

  1. Calculate the average of parameters and print the restart file
  2. Store the trajectory file
  3. Sample


IBIsCO units

Parameter

Unit

Length

Nanometre (nm)

Volume

Cubic nanometre (nm3)

Angle

Degree

Time

Pico-second (ps)

Energy

kJ/mol

Force

kJ/mol.nm

Mass

g/mol

Velocity

nm/ps

Temperature

kelvin (K)

Pressure

kPa

Mass density

kg/m3

The basic choice for the reduced units (indicated by an asterisk) is the following:

  1. Unit of length, 1 nanometre (Rscale)
  2. Unit of energy, kB300 (Escale)
  3. Unit of mass, atomic mass of carbon (Mscale)

and from these basic units, all other units follow (time, pressure and temperature).

t* = t kB / Escale   ,   P* = 1000 Rscale3 P / Escale   ,   T* = T / Rscale (Escale / Mscale)1/2

Boltzmann Inversion

In IBIsCO, there are two options to define the bond and bend interactions. You can prepare the Gaussian file and specify the Gaussian functions that have been fitted to the bond lengths and angle distributions, or you can prepare potential tables for bond and bend interactions. The strategy to calculate the potential and force during the simulation is the same for both options. We are using liner interpolation to find the forces and potentials. In the case of GAUSSIAN we need to make the force and potential tables at the start of the simulation.
In the Gaussian file, a multi peaked distribution of a structural parameter
θ (bond or angle) is approximated by a sum of n Gaussian functions characterized by their centers (θci), total area (Ai), and width (wi):

Given a distribution P(θ) of some structural parameter θ, the corresponding potential is derived by a simple Boltzmann inversion of P(θ). The corresponding potential, V(θ), obtained by Boltzmann inversion is written as follows:

where kB is the Boltzmann constant and T is the temperature. The corresponding expression for the force is:

If we define

,

the potential and the corresponding expression for the force can be written in a more compact form:




For the bond potentials, it is possible to derive similar equations for potential and forces:




where
lci is the center of the Gaussian describing bond length distributions.
Attention: the potential is made from a distribution by Boltzmann inversion. In the inversion, the temperature (as specified in the control file) is used. This means that, with the GAUSSIAN option, simulation at different temperature have different potential energy functions for the bonded interactions, if the same distribution is used as input. To avoid this, use tabulated potentials.

References

This is just filler text so you can see what the page looks like full of text.

  © 2008 Theoretical Physical Chemistry TU Darmstadt - All rights reserved.