Users Manual

Introduction

Input files

1. control

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

2. coordinate

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

3. 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.

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.

2. 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.

3. 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 f_ij(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 R_ij = R_i - R_j 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 f_nb(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 f_i(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.

4. 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, T₀. λ 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

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

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

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

Now, we can calculate the total pressure tensor:

5. Reset the net momentum to zero.
6. 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:

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

7. Calculate the average of parameters and print the restart file
8. Store the trajectory file
9. Sample

IBIsCO units

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, k_B×300 (Escale)
3. Unit of mass, atomic mass of carbon (Mscale)

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

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 (A_i), and width (w_i):

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 k_B 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 l_ci 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.