Monte-Carlo sampling

In order to sample the Boltzmann distribution in presence of excluded volume effects, we make use of an algorithm which is called theMetropolis Monte-Carlo,Markov-Chain-Monte-Carloor simply Monte-Carloalgorithm [126]. We simply refer to it asMC sampling. MC sampling constitutes an efficient method to sample a given probability distribution defined on a high-dimensional space Ω. It only requires that the values of the (not necessarily normalized) probability distribution can be numerically evaluated. For the Boltzmann distribution, this condition breaks down to the energy functionEbeing evaluable.

The MC sampling algorithm is shown in figure 4.2 as a flowchart and works as follows: A discrete random walkL:={x_i|i=1, 2, 3, . . . } withinΩis designed in such a way that the probability to go from a state Ato a stateBis the same as going fromBto A. The index iis a step index, not time, and xparametrizes the model chain. The random walk must satisfy the detailed balance condition: The transition ratekis then written as the product

(4.14) k(A,B)=k^trial(A,B)·p^accept(A,B)

wherek^trial(A,B) is isotropic and uniformly distributed around A, describing trial stepsfrom point Ato a pointB in its vicinity. It is crucial that these trial steps are unbiased and allow all points withinΩwith non-zero probability to be reached by the algorithm. In particular, this means thatk^trial(A,B)=k^trial(B,A), and thus the detailed balance condition above reduces to

(4.15) p^accept(B,A)

After a trial step has been chosen, it is accepted with probability p^accept(A,B), or otherwise rejected and a new trial step starting from Ais attempted. The core idea of MC sampling is to implement unbiased trial steps and then cleverly choose the functional form ofP^accept(A,B) such that detailed balance is obeyed while the ratio p^accept(A,B)/p^accept(B,A) equals the Boltzmann factor. This condition is fulfilled for

(4.16) P^accept(B,A)=min

which can be easily checked for bothE(A)≤E(B) and forE(A)>E(B). Therefore, it suffices that any chain conformation’s energy can be computed in order to sample the Boltzmann distribution.

draw random number calculate

no yes

make random step update list choose initial state start with empty list

update state

Figure 4.2: This flow chart shows the algorithm called Monte-Carlo (MC) sampling. Starting at the top, an arbitrary state of the system (e.g. a straight chain) is initiated. The state’s parametrization~xis then repeatedly varied in a random manner (trial steps) leading to a new state. If the new state’s energy is lower than the energy of the current state, the trial step is accepted and the current state of the system is updated to this new state. If the energy is higher than before, the trial step is only accepted with probability exp (−∆E/T) and otherwise rejected. Before each trial step, the current state is saved to a list, regardless if the preceding trial step was successful or not. After a sufficient number of trial steps, this list is saved to an external file. Its elements are distributed according to the Boltzmann distribution.

In this thesis, the parametrization used to represent chain conformations withinΩis the following:

A chain consisting ofN beads is constructed via itsN−1 tangent vectors, where the first tangent vector is held fixed because the chain’s global orientation in space is irrelevant for its energy function. The remainingN−2 tangent vectors are represented in spherical coordinates viaN−2 torsional (azimuthal) angles_φandN−2 cosines of bending (polar) angles cos_θrespectively. The polar axis for each bond is the previous bond’s tangent vector, and the azimuthal angle of tangent vector~t_i is its angle to the plane spanned by~t_i−1 and~t_i−2. Trial moves consist of changing these 2N−4 values by independently and simultaneously adding random numbers to each of them. These random numbers are uniformly distributed on the interval [−∆,∆], and we employ reflecting boundary conditions for cos_θin order to stay in the interval [−1, 1]. Note that, without further discrimination, these trial moves are unbiased and lead to a uniform distribution of cos_θ and_φ, which is the correct equilibrium distribution on a sphere’s surface (see the footnote on page 28) in the absence of bending energy and excluded volume effects. The trial moves are then accepted or rejected using the Boltzmann factor, and thus detailed balance is achieved. For calculating the Boltzmann factor, the energy of chain conformations where any two beads overlap is set to infinity such that corresponding trial moves are automatically rejected:

(4.17) p³

0 if any two beads overlap 1

Whenever MC sampling is used, the value of∆is chosen anew such that the average acceptance rate of trial moves is around 30%. The average acceptance rate of the MC sampling algorithm does not influence the sampled distribution, but does determine the efficiency (speed) with with different regions ofΩare sampled. It is typically larger for small step sizes∆because in that case, the energy differenceE(B)−E(A) is small and the Boltzmann factor close to 1. We remark that a Langevin simulation samples the Boltzmann ensemble as well, but is a computational overkill when used for this purpose alone. That is because MC sampling allows for much larger “time steps” and does not require the calculation of forces, diffusion matrices, Choleski decompositions and Lagrangian multipliers. Note also that the MC sampling algorithm is not altered by the presence or absence of the bound-end constraint, demonstrating that it results in correct statistics, notcorrect dynamics! To give the reader an intuition, in practice (using my code) a sufficiently long Langevin simulation takes a few days, and sufficiently long MC sampling a few minutes.