clear; format('e',12); // ************************************ // // 入力パラメータ // // ************************************ // アルミニウムの格子定数: a (bohr) a = 7.6517; // アルミニウムの価電子数 n = 3; // ************************************ // // 物理定数 // // ************************************ // リュードベリ原子単位系でのボルツマン定数 // Boltzmann constant kB = 1.3806488E-23 (J/K) // cf. 1 (Ry) = 2.179872E-18 (J) kB = 1.3806488E-23 / 2.179872E-18; // (Ry/K) // ************************************ // // 化学ポテンシャルの計算 // // ************************************ // エネルギーベクトルの生成 E = linspace(0,2,2000)'; // 温度ベクトルの生成 T = linspace(10,20000,1000); // アルミニウムの原子体積(bohr3/atom) // fccは単位格子に原子4個 vat = (a .^ 3) / 4; // フェルミエネルギー(Ry) ef = (3 * %pi * %pi * n / vat) ^ (2/3); // 自由電子モデルでの状態密度 (state/Ry/atom) function Dfree = Dfree(energy); Dfree = vat / 2 / %pi / %pi .* (energy ^ (1/2)) endfunction // Sommerfeld展開による化学ポテンシャル Musomm = ef - %pi * %pi * kB * kB / 12 / ef .* T .* T; // フェルミ分布関数の定義 function fermi = fermi(mu, energy, temp); fermi = 1 ./ (exp((energy - mu) ./ (kB * temp)) + 1) endfunction // 化学ポテンシャルの計算 deff('[y] = f1(x,temp)','y = intsplin(E,Dfree(E) .* fermi(x,E,temp)) - 3'); // 化学ポテンシャルベクトルの初期化 Munum = ones(T); for i = 1:length(T) do temp = T(i); // fsolve関数を使って化学ポテンシャルμを計算 Munum(i) = fsolve(0,f1); end // ************************************ // // グラフ描画 // // ************************************ // 化学ポテンシャルの描画 xsetech([0,0,0.95,0.95]); // Sommerfeld展開による化学ポテンシャル plot(T,Musomm - ef,'-b'); // 数値積分による化学ポテンシャル plot(T,Munum - ef,'--r'); legend(['Sommerfeld expansion';'Numerical calculation'],1); xlabel('Temperature (K)'); ylabel('Chemical potential (Ry)'); // 状態密度の描画(インセット) xsetech([0.15,0.38,0.4,0.4]); // 状態密度 plot(E-ef,Dfree(E)); // 最低温度(10K) plot(E-ef,Dfree(E).*fermi(Munum(1),E,T(1)),'-g'); // 最高温度(20,000K)での電子の分布 plot(E-ef,Dfree(E).*fermi(Munum(length(T)),E,T(length(T))),'-r'); legend(['DOS';'T = 10 K';'T = 20,000 K'],2); xlabel('Energy (Ry)'); ylabel('Densitiy of states (state/Ry/atom)');