% Find a0, B, Ecoh of Si, SiGe, Ge from VASP calculations
% a0 : eql. latt. const.
% B  : Bulk modulus (GPa)
% Ecoh: cohesive energy eV
%
% compare different crystal structures

generic = struct('N', 0, 'C1', 0, 'C2', 0, 'C3', 0, 'data', 0);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Si(1-x)Ge(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%MEAM Si(1-x)Ge(x) zinc-blend structure data with x = 1/512
% from ~/Codes/MD++/sige-meam/ELatt.dat
MEAM_SiGe = generic;
MEAM_SiGe.N = 512;
MEAM_SiGe.C1 = [4 0 0];
MEAM_SiGe.C2 = [0 4 0];
MEAM_SiGe.C3 = [0 0 4];
%MEAM_SiGe.data = [  5.440:0.001:5.449 ]'; % latt_const.(A)   Epot(eV) (32)
%MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
%MEAM_SiGe.data(:,3) = [
% -2.34405440143234e+03
% -2.34406174928596e+03
% -2.34406720346273e+03
% -2.34407076735192e+03
% -2.34407244433792e+03
% -2.34407223780089e+03
% -2.34407015111614e+03
% -2.34406618765477e+03
% -2.34406035078301e+03
% -2.34405264386292e+03 
%];
%MEAM_SiGe.data = [  5.450:0.001:5.466 ]'; % latt_const.(A)   Epot(eV) (64)
%MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
%MEAM_SiGe.data(:,3) = [
% -2.31773141774331e+03
% -2.31774521039280e+03
% -2.31775712179189e+03
% -2.31776715530201e+03
% -2.31777531428004e+03
% -2.31778160207838e+03
% -2.31778602204490e+03
% -2.31778857752300e+03
% -2.31778927185160e+03
% -2.31778810836514e+03
% -2.31778509039359e+03
% -2.31778022126244e+03
% -2.31777350429273e+03
% -2.31776494280100e+03
% -2.31775454009953e+03
% -2.31774229949584e+03
% -2.31772822429330e+03
%];
%MEAM_SiGe.data = [  5.461:0.001:5.482 ]'; % latt_const.(A)   Epot(eV) (96)
%MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
%MEAM_SiGe.data(:,3) = [
% -2.29158442898549e+03
% -2.29160282165048e+03
% -2.29161934879804e+03
% -2.29163401375692e+03
% -2.29164681985147e+03
% -2.29165777040154e+03
% -2.29166686872252e+03
% -2.29167411812539e+03
% -2.29167952191667e+03
% -2.29168308339849e+03
% -2.29168480586851e+03
% -2.29168469262002e+03
% -2.29168274694186e+03
% -2.29167897211850e+03
% -2.29167337142998e+03
% -2.29166594815196e+03
% -2.29165670555571e+03
% -2.29164564690806e+03
% -2.29163277547159e+03
% -2.29161809450451e+03
% -2.29160160726042e+03
% -2.29158331698870e+03 
%];
%MEAM_SiGe.data = [  5.475:0.001:5.496 ]'; % latt_const.(A)   Epot(eV) (128)
%MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
%MEAM_SiGe.data(:,3) = [
% -2.26564169154174e+03
% -2.26565912522176e+03
% -2.26567471886512e+03
% -2.26568847575495e+03
% -2.26570039916976e+03
% -2.26571049238390e+03
% -2.26571875866717e+03
% -2.26572520128517e+03
% -2.26572982349893e+03
% -2.26573262856533e+03
% -2.26573361973671e+03
% -2.26573280026122e+03
% -2.26573017338251e+03
% -2.26572574234006e+03
% -2.26571951036883e+03
% -2.26571148069965e+03
% -2.26570165655881e+03
% -2.26569004116851e+03
% -2.26567663774641e+03
% -2.26566144950608e+03
% -2.26564447965658e+03
% -2.26562573140287e+03
%];
%MEAM_SiGe.data = [  5.491:0.001:5.507 ]'; % latt_const.(A)   Epot(eV) (160)
%MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
%MEAM_SiGe.data(:,3) = [
% -2.24002872831514e+03
% -2.24004192804737e+03
% -2.24005331936961e+03
% -2.24006290550952e+03
% -2.24007068969042e+03
% -2.24007667513133e+03
% -2.24008086504700e+03
% -2.24008326264786e+03
% -2.24008387114005e+03
% -2.24008269372544e+03
% -2.24007973360161e+03
% -2.24007499396189e+03
% -2.24006847799531e+03
% -2.24006018888666e+03
% -2.24005012981646e+03
% -2.24003830396098e+03
% -2.24002471449232e+03
%];
%MEAM_SiGe.data = [  5.505:0.001:5.520 ]'; % latt_const.(A)   Epot(eV) (192)
%MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
%MEAM_SiGe.data(:,3) = [
% -2.21446842728671e+03
% -2.21448104612199e+03
% -2.21449188115643e+03
% -2.21450093557197e+03
% -2.21450821254633e+03
% -2.21451371525296e+03
% -2.21451744686109e+03
% -2.21451941053573e+03
% -2.21451960943764e+03
% -2.21451804672340e+03
% -2.21451472554534e+03
% -2.21450964905161e+03
% -2.21450282038614e+03
% -2.21449424268866e+03
% -2.21448391909474e+03
% -2.21447185273571e+03
%];
%MEAM_SiGe.data = [  5.520:0.001:5.533 ]'; % latt_const.(A)   Epot(eV) (224)
%MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
%MEAM_SiGe.data(:,3) = [
% -2.18915242979307e+03
% -2.18916299517619e+03
% -2.18917180420555e+03
% -2.18917886001304e+03
% -2.18918416572638e+03
% -2.18918772446913e+03
% -2.18918953936067e+03
% -2.18918961351625e+03
% -2.18918795004695e+03
% -2.18918455205970e+03
% -2.18917942265729e+03
% -2.18917256493839e+03
% -2.18916398199753e+03
% -2.18915367692510e+03
%];
%MEAM_SiGe.data = [  5.534:0.001:5.547 ]'; % latt_const.(A)   Epot(eV) (256)
%MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
%MEAM_SiGe.data(:,3) = [
% -2.16413667801199e+03
% -2.16414739375472e+03
% -2.16415637710621e+03
% -2.16416363115174e+03
% -2.16416915897239e+03
% -2.16417296364520e+03
% -2.16417504824308e+03
% -2.16417541583486e+03
% -2.16417406948529e+03
% -2.16417101225505e+03
% -2.16416624720074e+03
% -2.16415977737487e+03
% -2.16415160582591e+03
% -2.16414173559826e+03
%];
MEAM_SiGe.data = [  5.548:0.001:5.562 ]'; % latt_const.(A)   Epot(eV) (288)
MEAM_SiGe.data = [ MEAM_SiGe.data, ones(length(MEAM_SiGe.data(:,1)),2) ];
MEAM_SiGe.data(:,3) = [
 -2.13929824528630e+03
 -2.13930918618158e+03
 -2.13931841824069e+03
 -2.13932594450304e+03
 -2.13933176800393e+03
 -2.13933589177467e+03
 -2.13933831884248e+03
 -2.13933905223059e+03
 -2.13933809495823e+03
 -2.13933545004058e+03
 -2.13933112048885e+03
 -2.13932510931021e+03
 -2.13931741950787e+03
 -2.13930805408102e+03
 -2.13929701602487e+03
];
%
MEAM_SiGe = fit_aEB(MEAM_SiGe);


fs = 16;
figure(1);
plot(MEAM_SiGe.a,MEAM_SiGe.Epot/MEAM_SiGe.N - MEAM_SiGe.Ecoh,'b.', ...
     MEAM_SiGe.aaxis,MEAM_SiGe.Efit/MEAM_SiGe.N - MEAM_SiGe.Ecoh,'r-' );
 
set(gca,'FontSize',fs); 
%xlabel('a  (angstrom)');
%text('Interpreter','latex','String','$a \ \ (\AA)$',...
%     'Position',[5.52,-7e-5],'FontSize',fs);
ylabel('E - E_{coh} (eV)','FontName','Times New Roman');
%title('LDA + MEAM');

t1=text(LDA_SiGe.a0  -0.001,-2e-5,'LDA'); set(t1,'FontSize',fs-2);
t1a=text(MEAM_SiGe.a0-0.002,-2e-5,'MEAM');set(t1a,'FontSize',fs-2,'Color','blue');
t2=text(5.502,7e-5,'ZB');  set(t2,'FontSize',fs-2);
t2a=text(5.510,7e-5,'RH');  set(t2a,'FontSize',fs-2);
t3=text(5.529,7e-5,'ZB');  set(t3,'FontSize',fs-2);
t3a=text(5.534,7e-5,'RH');  set(t3a,'FontSize',fs-2);
%xlim([5.3 5.8]);
%ylim([-5 15]*1e-5);
title('SiGe');


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

x = [0:32:512]/512;
a0x = inline('5.431 + 0.20*x + 0.027*x.^2');
Bx  = inline('97.9 - 22.8*x');

disp(sprintf('experiment: Si: a0 = %.3f  SiGe: a0 = %.3f  Ge: a0 = %.3f', ...
              a0x(0),a0x(0.5),a0x(1)));
a0_EXP = a0x(x);
B_EXP  = Bx (x);

a0_MEAM = a0_EXP;
a0_MEAM(1:10) = [5.431,5.444395,5.457886,5.471460,5.485069,5.498856,5.512624, ...
    5.526551,5.540723,5.554943]; 
a0_MEAM(end)=5.6575;

B_MEAM = B_EXP;
B_MEAM(1:10) = [97.6141,96.121866,94.659885,93.226247,91.685625,90.575934,88.969749, ...
    87.529811,86.131865,84.643183];
B_MEAM(end) = 75.0858;
              
figure(2);
title('lattice constant of Si(1-x)Ge(x)');
plot(x,a0_EXP,'b-',x,a0_MEAM,'r.');
set(gca,'FontSize',fs);
xlabel('x','FontName','Times');
ylabel('a_0 (angstrom)','FontName','Times');
legend('Expt','MEAM',2);

figure(3);
title('bulk modulus of Si(1-x)Ge(x)');
plot(x,B_EXP,'b-',x,B_MEAM,'r.');
set(gca,'FontSize',fs);
xlabel('x','FontName','Times');
ylabel('B (GPa)','FontName','Times');
legend('Expt','MEAM');
