#72 Computer simulation of Langevin equation 朗之萬方程的電腦模擬

Today we will discuss the computer simulation of Langevin equation, implement it to a model system oscillating between 2 potential well, and verify the Arrhenius law for reaction rate constant.

Let's start from the overdamped Langevin equation we mentioned previously:

The equation could be written in a more general form:

in which $\eta(t)$ represents the white noise with magnitude 1:

Now let's consider the impulse it produces within a small time interval

The mean of this tiny impulse is obviously 0. Now let's consider its variance. From the definition of variance, since the mean of the total impulse is 0, we will get

Using the properties of delta function, we will get

which means the impulse in a short time interval is a Gaussian noise. This allows us to write the above stochastic differential equation in discrete form:

which could be readily simulated by any programming language.

Computer simulation in MATLAB   MATLAB電腦模擬

Now we will show how to implement the above discussion in computer simulation. Assume we have a system with a potential well

The system should oscillate between the 2 potential well. The activation energy is V0. By recording the oscillation pattern under different temperature, we could calculate the rate constant under different temperature and verify the Arrhenius law.

A test program which simulates the system under temperature 150K would looks like this:

%% basic setting

r = 5E-9; eta = 8.90E-4; zeta = 6*pi*eta*r;

kB = 1.38E-23;

syms x;

x0 = 1E-9;

V0 = 1E-20;

V = V0*(1-(x/x0).^2).^2;

F = matlabFunction(diff(V,x));

V_fun = matlabFunction(V);

%% test simulation

n_sim = 1E7;

x_total = zeros(1, n_sim);

x_total(1)=rand(1,1)*(1E-9);

T =150; A = sqrt(2*zeta*kB*T*dt);

for i = 2:1:n_sim

x_old = x_total(i-1);

x_total(i) = x_old + 1/zeta*(-dt*F(x_old)+A*randn(1,1));

end

figure;

plot(dt:dt:(n_sim*dt),x_total);

xlabel('time(s)');

ylabel('reaction coordinate');

This would yield a figure like this:

By recording the dwelling time in 1 potential well, we could calculate the reaction rate constant easily. However, to reliably calculate the dwelling time, we must remove some noise from the data.

HPD = 5000;
coeffMA = ones(1, HPD)/HPD;
avg24hx_total = filter(coeff24hMA, 1, x_total);
avg_total = filter(coeffMA, 1, x_total);
x_01 = (avg_total>0);
cal_temp = x_01(2:end) - x_01(1:(end-1));
trans_time_pt = find(cal_temp);
duration = (trans_time_pt(2:end) - trans_time_pt(1:(end-1)))*dt;
k_est = length(duration)/sum(duration);

The reaction rate constant estimated is 1.121E+06/s. Now let's try to simulate under different temperature:

k_est_total = zeros(1,7);
HPD = 5000;
coeffMA = ones(1, HPD)/HPD;
for T = 900:300:3000
n_sim = 2E7;
n_rec_per = 1;
n_rec = n_sim/n_rec_per;
x_total = zeros(1, n_rec);
x_total(1)=rand(1,1)*(1E-9);
x1 = x_total(1);
A = sqrt(2*zeta*kB*T*dt);
for i = 1:1:n_sim
x_old = x1;
x1 = x_old + 1/zeta*(-dt*F(x_old)+A*randn(1,1));
if mod(i,n_rec_per)==0
x_total(i/n_rec_per + 1) = x1;
end
end
avg_total = filter(coeffMA, 1, x_total);
x_01 = (avg_total>0);
cal_temp = x_01(2:end) - x_01(1:(end-1));
trans_time_pt = find(cal_temp);
duration = (trans_time_pt(2:end) - trans_time_pt(1:(end-1)))*n_rec_per*dt;
k_est_total((T-600)/300) = length(duration)/sum(duration);
display(T);
end

The above program would take about 30 minutes to run. After a long wait, we finally get our reaction rate constant under different temperature. The activation energy could be estimated by linear regression:

log_k_est_total = log(k_est_total');
T_total = 1./(900:300:3000);
T_total = T_total';
X = [ones(8,1) T_total ];
format long
b1 = X\log_k_est_total;
yCalc1 = X*b1;
scatter(900:300:3000,k_est_total);
hold on
plot(900:300:3000,exp(yCalc1),'LineWidth',2);
set(gca,'FontSize',14)
xlabel('Temperature(K)','FontWeight','bold','FontSize',14)
ylabel('reaction rate constant, k(1/s)','FontWeight','bold','FontSize',14)
grid on
display(b1(2)*kB);

And the result looks like this:

The result is compatible to Arrhenius law and the activation energy is estimated to be 1.267E-20, which is close to our setting 1.0E-20. Isn't it amazing?