脚本分享
之前我们已经详细介绍了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¶meters 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 \¶meters >> 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 #结束任务
脚本运行结束后便可得到所设定温度点内的热运输性质，并可进一步处理数据。
