MATLAB做晶體結(jié)構(gòu)圖(固體物理)

MATLAB做晶體結(jié)構(gòu)圖(固體物理).md

寫在前面

最近在復(fù)習(xí)考研復(fù)試《固體物理》這一門課，去年學(xué)的內(nèi)容已經(jīng)忘干凈了，所以就翻開前幾頁。突然看到了面心立方和體心立方結(jié)構(gòu)圖，想到了去年室友用Mathematica做了晶胞的結(jié)構(gòu)圖，于是就手癢癢自己也想來做一個。

具體物理內(nèi)容不會涉及到多少，但還是要求大家能對“簡單立方結(jié)構(gòu)”、“體心立方結(jié)構(gòu)”、“面心立方結(jié)構(gòu)”有一個簡單的理解，因為我比較懶，所以我就不放這些基礎(chǔ)內(nèi)容湊字?jǐn)?shù)了。

用MATLAB跑出來的效果圖為

面心立方-MATLAB跑出來的圖

我的所有程序都放在我的Github: https://github.com/HanpuLiang/Something-Small中，點(diǎn)擊即可查看。

基本的思考

如果我們想要做一個類似于這樣子的結(jié)構(gòu)圖的話，我們需要知道些什么?

想要達(dá)成的效果(圖源百度)

如果我們想要做成這么大個超胞的話，又需要在上面的基礎(chǔ)上怎么做?

超胞示意圖(圖源我自己)

有過做計算的同學(xué)已經(jīng)想到了，需要計算物理中，描述晶胞中原子位置的文件POSCAR，然后把POSCAR拖到VESTA中就可以上面這幅圖了。

所以我們這里模擬POSCAR來輸入?yún)?shù)：

  • 晶胞參數(shù)：描述晶胞大小的參數(shù)。
  • 超胞大小：當(dāng)晶胞數(shù)量大于1個并且周期性變化時，用三個數(shù)字描述其在三個方向上的晶胞疊加數(shù)目。
  • 各個原子的坐標(biāo)：沒有這個還怎么畫出來原子啊。

基本參數(shù)

晶胞參數(shù)

晶胞參數(shù)包括6個值：三條邊a1, a2, a3和三個角\alpha, \beta, \gamma，如下圖所示

晶胞參數(shù)

因為編寫程序的復(fù)雜程度問題，我們這里只考慮\alpha\neq 90^o, \beta=\gamma=90^o這樣的情況。如果不這樣的話，那就有點(diǎn)難了，需要考慮這個晶胞斜向的角度。

超胞大小

超胞大小是用來形容我們這個超胞到底由幾個晶胞組成，以及他們的排列方式是怎么樣子的。就比如下圖中，沿著a1方向的層數(shù)為1個，沿著a2方向的層數(shù)為3，沿著a3方向的層數(shù)為2.所以我們就可以設(shè)定為[1, 3, 2]。

值得注意的是，這里我用的并不是x, y, z軸的方向?。。?/strong>因為如果\alpha\neq 90^o的話，那么a2方向就不和y軸平齊了，所以為了保證兩個晶胞相連，就必須要沿著晶軸方向拓展。

超胞大小

原子位置

這個沒什么好說的，如果我們建立好了超胞，那么直接在對應(yīng)的坐標(biāo)畫上原子就好。

開始寫程序

設(shè)定參數(shù)

首先我們來設(shè)定好上面的參數(shù)

%% 參數(shù)設(shè)定
global cell_size a1 a2 a3 alpha beta gamma
% 晶胞參數(shù)
a = [1, 1, 1];
% 三個角度
angle = [pi/2, pi/2, pi/2];
% 超胞大小
cell_size = [2, 2, 2];
% 簡單立方
position1 = [0, 0, 0; ...
           1, 1, 0; ...
           1, 0, 0; ...
           0, 1, 0; ...
           0, 0, 1; ...
           1, 1, 1; ...
           1, 0, 1; ...
           0, 1, 1];
% 體心
position2 = [0.5, 0.5, 0.5];
% 面心
position3 = [0, 0.5, 0.5; ...
             0.5, 0, 0.5; ...
             0.5, 0.5, 0; ...
             1, 0.5, 0.5; ...
             0.5, 1, 0.5; ...
             0.5, 0.5, 1];

[a1, a2, a3] = deal(a(1), a(2), a(3));
[alpha, beta, gamma] = deal(angle(1), angle(2), angle(3));

這里將體心和面心的原子分別提取出來，作圖的時候再放上去和簡單立方的一起做就好了。

作圖

我們?yōu)榱藞D像美觀就得做一些處理坐標(biāo)軸的事情

%% 作圖
figure
% 作圖設(shè)置
hold on, axis equal
axis image off
view(-37.5, 30)

然后我們就可以愉快的畫超胞的框架和各個原子啦，下面的兩個函數(shù)是我自定義的兩個函數(shù)

% 做超胞框架
plotBox();
% 做各個原子
plotAtoms(position2, [244, 13, 100]/255, 40);
plotAtoms(position1, [29, 191, 151]/255, 50);   %簡單立方的一定要放在最后面

簡單立方因為最大，所以我在里面設(shè)定了一些自動變坐標(biāo)軸大小的內(nèi)容，所以要放在最后面。

這樣子主程序就完成了。后面詳細(xì)解釋這兩個函數(shù)的內(nèi)容。

做超胞的框架plotBox

我們需要定出立方體的8個頂點(diǎn)，然后做出12條邊，這一部分很簡單，根據(jù)簡單的數(shù)學(xué)就可以推導(dǎo)出公式，然后寫出程序來。

function plotBox()
% 做邊框
    global a1 a2 a3 alpha beta gamma cell_size 
    % 超胞的邊長
    [A1, A2, A3] = deal(a1*cell_size(1), a2*cell_size(2), a3*cell_size(3));
    % 8個頂點(diǎn)
    vertex = [0, 0, 0;...
              A1, 0, 0;...
              A2*cos(alpha), A2*sin(alpha), 0;...
              A2*cos(alpha)+A1, A2*sin(alpha), 0;...
              0, 0, A3;...
              A1, 0, A3;...
              A2*cos(alpha), A2*sin(alpha), A3;...
              A2*cos(alpha)+A1, A2*sin(alpha), A3];
    % 12個邊
    plotLine(vertex(1,:), vertex(2,:))
    plotLine(vertex(1,:), vertex(3,:))
    plotLine(vertex(2,:), vertex(4,:))
    plotLine(vertex(3,:), vertex(4,:))
    plotLine(vertex(5,:), vertex(6,:))
    plotLine(vertex(5,:), vertex(7,:))
    plotLine(vertex(6,:), vertex(8,:))
    plotLine(vertex(7,:), vertex(8,:))
    plotLine(vertex(1,:), vertex(5,:))
    plotLine(vertex(2,:), vertex(6,:))
    plotLine(vertex(3,:), vertex(7,:))
    plotLine(vertex(4,:), vertex(8,:))
end

function plotLine(x1, x2)
% 做兩個點(diǎn)之間的框架線
    plot3([x1(1) x2(1)], [x1(2), x2(2)], [x1(3), x2(3)], 'k', 'linewidth', 1.3)
end

做各個原子的圖

這個才是重頭戲。

我們首先要確定出在超胞內(nèi)，一共有多少個原子。我們之前設(shè)置的超胞大小就派上用場了，我們通過三個循環(huán)嵌套在一起，遍歷出超胞內(nèi)的所有晶胞，然后將其原子位置加到矩陣中，最后統(tǒng)一作圖。

function plotAtoms(position, markercolor, markersize)
% 做各原子圖像
    global cell_size a1 a2 a3 alpha beta gamma
    % 原始晶胞
    % plot3(position(:,1), position(:,2), position(:,3), 'ok', 'linewidth', 1.5, 'markersize', 50, 'markerfacecolor', [29,191,151]/255)
    % 超胞
    % 遍歷得到超胞所有原子
    cur_point = position;
    for i1 = 1:cell_size(1)
        for i2=1:cell_size(2)
            for i3=1:cell_size(3)
                x_plus = a1*(i1-1) + a2*cos(alpha)*(i1-1);
                y_plus = a2*sin(alpha)*(i2-1);
                z_plus = a3*(i3-1);
                cur_point = [cur_point; [position(:,1)+x_plus position(:,2)+y_plus position(:,3)+z_plus]];
            end
        end
    end
    plot3(cur_point(:,1), cur_point(:,2), cur_point(:,3), 'ok', 'linewidth', 1.5, 'markersize', markersize, 'markerfacecolor', markercolor)
    % 設(shè)置坐標(biāo)軸大小
    [x_min, x_max, y_min, y_max, z_min, z_max] = deal(min(cur_point(:,1)), max(cur_point(:,1)), ...
                                                      min(cur_point(:,2)), max(cur_point(:,2)), ...
                                                      min(cur_point(:,3)), max(cur_point(:,3)));
    x_len = (x_max - x_min)/6;
    y_len = (y_max - y_min)/6;
    z_len = (z_max - z_min)/6;
    axis([x_min-x_len x_max+x_len y_min-y_len y_max+y_len z_min-z_len z_max+z_len])
end

然后這樣就完事了。是不是很簡單的。

寫在后面

以上所有代碼我都放在了我的Github中，可以通過點(diǎn)擊我的Github: https://github.com/HanpuLiang/Something-Small去查看。代碼下載后可直接運(yùn)行。

如果喜歡的話，麻煩點(diǎn)個關(guān)注，給個贊，加個收藏噢。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

  • 本文轉(zhuǎn)載自博主一個人就是一個疊加態(tài)，有部分刪減修改，文中對相關(guān)概念方法做了詳細(xì)的總結(jié)，留坑待填... 1. 第一原...
    chempeng閱讀 44,555評論 1贊 33
  • ****1、****做 ****XRD ****有什么用途啊，能看出其純度****?****還是能看出其中含有某種...
    鴨梨山大哎閱讀 12,617評論 0贊 5
  • I'm so angry at Bob for buying a new bike.我因鮑勃買了一輛新自行車而對他...
    bayue_9閱讀 129評論 0贊 0
  • 來平遙之前二哥說“去看又見平遙” 山西人就是這么簡單 直 直的讓你無話可說 同行的老鍋在大街上歇斯底里的和我大聲說...
    哈哈同學(xué)閱讀 2,222評論 0贊 0
  • 小城樓外相思雨。風(fēng)滿相思，雨滿相思。翠壺銀盤冷荔枝。 覷得浮云分半日。柳也入詩，人也入詩。古巷...
    麥片加冰閱讀 412評論 0贊 0

友情鏈接更多精彩內(nèi)容