一种基于晶格重定向的晶体解离和滑移能垒自动计算方法

1.本发明属于材料计算领域，具体是一种基于晶格重定向的晶体解离和滑移能垒自动计算方法。


背景技术：

2.解离应变下的解离能垒和临界应力决定了：劈开给定晶面之间的化学键所需的能垒和临界应力；滑移应变下的滑移能垒和临界应力代表了：给定晶面沿特定滑移方向的能垒和临界应力。而且，基于滑移应变计算的层错能曲线决定了位错的错配能，是位错形核和滑移的主要决定参数。
3.然而截止到目前为止，现有技术仅采用解离和滑移能垒以及相应临界应力研究特定的材料和结构[physics reports 826,1
‑
49(2019)]，缺乏高通量计算的研究工作。这不仅严重限制了解离和滑移能垒以及相应临界应力在研究材料位错性质和强韧化机理中的应用，还难以满足“材料基因组计划”中对于解离和滑移能垒以及相应临界应力大数据的需求。


技术实现要素：

[0004]
针对上述问题，从实际应用需求出发，本发明提出了一种基于晶格重定向的晶体解离和滑移能垒自动计算方法，不仅为求解晶体材料解离和滑移能垒以及相应临界应力提供了高效、低成本的途径，还可以提高高强度/高韧性材料筛选和设计的效率，缩短新材料研发应用周期。
[0005]
所述的基于晶格重定向的晶体解离和滑移能垒以及相应临界应力自动计算方法，具体步骤如下：
[0006]
步骤一、针对某晶体，输入该晶体的结构文件，读取晶体学结构信息；
[0007]
结构信息包括标记名字、晶格基矢向量、元素类型、元素类型的个数、总原子数、不同元素的原子个数、是否选择selective dynamics、原子坐标的类型、原子坐标和原子位置约束等。
[0008]
步骤二、根据读取的晶体学结构信息确定该晶体结构的解离面或滑移面滑移方向；
[0009]
确定方法为：依据晶体拉伸应变和剪切应变的最小本征强度，获取晶体结构在沿最弱拉伸和剪切方向应变下的失稳结构和化学键断位置，从而判定该晶体结构中的最弱键合位置，选择最弱的键合位置所在的晶面晶向作为解离面或滑移面滑移方向；或自动解析该晶体结构中最大面间距的晶面，作为解离面或滑移面；
[0010]
步骤三、根据解离面或滑移面滑移方向，自动进行晶格基矢的重定向。
[0011]
对于解离应变，晶格基矢的重定向为：
[0012]
首先，自动设定两个位于解离面(h
c
k
c
l
c
)的不共线晶向[u
c1
v
c1
w
c1
]和[u
c2
v
c2
w
c2
]，并满足条件：
[0013]
u
c1
h
c
+v
c1
k
c
+w
c1
l
c
＝0
[0014]
u
c2
h
c
+v
c2
k
c
+w
c2
l
c
＝0
[0015]
然后，将晶格基矢重定向为[u
c1
v
c1
w
c1
]、[u
c2
v
c2
w
c2
]和[h
c
k
c
l
c
]，并将晶向[u
c1
v
c1
w
c1
]置于x轴方向，将晶向[u
c2
v
c2
w
c2
]置于x
‑
y面上，则解离面(h
c
k
c
l
c
)平行于x
‑
y面。
[0016]
对于滑移应变，晶格基矢的重定向为：
[0017]
首先，自动设定一个位于滑移面(h
s
k
s
l
s
)，且与滑移方向[u
s
v
s
w
s
]不共线的晶向[u
s1
v
s1
w
s1
]，满足条件：
[0018]
u
s
h
s
+v
s
k
s
+w
s
l
s
＝0
[0019]
u
s1
h
s
+v
s1
k
s
+w
s1
l
s
＝0；
[0020]
然后，将晶格基矢自动重定向为[u
s
v
s
w
s
]、[u
s1
v
s1
w
s1
]和[h
s
k
s
l
s
]，并将晶向[u
s
v
s
w
s
]置于x轴方向，将晶向[u
s1
v
s1
w
s1
]置于x
‑
y面上，则滑移面(h
s
k
s
l
s
)平行于x
‑
y面。
[0021]
步骤四、对重定向后的晶体结构施加解离应变或滑移应变，生成应变结构文件；
[0022]
应变过程分为两步，具体如下：
[0023]
第一步为：依据平行于x
‑
y面的解离面或滑移面相对于笛卡尔z轴的位置p
z
，将晶体结构中的原子分为两部分，即z轴坐标小于p
z
的原子为“下半部分”，z轴坐标大于p
z
的原子为“上半部分”。
[0024]
对于解离应变，将“上半部分”原子施加解离位移d，即将p
z
改为p
z
+d；
[0025]
对于滑移应变，将“上半部分”原子施加滑移位移u，即将p
x
和p
y
分别改为p
x
+u
x
和p
y
+u
y
，其中u
x
为滑移距离u在笛卡尔x的分量，u
y
为滑移距离u在笛卡尔y轴的分量；从而获得应变后的原子坐标。
[0026]
第二步为：将应变矩阵ε施加在晶格基矢向量a
i
上，i＝1,2,3，得到应变后的晶格基矢向量a
′
i
：
[0027][0028]
其中i为3
×
3单位矩阵，以及对于解离应变，应变矩阵ε为：其中ε
zz
＝d/l
cz
，l
cz
为重定向的晶格基矢c在笛卡尔z轴的分量；
[0029]
对于滑移应变，应变矩阵ε为：其中ε
zx
＝u
x
/l
cz
和ε
zy
＝u
y
/l
cz
。
[0030]
将应变后的原子坐标和晶格基矢向量共同输出作为应变结构文件，生成一系列不同解离应变值或滑移应变值的应变结构集合，分别命名为pos_[d]和pos_[u
x
]_[u
y
]
[0031]
步骤五、通过第一性原理法并行对应变结构集合pos_[d]和pos_[u
x
]_[u
y
]进行结构弛豫和静态计算，得到解离能χ(d)
‑
应变d曲线或层错能γ(u)
‑
应变u曲线；
[0032]
对于非弛豫模式，则始终保持应变结构中的晶格基矢和原子位置不变。
[0033]
对于弛豫模式，在解离应变下，只弛豫晶格基矢和解离面附近的原子坐标；在滑移
应变下，只弛豫晶胞基矢和垂直于滑移面方向的原子坐标。
[0034]
步骤六、对解离能χ(d)
‑
应变d曲线和层错能γ(u)
‑
应变u曲线进行拟合，求得解离和滑移能垒及相应临界应力；
[0035]
首先，对于解离能χ(d)
‑
应变d曲线，采用如下公式拟合：
[0036][0037]
其中m
i
和n为未知参数；
[0038]
对于层错能γ(u)
‑
应变u曲线，采用如下公式拟合：
[0039][0040]
且满足γ(0)＝0和dγ/du
u＝0
＝0；其中r0、r
i
和i
i
为未知参数。
[0041]
然后，对拟合解离能曲线或层错能曲线进行数值求导，获得解离应力曲线和滑移应力曲线
[0042]
最后，利用解离应力曲线和滑移应力曲线计算解离能垒和临界应力，以及滑移能垒和临界应力；
[0043]
解离能垒和临界应力分别为：
[0044]
滑移能垒和临界应力分别为γ
us
＝max{γ(u)}和
[0045]
本发明的优点在于：
[0046]
1)、一种基于晶格重定向的晶体解离和滑移能垒自动计算方法，可以计算任意对称性三维晶体的晶体解离和滑移能垒以及相应临界应力。
[0047]
2)、一种基于晶格重定向的晶体解离和滑移能垒自动计算方法，不仅可以实现材料力学性能的高通量计算，还可以满足“材料基因组计划”中对于晶体解离和滑移能垒以及相应临界应力大数据的需求。
[0048]
3)、一种基于晶格重定向的晶体解离和滑移能垒自动计算方法，能够用于快速筛选具有高强度/高硬度的三维材料，对超硬材料和高性能结构材料设计具有指导意义。
附图说明
[0049]
图1为本发明一种基于晶格重定向的晶体解离和滑移能垒自动计算方法的原理图；
[0050]
图2为本发明一种基于晶格重定向的晶体解离和滑移能垒自动计算方法的流程图；
[0051]
图3为本发明自动进行晶格基矢的重定向的示意图。
[0052]
图4为本发明解离应变和滑移应变的示意图。
具体实施方式
[0053]
下面结合实施例和附图对本发明做进一步详细说明；
[0054]
本发明公开了一种基于晶格重定向的晶体解离和滑移能垒自动计算方法，如图1
所示，包含了晶体学结构信息的读取、晶体对称性分析、晶格自动重定向、系列应变结构的生成、并行弛豫模式的选择以及第一性原理计算收敛、解离能或层错能曲线拟合求导、以及解离和滑移能垒和临界应力自动求解。不仅可以为求解晶体材料解离和滑移能垒和临界应力提供了高效、低成本的途径，还可以提高高强度/高韧性材料筛选和设计的效率，缩短新材料研发应用周期。
[0055]
如图2所示，具体步骤如下：
[0056]
步骤一、针对某晶体，输入该晶体的结构文件，读取晶体学结构信息；
[0057]
结构信息包括标记名字、晶格基矢向量、元素类型、元素类型的个数、总原子数、不同元素的原子个数、是否选择selective dynamics、原子坐标的类型(笛卡尔坐标或分数坐标)、原子坐标和原子位置约束(决定松弛过程中是否允许原子移动)等。
[0058]
如果输入结构原子坐标使用了分数坐标，则自动转换为笛卡尔坐标。若选择selective dynamics，则将之去掉。
[0059]
步骤二、晶体对称性解析，根据读取的晶体学结构信息确定该晶体结构的解离面或滑移面滑移方向；
[0060]
根据读取的晶体学结构信息与设定的对称性解析精度，输入结构的对称性由spglib接口程序确定。确定方法为：依据晶体拉伸应变和剪切应变的最小本征强度，获取晶体结构在沿最弱拉伸和剪切方向应变下的失稳结构和化学键断位置，从而判定该晶体结构中的最弱键合位置，选择最弱的键合位置所在的晶面晶向作为解离面或滑移面滑移方向；或自动解析该晶体结构中最大面间距的晶面，作为解离面或滑移面；
[0061]
步骤三、根据解离面或滑移面滑移方向，自动进行晶格基矢的重定向。
[0062]
如图3所示，对于解离应变，晶格基矢的重定向为：
[0063]
首先，自动设定两个位于解离面(h
c
k
c
l
c
)的不共线晶向[u
c1
v
c1
w
c1
]和[u
c2
v
c2
w
c2
]，并满足条件：
[0064]
u
c1
h
c
+v
c1
k
c
+w
c1
l
c
＝0
[0065]
u
c2
h
c
+v
c2
k
c
+w
c2
l
c
＝0
[0066]
然后，将晶格基矢重定向为[u
c1
v
c1
w
c1
]、[u
c2
v
c2
w
c2
]和[h
c
k
c
l
c
]，并将晶向[u
c1
v
c1
w
c1
]置于x轴方向，将晶向[u
c2
v
c2
w
c2
]置于x
‑
y面上，则解离面(h
c
k
c
l
c
)平行于x
‑
y面。
[0067]
对于滑移应变，晶格基矢的重定向为：
[0068]
首先，自动设定一个位于滑移面(h
s
k
s
l
s
)，且与滑移方向[u
s
v
s
w
s
]不共线的晶向[u
s1
v
s1
w
s1
]，满足条件：
[0069]
u
s
h
s
+v
s
k
s
+w
s
l
s
＝0
[0070]
u
s1
h
s
+v
s1
k
s
+w
s1
l
s
＝0；
[0071]
然后，将晶格基矢自动重定向为[u
s
v
s
w
s
]、[u
s1
v
s1
w
s1
]和[h
s
k
s
l
s
]，并将晶向[u
s
v
s
w
s
]置于x轴方向，将晶向[u
s1
v
s1
w
s1
]置于x
‑
y面上，则滑移面(h
s
k
s
l
s
)平行于x
‑
y面。
[0072]
为了尽可能减小重定向超胞的大小，提高计算效率，在晶体对称性解析的基础上，对于特定结构和晶向，允许晶向指数为分数(如1/2等)。
[0073]
步骤四、对重定向后的晶体结构施加系列解离应变或滑移应变，生成应变结构文件；
[0074]
应变过程分为两步，具体如下：
[0075]
第一步为：依据平行于x
‑
y面的解离面或滑移面相对于笛卡尔z轴的位置p
z
，将晶体结构中的原子分为两部分，即z轴坐标小于p
z
的原子为“下半部分”，z轴坐标大于p
z
的原子为“上半部分”。
[0076]
施加解离应变则将两个相邻晶面沿解离方向分开一定距离d，而施加滑移应变则将两个相邻晶面(滑移面)沿滑移方向滑移一定距离u，并同时保持其他原子层的相对位置不变。
[0077]
具体为：对于解离应变，将“上半部分”原子施加解离位移d，即将p
z
改为p
z
+d；
[0078]
对于滑移应变，将“上半部分”原子施加滑移位移u，即将p
x
和p
y
分别改为p
x
+u
x
和p
y
+u
y
，其中u
x
为滑移距离u在笛卡尔x的分量；u
y
为滑移距离u在笛卡尔y轴的分量，从而获得应变后的原子坐标。
[0079]
第二步为：通过矩阵运算对晶胞结构施加应变，即将应变矩阵ε施加在晶格基矢向量a
i
上，i＝1,2,3，得到应变后的晶格基矢向量a
′
i
：
[0080][0081]
其中i为3
×
3单位矩阵，以及对于解离应变，应变矩阵ε为：其中ε
zz
＝d/l
cz
，l
cz
为重定向的晶格基矢c在笛卡尔z轴的分量；
[0082]
对于滑移应变，应变矩阵ε为：其中ε
zx
＝u
x
/l
cz
和ε
zy
＝u
y
/l
cz
。
[0083]
将应变后的原子坐标和晶格基矢向量共同输出作为应变结构文件，生成一系列不同解离应变值或滑移应变值的应变结构集合，分别命名为pos_[d]和pos_[u
x
]_[u
y
]，如图4所示。
[0084]
步骤五、通过第一性原理法并行对应变结构集合pos_[d]和pos_[u
x
]_[u
y
]进行结构弛豫和静态计算，得到解离能χ(d)
‑
应变d曲线或层错能γ(u)
‑
应变u曲线；
[0085]
对于晶体解离和滑移能垒和临界应力的计算，可选择非弛豫模式和弛豫模式。对于非弛豫模式，则始终保持应变结构中的晶格基矢和原子位置不变。
[0086]
对于弛豫模式，在解离应变下，只弛豫晶格基矢和解离面附近的原子坐标；在滑移应变下，只弛豫晶胞基矢和垂直于滑移面方向的原子坐标。
[0087]
通过第一性原理(从头算)计算软件并行对应变结构pos_[d]和pos_[u
x
]_[u
y
]进行结构弛豫(只在弛豫模式下)和静态计算，并获取应变结构的能量和应力值，并判断结构是否收敛。如果是，则提取弛豫好结构的能量，否则，将该结构复制为第一性原理结构文件，重新进行结构弛豫。
[0088]
步骤六、对解离能χ(d)
‑
应变d曲线和层错能γ(u)
‑
应变u曲线进行拟合，求得解离和滑移能垒及相应临界应力；
[0089]
首先，基于不同应变值结构的第一性原理计算结果，自动获取解离能χ(d)
‑
应变曲
线或层错能γ(u)
‑
应变曲线。
[0090]
对于解离能χ(d)
‑
应变d曲线，采用如下公式拟合：
[0091][0092]
其中m
i
和n为未知参数；n的大小决定了拟合精度。
[0093]
对于层错能γ(u)
‑
应变u曲线，采用如下公式拟合：
[0094][0095]
且满足γ(0)＝0和dγ/du|
u＝0
＝0；其中r0、r
i
和i
i
为未知参数。
[0096]
然后，对拟合解离能曲线或层错能曲线进行数值求导，获得解离应力曲线和滑移应力曲线
[0097]
最后，利用解离应力曲线和滑移应力曲线计算解离能垒和临界应力，以及滑移能垒和临界应力；
[0098]
解离能垒和临界应力分别为：χ
c
＝lim
d
→
∞
χ(d)和
[0099]
滑移能垒和临界应力分别为γ
us
＝max{γ(u)}和
[0100]
实施例：
[0101]
(1)预处理：准备晶体结构文件以及结构完全弛豫等。
[0102]
(2)读取输入结构文件的晶体学结构信息，包括标记名字、晶格基矢向量、元素类型、元素种类的个数、总原子数、不同元素的原子个数、是否选择selective dynamics、原子坐标的类型、原子坐标以及原子位置约束等。
[0103]
(3)采用spglib接口程序确定晶胞的空间群。
[0104]
(4)根据输入的解离面(hkl)和滑移路径(hkl)[uvw]，自动进行晶格重定向，并将重定向后的结构复制为alipos0文件，作为解离或滑移应变的初始结构。
[0105]
(5)自动生成一系列解离或滑移应变的应变结构，并分别命名为pos_[d]和pos_[u
x
]_[u
y
]。
[0106]
(6)并行计算系列解离或滑移应变的应变结构的能量，对每个应变值，构建单独的文件夹，并在每个文件夹下，并行进行步骤(7)
‑
(9)。
[0107]
(7)选择弛豫模式。
[0108]
(8)通过第一性原理(从头算)计算软件，对应变晶胞进行结构弛豫并计算应变晶胞的能量和应力值。
[0109]
(9)判断应变结构是否收敛。如果是，则提取弛豫好的结构的能量；否则，将该结构复制为第一性原理结构文件，重新进行结构弛豫。
[0110]
(10)待所有应变值都计算完后，输出解离能
‑
应变曲线或层错能
‑
应变曲线。对解离能曲线或层错能曲线进行拟合并求导，自动获得解离应力曲线或滑移应力曲线。
[0111]
(11)基于解离能和解离应力曲线或层错能和滑移应力曲线，自动获取晶体解离和滑移能垒和临界应力。