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

% revision history:
% 7/22/2010 Mark D. Shattuck <mds> md.m
%           MD Demo for HandsOn 2010
%           000 mds Initial conditions and visualization
% 7/24/2010 001 mds Add Euler
% 7/24/2010 002 mds Add Interaction detection and Force Law
% 7/24/2010 003 mds Add KE, PE, Nplotskip, and plotit
% 1/21/2011 010 mds Add Velocity Verlet

%% Experimental Parmeters
N=10;  % number of particles
D=2;   % diameter of particles
K=100; % spring constant for harmonic force law
M=3;   % mass of particles

Lx=10*D;  % size of box
Ly=10*D;

TT=100; % total simulation time

%% Physical Parameters
g=0;

%% 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

%% Initial Conditions
[x y]=ndgrid(D/2:D:Lx-D/2,D/2:D:Ly-D/2);
[junk ii]=sort(rand(1,numel(x)));
x=x(ii(1:N));
y=y(ii(1:N));

vx=randn(1,N)/3;
vy=randn(1,N)/3;

ax_old=0*x;    % inital condition Fnm=0;
ay_old=0*y-g;

%% Output
Ek=zeros(Nt,1);    % Kinetic Energy
Ep=zeros(Nt,1);    % particle-particle potential
Ewp=zeros(Nt,4);   % wall-particle potential (1234)=>(LBRT)

%% Setup Plotting
if(plotit)
  clf;
  h=zeros(1,N);
  for np=1:N
    h(np)=rectangle('Position',[x(np)-.5*D y(np)-.5*D D D],'Curvature',[1 1],'edgecolor','b');
  end
  axis('equal');
  axis([0 Lx 0 Ly]);
end

%% Main Loop

for nt=1:Nt
  
  % plot particles
  if(plotit && rem(nt-1,Nplotskip)==0)
    for np=1:N
      set(h(np),'Position',[x(np)-.5*D y(np)-.5*D D D]);
    end
    drawnow;
  end
  
  x=x+vx*dt+ax_old.*dt.^2/2;  % first step in Verlet integration
  y=y+vy*dt+ay_old.*dt.^2/2;
  
  % Interaction detector and Force Law
  Fx=zeros(1,N);
  Fy=zeros(1,N);
  
  for nn=1:N
    for mm=nn+1:N
      dy=y(mm)-y(nn);
      if(abs(dy)<=D)      % check dy for speedup
        dx=x(mm)-x(nn);
        dnm=dx.^2+dy.^2;
        if(dnm<D^2)       % do not take sqrt unless needed
          dnm=sqrt(dnm);
          F=-K*(D/dnm-1);
          Ep(nt)=Ep(nt)+(D-dnm).^2;  % particle-particle PE
          Fx(nn)=Fx(nn)+F.*dx;       % particle-particle Force Law
          Fx(mm)=Fx(mm)-F.*dx;
          Fy(nn)=Fy(nn)+F.*dy;       % Newton's 3rd Law
          Fy(mm)=Fy(mm)-F.*dy;
        end
      end
    end
  end
  
  Ep(nt)=K/2*Ep(nt); % scale PE
  
  ii=x<D/2;
  dw=x(ii)-D/2; % Left wall
  Fx(ii)=Fx(ii)-K*dw;
  Ewp(nt,1)=K*sum(dw.^2)/2; %PE
  
  ii=y<D/2;
  dw=y(ii)-D/2;  % Bottom wall
  Fy(ii)=Fy(ii)-K*dw;
  Ewp(nt,2)=K*sum(dw.^2)/2; %PE
  
  ii=x>Lx-D/2;
  dw=x(ii)-(Lx-D/2);  % Right wall
  Fx(ii)=Fx(ii)-K*dw;
  Ewp(nt,3)=K*sum(dw.^2)/2; %PE
  
  ii=y>Ly-D/2;
  dw=y(ii)-(Ly-D/2);  % Top wall
  Fy(ii)=Fy(ii)-K*dw;
  Ewp(nt,4)=K*sum(dw.^2)/2; %PE
  
  ax=Fx./M;       % Newton's 2nd Law
  ay=Fy./M-g;
  
  vx=vx+(ax_old+ax).*dt/2;  % second step in Verlet integration
  vy=vy+(ay_old+ay).*dt/2;
  
  % Kinetic energy
  Ek(nt)=M*sum((vx.^2+vy.^2))/2;
  
  ax_old=ax;
  ay_old=ay;
end
