HA2EV = 27.21138602;
BOHR2ANG = 0.52917721067;
dlst = dir('results/scf_f*.mat');
npt = numel(dlst);
avol = []; anrg = [];
for ii = 1:numel(dlst)
   mrc = load(['results/',dlst(ii).name]);
   avec = mrc.domain.latvec;
   natom = size(mrc.atom.xyz,1);
   vol = abs(det(avec)) * BOHR2ANG^3 / natom;
   etot = mrc.energy.Etot * HA2EV / natom;
   avol = [avol;vol];
   anrg = [anrg;etot];
end
[avol,ivol] = sort(avol);
anrg = anrg(ivol);
[v0, bmod0, bmodder0, res] = eosfit(avol,anrg);
fprintf('v_eq = %f ang^3/atom\n',v0)

% figure(); hold on
% plot(avol,anrg)
% xlabel('volume (A^3/atom)')
% ylabel('energy (eV/atom)')
% saveas(gcf,'eos_diamond.png')
