% clusterMotion.m
% Dynamics of a cluster of n noble gas atoms
% Solution to problem 1 of control 2
% JMS, Feb.2019

% Set constants
clear
mass = 40;           % atomic mass (Da) (argon)
eps = 120;           % LJ energy parameter (K)
sigma = 3.79;        % LJ distance parameter (Ang)
T0 = 50;             % initial temperature (K)
dt = 5e-3;           % integration time step (time units)
tmax = 10;           % integration time (time units)
ntplot = 100;        % number of time steps between plots
ntwrite = 10;        % numer of time steps between writes
rnoise = 1.e-2;      % random noise in positions (Ang)

% Set initial positions and velocities
% l=sigma/2/sqrt(2); ra0=[l,l,l; -l,-l,l; ...
%                         l,-l,-l; -l,l,-l]';         % tetrahedron, na=4
% l=sigma/sqrt(2); ra0=[l,0,0; -l,0,0; 0,l,0; ...
%                       0,-l,0; 0,0,l; 0,0,-l]';      % octahedron, na=6
l=sigma/2; ra0=[l,l,l;  l,-l,l;  l,l,-l;  l,-l,-l; ...
               -l,l,l; -l,-l,l; -l,l,-l; -l,-l,-l]';  % cube, na=8
na = size(ra0,2);                  % number of atoms
va0 = randn(3,na).*sqrt(T0/mass);  % Boltzman distribution at temp. T0

% Subtract center of mass position and velocity
rcm = sum(ra0,2)/na;               % center of mass position
vcm = sum(va0,2)/na;               % center of mass velocity
ra0 = ra0-repmat(rcm,1,na);        % positions in CM system
va0 = va0-repmat(vcm,1,na);        % velocities in CM system

% Integrate equations of motion
myforce = @(r)forceLJ(r,eps,sigma);  % force function, as needed by verlet
ma = mass*ones(1,na);                % masses of each atom
[ra,va] = verlet( myforce, ma, ra0, va0, dt, tmax );
t = 0:dt:tmax;                       % trajectory times
nt = numel(t);                       % number of times

% Plot energy and temperature as a function of time
for it = 1:nt                      % loop on times
  ek(it) = 0.5*sum(mass*sum(va(:,:,it).^2)); % kinetic energy at time it
  [~,u(it)] = forceLJ(ra(:,:,it),eps,sigma); % potential energy
end
figure(1)
subplot(2,1,1)
plot(t,ek+u)               % plot total energy as a function of time
xlabel('time (Ut)')
ylabel('energy/kB (K)')
subplot(2,1,2)
plot(t,ek/(3*na/2))        % plot temperature as a function of time
xlabel('time (Ut)')
ylabel('temperature (K)')

% Plot motion
figure(2)
[xs,ys,zs] = sphere();     % sphere points, of radius one and center zero
xs = xs*sigma/2;           % scale sphere radius to sigma/2
ys = ys*sigma/2;
zs = zs*sigma/2;
for it = 1:ntplot:nt       % loop on selected times
  for ia = 1:na            % loop on atoms
    % Plot a sphere centered at ra(ia,it)
    surf(ra(1,ia,it)+xs, ra(2,ia,it)+ys, ra(3,ia,it)+zs) 
%    % Alterantively just plot a circle centered at ra(ia,it)
%    plot3(ra(1,ia,it), ra(2,ia,it), ra(3,ia,it), 'o', ...
%          'MarkerSize',30)
    hold on                % not to erase previous atoms
  end
  hold off                 % to allow a clean figure at next time
  axis equal               % to make x,y, and z distances equal
  axis(3*[-l,l,-l,l,-l,l]) % to keep plotted box fixed
  grid on                  % to have a positions reference
  xlabel('x (Ang)')
  ylabel('y (Ang)')
  zlabel('z (Ang)')
  pause(0.1)               % show each step during 0.1 seconds
end

% Add random noise to all positions
ra = ra+rnoise*randn(3,na,nt);

% Write time and coordinates every ntwrite steps, as a text file
file = fopen('coords.dat','w');  % open file for writting
for it = 1:ntwrite:nt
  fprintf(file,' %6.3f',t(it),ra(:,:,it));  % each time in a line
  fprintf(file,'\n');      % finish line with carriage return
end
fclose(file);              % close file

% Save selected times and coords as a .mat file
it = 1:ntwrite:nt;   % every ntwrite time steps
t = t(it);
ra = ra(:,:,it);
save('coords.mat','t','ra')
