function forces = LJ_Force(coords,L)

        % Initialize all forces to 0
        forces = zeros(size(coords));

        % Get the number of particles
        nPart = size(coords,2);

        % take the x, y, z coordinators, respectively
        xcoords = coords(1,:);
        ycoords = coords(2,:);
        zcoords = coords(3,:);

        % do xi-xj, yi-yj, zi-zj and take closest image
        dx = repmat(xcoords',[1,nPart])-repmat(xcoords,[nPart,1]);
        dy = repmat(ycoords',[1,nPart])-repmat(ycoords,[nPart,1]);
        dz = repmat(zcoords',[1,nPart])-repmat(zcoords,[nPart,1]);
        dx = dx - L*round(dx/L);
        dy = dy - L*round(dy/L);
        dz = dz - L*round(dz/L);

        % tak r^2, dx, dy, dz in vector form for i<j
        dr2 = dx.^2 + dy.^2 + dz.^2;
        [row,col,DR2] = find(triu(dr2,1));
        linearIndex = (col-1)*nPart+row;
        dx = dx(linearIndex);
        dy = dy(linearIndex);
        dz = dz(linearIndex);

        % calculate the pairwise force in vector form
        invDR2    = zeros(size(DR2));
        cfIndex = (DR2<0.25*L^2);
        invDR2(cfIndex) = 1./DR2(cfIndex);
        forceFact = 48*invDR2.^4.*(invDR2.^3 - 0.5);

        xForce = dx.*forceFact;
        yForce = dy.*forceFact;
        zForce = dz.*forceFact;

        % change the pairwise force in matrix form
        xForce = full(sparse(row,col,xForce,nPart,nPart));
        yForce = full(sparse(row,col,yForce,nPart,nPart));
        zForce = full(sparse(row,col,zForce,nPart,nPart));

        % sum up the pairwise force
        xForce = sum(-xForce'+xForce,2);
        yForce = sum(-yForce'+yForce,2);
        zForce = sum(-zForce'+zForce,2);

        forces = [xForce';yForce';zForce'];

end

Published with MATLAB® R2012b