Brownian Dynamics Simulation

Now let us introduce and compare several ways of of integrating the coupled stochastic equations of motion (9.65) in a so-calledBrownian dynamics(BD)simulation. We will close this section making some remarks on the generation of Gaussian random numbers.

9.2.1 Conventional Brownian Dynamics

Theconventional Brownian dynamics algorithm(CBD) of Ermak [Erm75,All87] uses a first-order Euler method to give a solution for equation (9.65) for the many body system which takes the form r_i(t+∆t) =r_i(t) +∆r_i(∆t) (9.72) with the displacement

∆r_i(∆t) =βD₀F_i∆t+W_i (9.73) after the time step ∆t.³⁷ The vector Wi of random displacements of particle iis sampled from a Gaussian distribution with zero mean and variance hW_i·W_ji=6D₀∆tδ_{i j}³⁸. This approach is based on the assumption that the systematic forces F_i remain essentially constant during the time step ∆t. The time scales with the friction coefficientξ, as∆t always occurs divided by ξ (D₀=k_BT/ξ) in the discretized version of equation (9.73) also.

For efficient sampling of the configuration space the basic time-step∆tshould be chosen as large as possible, but small enough so that the forces remain essentially constant. Typical choices are in the range∆t/τ_B≈0.5·10⁻⁴−0.5·10⁻². Simulations of dense systems with steeply repulsive pair interactions require smaller time-steps, whereas longer time-steps can be used for less strongly correlated systems.

37In the program this integration method will be used as default for given compiler directiveBD OVERDAMPEDduring the compilation process. Other integration methods may be specified by further compiler directives (cf. footnote41).

38In index notation:

In the simulation for all particles and each coordinate direction we can draw Gaussian random numbers from a single Gaussian distribution, because the correlation contains the productδ_{i j}δ_αβand thus no cross correlations do exist.

From equation (9.73) we immediately get the mean square displacement in a time step h∆r²_iαi=6D₀t+β²D²₀hF_iα²i(∆t)²+O (∆t)³

. (9.74)

For many physical realizations the following equation holds³⁹ hF_iα²i −kBT

∂²U

∂r_iα²

=0. (9.75)

Thus, comparison of equation (9.74) with equation (9.71) shows, that the CBD algorithm always overestimates the exact result by 2β²D²₀hF_iα²i(∆t)², which is positive irrespective of the form of the interaction potential. Due to thehF_iα²ifactor the error depends on the density and the steepness of the interaction potential. [Hey00]

The Euler integration scheme is the simplest choice, but besides the accuracy it is not very sta-ble [Pre92]. This is due to the fact that the formula (9.73) is not symmetric, i.e. it only uses derivative information at the beginning of the time step.

Usually the Euler scheme gives good numerical results when the drift and the diffusion coefficients are constant, as it is the case for our situation.⁴⁰

9.2.2 Stochastic Runge-Kutta Method

Bra´nka and Heyes [Bra99,Hey00] proposed an extension of the conventional Brownian dynamics algorithm of Ermak [Erm75], which gives more accurate results. The integration of equation (9.65) is carried out using a second-orderstochastic Runge-Kutta (SRK) algorithm, which updates the particle positions according to

ri(t+∆t) = ri(t) +∆ri(∆t)

= r_i(t) +βD₀ 2

F^a_i +F^b_i

∆t+W_i, (9.76)

where the forces are calculated at two points in time: F^a_i ≡F(r^N)at timet, and then F^b_i(R^N)at timet+∆twith

Ri=ri+βD0Fi∆t+Wi. (9.77)

This update scheme is also known asHeun predictor-corrector methodwhere the CBD displace-ment is used as the “predictor” step. TheWiin equations (9.76) and (9.77) are identical.⁴¹

39Unlike in molecular dynamics simulations, there exists no conserved quantity in a BD algorithm that can be used to check the correctness of the time stepping algorithm. But for many physical systems the relation (9.75) can serve as a cross check of the code.

40The more general Euler algorithm, which corresponds to the Itˆo stochastic differential equation dr(t) =A(t,r(t))dt+B(t,r(t))·dW(t),

reads

r(t_j+1) =r(t_j) +A(t_j,r(t_j))(t_j+1−t_j) +B(t_j,r(t_j))·(W_j+1−W_j)

forj=0,1,2...,N_step. It can be shown to converge strongly with order 1/2 (instead of order 1 as for the determin-istic case) [¨Ot96]. The source of the convergence order 1/2 are the configuration dependent coefficientsB, which in the Euler scheme are only evaluated at the initial configuration. Thus higher order of convergence can be expected for additive noise.

41The stochastic Runge-Kutta method will be used in the executable simulation program, when the compiler directives BD OVERDAMPEDandHEYES SRKare both set.

This predictor-corrector method can be shown to be weakly second-order convergent [¨Ot96,Klo99].

This means that, ifr^exact_i is the exact trajectory, then the error of any polynomial function ofr_iis bounded by

g

r^exact_i (t)

−g[ri(t)]≤c(∆t)², (9.78) wherecis a positive constant, which does not depend on∆t.

The mean square displacement in a time step is up to the second order in ∆t identical with the exact result of equation (9.71). Therefore the leading error term is of orderO((∆t)³)compared to O((∆t)²)of the CBD algorithm.

9.2.3 Smart Monte-Carlo Method

The Smart Monte-Carlo methodcan be considered as an extension of the CBD method [Ros78, Hey98,Hey00,Jar99,All87]. To ensure that the state distribution at each Monte-Carlo-time step conforms to the Boltzmann distribution

Pk= 1

Ze^−βUk, (9.79)

a Boltzmann weighted importance sampling is introduced into the generation of tentative trajecto-ries. HereZis the partition function andU_kis the potential energy of thekth configuration.

The transition probability from statekto statel,

Pkl=TklAkl, (9.80)

can be written as the product of the acceptance probabilityA_kltimes the state-to-state displacement transition vectorT_kl. To ensure the condition of detailed balanceA_klis chosen as

A_kl=min

Trajectories are generated using the algorithm of Ermak where now the particle indices are re-placed by a state indexkorl,i.e.the trial new statelis calculated from statekusing

R_l=R_k+∆R_k (9.82)

and

∆Rk=βD0Fk∆t+Wk. (9.83)

These equations of motion result in the underlying stochastic matrix T_kl=

which is needed for the evaluation ofAkl. Finally, in a Metropolis Monte-Carlo scheme one decides whether the new trial stateR_l gets accepted.

The smart Monte-Carlo method can be used to deal with hard spheres and particles with steep short-range repulsions [Cic90].

9.2.4 Gaussian Random Numbers

Any quantityAis called aGaussian random numberif its probability distributionρ(a), defined by

The generation of Gaussian random numbers is based on uniform distributed random numbers of the interval(0,1). As a direct application of the central limit theorem Gaussian random numbers can be generated by simple addition of enough uniformly distributed random numberszi∈(0,1):

r12

For most applicationsn=12 gives good results. But for the generation of every Gaussian random number one has to draw 12 uniformly distributed random numbers.⁴² In our simulation we use the more efficientBox-Muller method[Box58,Pre92,All87], which works as follows. Two uniform deviates on(0,1),x1,x2are transformed to the two quantitiesy1andy2according

y₁ = p

Evaluation of the general transformation law for the joint probability distributions [Pre92]

p(y₁,y₂,...) =p(x₁,x₂,...)

42On the level of the first-order Ermak algorithm only the first and second moments of the Gaussian distribution are needed, which can equally well obtained by using a random number constructed from the uniformly distributed random numberzaccording to

This random numberRhas the same first two moments as a Gaussian random number. [Gre88]

So the joint probability distribution is a product of two independent normal distributions and thus eachyis independently distributed according to the normal distribution (9.86), which is centered around zero and has the varianceσ=1.0.

The efficiency of the Box-Muller transformation (9.88) can be increased by applying the follow-ing trick [Pre92]. Instead of drawfollow-ing uniform deviates x₁ and x₂ in the unit square, we pick a random point with coordinates v₁ andv₂ inside the unit circle around the origin. The radial coordinate R² =v²₁+v²₂ is a uniform deviate, which can be used for x1, while the polar angle φ=arctan(v₂/v₁)can serve as the random angle 2πx₂. This allows us to avoid the numeric evalu-ation of the trigonometric function calls in (9.88), because the cosine and sine can now be written asv1/√

R²andv2/√

R²respectively.

This algorithm is implemented in the template classRandomGauss⁴³, which can be used in com-bination with any uniform distributed random number generator, for example the random number generators of the “GNU Scientific Library” [Gal06] or of the “blitz++ library” for scientific com-puting [Vel05].