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?
更多推荐
Matlab exercise07
发布评论