這幾天原本做電路實驗的時候看到FFT的分析，然后查了一下看到了有用matlab來模擬多少個波的疊加可以得到比較好看的三角波、方波之類的，一時心血來潮，就想著要不我也來試試用matlab畫出能帶圖，于是渣渣我又來了，順便復習一波matlab畫三維圖的知識，手冊走起。

---

首先，在利用MATLAB進行畫圖前，要先對公式進行處理，順便復習就從頭再推一次吧。

以下推導以簡單立方晶格s態電子為例，后面有MATLAB代碼附上（因為我自己學的不是很扎實，所以有很多小白注釋，大佬可以跳過）

簡單晶格中電子滿足薛定諤方程：

根據緊束縛近似，

可以寫成布洛赫和的形式：

其中

表示孤立原子波動方程的本征態，代入薛定諤方程中：

其中的周期勢

為孤立原子勢的疊加，利用：

可以得到：

左乘

并積分，利用正交關系，得到：

令

，則得到：

結合上面兩式，即可解出電子能譜為：

由于在緊束縛近似下，各格位上孤立原子的波函數之間交疊很少，我們保留到近鄰項，其他項忽略，得到：

對于簡單立方晶格而言，s態波函數是球對稱的，因此各個方向的交疊積分是相同的，具有相同的值，于是我們可以提出來，同時由于簡單立方晶格六個近鄰格點為：

因此代入得到：

以上理論推導已經完成，最終的式子也已經出來了，后面需要進行MATLAB的模擬。

---

首先我們構建三維數組，生成三維網格。

figure; [x, y, z] = meshgrid(lnspace(-2, 2, 101)); % lnspace是用于生成-2～2中間101個點，meshgrid可用于生成三維網格，x、y、z都是101x101x101的三維矩陣

2. 然后我們需要對圖形范圍進行限制，如果不限制，則等能面將在倒格子單胞里繪制，故要把簡約布里淵區以外的圖形切掉。

x1 = x(1, :, 1); y1 = y(:, 1, 1); z1 = z(1, 1, :); % 此處就是選擇了三個方向的坐標來尋找，具體為什么這樣選可以看生成meshgrid后的情況

3. 先寫出等能面的方程

首先又回到我們的公式，我們需要對能帶表達式無量綱化處理，

由此無量綱處理完畢，轉換成代碼即為（能帶方程），

val=1-(cos(pi*x)+cos(pi*y)+cos(pi*z));

4. 將布里淵區外的能帶取為非數，方便畫圖與看，此處利用的是如果超過了設定的布里淵區（1中已設定）的面則為非數，對于簡單立方而言很簡單，布里淵區是一個立方體，每個面都平行于兩個軸，十分容易操作，

for i = 1:101for j = 1:101for k = 1:101if x1(i)-1>0 | -x1(i)-1>0 | y1(j)-1>0 | -y1(j)-1>0 | z1(k)-1>0 | -z1(k)-1>0;val(i,j,k) = NaN;endendend end

5. 繪制面心立方等能面

p = patch(isosurface(x, y, z ,val, -1)); % patch函數就是布丁函數，在布丁的位置上色，isosurface則是對幾個坐標的體現，val是它們的體數據，-1是 % 對體數據的篩選，代表采用的等能面的補片 set(p, 'FaceColor', 'yellow', 'EdgeColor', 'none');

isovalue為-1時選取的等能面截面

isovalue為-1時等能面的立體圖

6. 繪制簡單立方的布里淵區：立方體

這個部分需要一定的幾何知識，可以直接引用現成的也可以自己推導。

各個坐標點的選擇順序% 各個點的順序選取 Vm1=[1 1 1;1 -1 1;-1 -1 1;-1 1 1;1 1 -1;1 -1 -1;-1 -1 -1;-1 1 -1]; % 各個面的坐標順序（連線順序） Fm1=[1 2 3 4;5 6 7 8;3 4 8 7;2 1 5 6;1 4 8 5;2 3 7 6]; patch('Vertices', Vm1, 'Faces', Fm1, 'Linestyle', '-', 'LineWidth', 2, 'FaceColor', 'none'); % Vertices為連接的點，Faces為點連接順序，此處相當于作出一個框架

出圖為，

邊邊的黑框就是所畫的邊界

利用以下代碼進行三維及光源的轉化，

axis vis3d; % 凍結屏幕高寬比，使得一個三維對象的旋轉不會改變坐標軸的刻度顯示，不凍結會變形 axis off; view(3); % 三維視角看圖 camlight % 光源 lighting gouraud; % 著色的模式

最后生成圖片

目前我還沒搞清楚為什么框架顯示是不完整的...

注：對于isovalue就是選擇一個基于其水平之下val的體數據，也就是U=-1，我們可以嘗試不同的J1的大小，選用不同的isovalue得到不同的能帶圖。

isovalue=1

還有比如fcc，fcc會比較復雜，我把代碼貼在下面，主要是要構造一個截角八面體的邊框，還有對于布里淵區的面方程的掌握，其他是一樣的。

% 面心立方等能面（定義倒格子常數為4） figure; [x,y,z]=meshgrid(linspace(-2, 2, 101)); x1=x(1,:,1); y1=y(:,1,1); z1=z(1,1,:); val=(cos(0.5*pi*x).*cos(0.5*pi*y)+cos(0.5*pi*y).*cos(0.5*pi*z)+cos(0.5*pi*z).*cos(0.5*pi*x)); % 在布里淵區之外的等能面設為非數 for i = 1:101for j = 1:101for k = 1:101if x1(i)+y1(j)+z1(k)-3>0 | x1(i)-y1(j)+z1(k)-3>0 | -x1(i)-y1(j)+z1(k)-3>0 | -x1(i)+y1(j)+z1(k)-3>0 | x1(i)+y1(j)-z1(k)-3>0 | x1(i)-y1(j)-z1(k)-3>0 | -x1(i)-y1(j)-z1(k)-3>0 | -x1(i)+y1(j)-z1(k)-3>0;val(i,j,k) = NaN;endendend end p=patch(isosurface(x,y,z,val,0.2)); set(p, 'FaceColor', 'yellow', 'EdgeColor', 'none');% 繪制面心立方的布里淵區：截角八面體 Vm1 = [1 0 2;0 -1 2;-1 0 2;0 1 2;1 0 -2;0 -1 -2;-1 0 -2;0 1 -2;1 2 0;0 2 1;-1 2 0;0 2 -1;1 -2 0;0 -2 1;-1 -2 0;0 -2 -1;-2 -1 0;-2 0 1;-2 1 0;-2 0 -1;2 -1 0;2 0 1;2 1 0;2 0 -1]; % 截角八面體的六個正方形頂點索引 k=1; Fm1 = zeros(6,4); for i = 1:6for j = 1:4Fm1(i, j) = k;k = k+1;end end % 截角八面體的八個正六邊形的頂點索引 Fm2 = [2 1 22 21 13 14;3 2 14 15 17 18;4 3 18 19 11 10;1 4 10 9 23 22;6 5 24 21 13 16;7 6 16 15 17 20;8 7 20 19 11 12;5 8 12 9 23 24]; % 繪制截角八面體的四個正方形 patch('Vertices', Vm1, 'Faces', Fm1, 'Linestyle', '-', 'LineWidth', 2, 'FaceVertexCData', rand(size(Fm1,1),1), 'FaceColor', 'none'); % 繪制截角八面體的八個正六邊形 patch('Vertices', Vm1, 'Faces', Fm2, 'Linestyle', '-', 'LineWidth', 2, 'FaceColor', 'none'); axis vis3d; axis off; view(3); camlight lighting gouraud;

isovalue=-0.2 fcc能帶圖

注：另外對于patch函數我也是一知半解，有時間我會去研究研究的！

參考文獻：

[1] 劉建軍,華厚玉.簡立方晶格s態電子等能面的計算機模擬[J].淮北煤炭師范學院學報,2005,9,26(3), 23-25

[2] 肖瑞春,陶松濤.緊束縛近似下立方晶系s態原子等能面的MATLAB模擬[J].上海工程技術大學學報.2012,9,26(3).

總結

以上是生活随笔為你收集整理的matlab patch函数_MATLAB实现紧束缚近似能带结构画图的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。