Back to Gallery
  • Remix
  • Share
  • New Entry
on 1 Dec 2023
  • 9
  • 22
  • 0
  • 0
  • 1996
drawframe(1);
Write your drawframe function below
%Molecular dynamics code with:
%-periodic boundary conditions
%-neighbour list
%-velocity Verlet
function drawframe(f)
global gg hg
if f==1
rng('default')
gg.boxL=1/2;
boxL2 = gg.boxL/2;
%number of neighbours estimate (algorithmic parameters)
gg.NBestimate=10;
%particle size
r1=0.2;
r2=0.05;
N1 = 10;
N2 = 30;
gg.N = N1+N2;
%particle radius here uniform berween rmin and rmax
gg.r = [r1*ones(N1,1);r2*ones(N2,1)];
rmax = max(gg.r);
%Cutoff distance for the neighbour list
gg.CutNbList = 2*rmax*2;
%criterion for updating the neighbour list
gg.criterionNBLisdt = 0.9*(gg.CutNbList-(2*rmax));
%Parameters for pseudo dynamics
gg.parameters.mass = 0.5;
gg.parameters.dt = 2e-2;
gg.parameters.viscosity=1e-3;
%system initialization
%positions
gg.X = gg.boxL*rand(gg.N,2)-boxL2;
%velocities
gg.V = zeros(gg.N,2);
%NB list has to be built at the first iteration
gg.updateNB = true;
figure;
hg = scatter(gg.X(:,1),gg.X(:,2),1000*gg.r,gg.r,'o','filled');
hold on
plot(boxL2*[-1 1 1 -1 -1],boxL2*[-1 -1 1 1 -1],'k-.','linewidth',2);
axis square off
for i=1:0
drawframe(2)
end
end
%neighbour list update
if gg.updateNB
%NB list
[gg.i_NB,gg.j_NB,gg.NBestimate,gg.ACC] = compute_NB(gg.X,gg.boxL,gg.CutNbList,gg.NBestimate);
%interparticle distances for all pairs in the NB list
gg.rij=gg.r(gg.i_NB)+gg.r(gg.j_NB);
%upper bound on travel distance since last NB list update
gg.DistNB = zeros(gg.N,1);
%no need to update next time
gg.updateNB = false;
end
%explicit time stepping verlet velocity
%force calculation using neighbour list
% see Allen and Tildesley for references
[F] = computeForces(gg.X,gg.V,gg.rij,gg.ACC,gg.boxL,gg.parameters);
X_new = gg.X + gg.parameters.dt*gg.V + (0.5*gg.parameters.dt^2/gg.parameters.mass)*F;
gg.V = gg.V + (0.5*gg.parameters.dt/gg.parameters.mass)*F;
[F] = computeForces(X_new,gg.V,gg.rij,gg.ACC,gg.boxL,gg.parameters);
gg.V = gg.V + (0.5*gg.parameters.dt/gg.parameters.mass)*(F);
%update bound on displacement to update neighbour list
gg.DistNB = gg.DistNB + vecnorm(X_new-gg.X,2,2);
%update position
gg.X = X_new;
%test for neighbour list update
if any(gg.DistNB>gg.criterionNBLisdt)
gg.updateNB=true;
end
gg.X = gg.X-gg.boxL*round(gg.X/gg.boxL);
hg.XData = gg.X(:,1);
hg.YData = gg.X(:,2);
end
function [F] = computeForces(X,V,rij,ACC,boxL,parameters)
%compute repulsive forces between spheres
%spheres pairs that aree further appart than the sum of their radii experience no
%force
%if they are closer, a simple linear spring is considered
%viscosity is optional
%compute all distances
n = ACC'*X;
%apply minimum image convention for interaction
n = n-boxL*round(n/boxL);
%compute distances
d = sqrt(sum(n.^2,2));
%select only distances smaller than the sum of radii
mask = d<=rij;
d(~mask) = rij(~mask);
%compute potential energy
%E = sum((d-rij).^2);
%total energy (Kinetic + potential)
%E = E + 0.5*parameters.mass*sum(V(:).^2);
%compute force
f = 2*((rij-d).^(2-1)).*(n./d);
%net force on particles
F = ACC*f;
F = F-parameters.viscosity*V;
end
function [i_NB,j_NB,new_estimate,ACC] = compute_NB(X,boxL,cutoff,estimate)
%compute the list of neighbours closer than cutoff
%estimate is the average number of interactions. if 10 it means than on
%average each particle has 10 neighbours
%ACC is a sparse matrix with 1 and -1, used to easily compute things
N = size(X,1);
%initialization using size estimates
sz = ceil(N*estimate/2);
i_NB = zeros(sz,1);
j_NB = zeros(sz,1);
%totam number on NB interactions
NNB = 0;
%loop over all particles
for i=1:N
%compute only interactions with remaining particles (i>j)
%connecting vector
dX = X(i,:)-X((i+1):end,:);
%apply minimum image convention
dX = dX -boxL*round(dX/boxL);
%compute distance
D = sum(dX.^2,2);
%apply cutoff
mask = D<(cutoff^2);
%numner of new interactions
nloc = nnz(mask);
%accumulate nex interactions
i_NB(NNB+(1:nloc))=i;
j_NB(NNB+(1:nloc)) = find(mask)+i;
D_NB(NNB+(1:nloc)) = D(mask);
NNB = NNB+nloc;
end
%finalize sizes, trim unused memory
i_NB = i_NB(1:NNB);
j_NB = j_NB(1:NNB);
new_estimate = ceil(2*NNB/N);
%compute ACC matrix easy on Matlab
ACC = sparse([i_NB;j_NB],[(1:NNB) (1:NNB)]',[-ones(NNB,1);ones(NNB,1)],N,NNB);
end
Animation
Remix Tree