% fitForces.m
% Find Lennard-Jones parameters from a fit of a cluster trajectory
% Solution to problem 2 of control 2
% JMS, Feb.2019

% Set constants
clear
mass = 40;                  % argon atomic mass (Da)

% Read data from text file
dat = load('coords.dat');   % get data as a matrix
t = dat(:,1);               % times in first column
ra = dat(:,2:end)';         % coordinates in rest of columns
na = size(ra,1)/3;          % number of atoms
nt = numel(t);              % number of times
ra = reshape(ra,3,na,nt);   % reshape coordinates as a three-index array

% Alternatively, read coords and times from .mat file
% load('coords.mat')
% na = size(ra,2);          % number of atoms
% nt = numel(t);            % number of times

% Find forces from trajectories, discarding initial and final times
dt = t(2)-t(1);             % time step
it = (2:nt-1);              % all times except first and last
fa = mass*(ra(:,:,it+1)-2*ra(:,:,it)+ra(:,:,it-1))/dt^2;  % f=m*a
ra = ra(:,:,it);            % discard first and last times
t = t(it);  
nt = nt-2;

% Set inital parameters for nlinfit
eps0 = 100;                 % LJ energy parameter (K)
sigma0 = 5;                 % LJ distance parameter (Ang)
par0 = [eps0,sigma0];

% Find LJ parameters that fit forces
[par,~,~,cov] = nlinfit(ra, reshape(fa,3*na*nt,1), @forceModel, par0);
dpar = sqrt(diag(cov));     % statistical parameter errors
eps   = par(1);             % get parameters
sigma = par(2);
deps   = dpar(1);           % get errors
dsigma = dpar(2);
fprintf('  eps = %6.3f +/- %6.3f K\n',eps,deps)
fprintf('sigma = %6.3f +/- %6.3f Ang\n',sigma,dsigma)
