function Figure8_5

% Matlab code for Figures 8.5 in book:
% Radek Erban and Jon Chapman, Stochastic Modelling of Reaction-Diffusion
% Processes, Cambridge University Press, 2019

close all;

rr=[0:0.01:4];
lj=4*((1./rr).^12-(1./rr).^6);
ljp=4*(-12*(1./rr).^13+6*(1./rr).^7);
E=-1/3;
cec=E-((E+1)/2)*2^(1/3)./(rr.*rr);

figure(1);
set(gca,'Fontsize',18);
plot(rr,lj,'b','Linewidth',4);
hold on;
plot(rr,0*rr,'k','Linewidth',1);
plot([2^(1/6) 2^(1/6)],[-1 0],'k:');
plot([0 2^(1/6)],[-1 -1],'k:');
plot([0 3.1],[E E],'k:');
plot(rr,cec,'g','Linewidth',2);
plot(rr,lj,'b','Linewidth',4);
axis([0 3.1 -1.1 1.1]);
xlabel('$r$','interpreter','latex');
ylabel('$\Phi_{\mathrm L\mathrm J}(r)$','interpreter','latex');
set(gca,'Xtick',[0 0.5 1 1.5 2 2.5 3]);
text(1.07,0.09,'$2^{1/6}\sigma$','fontsize',18,'interpreter','latex');
text(0.856,-0.65,'$r_{1}$','fontsize',18,'interpreter','latex');
text(1.375,-0.625,'$r_{2}$','fontsize',18,'interpreter','latex');
text(-0.15,1,'$\varepsilon$','fontsize',18,'interpreter','latex');
text(-0.22,-1,'-$\varepsilon$','fontsize',18,'interpreter','latex');
text(-0.23,E,'$E$','fontsize',18,'interpreter','latex');
text(0.96,-1.17,'$\sigma$','fontsize',18,'interpreter','latex');
text(1.93,-1.2,'$2\sigma$','fontsize',18,'interpreter','latex');
text(2.93,-1.2,'$3\sigma$','fontsize',18,'interpreter','latex');
set(gca,'XtickLabel',{'','','','','','',''});
set(gca,'Ytick',[-1 -0.5 0 0.5 1]);
set(gca,'YtickLabel',{'','','0','',''});
set(gca,'Fontsize',18);

figure(2);
set(gca,'fontsize',18);  
plot(rr,-ljp,'b','Linewidth',4);
hold on;
plot(rr,0*rr,'k','Linewidth',1);
axis([0.7 2.1 -2.5 2.5]);
xlabel('$r$','interpreter','latex');
ylabel('$-\Phi_{\mathrm L\mathrm J}^\prime(r)$','interpreter','latex');
set(gca,'Xtick',[0 0.5 1 1.5 2]);
text(1.05,-2.75,'$2^{1/6}\sigma$','fontsize',18,'interpreter','latex');
text(1.15,0.2,'$(26/7)^{1/6}\sigma$','fontsize',18,'interpreter','latex');
text(0.64,1,'$\frac{\varepsilon}{\sigma}$','fontsize',18,'interpreter','latex');
text(0.616,-1,'-$\frac{\varepsilon}{\sigma}$','fontsize',18,'interpreter','latex');
text(0.6035,2,'$\frac{2\varepsilon}{\sigma}$','fontsize',18,'interpreter','latex');
text(0.58,-2,'-$\frac{2\varepsilon}{\sigma}$','fontsize',18,'interpreter','latex');
text(0.98,-2.68,'$\sigma$','fontsize',18,'interpreter','latex');
text(1.965,-2.7,'$2\sigma$','fontsize',18,'interpreter','latex');
set(gca,'XtickLabel',{'','','','',''});
set(gca,'Ytick',[-2 -1.5 -1 -0.5 0 0.5 1 1.5 2]);
set(gca,'YtickLabel',{'','','','','0','','','',''});
plot([2^(1/6) 2^(1/6)],[-3 0],'k--');
plot([(26/7)^(1/6) (26/7)^(1/6)],[-(4*126/169)*((26/7)^(-1/6)) 0],'k--');
set(gca,'fontsize',18);

fs=50;
epsilon=1;
E=-1/50;
sigma=1;
inpos=(2^(1/6))/2;
inimp=sqrt((1+E)/2);
dt=0.001;
finaltime=1000;
finaltimedis=round(finaltime/dt);

initpos=[-inpos 0 0 inpos 0 0];
initvel=[inimp inimp 0 -inimp -inimp 0];
time=zeros(finaltimedis+1,1);
pos=zeros(finaltimedis+1,6);
vel=zeros(finaltimedis+1,6);
pos(1,:)=initpos;
vel(1,:)=initvel;

ALJ=4*epsilon*sigma^12;
BLJ=4*epsilon*sigma^6;

t=1;
sqdist=(pos(t,1)-pos(t,4))*(pos(t,1)-pos(t,4))+(pos(t,2)-pos(t,5))*(pos(t,2)-pos(t,5))+(pos(t,3)-pos(t,6))*(pos(t,3)-pos(t,6));
propconst=0.5*dt*(-(12*ALJ/(sqdist^7))+(6*BLJ/(sqdist^4)));
frc=zeros(1,6);
frc(1)=propconst*(pos(t,4)-pos(t,1));
frc(2)=propconst*(pos(t,5)-pos(t,2));
frc(3)=propconst*(pos(t,6)-pos(t,3));
frc(4)=propconst*(pos(t,1)-pos(t,4));
frc(5)=propconst*(pos(t,2)-pos(t,5));
frc(6)=propconst*(pos(t,3)-pos(t,6));

for t=1:finaltimedis
    time(t+1)=dt*t;
    vel(t+1,:)=vel(t,:)+frc(1,:);
    pos(t+1,:)=pos(t,:)+vel(t+1,:)*dt;
    sqdist=(pos(t+1,1)-pos(t+1,4))*(pos(t+1,1)-pos(t+1,4))+(pos(t+1,2)-pos(t+1,5))*(pos(t+1,2)-pos(t+1,5))+(pos(t+1,3)-pos(t+1,6))*(pos(t+1,3)-pos(t+1,6));
    propconst=0.5*dt*(-(12*ALJ/(sqdist^7))+(6*BLJ/(sqdist^4)));
    frc=zeros(1,6);
    frc(1)=propconst*(pos(t+1,4)-pos(t+1,1));
    frc(2)=propconst*(pos(t+1,5)-pos(t+1,2));
    frc(3)=propconst*(pos(t+1,6)-pos(t+1,3));
    frc(4)=propconst*(pos(t+1,1)-pos(t+1,4));
    frc(5)=propconst*(pos(t+1,2)-pos(t+1,5));
    frc(6)=propconst*(pos(t+1,3)-pos(t+1,6));
    vel(t+1,:)=vel(t+1,:)+frc(1,:);
end   
    
figure(3);
set(gca,'fontsize',18);
plot(pos(1:7800,1),pos(1:7800,2),'r','Linewidth',3);
hold on;
plot(pos(1:7800,4),pos(1:7800,5),'b','Linewidth',3);
plot(pos(7800,1),pos(7800,2),'or','Markersize',8,'Linewidth',10);
plot(pos(7800,4),pos(7800,5),'ob','Markersize',8,'Linewidth',10);
plot(pos(1,1),pos(1,2),'xr','Markersize',12,'Linewidth',4);
plot(pos(1,4),pos(1,5),'xb','Markersize',12,'Linewidth',4);
set(gca,'Xtick',[-1 -0.5 0 0.5 1]);
set(gca,'XtickLabel',{'','','','',''});
text(-1.07,-1.1,'$-\sigma$','fontsize',18,'interpreter','latex');
text(-0.66,-1.1,'$-\sigma/2$','fontsize',18,'interpreter','latex');
text(-0.02,-1.1,'$0$','fontsize',18,'interpreter','latex');
text(0.42,-1.1,'$\sigma/2$','fontsize',18,'interpreter','latex');
text(0.96,-1.1,'$\sigma$','fontsize',18,'interpreter','latex');
text(-0.02,-1.2,'$x$','fontsize',18,'interpreter','latex');
set(gca,'Ytick',[-1 -0.5 0 0.5 1]);
set(gca,'YtickLabel',{'','','','',''});
text(-1.2,-0.98,'$-\sigma$','fontsize',18,'interpreter','latex');
text(-1.2,-0.5,'$-\frac{\sigma}{2}$','fontsize',18,'interpreter','latex');
text(-1.1,0,'$0$','fontsize',18,'interpreter','latex');
text(-1.1,0.52,'$\frac{\sigma}{2}$','fontsize',18,'interpreter','latex');
text(-1.1,1.02,'$\sigma$','fontsize',18,'interpreter','latex');
text(-1.22,0.3,'$y$','fontsize',18,'interpreter','latex');
axis([-1 1 -1 1]);
set(gca,'fontsize',18);

figure(4);
set(gca,'fontsize',18);
plot(pos(:,4),pos(:,5),'b','Linewidth',3);
box on;
set(gca,'Xtick',[-1 -0.5 0 0.5 1]);
set(gca,'XtickLabel',{'','','','',''});
text(-1.07,-1.1,'$-\sigma$','fontsize',18,'interpreter','latex');
text(-0.66,-1.1,'$-\sigma/2$','fontsize',18,'interpreter','latex');
text(-0.02,-1.1,'$0$','fontsize',18,'interpreter','latex');
text(0.42,-1.1,'$\sigma/2$','fontsize',18,'interpreter','latex');
text(0.96,-1.1,'$\sigma$','fontsize',18,'interpreter','latex');
text(-0.02,-1.2,'$x$','fontsize',18,'interpreter','latex');
set(gca,'Ytick',[-1 -0.5 0 0.5 1]);
set(gca,'YtickLabel',{'','','','',''});
text(-1.2,-0.98,'$-\sigma$','fontsize',18,'interpreter','latex');
text(-1.2,-0.5,'$-\frac{\sigma}{2}$','fontsize',18,'interpreter','latex');
text(-1.1,0,'$0$','fontsize',18,'interpreter','latex');
text(-1.1,0.52,'$\frac{\sigma}{2}$','fontsize',18,'interpreter','latex');
text(-1.1,1.02,'$\sigma$','fontsize',18,'interpreter','latex');
text(-1.22,0.3,'$y$','fontsize',18,'interpreter','latex');
axis([-1 1 -1 1]);
set(gca,'fontsize',18);

