Matlab exercise07

Using potential.txt, where the first column is distance ® in Angstrom and second column is energy (E) in atomic units, find a potential representation in the following forms
(a) Lennard-Jones potential - find parameters epsilon and sigma
(b) Buckingham potential - find parameters A, B, C.

使用 potential.txt，其中第一列是以埃为单位的距离 ®，第二列是以原子单位为单位的能量 (E)，找到以下形式的潜在表示
(a) Lennard-Jones 势——找到参数 epsilon 和 sigma
(b) Buckingham potential - 找到参数 A、B、C。

% write solution here 
D = load('potential.txt'); 
r = D(7:end,1); 
e = D(7:end,2); 
% simplify Vlj(r) = a*r^-12 - b*r^-6
% linear parameters solution via set of equations
eq = [r.^-12  -r.^-6];
p = eq\e; 
elj = p(1).*r.^-12 - p(2).*r.^-6; 

plot(r,e,'rs'); hold on
plot(r,elj,'g-');
% NOTICE that elj > 0 p(2) < 0, thus fit provides a repulsive curve
% if you want to improve LJ-fit you need to remove some points from the
% repulsion region - we are fitting 2 parameters, thus there is enough
% points in the input data. 

% V(r) = A*exp(-Bx) - C/r^6
bug = @(p,r) p(1).*exp(-p(2).*r) - p(3).*r.^-6;
p0 = [1 1 1]; 
p = lsqcurvefit(bug,p0,r,e);

plot(r,bug(p,r),'b-');
legend('DATA','LJ','Buck'); 

% check the parameter stability for both p0 vectors? 
% Any ideas - how to improve on numerical stability?