How to include random force in the simulation (Classical Molecular Dynamics)

Question

I need to implement a random force in my code according to the fluctuation dissipation theorem. I have a Gaussian distribution function ready width average 0 and distribution 1 and I know I need to multiply it by something but I'm not sure what. The fluctuation dissipation theorem gives (in one dimension):

$\left \langle A (t_1) A (t_2) \right \rangle = 2 m \gamma k_B T \delta (t_1 - t_2 )$

I can't convince my self if I should multiply my Gaussian with my $\sqrt{dt}$ or divide by it. I can find justification for both direction:

1. The force should act like delta function (no correlation between times) hence, I should divide by it.
2. The units are ok when multiplying by $\sqrt{dt}$.
3. The larger the time step, the more collision happened hence the force should be larger.

I remember that in some place I read or heared I should divide by the $\sqrt{dt}$ but I couldn't find where.

I should note that I tried both approach, when I divide I suspect the forces to be too large (I only see noise) and when I multiply, I almost can't see any random effect.

Answer 1

I just happened to simulate the exact problem. Firstly, this falls under the topic of stochastic differential equations, as it will help us choose the right mathematics to work with.

For every physical process with dissipation, we can in-general associate a noise force to explain the physics of diffusion and noise dynamics. This is captured under the fluctuation dissipation theorem. There is a more deeper reason for this when treated at the quantum mechanical level (to maintain commutation relationship of quantized fields under going dissipation for example in a cavity). http://en.wikipedia.org/wiki/Fluctuation-dissipation_theorem

Now, to answer your question, the square root term appears when discretizing the differential equation.We need to write down the derivative as (F(n+1)-F(n))/DT and the variance equation in discrete form by replacing the continuous time delta function (which has units of inverse time) with Kronecker delta function divided by the time step. The variance of the noise for fixed time step then has a sqrt(DT) in the denominator and will produce the form as you show in the question.

Note: that the delta function has units of inverse time. so there is no discrepancy in units.

There are a few key points to keep in mind.

1. use a first order integration method such as forward euler or midpoint method. http://en.wikipedia.org/wiki/Backward_Euler_method#The_backward_Euler_method . We can easily show that, higher order methods converge to wrong noise variances. In general, stochastic differential equations do not admit the higher order integration methods like higher order runge kutta. see Phys. Rev. E 74, 068701 (2006) Comment on “Numerical methods for stochastic differential equations” http://pre.aps.org/abstract/PRE/v74/i6/e068701

2. If the noise is additive any first order integration will give a correct answer. But, if the noise is multiplicative, appropriate calculus must be employed (Ito or stratonovich). http://en.wikipedia.org/wiki/Stratonovich_integral http://en.wikipedia.org/wiki/Ito_Calculus

3. A quick check for the SDE will be to implement Brownian motion of a particle at temperature T and verify the average energy to be 0.5*kT.

Answer 2

Think about how your discrete time variables are being generated from continuous time variables.

If $A_i =\int_{t_i}^{t_i+dt}A(t)$ then $<A_i A_j> \propto dt * \delta_{ij}$. On the other hand if $A_i = {1 \over{dt}}\int_{t_i}^{t_i+dt}A(t)$ then $<A_i A_j> \propto {1\over{dt}} \delta_{ij}$

Answer 3

So... the delta function in time is defined as:

$\delta (t - t_0) = \int _0^\infty \delta (t) dt = \int_t^{t_0} \delta (t) dt \equiv 1$

which means that the delta function has the unit of $sec^{-1}$ hence I should divide by the time step.