new;

AVO = 6.0221367E+23;
EV_IN_J = 1.60217733E-19;
Zr_Z = 40;  Zr_A = 91.224;
C_Z = 6;    C_A = 12.0107;
ZrC_A = Zr_A + C_A;
ZrC_weight = ZrC_A / AVO;     
% one Zr atom + one C atom in g
 
% density [g/cm^3] vs energy/unit cell [eV]
NaCl = [
4.307770743806529 -1550.880075597751
4.44851538997441 -1551.142255315254
4.589139007123314 -1551.3782056237646
4.763166633509758 -1551.6534728142349
4.990531748559187 -1551.968040547747
5.160949684454654 -1552.1810156150127
5.341258716925984 -1552.3775891321993
5.5247992242039885 -1552.5544878694884
5.658188331743567 -1552.6625752554658
5.801513921741123 -1552.7640971197363
5.938073989577889 -1552.83939502132
6.084570539872632 -1552.9081274011967
6.277533156277965 -1552.9669856285605
6.46703736846604 -1552.9963384939601
6.765963916710855 -1552.9797000295764
6.9318947017277655 -1552.9402197584411
7.074509246238835 -1552.8876438448801
7.217033018985671 -1552.8153958745745
7.329671701221713 -1552.7464510887694
7.4290213971741235 -1552.677517195576
7.534955071758862 -1552.5954631515806
7.624292642116434 -1552.5167013994737
7.7103079659031 -1552.4379423705197
7.79629303243502 -1552.3526259893176
7.882247841712198 -1552.2607522558674
7.945037696757255 -1552.188569641232
4.381510699343493 -1551.0209987617218
4.532237213851468 -1551.2864489858903
4.689547706991771 -1551.5387790592565
4.893581452907987 -1551.836972474219
5.074087157535153 -1552.0761687810184
5.264468830110809 -1552.2956848616147
5.444732476700023 -1552.4824223504288
6.6132464748407145 -1553.0027760274795
7.154061620412149 -1552.848234375721
];


ZnS = [
4.133416339068848 -1551.2539890030007
4.2302456057010716 -1551.3588276675364
4.377059857170629 -1551.4964122460187
4.493777217346182 -1551.5913985432644
4.603759312615686 -1551.6667182300716
4.776954864496667 -1551.7616582337175
4.9134241605692 -1551.817284078557
5.069751292930317 -1551.8565002038583
5.225957396272458 -1551.869486920167
5.3820575992229935 -1551.8595229036073
5.5513106308108116 -1551.820039909319
5.760203691838938 -1551.7313440953346
5.833096444243066 -1551.6886613963572
5.912603432534265 -1551.6394158988257
6.0152299891754675 -1551.5606432541072
6.091369345013643 -1551.5015644513564
6.167463314969702 -1551.4326496202336
6.236912791783945 -1551.3637402354166
6.299717775456374 -1551.2948362969053
6.365829877072338 -1551.2226509591171
6.431896592806186 -1551.1406295929567
6.4880570833368 -1551.0717311007513
6.752097016861612 -1550.6944654942486
6.844424003987906 -1550.543570144259
6.936735862486829 -1550.3893961181454
7.035571185108589 -1550.2089872367335
7.101471485941345 -1550.0909004332084
7.167296143637239 -1549.9564202490626
7.216675983379689 -1549.8580191180465
];

CsCl = [
5.335709536405907 -1549.7349707298883
5.399134793800592 -1549.8004925124642
5.455870172171344 -1549.8561837129737
5.509253046716445 -1549.905320284388
5.565973296459825 -1549.9577328087735
5.612666291980995 -1549.9970387981214
5.666018909271353 -1550.0396180172875
5.816019249665468 -1550.1476917875002
5.879338606668548 -1550.1902628372077
5.945980210242535 -1550.232831163762
6.019266306958336 -1550.2753940440107
6.08587765327758 -1550.311405018317
6.189094226386294 -1550.3605007424374
6.23903382534151 -1550.3834106280121
6.2822986739001685 -1550.3997686076443
6.378779982102834 -1550.429197721326
6.435348945572494 -1550.448823484471
6.568374966055142 -1550.4782226434704
6.618254050500871 -1550.4880178245487
6.66479575974832 -1550.494537052656
6.761186296186754 -1550.5042941095933
6.821001863090461 -1550.5075237689648
6.963843337012111 -1550.5041279972647
7.033565129119053 -1550.4942347826811
7.090013063569735 -1550.4876311368332
7.202863546588986 -1550.464587816765
7.279184445955627 -1550.4448531275034
7.401926025550739 -1550.4054082573562
7.468249927949913 -1550.379124385305
7.547832559377971 -1550.346272268394
7.620755569036843 -1550.3101469216647
7.690341203497436 -1550.2707456219644
7.756604591387122 -1550.231347045417
7.82615996859297 -1550.1853883934684
7.935461275630732 -1550.1131676546922
  ];

xmin = 4; xmax = 8; xmesh = 100;
x = xmin:(xmax-xmin)/xmesh:xmax;
ymin = -1553.3; ymax = -1549.7;

clf;
plot (NaCl(:,1), NaCl(:,2), 'go');
hold on;
[NaCl_P,NaCl_S] = polyfit(NaCl(:,1),NaCl(:,2),2);
NaCl_x = x;  NaCl_y = polyval (NaCl_P, NaCl_x);
NaCl_xmin = -NaCl_P(2)/2/NaCl_P(1);
NaCl_ymin = polyval (NaCl_P, NaCl_xmin);
NaCl_volume = (ZrC_weight / NaCl_xmin) * (1E8)^3;
NaCl_constant = (4 * NaCl_volume) ^ (1./3);
NaCl_B = 2 * NaCl_P(1) * NaCl_xmin^2 * EV_IN_J / (NaCl_volume*1E-30);

% purple
COLOR_PURPLE = [0.5742    0.4375    0.8555];
h = plot (ZnS(:,1), ZnS(:,2), 'o');
set(h, 'Color', COLOR_PURPLE);
[ZnS_P,ZnS_S] = polyfit(ZnS(:,1),ZnS(:,2),2);
ZnS_x = x;  ZnS_y = polyval (ZnS_P, ZnS_x);
ZnS_xmin = -ZnS_P(2)/2/ZnS_P(1);
ZnS_ymin = polyval (ZnS_P, ZnS_xmin);
ZnS_volume = (ZrC_weight / ZnS_xmin) * (1E8)^3;
ZnS_constant = (4 * ZnS_volume) ^ (1./3);
ZnS_B = 2 * ZnS_P(1) * ZnS_xmin^2 * EV_IN_J / (ZnS_volume*1E-30);

plot (CsCl(:,1), CsCl(:,2), 'ro');
[CsCl_P,CsCl_S] = polyfit(CsCl(:,1),CsCl(:,2),2);
CsCl_x = x; CsCl_y = polyval (CsCl_P, CsCl_x);
CsCl_xmin = -CsCl_P(2)/2/CsCl_P(1);
CsCl_ymin = polyval (CsCl_P, CsCl_xmin);
CsCl_volume = (ZrC_weight / CsCl_xmin) * (1E8)^3;
CsCl_constant = (1 * CsCl_volume) ^ (1./3);
CsCl_B = 2 * CsCl_P(1) * CsCl_xmin^2 * EV_IN_J / (CsCl_volume*1E-30);

grid on;
xlabel ('Density [g/cm^3]');
ylabel ('Energy/unit cell [eV]');
axis ([xmin xmax ymin ymax]);
title ...
('Cohesive energies of ZrC in NaCl, ZnS, CsCl structures with quadratic fit');

plot (NaCl_x, NaCl_y, 'g');
h = plot (ZnS_x, ZnS_y);
set(h, 'Color', COLOR_PURPLE);
plot (CsCl_x, CsCl_y, 'r');

fprintf ('NaCl = (%f) * x.^2 + (%f) * x + (%f);\n',  NaCl_P);
fprintf ('xmin = %f g/cm^3; volume = %f A^3/cell; ymin = %f eV/cell;\n', ...
    NaCl_xmin, NaCl_volume, NaCl_ymin);
fprintf ('lattice = %f A; B =%e Pa;\n\n', NaCl_constant, NaCl_B);

fprintf ('ZnS = (%f) * x.^2 + (%f) * x + (%f);\n',  ZnS_P);
fprintf ('xmin = %f g/cm^3; volume = %f A^3/cell; ymin = %f eV/cell;\n', ...
    ZnS_xmin, ZnS_volume, ZnS_ymin);
fprintf ('lattice = %f A; B =%e Pa;\n\n', ZnS_constant, ZnS_B);

fprintf ('CsCl = (%f) * x.^2 + (%f) * x + (%f);\n',  CsCl_P);
fprintf ('xmin = %f g/cm^3; volume = %f A^3/cell; ymin = %f eV/cell;\n', ...
    CsCl_xmin, CsCl_volume, CsCl_ymin);
fprintf ('lattice = %f A; B =%e Pa;\n\n', CsCl_constant, CsCl_B);

NaCl_s = sprintf ('NaCl: %fx^2%fx%f', NaCl_P);
ZnS_s = sprintf ('ZnS:  %fx^2%fx%f', ZnS_P);
CsCl_s = sprintf ('CsCl: %fx^2%fx%f', CsCl_P);

legend (NaCl_s, ZnS_s, CsCl_s, 2);
print -depsc cohesive_energy.eps;
