单层BiBr的结构优化踩坑（一）

发表于： 2023-10-13 分类于：计算笔记

字数： 1785 阅读：≈ 4分钟浏览：

写在前面

初始结构的合理性对于结构弛豫来说非常重要，尽管设置ISIF=3，ISIF=4时对晶格参数的要求没有那么敏感，但是这种情况下研究二维材料时，弛豫结果不一定是我们需要的（晶胞可能会变形），这时就只能选择设置ISIF=2（保持晶格形状、体积不变，只改变原子坐标）。一般来说，只要赝势文件里面的原子顺序没搞错，初始结构不是非常非常糟糕的话，弛豫后的结构肉眼看上去不会跟弛豫前有很大区别，因此弛豫后的结构是否合理反而由初始结构决定了。我个人的理解是ISIF=2适合在初始结构附近找一个计算上能够收敛的结果，但这个结果不一定合理。此外，在复现别人结果的时候一定要仔细，赝势、交换关联泛函、K点选取，vasp版本号不同可能对计算结果也有一定影响（有一个计算结果是用5.3.3版本算出来的，我用5.4.4, 6.3, 6.4都算不出，后三个的结果倒是一样。）

问题描述

我遇到的问题是按照实验给的大致数据建立Bi4Br4的晶胞并使用ISIF=2优化后，4个Bi-Br键，其中有一个键长的参数不太合理（下面的键长会比上面长一些。）

这是一个新型的Bi4Br4单层结构，相邻的两条链在b方向上相差大约b/3。这个工作计算了bulk Bi4Br4中的单层结构，相邻的两条链在b方向上正好相差b/2。晶格常数来源于实验数据。

如果我使用ISIF=4优化新型的结构，那么他将变成类似上图中的结构，不是我想研究的对象。可是如果使用ISIF=2，又没有办法优化晶格参数，此时应该怎么办？

解决思路

老师提示我可以使用shell脚本修改晶格参数后再使用ISIF=2批量提交优化任务，最后比较弛豫结果的总能，看看总能最低的晶格参数优化出来的结构是否合理。二维材料的晶格参数只有两个自由度，相当于建立初始结构后对这个结构进行缩放，剩下的就交给vasp自己搞定。因此，为了直观，方便的改变晶格参数，我建立了一个正交的conventional cell。我们来看一下POSCAR文件的前几行：

Generated by cif2cell 2.0.0.  :  Failed to get author information, No journal information. Species order: Bi Br
  39.192000
  1.000000000000000   0.000000000000000   0.000000000000000
  0.000000000000000   0.110685854255971   0.000000000000000
  0.000000000000000   0.000000000000000   0.637885282710757
  Bi  Br
  24  24

可以看到，正交晶系的好处就是修改晶格参数非常简单，只需要改每一列的一个数就可以了（每列代表一个晶格矢量，注意不是行）。因为这个晶胞本身在a方向上不具备周期性，所以周期性边界条件可能会导致结构边界处的对称性变差，为了保证不影响弛豫结果，我首先建立了一个包含6条链的超胞，然后使用materials studio调整每条链的位置，使相邻两条链在b方向上相差b/3。接着写好基本的输入文件，就可以让我们的脚本登场了。

#!/bin/bash
file=../relax-ref/POSCAR #以此为基础进行修改
cell_len=$(awk 'NR==2 {print $1}' "$file") #获取POSCAR的内容
cell_a=$(awk 'NR==3 {print $1}' "$file")
cell_b=$(awk 'NR==4 {print $2}' "$file")
a_min=7.05 #自定义晶格参数
b_min=4.22
grid_step=0.005 #步长
grid_stepb=0.001
for i in {0..10}
do
        for j in {0..10}
        do  
                a=$(echo "scale=15;($a_min+$i*$grid_step)*6/$cell_len"|bc|awk '{printf "%f", $0}') #换算成POSCAR的内容
                b=$(echo "scale=15;($b_min+$j*$grid_stepb)/$cell_len"|bc|awk '{printf "%f", $0}')
                cp -r ../relax-ref/ relax-${i}-${j}
                newfile=relax-${i}-${j}/POSCAR
                awk -v a="$a" -v b="$b" '
                    NR==3 {$1="  "a"  "} 
                    NR==4 {$2="  "b"  "} 
                    {print}
                    ' "$file" > "$newfile"
                cd relax-${i}-${j}
                sbatch vasp-cqum.sh #提交任务
                cd ..
        done
done

其实如果提取，计算全部使用awk来完成，可能会好一点，shell里面的bc还是不太好用，因为不熟悉awk，再加上很多都是chatgpt帮我写的，所以实现方式非常丑陋，有时间的话可能会考虑改进一下。写好了之后，老师批评我说这种地毯式搜索实在太笨了，非常浪费计算资源，但是我实在是想偷懒（）。

等任务都计算完毕后，再执行另外一个脚本将OSZICAR中的总能提取出来。

#!/bin/bash
for dir in */
do
        out_file=${dir}slurm-* #如果没有标准输出文件，则需要去OUTCAR确认计算是否收敛，去OSZICAR提取总能结果
        in_file=${dir}POSCAR
        i=$(echo $dir | awk -F- '{print $2}')
        j=$(echo $dir | awk -F- '{print substr($3,1,length($3)-1);}')
        a=$(awk 'NR==2 {val_l=$1} NR==3 {val_a=$1; print val_l*val_a/6}' "$in_file") #换算成单链的参数
        b=$(awk 'NR==2 {val_l=$1} NR==4 {val_b=$2; print val_l*val_b}' "$in_file")
        if grep -q "reached required accuracy" $out_file; then  #检测计算结果是否收敛
                F_value=$(grep "F=" $out_file | tail -1 | awk '{print $3}')
        else
                F_value=NaN
        fi  
        echo "$i  $j  $a  $b  $F_value" >> tot_energy.dat
done

最后写一个julia脚本绘制图像。

using DelimitedFiles #读取文件所需的库
using CairoMakie #绘图所需的库
data=readdlm("tot_energy.dat"); #1-5列分别为i j a b F
idx=broadcast(Int,data[:,1:2]).+1; #julia的数组从1开始索引
idx=[tuple(row...) for row in eachrow(idx)] #chatgpt写的，idx是n行2列的矩阵，row就是矩阵的每一行，...表示某种展开
map=zeros(17,9);
for (i, position) in enumerate(idx) #i是循环变量，position是i对应的idx值
    map[position...] = data[i, 5]
end
fig, ax, hm = heatmap(map) #热力图
ax.aspect=2 #坐标轴比例(x/y)
ax.xticks=(1:2:17,broadcast(string,sort(unique(data[:,3])))[1:2:17]);
ax.yticks=(1:9,broadcast(string,sort(unique(data[:,4]))));
Colorbar(fig[1, 2], hm) 
save("tot_energy_2b_wide.png",fig)

为了验证这个方法的可行性，我首先研究了两条链在b方向上正好相差b/2的结构，结果如下

横轴是换算成单链的晶格参数，纵轴是b方向的晶格参数，颜色越深代表总能越低，这与ISIF=4得到的晶格参数相差不大，因此可以尝试用这种方法研究新的结构。