Fit a0EB.m

function [a0, Ecoh, B] = fit_a0EB( datafile )

%data = load('Elatt.B.dat'); data = load( datafile );

a = data(:,1); E = data(:,3); V = a.^3; n = 11; % number of atoms in computational cell N = 11; % number of atoms in an FCC unit cell

[P1, S1] = polyfit(a, E, 2); [P2, S2] = polyfit(V, E/n*N, 2);

a0 = -P1(2)/(2*P1(1)); V0 = a0^3; Ecoh = polyval(P1,a0)/n; B = 2*V0*P2(1) * 160.2; % conversion from eV/A^3 to GPa

disp(sprintf('a0 = %.6f (A)',a0)); disp(sprintf('Ecoh = %.6f (eV)',Ecoh)); disp(sprintf('B = %.6f (GPa)',B));

aaxis = [0:0.01:1]*(max(a)-min(a))+min(a); Efit = polyval(P1,aaxis); vaxis = aaxis.^3;

figure(1) plot(a,E,'o',aaxis,Efit,'-');