HA2EV = 27.21138602;
BOHR2ANG = 0.52917721067;
ve = [];
vc = [];
for ii = [10,100,1000,10000]
   tmp = sprintf('results/scf_f%d.mat',ii);
   mrc = load(tmp,'energy','atom');
   ve(end+1) = mrc.energy.Etot;
   vc(end+1) = mrc.atom.nvalence;
end

x = diff(vc);
y = diff(ve);
vbm = y./x * HA2EV;
fprintf('vbm = %f eV\n', vbm(end))

% figure()
% semilogx(x,vbm); hold on
% xlabel('\Delta N')
% ylabel('\Delta E / \Delta N (eV)')
% saveas(gcf,'eos_vbm.png')