• 淡江時報報導 Chi-Cheng 與合作者發現水在室溫常壓下也有結構，不是一般所認知無結構的液態。posted on Jul 25, 2025.

• Chi-Cheng was awarded IOP Trusted Reviewer status. posted on Oct 1, 2024.

• Up to 50 free e-prints of our recent work on NbSe2 can be downloaded via this link. posted on Jun 12, 2024.

• Chi-Cheng is selected as an outstanding referee by the American Physical Society (APS) in 2024. [淡江時報報導] posted on Mar 11, 2024.

• 賀Chin-En參加淡江大學理學院學生論文展獲得院長獎。posted on Jun 7, 2023.

• The 21th Workshop on First-Principles Computational Materials Physics will be held in National Sun Yat-Sen University, Kaohsiung, Taiwan on June 29-30, 2023. posted on Jun 7, 2023.

• 2023 NCTS-KIAS workshop on ab initio approaches to quantum materials wil be held in National Taiwan University, Taipei, Taiwan on April 27-29, 2023. posted on Mar 28, 2023.

• The 21st Taiwan-Japan-Korea symposium on strongly correlated electron systems (TJK21) will be held in NSRRC, Hsinchu, Taiwan on April 6-7, 2023. posted on Mar 28, 2023.

• The work on "Atomically thin metallic Si and Ge allotropes with high Fermi velocities" has been accepted for publication in Physical Review B. posted on Feb 24, 2023.

• Chi-Cheng 被選入全球前 2% 頂尖科學家之2021年度科學影響力排行榜。[淡江時報報導]. posted on Nov 12, 2022.

• The 20th workshop on first-principles computational materials physics will be held on 9/1 and 9/2. posted on Aug 10, 2022.

• 賀Chieh-Chun申請科技部大專學生研究計畫通過。posted on Jun 26, 2022.

• Chi-Cheng 被選入全球前 2% 頂尖科學家之2020年度科學影響力排行榜。[淡江時報報導]. posted on Dec 4, 2021.

• Chieh-Chun, Chang-Yu, Jhe-Hao, and Chin-En will present posters discussing optical properties of diamond, silicon, and cobalt and Raman spectrum in graphene/h-BN heterostructures, respectively, during the Physics Week held in Science Hall, Tamkang University between Nov 29 and and Dec 3, 2021. posted on Nov 28, 2021. 

• Chi-Cheng will present a talk, "固態材料裡的世界：磁單極", in S215, Science Hall, Tamkang University between 12:40 and 13:00 on Nov 29, 2021. posted on Nov 26, 2021. 

• 2021物理週-學術發表@淡江大學科學館11/29~12/03 (12:00~13:00 everyday)。歡迎各位同學前來共襄盛舉。posted on Nov 24, 2021.  

• The work on "Photocurrent-driven transient symmetry breaking in the Weyl semimetal TaAs" is just published online in Nature Materials (2021). [A view-only version is available here]. posted on Nov 9, 2021.  

• The work on "Modulating Chemical Composition and Work Function of Suspended Reduced Graphene Oxide Membranes Through Electrochemical Reduction" has been accepted for publication in Carbon. posted on Sep 15, 2021. [free access before Nov. 14, 2021]

• NTNU & Tamkang Physics Online Workshop will be held today. posted on Sep 7, 2021.

• The 19th workshop on first-principles computational materials physics (online workshop) will be held on Sep 2 and Sep 3, 2021. Most researchers working in the related fields in Taiwan will join this workshop. posted on Aug 31, 2021.

• The work on "The hidden competing phase revealed by first-principles calculations of phonon instability in the nearly optimally doped cuprate La1.875Sr0.125CuO4" has been accepted for publication in Physical Review B. posted on Aug 17, 2021.

• Chi-Cheng will present a talk, "Fermi arc, pseudo gap, and the LTLO Phase in Sr-doped La2CuO4", in an online seminar, NCTS beween 14:00 and 15:00 on Aug 12, 2021. posted on Jul 23, 2021.

• Chi-Cheng will present a talk, "First-principles calculations of phonon instability in the nearly optimally doped cuprate La1.875Sr0.125CuO4", in S215, Science Hall, Tamkang University between 14:10 and 15:00 on Apr 13, 2021. posted on Apr 7, 2021.         

• Chi-Cheng will present a talk, "Applications of density functional theory: from single-particle picture to many-body physics", in the 1st meeting room, 5th floor, Institute of Physics, Academia Sinica between 11:00 and 12:00 on Dec 17, 2020. posted on Dec 11, 2020.

• The news for our work on the flat band in 2-dimensional Ge layer can be found here. posted on Dec 11, 2020.

• The work on "Emergence of Nearly Flat Bands through a Kagome Lattice embedded in an Epitaxial Two-dimensional Ge layer with a Bitriangular Structure" has been accepted for publication as a Rapid Communication in Physical Review B. posted on Oct. 25, 2020.

• The work on "Partitioning interatomic force constants for first-principles phonon calculations: Applications to NaCl, PbTiO3, monolayer CrI3, and twisted bilayer graphene" has been accepted for publication in Journal of Physics: Condensed Matter. posted on Oct. 25, 2020. 

We are interested in the electronic-structure calculations using first-principles approaches and tight-binding model Hamiltonians. We are also interested in developing new methods that can be applied to the research related to condensed matter physics and materials science based on quantum mechanics. 

Focused topics are

‧ Electronic structures of real materials

‧ Lattice dynamics in real materials

‧ Phase transitions under high pressure 

‧ Magnetic, orbital, charge, and lattice structures in correlated materials 

‧ Local excitations of correlated materials

‧ Two-dimensional materials: Graphene, silicene, borophene, etc.

‧ Topological insulators and Weyl semimetals

‧ Exciton dynamics in weakly interacting systems and correlated Mott insulators

對Linux作業系統有基本認識與會撰寫簡單程式(或有高度興趣)的同學，歡迎加入我們的研究行列。

你可能會學習到的內容：

1. 如何描述金屬、絕緣體與半導體材料的幾何結構。

2. 如何描述金屬、絕緣體與半導體材料的電子能帶結構。

3. 了解新穎材料中的原子如何振動。

4. 探究較大尺度下凝態物質的幾何結構與電子結構。

5. 學習如何計算新穎材料中電子的激發態。

6. 你有聽過第一原理計算與密度泛函理論嗎? 想要學習更多與密度泛函理論相關的知識嗎?

7. 你有聽過目前科學家還沒找到磁單極嗎? 在固態物質中，物理學家發現可以找到與磁場行為一樣的物理量，而在動量空間中，磁單極是存在的!? 這與電子波函數的拓墣性質有關。想要多了解這方面的物理嗎?

8. 你有聽過費曼圖嗎? 量子力學中的微擾理論所給出的數學式子很複雜，但若用畫圖的方式來表達格林函數，數學與物理都變得有趣了。我們在了解固態系統的多體行為時，常常會用到費曼圖，想了解多體物理與費曼圖嗎?

我們主要在做些什麼? (orientation)

有聽過薛丁格方程式嗎？在學量子物理時，我們會給定一已知的位能來找出受到此位能影響的粒子所對應的能量本徵值與本徵函數。要了解固態系統中電子的能階，我們也是要寫下此系統的薛丁格方程式，然後對其求解。我們通常解的是與時間無關的薛丁格方程式，也就是能量本徵值問題。若此固態系統有週期性，其解就是你可能常聽到的能帶結構。固態物理告訴我們，大部分材料的性質都與能帶結構有關。因此，我們可以利用電腦，把想要研究的材料創建出來，也就是給定此材料是由哪些原子組成與原子位置，然後利用第一原理軟體讀取我們所給定的材料結構，解出該材料的能帶結構。以操作層面而言，我們需要在Linux作業系統中編輯一個文字檔案(input檔)，主要內容就是該材料的幾何結構。第一原理軟體不只能給我們能帶結構，還會輸出其他有用的資訊，例如：總能量、電荷密度、波函數、壓力等等。可以補充說明的是，有些同學可能會質疑固態材料是一個多體的系統，為何我們要解的還是像是本網頁首頁中看到的單電子薛丁格方程式? 這個答案是：其實我們解的是Kohn-Sham方程式，在利用密度泛函理論與一些近似下，多體的薛丁格方程式可以寫成單體的Kohn-Sham方程式，其中的對應關係，也還是目前值得研究的課題。有趣的是，教科書中的位能形式，例如：有限深位能阱，是固定不變的。但是在Kohn-Sham方程式中的位能是與解有關的，也就是說在求解的過程中，方程式的位能會一直被更新，直到新的解與位能項達到一致的情況時，計算就會停止，這樣求出來的解我們稱為自洽解。注意，該自洽解是針對某一材料幾何結構的自洽解，得到該解後，若想優化原子位置找出新結構，可以利用原子受力(量子力學中沒有力，你可以想像是能量對原子位移的微分或是期望值)給出一個新結構，然後再針對新的結構求出電子波函數的自洽解。目前的第一原理軟體，都可以完成這些操作，那我們為何還需要寫程式? 簡單回答：當所採用的第一原理軟體沒有辦法進行你想要的計算時，我們就需要自己發展新程式，自己的程式也可以幫你更有效率地準備input檔與讀取output檔的資訊。

1. 進入 https://www.openmx-square.org/viewer/index.html   

2. 點選 Examples->select V2O5-H2O.dat

3. 點選 Save->select OMX (frac) 後可下載範本檔案 abc.dat

4. 範本檔案中 # 開頭的是註記，不會被OpenMX讀取。

   我們可將內容修改如下：

# 因為 graphene 只有一種種類的原子，我們可將 Species.Number 改為 1。

# 我們也可以增加 C 進去，但是不刪除 H, O, V，這種作法需將 Species.Number 改為 4。

# 我們將採用只給 C 的作法。

# 參考 https://www.openmx-square.org/openmx_man3.9/node27.html 可設定如下：

Species.Number                     1

<Definition.of.Atomic.Species

  C     C6.0-s2p2         C_PBE19

Definition.of.Atomic.Species>

# 接下來要設定單位晶胞中的原子數目，與我們將選定用fractional coordinate表示原子位置。

Atoms.Number                       2

Atoms.SpeciesAndCoordinates.Unit   FRAC

# 接下來給定原子位置，最後還要給定你猜測每顆原子的 spin-up 與 spin-down 電子數目，

# 各原子電子數的總和應該要與 https://www.openmx-square.org/openmx_man3.9/node27.html 

# 中所給定的 Valence electrons 一致。沒有考慮磁性時，值各給一半即可。

<Atoms.SpeciesAndCoordinates

  1     C     0.3333333   0.3333333   0.5          2.0  2.0

  2     C     0.6666667   0.6666667   0.5          2.0  2.0

Atoms.SpeciesAndCoordinates>

# 接下來給定描述單位晶胞的單位，與晶格向量 (A,B,C) 在 x,y,z 的分量。

# Ax Ay Az

# Bx By Bz

# Cx Cy Cz

Atoms.UnitVectors.Unit             Ang

<Atoms.UnitVectors

    2.1304224933 -1.23 0.

    2.1304224933  1.23 0.

    0 0 10.

Atoms.UnitVectors>

# 各區塊的資訊在abc.dat中出現的順序並不重要，也可以先定義晶格向量後再給原子位置。

# 上面檔案的全部內容如下：

============ abc.dat ============

Species.Number                     1

<Definition.of.Atomic.Species

  C     C6.0-s2p2         C_PBE19

Definition.of.Atomic.Species>

Atoms.Number                       2

Atoms.SpeciesAndCoordinates.Unit   FRAC

<Atoms.SpeciesAndCoordinates

  1     C     0.3333333   0.3333333   0.5          2.0  2.0

  2     C     0.6666667   0.6666667   0.5          2.0  2.0

Atoms.SpeciesAndCoordinates>

Atoms.UnitVectors.Unit             Ang

<Atoms.UnitVectors

    2.1304224933 -1.23 0.

    2.1304224933  1.23 0.

    0 0 10.

Atoms.UnitVectors>

============ abc.dat ============

5. 將此檔案 abc.dat (從檔案總管、Finder、...)拉入你剛剛看到V2O5-H2O.dat的瀏覽器畫面中，

    你應該可以看到 OpenMX Viewer 的畫面已經變成剛剛我們所定義出的 graphene 結構。

    這時你可以試試下方的其他功能。 例如，將 Cells 打勾才可看到 unit cell。改變 Supercell 

    右方的數字，你應該可以看到用不同週期性呈現 graphene 結構的結果。若改變 Bond Factor

    會影響鍵結出現的數量(利用原子-原子間的距離判斷出是否畫出鍵結)。

1. 在練習產生 graphene 的結構時，我們建立了描述 graphene 的 abc.dat。

    我們可以加入一些新的參數。雖然這些參數出現的順序不重要，我們儘量採用與 

    OpenMX Viewer (Save->select OMX) 一致的順序。

============ abc.dat ============

System.CurrrentDirectory         ./

System.Name                     abc

level.of.stdout                   1

level.of.fileout                  1

DATA.PATH           /opt/openmx/openmx3.9/DFT_DATA19

Species.Number                     1

<Definition.of.Atomic.Species

  C     C6.0-s2p2         C_PBE19

Definition.of.Atomic.Species>

Atoms.Number                       2

Atoms.SpeciesAndCoordinates.Unit   FRAC

<Atoms.SpeciesAndCoordinates

  1     C     0.3333333   0.3333333   0.5          2.0  2.0

  2     C     0.6666667   0.6666667   0.5          2.0  2.0

Atoms.SpeciesAndCoordinates>

Atoms.UnitVectors.Unit             Ang

<Atoms.UnitVectors

    2.1304224933 -1.23 0.

    2.1304224933  1.23 0.

    0 0 10.

Atoms.UnitVectors>

scf.XcType                    GGA-PBE

scf.SpinPolarization          off

scf.ElectronicTemperature     300.0

scf.energycutoff              300.0

scf.maxIter                   300

scf.EigenvalueSolver          band

scf.Kgrid                     20 20 1

scf.Mixing.Type               rmm-diisk

scf.Init.Mixing.Weight        0.05

scf.Min.Mixing.Weight         0.01

scf.Max.Mixing.Weight         0.30

scf.Mixing.History            50

scf.Mixing.StartPulay         13

scf.criterion                 1.0e-8

scf.restart    on

scf.fixed.grid   0. 0. 0.

============ abc.dat ============

2. 假設你已經將OpenMX安裝在Linux作業系統下，且執行檔的完整路徑是

   /opt/openmx/openmx3.9/source/openmx 

   你可以直接執行下列指令進行計算

   /opt/openmx/openmx3.9/source/openmx abc.dat -nt 2 >& log &

   上面說明目前input檔的名稱是 abc.dat，我們採用2個核心來計算，

   OpenMX output檔案(像是 abc.out)以外的標準輸出會存入 log 檔。最後的 & 是指丟到背景算。

   如果你熟悉Linux作業系統，這些東西應該不陌生。你的input檔、log檔也可以取其他名稱。

   (如果你的執行檔執行時需要先載入一些環境設定，在執行     

    /opt/openmx/openmx3.9/source/openmx abc.dat 

    之前可能要先載入環境設定，例如

    source /opt/intel/oneapi/setvars.sh intel64

    /opt/openmx/openmx3.9/source/openmx abc.dat -nt 2 >& log &)

   當然，我們建議你用 PBS 等排隊系統丟工作到後端電腦做計算(如 pc cluster)，

   此時你需要學習如何利用 PBS script (test.pbs) 丟工作(qsub test.pbs)。

   test.pbs 可能包含：

   source /opt/intel/oneapi/setvars.sh intel64

   cd $PBS_O_WORKDIR

   mpirun /opt/openmx/openmx3.9/source/openmx abc.dat -nt 2 >& log

3. 計算完成後，可以在 abc.out 中搜尋Utot。Utot的值就是對應此結構的總能(Total Energy)，

   也就是熱力學中的內能U。如果你改變結構，例如，將晶格常數 Ax Ay Az 與 Bx By Bz 放大或

   縮小，則你可以比較不同幾何結構的總能量(Utot)。大自然喜歡能量低的結構，所以你所得到

   對應總能最低的結構就是你所預測出來的穩定結構。Utot 的單位是 Hartree。

   注意：在這個例子中我們不需要改變 Cx Cy Cz，因為採用有考慮週期性邊界條件的程式時，

   我們可以模擬沒有週期性的方法之一就是在沿著沒有週期性的方向上加一層真空層，目前的

   例子就是用一比較長的 c 軸來實現。通常真空層的厚度至少要 10 埃，甚至 20 埃以上。

   c 軸長度不可以過長，過長記憶體會不夠，通常是給 10~50 埃。

4. 你可以試著計算出FCC結構的矽或鑽石的總能與體積之關係(EV圖)。在計算之前，

   請先看一下我們對新增參數的說明。

System.CurrrentDirectory         ./   

# output 檔案會放在的資料夾，./ 表示當下的資料夾，一般不需要更改。

System.Name                     abc    

# output 檔的檔名會以這個名稱作開頭(abc.out)，如果你會在同一資料夾下跑很多計算，

# 則有需要改這名稱，不然先前同檔名的檔案內容會被新的計算所覆蓋。

level.of.stdout                   1  

level.of.fileout                  1

# 選定輸出(e.g. log and abc.out)內容的資訊多寡。

# 數字愈大輸出的資訊就愈多，但是檔案大小也就愈大，硬碟空間可能很快就不夠。

DATA.PATH           /opt/openmx/openmx3.9/DFT_DATA19

# 安裝OpenMX後，會有存放虛位勢與基底的資料夾，該路徑必須要讓OpenMX執行檔知道。

# 若路徑不正確，將無法執行。

Species.Number                     1

<Definition.of.Atomic.Species

  C     C6.0-s2p2         C_PBE19

Definition.of.Atomic.Species>

Atoms.Number                       2

Atoms.SpeciesAndCoordinates.Unit   FRAC

<Atoms.SpeciesAndCoordinates

  1     C     0.3333333   0.3333333   0.5          2.0  2.0

  2     C     0.6666667   0.6666667   0.5          2.0  2.0

Atoms.SpeciesAndCoordinates>

Atoms.UnitVectors.Unit             Ang

<Atoms.UnitVectors

    2.1304224933 -1.23 0.

    2.1304224933  1.23 0.

    0 0 10.

Atoms.UnitVectors>

# 以上已經在練習產生 graphene 的結構時說明過，當你給定原子位置時，

  其實也在建構 Kohn-Sham equation 的位能項。

scf.XcType                    GGA-PBE    

# 選定交換相干能的近似方法，此處採用 GGA 近似。通常我們會先採用 GGA 近似。

scf.SpinPolarization          off        

# 是否考慮有磁性? 若設為 on 則 spin-up 與 spin-down 的電子數目應該也要設成不一樣。

scf.ElectronicTemperature     300.0      

# 電子溫度(K)。因應k-points不夠密時，要做能量smearing時需要用到的參數，通常不用更改。

scf.energycutoff              300.0      

# 此數字(單位: Ryd)反映了描述實空間物理量之格點數目。數字愈大表示格子點愈密、數目愈多。

# 這個值通常會給 220 ~ 1000。300 通常足夠，1000 幾乎是上限。對應的 grids可在log中找到。

# 可用 scf.Ngrid 取代，scf.Ngrid 可直接給出沿 a, b, c 軸的格子點數目。

scf.maxIter                   300        

# 解 Kohn-Sham equation 時需要找自洽解，當找自洽解的步數超過此設定值時，程式會停止。

scf.EigenvalueSolver          band       # 執行對角化的方法，可不用更改。

scf.Kgrid                     20 20 1    #   number of grids along a*, b*, c*

# 每一個實空間的晶胞都有一對應的倒空間晶胞，此數字給出沿著倒晶格向量的格子點數目。

# 與scf.Ngrid 一樣，愈密會愈準，通常此參數與倒晶格向量長度成正比。最小為 1。

# 因此，模擬中若有沒有週期性的方向時，沿該方向的 k-point 應設為 1 。

# k-points的多寡需要測試，通常增加k-points若總能量變化小於 0.1 meV/atom時，不需再增加。

# 在能帶結構中的Fermi level附近有Dirac cone的情況下，可能需要更多的k-points。

# 倒空間晶格向量若用*表示：a* = 2 pi ( b x c ) / V, b* = 2 pi ( c x a ) / V, c* = 2 pi ( a x b ) / V。

# V = a dot ( b x c )。

# 下面 6 行是找自洽解會用到的參數，與如何混合新舊電荷密度有關，

# 有效調整這些參數可以加快收斂速度，但是一旦可收斂，其解應該是一致的，只影響時間。

# 除非計算無法收斂，通常不太需要更改，如何調整這些參數也與系統有關，需要經驗。

scf.Mixing.Type               rmm-diisk

scf.Init.Mixing.Weight        0.05

scf.Min.Mixing.Weight         0.01

scf.Max.Mixing.Weight         0.30

scf.Mixing.History            50

scf.Mixing.StartPulay         13

scf.criterion                 1.0e-8 

# 通常兩次計算出來的自洽解能量(eigenvalue 總和)差小於此值時(單位: Hartree)，計算會結束。

# 所以這個參數決定計算何時收斂，若步數超過 scf.maxIter (我們設為300) 還無法收斂，

# 有可能是給定的結構有錯誤或有需要調整上面的參數。

scf.restart    on  

# 若加上這一行，會先檢查有沒有先前計算的結果，如果有則會讀取先前的結果加快收斂的速度。

scf.fixed.grid   0. 0. 0. # 我們一般會加入這一行來固定實空間grids的參考原點。

1. 考慮有週期性的固態材料時，其電子能階與孤立原子中的電子能階不同。

    原子的能量是離散的，可用能階：En 表示。Bloch's theorem 證明在有

    週期性的位能下，描述電子的能量不只有量子數 n ，還多了波向量 k 。

    能帶結構的圖一般是指畫出 En(k)。通常我們會用橫軸表示 k ，縱軸表示

    能量來畫圖。不同的 n 可以用不同的 E(k) 曲線來表達。我們也習慣將

    OpenMX 輸出檔案 (abc.out) 中的 Chemical Potential 當作是參考的能量，

    也就是縱軸的零點 (須將能量平移)。也就是說，每個能階有多少電子佔據

    可由 Fermi-Dirac distribution 決定。溫度為零時，縱軸的正負值剛好告訴

    我們哪些能階是空的，哪些是填滿電子。若要計算 graphene 的能帶結構，

    可將上例 (Total energy) 中的 input 檔案再加入下面的參數。因為已經加入

    "scf.restart on"，再執行一次時，找自洽解的時間會變短，除非你刪除了

    先前的檔案。(先前的檔案放在資料夾 abc_rst 下，若不需要也可以刪除

    以增加硬碟空間：rm -rf ./abc_rst) 

Band.dispersion         on  # 內定值為off，改為on則會計算能帶結構

Band.Nkpath             3   

# 通常呈現能帶的方式是在倒空間中選定一些路徑，我們設定有 3 條路徑。

# 可以在 OpenMX Viewer 中點選 Brillouin zone 參考 OpenMX 所給出的

# 建議路徑。

<Band.kpath 

 26  0 0 0  0.5 0.5 0 G M   

# 決定第一條路徑要用幾個點畫圖，然後給定點的數目，之後給兩端點的座標。

# 此座標是以倒空間晶格向量為單位。最後，如果要用 OpenMX 提供的工具

# 畫圖，可給定這兩個端點在圖中出現的代號，也可隨意給定代號，如 A B 。

 15  0.5 0.5 0  0.333333333333 0.666666666666 0 M K     

# 設定第二段路徑所需的資訊。

 30  0.333333333333 0.666666666666 0  0 0 0 K G         

# 設定第三段路徑所需的資訊。

Band.kpath>

2. 執行下方完整的執行檔後，你會得到 abc.Band 的檔案。   

System.CurrrentDirectory         ./

System.Name                     abc

level.of.stdout                   1

level.of.fileout                  1

DATA.PATH           /opt/openmx/openmx3.9/DFT_DATA19 

# 請依你的情況更改 DATA.PATH 路徑

Species.Number                     1

<Definition.of.Atomic.Species

  C     C6.0-s2p2         C_PBE19

Definition.of.Atomic.Species>

Atoms.Number                       2

Atoms.SpeciesAndCoordinates.Unit   FRAC

<Atoms.SpeciesAndCoordinates

  1     C     0.3333333   0.3333333   0.5          2.0  2.0

  2     C     0.6666667   0.6666667   0.5          2.0  2.0

Atoms.SpeciesAndCoordinates>

Atoms.UnitVectors.Unit             Ang

<Atoms.UnitVectors

    2.1304224933 -1.23 0.

    2.1304224933  1.23 0.

    0 0 10.

Atoms.UnitVectors>

scf.XcType                    GGA-PBE

scf.SpinPolarization          off

scf.ElectronicTemperature     300.0

scf.energycutoff              300.0

scf.maxIter                   300

scf.EigenvalueSolver          band

scf.Kgrid                     20 20 1

scf.Mixing.Type               rmm-diisk

scf.Init.Mixing.Weight        0.05

scf.Min.Mixing.Weight         0.01

scf.Max.Mixing.Weight         0.30

scf.Mixing.History            50

scf.Mixing.StartPulay         13

scf.criterion                 1.0e-8

scf.restart    on

scf.fixed.grid   0. 0. 0.

Band.dispersion               on

Band.Nkpath                   3

<Band.kpath

 26  0 0 0  0.5 0.5 0 G M

 15  0.5 0.5 0  0.333333333333 0.666666666666 0 M K

 30  0.333333333333 0.666666666666 0  0 0 0 K G

Band.kpath>

3. abc.Band 內容可能如下，你可以利用 gnuplot 畫圖。參考

   https://www.openmx-square.org/openmx_man3.9/node68.html 

 16  0  -0.133956505561481

 0.780342567323407 -1.351592973906713  0.000000000000000  0.780342567323407  1.351592973906713 -0.000000000000000 -0.000000000000000  0.000000000000000  0.332491871581051 

3

26  0.000000000000000  0.000000000000000  0.000000000000000   0.500000000000000  0.500000000000000  0.000000000000000  A B

15  0.500000000000000  0.500000000000000  0.000000000000000   0.333333333333000  0.666666666666000  0.000000000000000  A B

30  0.333333333333000  0.666666666666000  0.000000000000000   0.000000000000000  0.000000000000000  0.000000000000000  A B

16  0.000000000000000  0.000000000000000  0.000000000000000

-0.865055214011247 -0.423041048682285 -0.254865102887435 -0.254819098403511  0.222664255807591  0.222767210253693  0.288037661924240  0.340759814014258  0.345500595480136  0.367329626479496  2.491311910537318  2.491333821829503  3.102589093474561  3.324052934461580  3.324324060007386  6.671316390468590 

16  0.020000000000000  0.020000000000000  0.000000000000000

-0.864708984046771 -0.422614583224164 -0.256306246595493 -0.255672495687420  0.222812818207331  0.224408717557488  0.287099770044477  0.341056315462985  0.346151316151122  0.367303127253048  2.413470627892193  2.476903864565122  3.090595831257491  3.285912608519409  3.325356807449333  6.681472728304298 

...

4. 上面這個檔案 (abc.Band) 包含了有沒有考慮磁性、化學勢、倒空間晶格

    向量、所計算的k-points與對應的能階等資訊。你也可以利用下面的 fortran

    code 將 abc.Band 轉換成可用 xmgrace 直接讀取的格式 (NXY)，但你需要

   自己修飾圖片，例如更改 x y 的範圍與自己標示 k 點的符號等等。你可以

   嘗試自己寫個程式產生資訊更完整的檔案或圖片。你也可以用EXCEL畫圖，

   轉檔後的格式為： 

   x1 y1 y2 y3 ...

   ...

   例如，

   gfortran getbands.f90

   ./a.out < abc.Band

   xmgrace -nxy band

============ getbands.f90 ============

Program change_format_OpenMX_Band_4_plot

integer::nband,nkpts,tmpi,nsec,spin

real(kind=8)::chem,pdis,dis

real(kind=8)::a(10000),rv(3,3),q1,q2,q3,p1,p2,p3

character ::tmp

read *,nband,spin,chem

if (spin.eq.0) then

open(11,file="band")

read *,rv(1,1),rv(1,2),rv(1,3),rv(2,1),rv(2,2),rv(2,3),rv(3,1),rv(3,2),rv(3,3)

read *,nsec

nkpts=0

do i=1,nsec

read *,tmpi

nkpts=nkpts+tmpi

enddo

p1=0.d0

p2=0.d0

p3=0.d0

read *,tmpi,q1,q2,q3

read *,(a(j),j=1,nband)

dis=0.d0

write(11,'(i5,10000f12.6)') 0,((a(j)-chem)*27.2113845d0,j=1,nband)

p1=q1

p2=q2

p3=q3

pdis=0.d0

do i=2,nkpts

read *,tmpi,q1,q2,q3

read *,(a(j),j=1,nband)

dis=dsqrt(((q1-p1)*rv(1,1)+(q2-p2)*rv(2,1)+(q3-p3)*rv(3,1))**2+((q1-p1)*rv(1,2)+(q2-p2)*rv(2,2)+(q3-p3)*rv(3,2))**2 &

+((q1-p1)*rv(1,3)+(q2-p2)*rv(2,3)+(q3-p3)*rv(3,3))**2)

write(11,'(f14.8,1x,10000f12.6)') dis+pdis,((a(j)-chem)*27.2113845d0,j=1,nband)

p1=q1

p2=q2

p3=q3

pdis=dis+pdis

enddo

close(11)

else

open(12,file="band_up")

open(13,file="band_dn")

read *,rv(1,1),rv(1,2),rv(1,3),rv(2,1),rv(2,2),rv(2,3),rv(3,1),rv(3,2),rv(3,3)

read *,nsec

nkpts=0

do i=1,nsec

read *,tmpi

nkpts=nkpts+tmpi

enddo

p1=0.d0

p2=0.d0

p3=0.d0

dis=0.d0

read *,tmpi,q1,q2,q3

read *,(a(j),j=1,nband)

write(12,'(i5,10000f12.6)') 0,((a(j)-chem)*27.2113845d0,j=1,nband)

read *,tmpi,q1,q2,q3

read *,(a(j),j=1,nband)

write(13,'(i5,10000f12.6)') 0,((a(j)-chem)*27.2113845d0,j=1,nband)

p1=q1

p2=q2

p3=q3

pdis=0.d0

do i=2,nkpts

read *,tmpi,q1,q2,q3

read *,(a(j),j=1,nband)

dis=dsqrt(((q1-p1)*rv(1,1)+(q2-p2)*rv(2,1)+(q3-p3)*rv(3,1))**2+((q1-p1)*rv(1,2)+(q2-p2)*rv(2,2)+(q3-p3)*rv(3,2))**2 &

+((q1-p1)*rv(1,3)+(q2-p2)*rv(2,3)+(q3-p3)*rv(3,3))**2)

write(12,'(f14.8,1x,10000f12.6)') dis+pdis,((a(j)-chem)*27.2113845d0,j=1,nband)

read *,tmpi,q1,q2,q3

read *,(a(j),j=1,nband)

write(13,'(f14.8,1x,10000f12.6)') dis+pdis,((a(j)-chem)*27.2113845d0,j=1,nband)

p1=q1

p2=q2

p3=q3

pdis=dis+pdis

enddo

close(12)

close(13)

endif

end

============ getbands.f90 ============

5. 比較你的能帶結構是否與網路上找到的類似？有在 K 點看到 Dirac point？

    getbands 已經把 Chemical Potential 平移到縱軸的零點。能量的單位是 eV。

1. 能帶結構圖一般只能呈現沿著倒空間中某些路徑的電子能階，無法展現所有

    k-points 的電子能量。二維系統可以採用 3d 繪圖的方式呈現所有 k-points 

    的能帶。例如，可用 xy 平面表示二維的倒空間，能量則用 z 軸呈現。此時， 

    能帶可用曲面的方式呈現出來。但是有時我們不關心個別 k-point 的能階，

    而是關心考慮所有 k-points 後，對應某個能量 E 的電子態數目，這就是

    態密度。態密度的單位是單位能量有多少電子態，若每個電子態都填滿著

    一顆電子，態密度反映了單位能量有多少顆電子。計算時，因為倒空間晶胞

    也有週期性，所以我們只需要考慮最小倒空間晶胞中的 k-points。

2. 要計算態密度可以在 abc.dat 中加入下面參數，然後再跑一次 OpenMX：

DosGauss.fileout             on           

DosGauss.Num.Mesh          1000         

DosGauss.Width             0.10         

Dos.Erange           -25.0  5.0        

Dos.Kgrid             200 200 1      

3. 說明如下：

DosGauss.fileout             on # 設定將執行態密度計算

DosGauss.Num.Mesh          1000 

# 給定態密度計算的能量範圍中要考慮的能量點數目，這會影響繪圖解析度。

DosGauss.Width             0.10 

# 我們不可能執行無窮密之 k-points 的計算，所以需要採用 Gauss smearing。

# 較密的 k-points 可選較小的 Gauss width，能量解析的程度也愈好，但耗時。

# 較稀疏的 k-points 需要用較大的 Gauss width 來美化曲線，避免鋸齒狀。

Dos.Erange           -25.0  5.0 # 設定計算態密度的能量範圍。單位是 eV。

                                              # 注意：Chemical Potential 為能量零點。

Dos.Kgrid             200 200 1 # k-points。通常與 scf.Kgrid 成正比。

# 態密度計算不需要重新找自洽解，可以給較密的 k-points。但是愈密愈耗時。

# Graphene 需要很密的 k-points 才可計算出平滑的態密度曲線。

4. 計算完畢後會產生 abc.Dos.val、abc.Dos.vec、abc.FermiSurf0.bxsf 檔案。

   abc.FermiSurf0.bxsf 可用 xcrysden 或其他可讀 bxsf 格式的軟體開啟，此檔案

   含有 Fermi surface 的資訊。若想畫解析度高的費米面，需要調密 k-points。

   若想要畫態密度，先檢查安裝完 OpenMX 的資料夾下是否有 DosMain。

   (ls /opt/openmx/openmx3.9/source/DosMain)

   若無，你可能需要在 /opt/openmx/openmx3.9/source 下執行下列指令產生 

   DosMain 執行檔：

   make DosMain

   然後你可執行下列指令：

   /opt/openmx/openmx3.9/source/DosMain abc.Dos.val abc.Dos.vec

   你可能會看到

   Max of Spe_Total_CNO = 8

   1 1 101 102 103 101 102 103 

   <abc.Dos.val>

   <abc>

   The input data caluculated by the Gaussian broadening method

   Do you want Dos(1) or PDos(2)? 

   (選 1 產生態密度，選 2 可以產生選定原子之各軌域的貢獻。)

   1 

   <Dos_Gaussian> start

   <Dos_Gaussian> make abc.DOS.Gaussian

5. abc.DOS.Gaussian 中存放著三個資訊：

   能量、態密度、從給定能量範圍起點開始積分之電子數目。

   單位是：eV、# of states/eV per unit cell)、electrons。

   你可以用 xmgrace 一次畫出三條曲線：xmgrace -nxy abc.DOS.Gaussian。 

1. 在"材料的總能量計算與結構預測"的練習中，我們需要手動調整原子位置與

    晶格常數。有時這進行這樣的計算是必要的，因為我們可以掌控材料的幾何

    結構。但是很多情況，這樣手動的調整原子位置與晶格常數是很費時的，

    尤其是如果我們只想求得穩定結構的時候 (或結構是未知的，沒有實驗數據)。

    我們之前的計算是要求結構不動，其實是採用預設值：MD.Type  Nomd。   

2. 晶格常數不變，原子位置可改變。

   若要靠 OpenMX 自動尋找原子受力為零的結構，可新增下列參數：

MD.Type                     Opt      

MD.Opt.DIIS.History           5       

MD.Opt.StartDIIS             30       

MD.Opt.EveryDIIS            200       

MD.maxIter                  300       

MD.Opt.criterion          0.0001      

3. 參數說明如下：

MD.Type                     Opt  # 將 Nomd 改為 Opt 可以讓原子自動移動

# 除了 Opt 外，還可以採用不同的演算法，請參考:

# https://www.openmx-square.org/openmx_man3.9/node49.html    

MD.Opt.DIIS.History           5        

MD.Opt.StartDIIS             30        

MD.Opt.EveryDIIS            200        

# 上面三個參數可以加速尋找原子受力為零的位置，但是一旦找到該位置，

# 原子受力就為零，所以除非太耗時，一般並不需要特別調整這些參數。

MD.maxIter                  300        

# 這和 scf.maxIter 類似，是指若尋找超過 300 個結構還找不到受力為零的

# 結構，計算會停止。注意：每一個結構都會進行找自洽解的計算，因為針對

# 每一個結構都要找到自洽解後才能計算每顆原子的受力。因此，MD.maxIter

# 與 scf.maxIter 不一定要設成一樣的值。

MD.Opt.criterion          0.0001 

# 通常會給 0.0001 ~ 0.0003，對於較難收斂的系統，有時會設為 0.0005。      

# 原子受力數值上不可能真的為零，因此，需要設定 MD.Opt.criterion，

# 當每顆原子受力都小於 MD.Opt.criterion Hartree/Bohr 時，計算就算完成。

4. 原子受力的資訊可在兩個地方找到。一個是 abc.out，例如我們之前計算 graphene 的 abc.out (Nomd)：

***********************************************************

***********************************************************

       xyz-coordinates (Ang) and forces (Hartree/Bohr)

***********************************************************

***********************************************************

<coordinates.forces

  2

    1     C     1.42028   0.00000   5.00000   0.000286842936  0.000000498979 -0.000000000009

    2     C     2.84056  -0.00000   5.00000  -0.000286842935 -0.000000498970  0.000000000013

coordinates.forces>

   但是 abc.out 只會紀錄最後結構中每顆原子沿著 x y z 的受力。若想知道所有

   計算過的結構與原子受力，可在 abc.md 中找到。值得一提的是，每一次更新

   的結構會存在 abc.dat# 中。OpenMX 不會更改你的 input 檔，因此 OpenMX

   自動產生的新檔案會加上 # 號以示區別。如果你最後不需要保存一開始的

   input 檔，你可以執行 cp abc.dat# abc.dat 覆蓋過去。

5. 若要限制某個原子沿某個方向是不可被移動的，可以加入下列參數：

<MD.Fixed.XYZ

    1  0 0 1

    2  0 0 1

MD.Fixed.XYZ>

# 最前面是原子的編號，後面是告訴 OpenMX 該原子是否可以沿著 x y z 移動。

#  0 是可移動，1 是不可移動。

6. 晶格常數與原子位置都可改變。

   有時我們需要固定體積 (Opt)，但有時我們只想找出總能量最低的結構。

   要讓晶格常數也可變動，只要將 Opt 改成 OptC5 ：

MD.Type                     OptC5      

MD.Opt.DIIS.History           5       

MD.Opt.StartDIIS             30       

MD.Opt.EveryDIIS            200

MD.maxIter                  300

MD.Opt.criterion          0.0001

   請參考：https://www.openmx-square.org/openmx_man3.9/node54.html    

   採用 OptC5 時，OpenMX 會尋找 Stress 也為零的結構。Stress 的資訊可以

   在 log 檔中找到。要注意的是，OptC5 不是萬能的，通常最後找到的結構與

   一開始給的結構有密切的關係，系統會被限制在能量的相對低點中。因此，

   一開始的 initial structure 通常決定了最後找到的結構。為了更有效的進行

   structural optimization，我們還可以加入對 unit cell 的 constraint：

<MD.Fixed.Cell.Vectors

     0 0 1

     0 0 0

     0 0 0

MD.Fixed.Cell.Vectors>

   這與限定原子是否可被移動類似，0 代表可動，1 代表固定。

   上面參數是指對晶格常數

   Ax Ay Az

   Bx By Bz

   Cx Cy Cz

   的限制。模擬二維材料時，我們可不改變 c 軸，因此可選擇固定 Cx Cy Cz。

7. 你也可以給定想要模擬的壓力(單位：GPa)。給定後，OptC5會考慮此壓力。

MD.applied.pressure         10.0       

    或是，你可以手動改變 unit cell 的大小，僅採用 opt，此時你可以加入：

scf.stress.tensor  on

    要求 OpenMX 即使在不動晶胞大小的情況下，也計算 Stress tensor。

8. 如果你是用 scf.energycutoff 來決定計算用的 grids 數目 (不是用 scf.Ngrid)，

    要注意 OptC5 最後產生的結構因為大小改變，對應相同 scf.energycutoff 的

    grids 數目可能也會改變。因此，為了讓 scf.energycutoff 的參數有意義，

    你需要再重新跑一次計算：首先，你可以用最後的結構覆蓋掉原始的結構

     cp abc.dat# abc.dat

     vi abc.dat

     qsub test.pbs (or /opt/openmx/openmx3.9/source/openmx abc.dat >& log)

     vi abc.dat 的部分是建議把 OpenMX 在檔案最後新增的幾行參數都刪除掉：

     例如，你可能發現檔案多了：

     geoopt.restart    on (刪除)

     MD.Current.Iter  27 (刪除)

     若沒刪除，OpenMX 會繼續之前的步數繼續跑下去，如果你是跑到 

     MD.maxIter 的上限 (例如，300) 而停止，若想要再從最新的結構繼續尋找

     新結構，可能會因為 MD.Current.Iter 300 讓程式只跑一步就停止。

     因為 grids 改變，OptC5 也需要多跑幾步找出更準確的結構，你可能需要

     重複上述步驟直到 cp abc.dat# abc.dat 後，只要再跑一次就能完成計算。

1. 計算聲子頻譜之前，你必須找到一個原子受力為零 (<0.0001 Hartree/Bohr) 的結構，並且要確認

   你的計算參數可讓總能量與力之值收斂，例如：kpoints (scf.Kgrid) 的數目必須夠多。

2. 因為聲子計算是假設原子間由彈簧相連接，所以需要力常數的資訊，而力常數可以利用在超晶胞中

   移動原子，然後將 abc.out 中所列出的原子受力 ( <coordinates.forces ... ) 除以位移大小而得到 

   (基本上就是 F = - k x，k = - F / x)。這樣的超晶胞計算需要對各種不同的原子進行位移，因此需要許多

   不同的 input 檔案。我們建議用 scf.Ngrid 來做聲子計算，若你是用 scf.energycutoff 來設定 grids 數目，

   你可以在 log 檔中查到 scf.Ngrid 的值，例如：Num. of grids of a-, b-, and c-axes = 21, 21, 105 。

   你可以執行

   /data2/LEE/shared/phonon_bin/geninput4openmx_grids phonon.input 21 21 105

   幫你產生所有需要的 input 檔案。你需要提供 phonon.input 檔案。21 21 105 是 scf.Ngrid 的值。

   phonon.input 的內容如下：

======== phonon.input ========

abc.dat # 提供受力為零之結構的 OpenMX input 檔案 (例如：cp abc.dat# abc.dat)

abcphonon # 提供你希望 geninput4openmx 幫你自動產生聲子計算的檔案名稱

2 2 1 # 超晶胞 (supercell) 大小，例如：2 x 2 x 1 。愈大愈好，但是算得愈慢。

0.01 # 原子位移大小 (in Ang)

6  # 提供 abc.dat 中 <Definition.of.Atomic.Species 列出之所有種類原子的原子序

4  #  要計算的 q points 數目

 0 0 0  # 列出 q points，和能帶結構一樣，是以倒晶格向量當作單位。

 0.5 0.5 0

 0.333333333333 0.666666666666 0

 0 0 0

100  # 在上述的總路徑中，你期望切割多少 q points (畫圖用)。程式會給你一個接近你期望的數字。

addbefore mpirun /opt/openmx/openmx3.9/source/openmx    # 參考下方說明

addafter -nt 4 >& log

======== phonon.input ========

   最後兩行 addbefore 與 addafter 會在螢幕輸出你執行 qsub test.pbs 時，可能需要在 test.pbs 檔案中

   插入的敘述。例如，你的 test.pbs 可能如下：

======== test.pbs ========

#!/bin/bash

#PBS -N jobname

#PBS -e para.err

#PBS -o para.log

#PBS -V

#PBS -q dell-3

#PBS -l select=1:ncpus=32:mpiprocs=8:ompthreads=4,walltime=9999:00:00

source /home/pkg/intelAPI_2021/setvars.sh --force

cd $PBS_O_WORKDIR

mpirun /opt/openmx/openmx3.9/source/openmx abc.dat -nt 4 >& log

======== test.pbs ========

   若你如上設定 addbefore 與 addafter 的敘述，geninput4openmx 可能會輸出：

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0000  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0001  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0002  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0003  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0004  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0005  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0006  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0007  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0008  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0009  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0010  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0011  -nt 4 >& log

 mpirun /opt/openmx/openmx3.9/source/openmx abcphonon_0012  -nt 4 >& log

 (會在自動產生的 input 檔名前後幫你加上你給的 addbefore 與 addafter 參數。)

   你需要自行複製這幾行或用 >> 將其導入至 test.pbs 中，然後再依你的需求修改 test.pbs。

   產生的 abcphonon 檔應該要重複你的 abc.dat 檔，abcphonon_0000 檔會是沒有移動原子的超晶胞檔案，

   剩下編號大於零的檔案是有移動過原子的 input 檔案。若 abcphonon_0000 所算出的原子受力很大，

   可能是你一開始計算所用的參數不夠好 (abc.dat 中的參數可能需要再調整)。例如，可能是 kpoints 數目

   用反比轉換給 supercell 時，因為不能整除而產生計算上的誤差，或是 scf.energycutoff (or scf.Ngrid) 

   的值可能要再增加。

3. 當你完成 步驟2 的計算後 (qsub test.pbs)，需要將 *.out 中的 force 資訊整合給 phonon 程式以建構

   動力學矩陣，接著就可以計算聲子頻譜。phonon 是一個獨立的程式，需要 cell.dat 與 force.dat 兩個 

   input 檔案。你可以執行：

   /data2/LEE/shared/phonon_bin/geninput4phonon phonon.input cell.dat force.dat

   geninput4phonon 會依據你先前的 phonon.input 幫你產生 cell.dat 與 force.dat 這兩個檔案。之後若需要

   更改 q-point 的路徑，可以直接修改 cell.dat 最後面的 q-point list (不需要再執行 geninput4phonon)。

4. 執行

   /data2/LEE/shared/phonon_bin/phonon cell.dat force.dat 3 1 1 1

   可以產生聲子頻譜的結果，與能帶結構一樣，輸出的格式為 x y1 y2 y3 ...。

   你可以將結果導入 plot.dat 然後畫圖 (例如，xmgrace -nxy plot.dat)。你可以試著計算 graphene 的 

   phonon dispersion，並比較你的結果是否與網路上找到的一致。頻率的單位是 (1/cm)。

   (/data2/LEE/shared/phonon_bin/phonon cell.dat force.dat 3 1 1 1 > plot.dat)

   上面的 3 1 1 1 (p1 p2 p3 p4) 是指：

   p1: 可設定力常數隨距離衰減的情況 ~ 1/r^(p1)，可慢慢增加，看看是否可以消去 Gamma 點附近的虛頻。

   p2: 1 代表會強加 sum rule，讓 Gamma 點會有三個值為零的頻率，0 則不考慮 sum rule。

   p3: <=1，設定上述 p1 會影響的範圍，通常不需要調整。

   p4: >=1，設定上述 p1 會影響的範圍，通常不需要調整。

   p1, p3, p4 不會改變相稱 (commensurate) q points 的頻率，只會改變本來就是內差出來的頻率值。

   相稱的 q points 是指其對應的波可以被所選定用來計算的超晶胞來描述 (週期性必須吻合)。

1. 若執行 (參考前一個練習)

   /data2/LEE/shared/phonon_bin/phonon cell.dat force.dat 3 1 1 1

   可計算出聲子在所給定 q paths 上的頻率。若你想分析某個 q point 上對應某個頻率的

   振動模式，可以執行

   /data2/LEE/shared/phonon_bin/phonon_one_q cell.dat force.dat 3 1 1 1 0.5 0.5 0.0

   最後三個參數 0.5 0.5 0.0 是指 q=(q1, q2, q3)，以倒晶格向量為單位。你可能會得到：

===========================

q vector: 0.5 0.5 0

===========================

Frequency in units of cm^-1

1 : 482.006

2 : 635.388

3 : 659.93

4 : 1375.55

5 : 1397.49

6 : 1460.72

   執行 phonon_one_q 的目的是獲得各個振動模式的編號 (目前的例子是 1~6 號)。

2. 若你想知道第 i 個振動模式的本徴位移(eigen-displacement)，你需要提供振動模式之編號 (i)、

   要呈現此振動模式的超晶胞大小 (超晶胞應該要能表現出 q vector 的週期性)、位移向量的長度 (值給的

   愈大，位移向量就愈長)。例如 i = 2，超晶胞大小選 2 x 2 x 1 (for q = 0.5 0.5 0.0)，你可執行：

   /data2/LEE/shared/phonon_bin/phonon_displacement_xsf cell.dat force.dat 3 1 1 1 0.5 0.5 0.0 2 2 2 1 0.1 > abc.xsf

   我們把輸出結果導入 abc.xsf 檔，abc.xsf 可以用 xcrysden (Display -> Forces)、VESTA 等軟體開啟。  

3. 另一種呈現振動模式的方式是執行：

   /data2/LEE/shared/phonon_bin/phonon_displacement_md cell.dat force.dat 3 1 1 1 0.5 0.5 0.0 2 2 2 1 10 20 0.2 > abc.md

   前面的參數與 phonon_displacement_xsf 一致，最後三個參數是：要用幾個週期呈現此振動模式 (10)、

   每個週期要分割成多少步來呈現 (20)、原子位移大小 ( 0.2，數字愈大原子位移就愈大)。

   此 10 x 20 = 200 個結構畫面可將 abc.md 檔拉至 OpenMX Viewer 以動畫的方式來呈現。