function [fa,u] = forceLJ(ra,eps,sigma)
% Lennard-Jones forces and energy for a cluster of noble gas atoms
%   uij = eps*((sigma/dij)^12-2*(sigma/dij)^6)
%   fij = -duij/drij 
%       = 12*eps*((sigma/dij)^12-(sigma/dij)^6)*rij/dij^2
%   u = sum_(i<j)(uij)
%   fi = sum_(j~=i)(fij)
%   rij=ri-rj, dij=|rij|
% Input:
%   ra(3,na) : atomic coordinates (Ang)
%   eps      : LJ energy parameter (K)
%   sigma    : LJ distance parameter (Ang)
% Output:
%   fa(3,na) : atomic forces (K/Ang)
%   u        : potential energy (K)
% JMS, Feb.2019

na = size(ra,2);                     % number of atoms
u = 0;
fa = zeros(3,na);
for ia = 1:na-1                      % double loop on pairs of atoms
  for ja = ia+1:na
    rij = ra(:,ja)-ra(:,ia);         % vector from ia to ja
    dij = norm(rij);                 % distance between ia and ja
    uij = eps*((sigma/dij)^12-2*(sigma/dij)^6);            % interaction energy
    fij = 12*eps*((sigma/dij)^12-(sigma/dij)^6)*rij/dij^2; % force of ia on ja
    u = u+uij;                       % add to total potential energy
    fa(:,ia) = fa(:,ia)-fij;         % add to force on ia
    fa(:,ja) = fa(:,ja)+fij;         % add to force on ja
  end
end

end
