%% Molecular Dynamics Simulator
% <md.m> Mark D. Shattuck 7/22/2010

% revision history:
% 07/22/2010 Mark D. Shattuck <mds> md.m
%            MD Demo for HandsOn 2010
%            000 mds Initial conditions and visualization
% 07/24/2010 001 mds Add Euler
% 07/24/2010 002 mds Add Interaction detection and Force Law
% 07/24/2010 003 mds Add KE, PE, Nplotskip, and plotit
% 01/21/2011 010 mds Add Velocity Verlet
% 03/13/2012 000 mds cell sim mdcNNN.m
% 04/22/2012 006 mds redo from mdc005
% 07/29/2016 000 mds update to current cells code
% 02/11/2019 000 mds Cell diffusion celldiff.m
% 02/08/2023 000 mds demo version dpmNNN.m

%% Experimental Parmeters
Nc=38; % number of cells
Ns=24; % number of sides

Gn=1+1*rand(1,Nc)/2; % diameter ratio
Dc=10*Gn;            % diameter of cells

L0=Dc*sin(pi/Ns);  % length of cell wall side as regular Ns-gon
D0=L0;             % Interaction distance

Alf=1.5;  % initial shape parameter P^2/(4*pi*A)
T=.001;   % "Temperature"
B=1;      % Drag coefficient;


KC=200;   % Cell-Cell interaction strength
Kp=50;    % Cell perimeter compressibility
KA=200;   % Cell Area compressibility
Kb=10;    % Bending strength

Lx=20;  % inital/current size of box
Ly=Lx;
Atot=Lx*Ly;

TT=100; % total simulation time
%% Calculated Parameters
N=Nc*Ns;
ift=[2:Ns 1];    % used to compute circular diff
jft=[Ns 1:Ns-1]; % used to compute circular diff
kft=[Ns-1 Ns 1:Ns-2]; % used to compute circular diff

%% Display Parameters
plotit=true;   % plot ?
Nplotskip=50;  % number of timesteps to skip before plotting

%% Simulation Parmeters
dt=1e-2;
Nt=fix(TT/dt); % Number of time steps
xi=sqrt(2*B*T/dt);  % Langevin random drive

%% Initial Conditions
[xc yc D Gn]=ppack(1,Ly/Lx,0,Dc+D0);    % use packing of disks for IC
dth=2*pi/Ns;                            % angle for Ns-gon
th=(0:Ns-1)'*dth.*(1+randn(Ns,1)/1000); % add some randomness

xc0=xc;   % save initial positions and sizes
yc0=yc;
Dp0=D;

d0=min(D*Gn);              % smallest diameter
Dc=D*Gn/(1+sin(pi/Ns))*Lx; % Ns-gon Diameters scale to Lx
A0=pi*Dc.^2/4;             % set area as regular Ns-gon
L0=2*sqrt(pi*Alf*A0)/Ns;   % length of cell wall side
D0=L0;                     % Interaction distance
L0=repmat(L0,Ns,1);        % Individual preferred

xc=(xc(1:Nc))*Lx;          % scale to Lx
yc=(yc(1:Nc))*Lx;

x=repmat(xc,Ns,1)+cos(th)*Dc/2;  % make each vertex
y=repmat(yc,Ns,1)+sin(th)*Dc/2;

vx=reshape(randn(Ns,Nc),Ns,Nc)/10;  % velocities for each vertex
vy=reshape(randn(Ns,Nc),Ns,Nc)/10;

vx=vx-mean(vx(:));  % move to COM
vy=vy-mean(vy(:));

ax_old=0*x;   % need for velocity verlet
ay_old=0*y;

th=2*pi*rand(1,Nc); % initial drive angle


%% Output
Ek=zeros(Nt,1);    % Kinetic Energy
Ep=zeros(Nt,1);    % cell-cell potential
Epl=zeros(Nt,1);   % total side length potential
EpA=zeros(Nt,1);   % Area potential
Epb=zeros(Nt,1);   % bending potential
As=zeros(Nt,Nc);   % area
Ps=zeros(Nt,Nc);   % perimeter

%% Setup Plotting
if(plotit)
  clf;
  xx=zeros(2*Ns,Nc);
  yy=zeros(2*Ns,Nc);
  % find outline of particle
  [xx,yy]=calcExtHull(x,y,.9*D0/2,Lx,Ly,ift,jft,xx,yy,Ns,Nc);
  pause(.01);
  cc=(jet(Nc)-.5)/2+.5;%cc=[1 .5 .5];
  xp=reshape([xx xx+Lx xx+Lx xx]-Lx/2,2*Ns,4*Nc);
  yp=reshape([yy yy yy+Ly yy+Ly]-Ly/2,2*Ns,4*Nc);
  h=patch(xp,yp,'r');
  set(h,'FaceVertexCData',repmat(cc,4,1),...
    'facecolor','flat');
  set(h,'linewidth',2,'edgecolor',0*[1 1 1]);
  
  axis('equal');
  axis([0 Lx 0 Ly]);
end

%% Main Loop
for nt=1:Nt
  % plot cells
  if(plotit && rem(nt-1,Nplotskip)==0)
    [xx,yy]=calcExtHull(x,y,.9*D0/2,Lx,Ly,ift,jft,xx,yy,Ns,Nc);
    xp=reshape([xx xx+Lx xx+Lx xx]-Lx/2,2*Ns,4*Nc);
    yp=reshape([yy yy yy+Ly yy+Ly]-Ly/2,2*Ns,4*Nc);
    set(h,'xdata',xp,'ydata',yp);
    drawnow;
  end
  if(rem(nt-1,fix(Nt/10))==0)
    fprintf('.');
  end
  
  x=x+vx*dt+ax_old.*dt.^2/2;  % first step in Verlet integration
  y=y+vy*dt+ay_old.*dt.^2/2;
  
  % Interactions
  
  % Perimeter Force
  lx=x(ift,:)-x;  % circular diff
  ly=y(ift,:)-y;
  lx=lx-Lx*round(lx/Lx); % periodic convention
  ly=ly-Ly*round(ly/Ly);
  lk=sqrt(lx.^2+ly.^2);
  P=sum(lk);
  F=-Kp*(L0./lk-1);
  Fx=F.*lx;
  Fy=F.*ly;
  Fx=Fx-Fx(jft,:);
  Fy=Fy-Fy(jft,:);
  Epl(nt)=Kp*sum((lk(:)-L0(:)).^2)/2;  % perimeter energy
  Ps(nt,:)=P;
  
  % Bending force
  Fcx=-Kb*(lx(ift,:)-lx);
  Fcy=-Kb*(ly(ift,:)-ly);
  Fx=Fx+Fcx-2*Fcx(jft,:)+Fcx(kft,:);
  Fy=Fy+Fcy-2*Fcy(jft,:)+Fcy(kft,:);
  Epb(nt)=sum(Fcx(:).^2+Fcy(:).^2)/(2*Kb); % Bending energy
  
  % Area Force
  A=sum(x.*y(ift,:)-y.*x(ift,:))/2;
  dx=x(ift,:)-x(jft,:);
  dy=y(jft,:)-y(ift,:);  % -dy
  dx=dx-Lx*round(dx/Lx); % periodic convention
  dy=dy-Ly*round(dy/Ly);
  dA=KA*(A-A0);
  for ns=1:Ns;
    Fx(ns,:)=Fx(ns,:)+dA.*dy(ns,:);
    Fy(ns,:)=Fy(ns,:)+dA.*dx(ns,:);
  end
  EpA(nt)=KA*sum((A-A0).^2);  % Area Energy
  As(nt,:)=A;
  
  % Cell-Cell Force
  for nc=1:Nc
    for mc=nc+1:Nc
      Dnm=(D0(mc)+D0(nc))/2;
      for ns=1:Ns
        for ms=1:Ns
          dx=x(ms,mc)-x(ns,nc);
          dx=dx-Lx*round(dx/Lx);
          if(abs(dx)<Dnm)
            dy=y(ms,mc)-y(ns,nc);
            dy=dy-Ly*round(dy/Ly);
            dnm=dx.^2+dy.^2;
            if(dnm<Dnm^2)
              dnm=sqrt(dnm);
              F=-KC*(Dnm/dnm-1);
              Ep(nt)=Ep(nt)+(Dnm-dnm).^2;  % cell-cell PE
              Fx(ns,nc)=Fx(ns,nc)+F.*dx;
              Fy(ns,nc)=Fy(ns,nc)+F.*dy;
              Fx(ms,mc)=Fx(ms,mc)-F.*dx; % 3rd law
              Fy(ms,mc)=Fy(ms,mc)-F.*dy;
            end
          end
        end
      end
    end
  end
  
  % intra-cell interaction
  for nc=1:Nc
    mc=nc;
    Dnm=D0(nc);
    for ns=1:Ns
      for ms=ns+2:Ns-1
        dx=x(ms,mc)-x(ns,nc);
        dx=dx-Lx*round(dx/Lx);
        if(abs(dx)<Dnm)
          dy=y(ms,mc)-y(ns,nc);
          dy=dy-Ly*round(dy/Ly);
          dnm=dx.^2+dy.^2;
          if(dnm<Dnm^2)
            dnm=sqrt(dnm);
            F=-KC*(Dnm/dnm-1);
            Ep(nt)=Ep(nt)+(Dnm-dnm).^2;  % cell-cell PE
            Fx(ns,nc)=Fx(ns,nc)+F.*dx;
            Fy(ns,nc)=Fy(ns,nc)+F.*dy;
            Fx(ms,mc)=Fx(ms,mc)-F.*dx; % 3rd law
            Fy(ms,mc)=Fy(ms,mc)-F.*dy;
          end
        end
      end
    end
  end
  
  Ep(nt)=KC/2*Ep(nt);
  
  th=th+randn(1,Nc)/20;
  
  Fx=(Fx-B*(vx+ax_old*dt/2))/(1+B*dt/2)+xi*repmat(cos(th),Ns,1);
  Fy=(Fy-B*(vy+ay_old*dt/2))/(1+B*dt/2)+xi*repmat(sin(th),Ns,1);
  
  ax=Fx;
  ay=Fy;
  
  vx=vx+(ax_old+ax).*dt/2;  % second step in Verlet integration
  vy=vy+(ay_old+ay).*dt/2;
  
  % Kinetic energy
  Ek(nt)=sum((vx(:).^2+vy(:).^2))/2;
  
  ax_old=ax;
  ay_old=ay;
end
fprintf('\n');
