% cohesive.m %

set(0,'DefaultAxesFontSize',18);
set(0,'DefaultTextFontSize',18);

np = 2;
load coh.out;

clf;
plot( coh(:,1), coh(:,2), 'ro' );
hold on;

Xmin = min(coh(:,1));
Xmax = max(coh(:,1));

xmin = Xmin - 0.05 * (Xmax - Xmin);
xmax = Xmax + 0.05 * (Xmax - Xmin);
xmesh = 1000;
xdel = (xmax - xmin) / xmesh;
x = xmin + xdel/2 : xdel : xmax;

[p,s] = polyfit( coh(:,1), coh(:,2), 4 );

y = polyval(p,x);

plot(x,y);
[ymin i] = min(y);
xminapprox = x(i);

p1 = polyder(p);
xmin = fsolve('cohesivepoly', xminapprox, optimset('Display','off'), p1);
ymin = polyval(p, xmin);

units; 

% volume0 is volume of the unit cell = atomic volume * np %
volume0 = xmin^3 / 4;
p2 = polyder(p1);
y2min = polyval(p2, xmin);
bulkmodulus = y2min * xmin^2 / volume0 / 9 * EV__A3_IN_GPA;

fprintf (1, 'lattice constant = %g A,\n', xmin);
fprintf (1, 'energy = %g eV/atom,\n', ymin / np);
fprintf (1, 'bulk modulus = %g GPa.\n', bulkmodulus);

clf;
plot( coh(:,1), coh(:,2), 'ro' ); hold on; plot(x,y);
legend ('DFT data', '4th-order polynomial fit');
xlabel('Lattice constant [Angstrom]');
ylabel('Energy per Si-C atom pair [eV]');
axis([3.4 5.2 -16 -5 ]);

% input('Press a key.. '); %
print -depsc cohesive.eps;
break;
