熱膨脹係數實驗數據分析（bashvaspShengBTE自動計算材料熱運輸性質腳本）

2023-09-21 09:55:25 1

腳本分享

之前我們已經詳細介紹了ShengBTE軟體與其安裝教程，詳細請查看ShengBTE簡單介紹、安裝與使用

由於在計算三階力學常數時，根據擴胞比例以及近鄰原子受力情況的不同，生成的3RD.POSCAR數量也盡不相同，筆者最近一次的計算總共計算了700 個任務，如果在此過程中手動布置任務則費時費力，效率低下。而在後續的ShengBTE計算時，通過bash腳本進行任務布置也可以實現多任務順序計算，節省精力與時間。

本文針對材料的熱輸運性質計算主要參考天璣算所提供的教程，詳細請見

https://phadwiki.com/page237.html?article_id=79&pagenum=all

這裡我們使用兩個原子的矽原胞做演示

Si 5.38930000000000 0.0000000000000000 0.5071343999939496 0.5071343999939496 0.5071343999939496 0.0000000000000000 0.5071343999939496 0.5071343999939496 0.5071343999939496 0.0000000000000000 2Direct 0.8750000000000000 0.8750000000000000 0.8750000000000000 0.1250000000000000 0.1250000000000000 0.1250000000000000

對其進行結構優化後擴胞並使用vasp結合phonopy進行聲子譜計算，這裡選擇擴胞比例為4 4 4，以得到二階力學常數備用。

下圖為計算所得聲子譜與聲子態密度，不存在虛頻。

有關phonopy的使用方法請參見官網

https://phonopy.github.io/phonopy/examples.html

以後如果有時間也會整理一些聲子譜或者熱膨脹係數、格林奈森常數等計算的案例。

在計算三階力學常數時，通過執行thirdorder命令生成擴胞文件（thirdorder已添加環境變量），這裡我們選擇擴胞比例為2 2 2，擴胞後原子總數為8 8=16，近鄰原子數選為12（這裡的計算選擇僅做參考，實際計算時需要自行調整參數等進行測試並選取可行方案）

#擴胞命令模版 thirdorder_vasp.py sow x y z -c (x,y,z是三個方向的擴胞倍數，c是考慮多少個臨近原子的受力來計算力常數矩陣) thirdorder_vasp.py sow 2 2 2 -12

根據顯示結果這裡有72個DFT任務，相對來說較少，但是如果手動一個任務一個任務提交也是非常麻煩且耗時，這裡我們通過一個腳本來整理文件，並一次性安排所有的計算任務。在執行腳本之前應當提前準備好INCAR、KPOINTS、POTCAR，並存放在當前目錄。其中INCAR為自洽運算，KPOINTS應當與計算聲子譜的任務相同。

mkdir 3RD && mv 3RD* 3RD #這裡將所有thirdorder生成的POSCAR文件統一存放在3RD文件夾內mkdir to-run #建立計算文件夾for i in {01..72} #總計72個計算任務，這裡對i賦值。如果總任務數是三位數的，開頭是｛001..do #開始for循環，總共會執行72個循環（這裡的循環次數由上一行｛｝內確定）mkdir to-run/$i #建立單一任務計算的文件夾 cp 3RD/3RD.POSCAR.$i to-run/$i/POSCAR #將每個thirdorder生成的任務文件複製到對應的計算文件夾並重命名為POSCAR cp INCAR KPOINTS POTCAR to-run/$i/ #準備其他計算文件，這裡已經默認準備好了。 cd to-run/$i/ #進入單一任務計算文件夾 mpirun -np 24 vasp |tee runlog #執行vasp計算，不同伺服器、集群、超算等執行任務方式不一樣，這裡只適用最簡單的情況 echo $i #返回當前任務編號 echo $i #返回當前任務編號 echo $i #返回當前任務編號 ，這裡多次返回任務編號主要是為了在顯示界面可以很清楚的看到當前計算的位置，以估算總體任務的時間與其他安排 cd ../../ #返回初始目錄 ，如果自身氣場特殊極易吸引bug，建議更改為直接返回到固定的指定目錄，在此固定目錄下進行下一步操作。 done #待for循環執行完72次後循環結束

在這個腳本中，主要是兩個地方的修改，一個是第三行的「｛｝」內的數字區間，這個要根據thirdorder實際生成的任務數量進行更改；一個是第九行的vasp執行命令，這裡是直接進行的24核並行運算，而在實際操作過程中要根據所在機器的情況改動，視情況而定。

在計算腳本執行完畢後，執行thirdorder的命令去獲取三階力學常數

#命令模版 find job* -name vasprun.xml|sort -n|/ thirdorder_vasp.py reap x y z -c find to-run/{01..72}/ -name vasprun.xml|sort -n| thirdorder_vasp.py reap 2 2 2 -12#其中job* 是指所有單一任務計算的文件夾，這裡我們也可以把第二行的相應部分改為 to-run/*/

執行完後會在當前文件夾得到文件FORCE_CONSTANTS_3RD

然後我們將前面計算聲子譜所得到的FORCE_CONSTANTS改名為FORCE_CONSTANTS_2ND

除此之外還需要準備的是一個CONTROL文件，這裡展示一個模版（來源：天璣算），具體如何編寫應閱讀ShengBTE的說明文件（在ShengBTE的安裝目錄裡，或在其官網都可以找到，在我們介紹ShengBTE的推送裡也有附錄上），並視具體情況改動。

&allocations nelements=2 #元素種類 natoms=16 #原子數量 ngrid(:)=21 21 1 #計算網格取點，需要進行嚴格測試&end&crystal lfactor=0.100000 #晶格縮放比例，ShenBTE採用nm作為單位，因此與vasp差一個數量級 lattvec(:,1)=10.1939042571769480 0.0000000000000000 -0.0005202257779721 lattvec(:,2)=0.0000000000000000 5.8835567629900831 0.0000000000000000 lattvec(:,3)=-0.0016161778037250 0.0000000000000000 24.2394661807105862 elements= "Au" "S" #元素符號 types= 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 #對應不同元素 positions(:,1)=0.0792887906745038 0.7628181349429907 0.5477417153337826 positions(:,2)=0.5792887900793176 0.2628181570448227 0.5477417107462345 positions(:,3)=0.3420845204951478 0.5000001343163234 0.5479707485969261 positions(:,4)=0.8420845388897342 0.0000001380045518 0.5479707590273382 positions(:,5)=0.0792887988980483 0.2371820787407889 0.5477416212454472 positions(:,6)=0.5792887903691906 0.7371820856979016 0.5477416090636752 positions(:,7)= 0.4207111952601045 0.7371818360856525 0.4522582828547815 positions(:,8)=0.9207112108774353 0.2371818574692190 0.4522582909979565 positions(:,9)=0.4207111998038618 0.2628179122832144 0.4522583851953175 positions(:,10)=0.9207111996088315 0.7628179378276924 0.4522583841178008 positions(:,11)=0.1579154828371944 0.9999998672025896 0.4520292560634112 positions(:,12)= 0.6579154646203824 0.4999998598272342 0.4520292358640028 positions(:,13)=0.3333765322103900 0.0000000837505425 0.3901738436731563 positions(:,14)=0.8333765550624702 0.5000000766771890 0.3901738421682481 positions(:,15)=0.1666234737495292 0.4999999208298789 0.6098261606528128 positions(:,16)=0.6666234565638564 0.9999999192994125 0.6098261543990704 scell(:)= 2 4 1 #二階力常數計算擴胞倍數&end&parameters T=300 #溫度控制 scalebroad=0.1 #控制高斯展寬，默認取1，手冊中說改小後對計算精度沒有影響，但是最近的一些研究發現對計算精度影響較大。收斂測試反應，取到4才收斂&end&flags nonanalytic=.TRUE. #控制參數具體參考手冊 nanowires=.FALSE. #控制參數具體參考手冊&end

特別注意，模版中含有中文注釋，在實際操作過程中應當刪除中文字符以防出錯。

由於ShengBTE可能沒有辦法同時設置多個溫度（存疑，目前權且認為如此），我們編寫自動計算腳本的時候要把上面31行到38行都刪除掉，再通過腳本補充進CONTROL文件內。

根據已有模型的參數，編寫的CONTROL文件如下

&allocations nelements=1 natoms=2 ngrid(:)=10 10 10 &end&crystal lfactor=0.100000 lattvec(:,1)= 0.0000000000000000 0.5070246140293163 0.5070246140293163 lattvec(:,2)= 0.5070246140293163 0.0000000000000000 0.5070246140293163 lattvec(:,3)= 0.5070246140293163 0.5070246140293163 0.0000000000000000 elements= "Si" types= 1 1 positions(:,1)= 0.8750000000000000 0.8750000000000000 0.8750000000000000 positions(:,2)= 0.1250000000000000 0.1250000000000000 0.1250000000000000 scell(:)= 4 4 4 &end

我們新建一個文件夾存放好ShengBTE計算所需的CONTROL， FORCE_CONSTANTS_2ND， FORCE_CONSTANTS_3RD文件，並編寫運行腳本

具體腳本如下：

for i in `seq 300 10 900` # 對i賦值，區間為300到900，每擱10取一個數，即（300 310 320 330 ....890 900），可以根據自身情況自行設置do #開始循環echo $i #返回任務號echo $i #返回任務號echo $i #返回任務號mkdir $i #建立單一任務文件夾cp CONTROL FORCE_CONSTANTS_2ND FORCE_CONSTANTS_3RD $i #將二階三階力學常數文件和CONTROL文件複製到單一任務文件夾#注意，這裡的CONTROL文件夾是刪減過的cd $i #進入單一文件夾內echo \&parameters >> CONTROL #補充CONTROL文件 echo T= $i >> CONTROL #補充CONTROL文件 #溫度設置 如果同時設置多個溫度，同一行只會算第一個並報錯。如果多次設定一個溫度到同一個CONTROL文件，只會運算最後一個。 echo scalebroad=0.1 >> CONTROL #補充CONTROL文件 echo \&end >> CONTROL #補充CONTROL文件echo \&flags >> CONTROL #補充CONTROL文件echo nonanalytic=.TRUE. >> CONTROL #補充CONTROL文件echo nanowires=.FALSE. >> CONTROL #補充CONTROL文件 echo \&end >> CONTROL #補充CONTROL文件 #這裡補充CONTROL文件用的是笨方法，如果有相應改進方法可以自行更改，簡潔代碼mpirun -np 24 ShengBTE | tee BTElog #執行ShengBTE計算，這裡的執行命令應視情況而定。cd .. #返回初始目錄done #結束任務

腳本運行結束後便可得到所設定溫度點內的熱運輸性質，並可進一步處理數據。