function [Ek,Ep,vxs,vys,xs,ys]=verletMDpbc(N,KEset,TT,plotit,D)
% verletMDpbc <Velocity Verlet MD in periodic boundary conditions.>
% Usage:: [Ek,Ep,vxs,vys,xs,ys]=verletMDpbc(N,KEset,plotit[true],D[2])
%

%% 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
% 04/18/2020 000 mds Add comments clean up for physics 471 introMDpbc.m
% 04/22/2020 000 mds euler version
% 04/22/2020 000 mds velocity verlet version
% 04/23/20 Mark D. Shattuck <mds> eulerMDpbc.m
%            000 mds function verstion  

%% Parse Input
if(~exist('D','var') || isempty(D))
	D=2; % diameter of particles
end
if(~exist('TT','var') || isempty(TT))
	TT=100; % Total simulation time
end
if(~exist('plotit','var') || isempty(plotit))
	plotit=true; % plot ?
end

%% Experimental Parameters
% N=10;  % number of particles   (passed in)
% D=2;   % diameter of particles (passed in)
K=100; % force vector n<-m Fnm=-K*(D/|dnm|-1)*dnm, where vect dnm=Xm-Xn
M=3;   % mass of particles

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

% TT=100;   % total simulation time (passed in)

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

%% Simulation Parmeters
dt=1e-2;       % integration time step
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); % place particles on grid 
ii=randperm(numel(x),N);                 % N random position on grid
x=x(ii);                                 % set x-positions
y=y(ii);                                 % set y-positions

vx=randn(1,N);  % normal (maxwell) distribution of velocities
vy=randn(1,N);

vx=vx-mean(vx);   % remove center of mass drift
vy=vy-mean(vy);

tmp=sqrt(2*KEset/M/sum((vx.^2+vy.^2)));  % set Kinetic energy
vx=vx*tmp;
vy=vy*tmp;

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

%% Save variables
Ek=zeros(Nt,1);    % Kinetic Energy
Ep=zeros(Nt,1);    % particle-particle potential
xs=zeros(Nt,N);    % x-position
ys=zeros(Nt,N);    % y-position
vxs=zeros(Nt,N);   % x-velocity
vys=zeros(Nt,N);   % y-velocity

%% Setup Plotting
% NB: periodic images not shown
if(plotit)                    % only if plotit true
  clf;                        % clear figure
  h=plotNCirc(x,y,D,'none');  % plot particles no color; store handle h
  axis('equal');              % square pixels
  axis([0 Lx 0 Ly]);          % ploting range
end

%% Main Loop
for nt=1:Nt
  
  % plot particles
  if(plotit && rem(nt-1,Nplotskip+1)==0)  % nt-1 divides Nplotskip+1
    xp=mod(x,Lx);
    yp=mod(y,Ly);
    for np=1:N
      set(h(np),'Position',[xp(np)-.5*D yp(np)-.5*D D D]); % update handle 
    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;
  
  % position dependent calculations
  xs(nt,:)=x;                         % save positions
  ys(nt,:)=y;

  % Force between particles
  Fx=zeros(1,N);                      % zero forces
  Fy=zeros(1,N);
  
  for nn=1:N                          % check particle n 
    for mm=nn+1:N                     % againt particle m from n+1->N
      dy=y(mm)-y(nn);                 % y-comp of dnm vector from from n->m
      dy=dy-round(dy/Ly)*Ly;          % closest periodic image
      if(abs(dy)<D)                   % y-distance close enough?    
        dx=x(mm)-x(nn);               % x-comp of dnm vector from from n->m   
        dx=dx-round(dx/Lx)*Lx;        % closest periodic image
        dnm=dx.^2+dy.^2;              % squared distance between n and m
        if(dnm<D^2)                   % overlaping?
          dnm=sqrt(dnm);              % distance
          F=-K*(D/dnm-1);             % force magnitude 
          Fx(nn)=Fx(nn)+F.*dx;        % accumulate x force on n
          Fx(mm)=Fx(mm)-F.*dx;        % x force on m (equal opposite) 
          Fy(nn)=Fy(nn)+F.*dy;        % accumulate y force on n
          Fy(mm)=Fy(mm)-F.*dy;        % y force on m (equal opposite)
          Ep(nt)=Ep(nt)+(D-dnm).^2;   % particle particle potential calc I
        end
      end
    end
  end
  Ep(nt)=K/2*Ep(nt);                  % particle particle potential calc II 
  
  ax=Fx./M;                 % calc a from F=Ma
  ay=Fy./M;
  
  vx=vx+(ax_old+ax)*dt/2;   % second step in Verlet integration
  vy=vy+(ay_old+ay)*dt/2;
   
  % velocity dependent calculations
  Ek(nt)=M*sum((vx.^2+vy.^2))/2; % Kinetic energy
  vxs(nt,:)=vx;                  % save velocities
  vys(nt,:)=vy;
  
  ax_old=ax;                % save for next step
  ay_old=ay;
end




