bash+vasp+ShengBTE自動計(jì)算材料熱運(yùn)輸性質(zhì)腳本
腳本分享
之前我們已經(jīng)詳細(xì)介紹了ShengBTE軟件與其安裝教程，詳細(xì)請查看ShengBTE簡單介紹、安裝與使用
由于在計(jì)算三階力學(xué)常數(shù)時(shí)，根據(jù)擴(kuò)胞比例以及近鄰原子受力情況的不同，生成的3RD.POSCAR數(shù)量也盡不相同，筆者最近一次的計(jì)算總共計(jì)算了700+個(gè)任務(wù)，如果在此過程中手動布置任務(wù)則費(fèi)時(shí)費(fèi)力，效率低下。而在后續(xù)的ShengBTE計(jì)算時(shí)，通過bash腳本進(jìn)行任務(wù)布置也可以實(shí)現(xiàn)多任務(wù)順序計(jì)算，節(jié)省精力與時(shí)間。
本文針對材料的熱輸運(yùn)性質(zhì)計(jì)算主要參考天璣算所提供的教程，詳細(xì)請見https://phadwiki.com/page237.html?article_id=79&pagenum=all
這里我們使用兩個(gè)原子的硅原胞做演示
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
對其進(jìn)行結(jié)構(gòu)優(yōu)化后擴(kuò)胞并使用vasp結(jié)合phonopy進(jìn)行聲子譜計(jì)算，這里選擇擴(kuò)胞比例為4 4 4，以得到二階力學(xué)常數(shù)備用。
下圖為計(jì)算所得聲子譜與聲子態(tài)密度，不存在虛頻。
有關(guān)phonopy的使用方法請參見官網(wǎng)https://phonopy.github.io/phonopy/examples.html
也可參考天璣算的教程：PHONOPY使用教程（上）PHONOPY使用教程（下）
以后如果有時(shí)間也會整理一些聲子譜或者熱膨脹系數(shù)、格林奈森常數(shù)等計(jì)算的案例。
在計(jì)算三階力學(xué)常數(shù)時(shí)，通過執(zhí)行thirdorder命令生成擴(kuò)胞文件（thirdorder已添加環(huán)境變量），這里我們選擇擴(kuò)胞比例為2 2 2，擴(kuò)胞后原子總數(shù)為8+8=16，近鄰原子數(shù)選為12（這里的計(jì)算選擇僅做參考，實(shí)際計(jì)算時(shí)需要自行調(diào)整參數(shù)等進(jìn)行測試并選取可行方案）
#擴(kuò)胞命令模版 thirdorder_vasp.py sow x y z -c (x,y,z是三個(gè)方向的擴(kuò)胞倍數(shù)，c是考慮多少個(gè)臨近原子的受力來計(jì)算力常數(shù)矩陣)thirdorder_vasp.py sow 2 2 2 -12
根據(jù)顯示結(jié)果這里有72個(gè)DFT任務(wù)，相對來說較少，但是如果手動一個(gè)任務(wù)一個(gè)任務(wù)提交也是非常麻煩且耗時(shí)，這里我們通過一個(gè)腳本來整理文件，并一次性安排所有的計(jì)算任務(wù)。在執(zhí)行腳本之前應(yīng)當(dāng)提前準(zhǔn)備好INCAR、KPOINTS、POTCAR，并存放在當(dāng)前目錄。其中INCAR為自洽運(yùn)算，KPOINTS應(yīng)當(dāng)與計(jì)算聲子譜的任務(wù)相同。
mkdir 3RD && mv 3RD* 3RD #這里將所有thirdorder生成的POSCAR文件統(tǒng)一存放在3RD文件夾內(nèi)mkdir to-run #建立計(jì)算文件夾for i in {01..72} #總計(jì)72個(gè)計(jì)算任務(wù)，這里對i賦值。如果總?cè)蝿?wù)數(shù)是三位數(shù)的，開頭是｛001..do ?#開始for循環(huán)，總共會執(zhí)行72個(gè)循環(huán)（這里的循環(huán)次數(shù)由上一行｛｝內(nèi)確定）mkdir to-run/$i #建立單一任務(wù)計(jì)算的文件夾cp 3RD/3RD.POSCAR.$i to-run/$i/POSCAR #將每個(gè)thirdorder生成的任務(wù)文件復(fù)制到對應(yīng)的計(jì)算文件夾并重命名為POSCARcp INCAR KPOINTS POTCAR to-run/$i/ #準(zhǔn)備其他計(jì)算文件，這里已經(jīng)默認(rèn)準(zhǔn)備好了。cd to-run/$i/ #進(jìn)入單一任務(wù)計(jì)算文件夾mpirun -np 24 vasp |tee runlog ?#執(zhí)行vasp計(jì)算，不同服務(wù)器、集群、超算等執(zhí)行任務(wù)方式不一樣，這里只適用最簡單的情況echo $i #返回當(dāng)前任務(wù)編號echo $i ?#返回當(dāng)前任務(wù)編號?echo $i ?#返回當(dāng)前任務(wù)編號 ，這里多次返回任務(wù)編號主要是為了在顯示界面可以很清楚的看到當(dāng)前計(jì)算的位置，以估算總體任務(wù)的時(shí)間與其他安排cd ../../ #返回初始目錄 ，如果自身氣場特殊極易吸引bug，建議更改為直接返回到固定的指定目錄，在此固定目錄下進(jìn)行下一步操作。done ?#待for循環(huán)執(zhí)行完72次后循環(huán)結(jié)束在這個(gè)腳本中，主要是兩個(gè)地方的修改，一個(gè)是第三行的“｛｝”內(nèi)的數(shù)字區(qū)間，這個(gè)要根據(jù)thirdorder實(shí)際生成的任務(wù)數(shù)量進(jìn)行更改；一個(gè)是第九行的vasp執(zhí)行命令，這里是直接進(jìn)行的24核并行運(yùn)算，而在實(shí)際操作過程中要根據(jù)所在機(jī)器的情況改動，視情況而定。
在計(jì)算腳本執(zhí)行完畢后，執(zhí)行thirdorder的命令去獲取三階力學(xué)常數(shù)
#命令模版 find job* -name vasprun.xml|sort -n|/ thirdorder_vasp.py reap x y z -cfind to-run/{01..72}/ -name vasprun.xml|sort -n| thirdorder_vasp.py reap 2 2 2 -12#其中job* 是指所有單一任務(wù)計(jì)算的文件夾，這里我們也可以把第二行的相應(yīng)部分改為 to-run/*/執(zhí)行完后會在當(dāng)前文件夾得到文件FORCE_CONSTANTS_3RD
然后我們將前面計(jì)算聲子譜所得到的FORCE_CONSTANTS改名為FORCE_CONSTANTS_2ND
除此之外還需要準(zhǔn)備的是一個(gè)CONTROL文件，這里展示一個(gè)模版（來源：天璣算），具體如何編寫應(yīng)閱讀ShengBTE的說明文件（在ShengBTE的安裝目錄里，或在其官網(wǎng)都可以找到，在我們介紹ShengBTE的推送里也有附錄上），并視具體情況改動。
&allocations? ?nelements=2 ? #元素種類? ?natoms=16 ? ? ?#原子數(shù)量? ?ngrid(:)=21 21 1 ?#計(jì)算網(wǎng)格取點(diǎn)，需要進(jìn)行嚴(yán)格測試&end&crystal? ?lfactor=0.100000 ?#晶格縮放比例，ShenBTE采用nm作為單位，因此與vasp差一個(gè)數(shù)量級? ?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 ? #對應(yīng)不同元素?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 ? ? #二階力常數(shù)計(jì)算擴(kuò)胞倍數(shù)&end¶meters? ?T=300 ? ?#溫度控制? ?scalebroad=0.1 ? #控制高斯展寬，默認(rèn)取1，手冊中說改小后對計(jì)算精度沒有影響，但是最近的一些研究發(fā)現(xiàn)對計(jì)算精度影響較大。收斂測試反應(yīng)，取到4才收斂&end&flags? ?nonanalytic=.TRUE. ?#控制參數(shù)具體參考手冊? ?nanowires=.FALSE. ?#控制參數(shù)具體參考手冊&end特別注意，模版中含有中文注釋，在實(shí)際操作過程中應(yīng)當(dāng)刪除中文字符以防出錯(cuò)。
由于ShengBTE可能沒有辦法同時(shí)設(shè)置多個(gè)溫度（存疑，目前權(quán)且認(rèn)為如此），我們編寫自動計(jì)算腳本的時(shí)候要把上面31行到38行都刪除掉，再通過腳本補(bǔ)充進(jìn)CONTROL文件內(nèi)。根據(jù)已有模型的參數(shù)，編寫的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
我們新建一個(gè)文件夾存放好ShengBTE計(jì)算所需的CONTROL， FORCE_CONSTANTS_2ND， FORCE_CONSTANTS_3RD文件，并編寫運(yùn)行腳本
具體腳本如下：
for i in `seq 300 10 900` # 對i賦值，區(qū)間為300到900，每擱10取一個(gè)數(shù)，即（300 310 320 330 ....890 900），可以根據(jù)自身情況自行設(shè)置do ?#開始循環(huán)echo $i ?#返回任務(wù)號echo $i ?#返回任務(wù)號echo $i ?#返回任務(wù)號mkdir $i #建立單一任務(wù)文件夾cp CONTROL FORCE_CONSTANTS_2ND FORCE_CONSTANTS_3RD $i ?#將二階三階力學(xué)常數(shù)文件和CONTROL文件復(fù)制到單一任務(wù)文件夾#注意，這里的CONTROL文件夾是刪減過的cd $i #進(jìn)入單一文件夾內(nèi)echo \¶meters >> CONTROL ? ? ?#補(bǔ)充CONTROL文件?echo ?T= $i ?>> CONTROL ? ?#補(bǔ)充CONTROL文件?#溫度設(shè)置 如果同時(shí)設(shè)置多個(gè)溫度，同一行只會算第一個(gè)并報(bào)錯(cuò)。如果多次設(shè)定一個(gè)溫度到同一個(gè)CONTROL文件，只會運(yùn)算最后一個(gè)。echo ? scalebroad=0.1 ? ?>> CONTROL ? ?#補(bǔ)充CONTROL文件echo \&end ? ?>> CONTROL ? ?#補(bǔ)充CONTROL文件echo \&flags ?>> CONTROL ? ?#補(bǔ)充CONTROL文件echo ? ?nonanalytic=.TRUE. ? >> CONTROL ? ?#補(bǔ)充CONTROL文件echo ? ?nanowires=.FALSE. ? >> CONTROL ? ?#補(bǔ)充CONTROL文件echo \&end ?>> CONTROL ? ?#補(bǔ)充CONTROL文件#這里補(bǔ)充CONTROL文件用的是笨方法，如果有相應(yīng)改進(jìn)方法可以自行更改，簡潔代碼mpirun -np 24 ShengBTE | tee BTElog ?#執(zhí)行ShengBTE計(jì)算，這里的執(zhí)行命令應(yīng)視情況而定。cd .. #返回初始目錄done ?#結(jié)束任務(wù)腳本運(yùn)行結(jié)束后便可得到所設(shè)定溫度點(diǎn)內(nèi)的熱運(yùn)輸性質(zhì)，并可進(jìn)一步處理數(shù)據(jù)。
愿有所成
IEchoQ
引喻失義??妄自菲薄
IEchoQ
