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):


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

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

T(t) = [

N ∑ i=1

mi v2i(t)]/(kBNf)

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):

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)

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:

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

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