% MurnaghanSi.m
% Calculate Murnaghan's equation of state for Si under hydrostatic pressure
% Prepared by Stefan Zollner, New Mexico State University, Las Cruces, NM.
% Date: 07 January 2024
% Required package: None.
%
% References:
% F. D. Murnaghan, The compressibility of media under extreme pressures,
% Proceedings of the National Academy of Sciences 30, 244 (1944).
% A. K. Singh and G. C. Kennedy, Compressions of Si, MgO, and ZrSiO4 to 8 GPa as measured
% with a WC-anvil x-ray apparatus and epoxy pressure medium, J. Appl. Phys. 48, 3362 (1977).
% S. Anzellini, M. T. Wharmby, F. Miozzi, A. Kleppe, D. Daisenberger, and H. Wilhelm,
% Quasi-hydrostatic equation of state of silicon up to 1 megabar at ambient temperature,
% Scientific Reports 9, 15537 (2019). 
%
% Material parameters
B0Si=97.89; % bulk modulus of Si at 300 K in units of GPa (Anzellini 2019)
B0pSi=4.24; % derivative of bulk modulus versus pressure (Anzellini 2019)
a0Si=5.4316; % lattice constant of Si at 300 K at ambient pressure (Anzellini 2019)
BSi=@(P)B0Si+P*B0pSi; 
BSiconstant=@(P)B0Si+P*0; 
% The multiplication by zero is needed to make this a function of P.
% Otherwise MATLAB won't take this. 
%
fig=figure('Name','Si substrate properties versus pressure','Color','white');
tiledlayout(2,2,'TileSpacing','compact','Padding','none');
% plot a figure of bulk modulus versus pressure using Murnaghan's equation of state.
% BvsP=figure('Name','Bulk modulus versus pressure','Color','white');
P=0:0.1:10; % create a pressure vector from 0 to 10 GPa with 0.1 GPa steps
nexttile;
plot(P,BSi(P),'k-',P,BSiconstant(P),'k:'); 
% xlabel('pressure P (GPa)'); 
xticks(0:2:10);
ylabel('bulk modulus B (GPa)');
ylim([90 140]);
legend('Murnaghan','constant B_0', ...
       'Location','northwest','FontSize',10); 
legend('boxoff'); 
text(0.8,0.5,'(a)','FontSize',10,'Units','normalized');
% 
% plot a figure of hydrostatic strain as a function of pressure
% epsvsP=figure('Name','hydrostatic strain versus pressure','Color','white');
nexttile;
epsSiMu=@(P)(1+B0pSi/B0Si*P).^(-1/3/B0pSi)-1; % Murnaghan's EOS
epsSi=@(P)-P/3./BSi(P);                       % small linear coefficient B0'
epsSiconst=@(P)-P/3/B0Si;                     % constant bulk modulus B0
plot(P,epsSiMu(P),'k-',P,epsSiconst(P),'k:',P,epsSi(P),'k--');
% xlabel('pressure P (GPa)'); 
xticks(0:2:10);
ylabel('hydrostatic strain \epsilon_H');
ytickformat('%.2f');
legend('Murnaghan','constant B_0',"small B_0'", ...
       'Location','southwest','FontSize',10); 
legend('boxoff'); 
text(0.8,0.9,'(b)','FontSize',10,'Units','normalized');
%
% plot a figure of Si lattice constant versus pressure
% avsP=figure('Name','Si lattice constant versus pressure','Color','white');
nexttile;
aSiMu=@(P)a0Si*(1+epsSiMu(P));
aSiconst=@(P)a0Si*(1+epsSiconst(P));
aSi=@(P)a0Si*(1+epsSi(P));
plot(P,aSiMu(P),'k-',P,aSiconst(P),'k:',P,aSi(P),'k--');
xlabel('pressure P (GPa)'); xticks(0:2:10);
ylabel('Si lattice constant (A)');
ytickformat('%.2f');
legend('Murnaghan','constant B_0',"small B_0'", ...
       'Location','southwest','FontSize',10); 
legend('boxoff'); 
text(0.8,0.9,'(c)','FontSize',10,'Units','normalized');
% plot a figure of volume per atom versus pressure
% VvsP=figure('Name','Si relative volume change versus pressure','Color','white');
% Compare Fig. 2 in Singh (1977) or Fig. 1 in Anzellini (2019).
nexttile;
VSiMu=@(P)(1+epsSiMu(P)).^3;
VSi=@(P)(1+epsSi(P)).^3;
VSiconst=@(P)(1+epsSiconst(P)).^3;
colororder({'k','k'});
yyaxis left;
plot(P,VSiMu(P),'k-',P,VSiconst(P),'k:',P,VSi(P),'k--');
xlabel('pressure P (GPa)'); xticks(0:2:10);
ylabel('relative volume V/V_0');
ylim([0.9 1]);
ytickformat('%.2f');
yyaxis right;
yticklabels({'0.1','0.08','0.06','0.04','0.02','0.00'});
ylabel('relative volume change \Delta{}V/V_0');
legend('Murnaghan','constant B_0',"small B_0'", ...
       'Location','southwest','FontSize',10); 
legend('boxoff'); 
text(0.8,0.9,'(d)','FontSize',10,'Units','normalized');
