本文版權歸ponychen所有，請關注他的GitHub的主頁獲取腳本！
在材料力學強度的表徵測試中，理想拉伸曲線是其中典型的一種方式。從拉伸曲線中我們可以獲得彈性模量，屈服強度等等有用的力學係數。使用第一性原理模擬理想拉伸(不局限於拉伸，但要注意，不要將宏觀應力應變與DFT計算的理想應力應變混淆)曲線當然是一種很fancy的方法，但是目前我沒有看到一個自動化實現這個過程的腳本。我用bash寫了一個idealdeform.sh腳本，可以自動化實現理想拉伸或者剪切過程，目前主要用於VASP。
Nanoscale, 2017, 9, 850–855
還是先講講原理吧。腳本原理相當簡單，我主要先講一下兩類應變和應力。我們知道，我們平時主要用的是工程應變和工程應力。這類應變應力定義如下
各項符號我就不過多解讀了，這個應該是基礎的東西。我們注意到，l_0 和A_0都是樣品在未拉伸時的初始狀態，這也就導致工程應變和工程應力並非材料實際所受的狀態，因此我們日常看到的實際樣品的工程應力應變曲線往往都是這種低碳鋼形狀的1。
材料在B點開始屈服，D點開始發生頸縮，E點發生斷裂。但實際上材料發生頸縮後，材料繼續加工硬化，真實應力實際上是繼續上升，直到斷裂。為了更加準確地描述這個過程，有時候大家也會使用真實應變和真實應力來繪製曲線。下面我們來簡單推導一下兩者得表達式。
考慮塑性變形的不可壓縮性，我們有：
我們將當前面積A的定義代入到真實應力的定義式中即得真實應力和工程應力的換算關係。
真實應變就更加簡單了，我們只需要對每個瞬時應變進行積分即可
但實際上這種真實應變與工程應變的換算關係僅存在於頸縮初期之前，因為後期的時候材料的不可壓縮的假設已經崩潰。因此實際材料變形曲線還需進一步校正，不過本文就不再過多討論了。再次強調，我只是介紹了兩類應力應變，請不要將宏觀應力應變和理想應力應變聯繫起來。
下面我開始講一下腳本的使用方法，原理就不細講了，根本在於假設應變過程足夠慢，可以看作準靜態應變，這樣我們只需按一定應變間隔施加應變矩陣，優化晶胞和原子即可。首先我們需要準備優化晶胞所需要的VASP四個輸入文件。注意，POSCAR的原子坐標必須是分數坐標 。INCAR中ISIF可以用4，但請在OPTCELL中關閉你變形方向的弛豫 (OPTCELL是什麼？自己去DET群文件搜搜吧)。最好不要用3，晶格的體積發生變化會導致曲線失真。然後將idealdeform.sh放到同一文件夾。我測試的體系是fcc Ni, a為3.5058。
現在我們來修改腳本中相關參數。
orientation表示施加應變的方式，你可以選擇XX, YY, ZZ, XY, XZ, YZ六種變形方式中的一種，前三種為拉伸，後三種為剪切。initial為施加的初始應變，step表示每一步應變的間隔，num表示從初始應變開始按照step間隔一共連續施加了多少次應變。圖中的例子就表示對fcc Ni沿著X軸拉伸，初始應變為0，每隔0.01施加應變進行一次拉伸，一共施加了100次應變，最後施加的應變為1.000(100%)。mpiexec表示你的運行語句，由你的系統決定。最後直接運行腳本即可。
腳本會告訴你目前進行的進度，我們只需靜待完成。計算完成後，程序默認會用gnuplot畫一個工程應力應變曲線。
單位為GPa。我已經重新修改過單位，證號就表示拉應力。如果你覺得這樣子看很辣眼睛，你可以將plot.sh也複製到當前文件夾，然後運行。
腳本運行完後，工程應力應變曲線的數據會被保存到engineeringstressstrain.all中，真實應力應變曲線的數據會被保存到truestressstrain.all中。如果某一步中出現無法收斂的情況，程序會停止。你應該在OUTCAR中查找原因，並從當前階段繼續運算。
下面說幾個注意事項，非常重要：
大部分時間你並不需要用到真實應變，尤其當你的體系存在表面的時候，你不用去管真實應力應變曲線DFT獲得的應力應變曲線跟宏觀應力應變曲線沒有任何可比性，使用者應該清楚自己使用這個腳本獲得的應力應變曲線到底應該如何理解。宏觀應力應變來源於位錯運動，跟DFT結果完全不是一個概念，我之所以在開頭描述了一堆主要是為了介紹真實應力應變