### 82 Computer simulation of entropy production 熵產生的電腦模擬

In the previous post we described how to derive the entropy production rate by Fokker-Planck equation. I think this may be a little bit abstract to our readers, so I decide to make this small supplement.

### Today's problem

Now let's consider a simple system:
a localized 3-dimensional particle. Its location is initially confined by an optical tweezer, so the probability distribution of its location follows a normal distribution in x direction, expressed in terms of Dirac delta series with n = 10^10. Its motion in y & z directions are completely restricted and ignorable. The particle has a radius of 1nm and is immersed in water at 350K. Now there is an engineering accident near by the lab so there is a power outage, and the optical tweezer shuts down. The particle will diffuse out and its entropy will increase. What is the entropy production rate of the system when the electricity just goes out (a.k.a. the initial entropy production rate)?

### Theoretical derivation

Okay this sounds a little bit complex and difficult. Please calm down and let's start from what we know. Since this is an overdamped system without no extra forces or potential, it's super-easy to write down its Langevin equation:
Its associated Fokker-Planck equation and probability flux would be:
According to our previous episode, the entropy production rate could be expressed as
If we substitute in the initial probability distribution (the Dirac delta series, this is nothing more than some Gaussian integrals), we will get:
Perfect, isn't it?

### Computer simulation -- Do they match!?

Now let's see does our computer simulation matches our prediction. Since this simulation involves lots of matrix manipulation, we will use MATLAB. The code looks something like this:

%% constants
tic;
n_list = 5*10^9:10^9:2*10^10; %n of delta series
nSim = 50000; dt = 3*10^(-13); nRun = 400;
kB = 1.38*10^(-23); T = 350;
r = 1*10^(-9); rho = 10^3; eta = 10^(-3); zeta = 6*pi*eta*r;
D = kB*T/zeta;%diffusion constant
A = sqrt(2*zeta*kB*T*dt);

%% simulation
S_record = zeros(length(n_list),nRun);
for j =1:length(n_list)
n = n_list(j);
x_data = zeros(nSim, nRun);
x_data(:,1) = randn(nSim,1)/(sqrt(2)*n);
S_record(j,1) = kB*EntropyCal(x_data(:,1));
for i = 2:nRun
x_data(:,i) = x_data(:,i-1) + 1/zeta*A*randn(nSim,1);
S_record(j,i) = kB*EntropyCal(x_data(:,i));
end
display(j);
end

%% entropy production rate
figure, plot(0:dt:(nRun-1)*dt,S_record(5,:));
xlabel('t','FontSize',14);ylabel('S (J/K)','FontSize',14);
S_prod_exp_initial = (S_record(:,2)-S_record(:,1))/dt;
S_prod_theo_initial = 2*kB*D*n_list.^2;
display([S_prod_exp_initial]); display([S_prod_theo_initial]);
figure, plot(n_list, S_prod_theo_initial);
hold on;
scatter(n_list, S_prod_exp_initial);
xlabel('n', 'FontSize', 14);
ylabel('dS(t=0)/dt (J/K*sec)','FontSize',14);
toc;

function S = EntropyCal(record)
[~, nLength] = size(record);
S = zeros(1, nLength);
for i = 1:1:nLength
h = histogram(record(:,i), 25,'Normalization','pdf');
pdf = h.Values; d = h.BinWidth; close all;
S(i) = - nansum(pdf.*log(pdf))*d;
end
end


It is basically the same as the code we introduced before. We run 50000 particles at the same time and use this to generate the evolution of probability distribution. We then use this probability distribution to monitor the change in entropy. The time evolution of entropy would look like this:
The curve is relatively smooth because we use many particles in parallel and have the simulation time interval small enough. We also change the n of Dirac delta series to adjust how much the particle was originally confined, and we see how much our predictions and experimental results match. Here is the results:
They match quite well in a certain range. However, since the accuracy of the experiment is highly dependent on the time scale of simulation, the number of simulations we run in parallel and the algorithm how we calculate entropy, certain errors are inevitable in this simple simulation program. (There are some ways to eliminate them though.)

### One extra problem if you like it

If the above problem is not exciting enough, try calculate the initial entropy production rate of the following system:
2 identically charged particles are confined in a potential well
and their locations are further restricted by 2 optical tweezers at -x0 & +x0, with probability distributions expressed in terms of Gaussian Dirac delta series with n = 10^10. Their motions in y & z directions are completely restricted and ignorable. The particles have a radius of 1nm and are immersed in water at 350K. Now there is a power outage again, and the optical tweezers shuts down simultaneously. Let k be the Coulomb's constant and each particle is charged with +q. What is the entropy production rate of the system when the electricity just goes out?

If I derive the answer correctly, it would be
Please tell me if I did it wrong XD. Many thanks.

In our next episode, we will move on to energy trade-off in sensory adaptation! Stay tuned!