Simulation algorithm

The simulation algorithm is shown as a flow chart in figure 4.1. It involves a method called Monte-Carlo sampling(MC sampling) which is explained in the section following this one. In order to start a simulation, one needs to specify six dimensionless parameters – the chain length N, the bending stiffness_σ, the hydrodynamic bead radiusa, the excluded volume radiusR_EV and the first bead’s hydrodynamic radiusa₁and the first bond lengthd₁(the latter 5 in units ofd). The algorithm itself requires three more parameters – the numerical time step size∆t, the constraint tolerancetol₁ and the material frame tolerancetol₂. If applicable, the twisting stiffness and reference curvatures (see equation 2.156, page 61) must also be specified.

Initialization MC sampling is used to draw N−2 bending and azimuthal angles from the Boltzmann distribution. Starting the chain from the origin into z-direction, these angles are then used to construct the initial chain. The initial chain is equipped with a Bishop frame which starts at the first chain edge as unit vectors inxand ydirections, and is then iteratively con-structed using parallel transport. To find the initial material frame, the Bishop frame is rotated by integrated twist angles which are also drawn from their Boltzmann distribution via MC sampling.

Integration step The integration step itself is a midpoint algorithm adapted from refer-ence [67]. It requires that we first calculate the diffusion matrix Dand the stochastic force p2 (B^T)⁻¹∆W/∆tby performing a Cholesky decomposition ofD, and thenpreprojectit. Preprojec-tion means that the force itself is modified using Lagrangian multipliers (as already explained) such that it fulfills the geometric constraints at the level of forces, not velocities. This stochastic force is held fixed throughout the integration step. Then the deterministic forces are evaluated at the chain’s state~rat the beginning of the integration step. Next, the Lagrangian multipliers are determined and the chain’s coordinates are updated with a halved time step size∆t/2, arriving at the midpoint~r^mid. At the midpoint, the diffusion matrix, deterministic forces and Lagrangian multipliers (but not the stochastic force) are reevaluated and the resulting velocity d_t~ris used to perform the full time step~r7→~r+d_t~r·∆t. Throughout this thesis, the value∆t=10⁻³t₀ is used. Simulated trajectory lengths range from at least 10⁴t₀to several 10⁵t₀, corresponding to

≈10⁷−10⁸or more integration steps. One such trajectory typically takes from several hours up to two days to simulate on a regular desktop PC, using four processor cores such as an Intel Core i5-6600 at 3.30 GHz.

draw initial twist angles using the stochastic force at the

initial coordinates, not the

Figure 4.1: This figure shows the numerical algorithm used to generate trajectories from the nondimensional, constrained Langevin equation (equation 4.2, page 86). Starting the algorithm at the top left corner, the initial state of the chain is drawn from the Boltzmann distribution using MC sampling. Afterwards we repeatedly perform integration steps, in between which the observables of interest are stored to an external file. Each integration step that results in any two beads overlapping gets rejected (and the chain before the step restored), and accepted if it does not. After successful integration steps, the geometric constraints and the orthonormality of the material frame are repaired if a predefined tolerance has been violated. For simulations where the material frame is not required, all orange boxes can be omitted.

Overhead Before each integration step, the chain’s current state is stored and observables of interest are saved to an external file. After the integration step, it is checked whether any two beads overlap sterically, i.e. min

i,j |~r_i−~r_j| <2R_EV. If this is the case, the previously stored state before the integration step is retrieved and the integration step repeated (with a new realization of the random force). In this case the retrieved state gets stored to the external file twice for two consecutive time steps, which is necessary for correct Boltzmann statistics¹. When a time step is rejected, one may therefore imagine the chain remaining in its state for a time interval of length ∆t. If the integration step was successful, the material frame is updated via paral-lel transport and it is checked whether either the constraints or the material frame require repair.

Repair step There are two independent repair steps. Firstly, whenever any bond length deviates from 1 (or d₁/d) by an amounttol₁=10⁻³ or more, the chain is separated into two parts: the part before and after that bond. Both parts are then rigidly shifted into the tangential direction of that bond, such that (i) the bond length after the shifting process equals 1, and that (ii) the center of mass of the entire chain is preserved. This “repair” process for the constraint lengths is required because the method of using Lagrangian multipliers does lead to time increments which are parallel to the constraints, but due to the finite numerical time step size, higher order errors slowly add up during a simulation. Secondly, the material vectorsm~⁽³⁾are known exactly in each time step because they equal the tangent vectors. The material vectorsm~⁽¹⁾, however, are correctly initialized before the simulation, and are then iteratively parallel transported along with the dynamics of the tangent vectors. While normalizing these vectors is trivial, we cannot be certain that they do not slowly acquires a tangential component over the course of the simulation due to numerical errors. Thus, whenever any|~m⁽¹⁾·m~⁽³⁾| >tol₂=10⁻⁵, we subtract the tangential component via

(4.12) m~⁽¹⁾7→m~⁽¹⁾−(m~⁽¹⁾·m~⁽³⁾)~m⁽³⁾

and normalize~m⁽¹⁾afterwards. The remaining material vector~m⁽²⁾is found as the cross product of the other two.

1According to the detailed balance condition, the probability outflux from a state close to a region of zero probability must be lowered.