// *** 初期設定 *** clear; format('v',15); // *** 定数 *** // エネルギー ry = 2.179872E-18; // 1 Ry = 2.179872E-18 J // 長さ bohr = 0.5291772E-10; // 1 Bohr = 0.5291772E-10 m // 圧力 au = ry / (bohr ^ 3); // 1 a.u. = 1.471E4 GPa // *** ファイルの読み込み *** X = fscanfMat("CuVolume.txt"); A = X(:,1); // 格子定数 a (bohr) V = A .^ 3; // 格子体積 v (bohr^3) // エネルギー E1 = X(:,2); // bzqlty = 1 E3 = X(:,3); // bzqlty = 3 E9 = X(:,4); // bzqlty = 9 E16 = X(:,5); // bzqlty = 16 Vplot = linspace(min(V),max(V)); // *** Murnaghan状態方程式の定義 *** // エネルギー function Etot = eos(V) Etot = (b / (dbdp * (dbdp - 1))) * V .* (dbdp .* (1 - v0 ./ V) + ((v0 ./ V) .^ dbdp) - 1) + e0; endfunction // *** フィッティング関数の定義 *** function err = f1(fp,m) v0 = fp(1); // ゼロ気圧の原子体積 e0 = fp(2); // ゼロ気圧の全エネルギー b = fp(3); // 体積弾性率 dbdp = fp(4); // 体積弾性率の圧力微分 err = E - eos(V); endfunction V0 = []; E0 = []; B = []; dBdP = []; Eplot = []; // *** フィッティング *** for i = [3:5] do E = 4 * X(:,i); // fcc単位格子に4個の原子 [emin,ind] = min(E); // フィッティング実行 [fp, v] = lsqrsolve([V(ind); E(ind); 0.01; 4],f1,size(V,1)); v0 = fp(1); // ゼロ気圧の原子体積 e0 = fp(2); // ゼロ気圧の全エネルギー b = fp(3); // 体積弾性率 dbdp = fp(4); // 体積弾性率の圧力微分 // パラメータの保存 V0 = [V0; fp(1)]; // ゼロ気圧の原子体積 E0 = [E0; fp(2)]; // ゼロ気圧の全エネルギー B = [B; fp(3)]; // 体積弾性率 dBdP = [dBdP; fp(4)]; // 体積弾性率の圧力微分 Eplot = [Eplot; eos(Vplot)]; end // *** グラフのプロット *** plot(Vplot,Eplot - min(Eplot)); plot(V,4 * X(:,3:5) - min(Eplot),'o'); legend(["bzqlty=3","bzqlty=9","bzqlty=16"]); xlabel("Volume (bohr^3)"); ylabel("Energy (Ry)"); bohr * 1E10 * (V0 ^ (1 / 3)) au * 1E-9 * B dBdP