晶格参数确定-BM状态方程方法

Birch-Murnaghan状态方程

Birch-Murnaghan状态方程描述了固体体积与其所受压力之间的关系。三阶的等温状态方程如下：

P(V)=\frac{3B_0}{2}\Big[{(\frac{V_0}{V})}^{7/3}-{(\frac{V_0}{V})}^{5/3}\Big]\Big\{1+\frac{3}{4}(B_0^{\prime}\Big[(\frac{V_0}{V}^{2/3}-1)\Big]\Big\}

B_0=-V\Big(\frac{\partial P}{\partial V}\Big)_{P=0}

B_0^{\prime}=\Big(\frac{\partial B}{\partial P}\Big)_{P=0}

其中， $P$ 为压强， $V_0$ 为参考体积， $V$ 为形变后体积， $B_0$ 为体积模量， $B_0^{\prime}$ 为体积模量相对于压强的导数。而系统内能可以通过积分获得：

E=E_0+\frac{9V_0B_0}{16}\Big\{\Big[(\frac{V_0}{V})^{2/3}-1\Big]^3B_0^{\prime}+\Big[(\frac{V_0}{V})^{2/3}-1\Big]^2\Big[6-4(\frac{V_0}{V})^{2/3}\Big]\Big\}

将原式中的体积转化为晶格参数的函数：

E(a)=E_0+\frac{9V_0B_0}{16}\Big\{\Big[(\frac{a_0}{a})^{2}-1\Big]^3B_0^{\prime}+\Big[(\frac{a_0}{a})^{2}-1\Big]^2\Big[6-4(\frac{a_0}{a})^{2}\Big]\Big\}

其中， $E(a)$ 和 $E_0$ 为晶格参数为 $a$ 和 $a_0$ 时的能量。令 $(1/a)^2=x$ 则上述方程能够写成：

y(x)=c_0+c_1x+c_2x^2+c_3x^3

为了得到最佳的晶格参数，需要拟合以上形式的曲线去寻找y(x)最小时x的取值。

拟合数据的获取

通常使用更改POSCAR中的缩放系数的方式来实现批量拟合数据的获取。需要注意的是，POSCAR文件中的坐标形式应当更改为分数坐标（Direct）的形式。

批量提交任务后提取数据的脚本（catBMData.sh)：

for dir in *; do
	if [ -d "$dir" ]; then
		echo -e $dir "\t" $(grep ' without' $dir/OUTCAR | tail -n 1| awk '{print$7}')
	fi
done; > data

脚本解析：

for dir in *; do
遍历当前目录下的所有文件和子目录，* 表示所有项。
if [ -d "$dir" ]; then
检查 dir 是否为一个目录，-d 是用来测试是否是目录。
echo -e $dir "\t" $(grep ' without' $dir/OUTCAR | tail -n 1| awk '{print$7}')
对每个目录：

※ 输出目录名 $dir。

※ 使用 grep ' without' $dir/OUTCAR 查找 OUTCAR 文件中包含 ' without' 的行。

※ 使用 tail -n 1 获取匹配行的最后一行（假设每个目录中可能有多行匹配的结果）。

※ 使用 awk '{print$7}' 提取该行的第7个字段。

※ echo -e 将输出目录名和第7个字段（用制表符 \t 分隔）。
done;
结束 for 循环。
> data
将所有输出重定向到 data 文件，覆盖之前的内容。

数据的拟合以及最佳晶格参数的获取

使用最小二乘法对三次函数的参数进行拟合，并按照最小值点公式获得最佳晶格参数。

#!/bin/env python

import math
import numpy as np

lattice_parameter = 2.8664

a, E = np.loadtxt('data', usecols=(0,1), delimiter='\t', unpack=True)
x = (a * lattice_parameter) ** (-2)
p = np.polyfit(x, E, 3)

c0 = p[3]
c1 = p[2]
c2 = p[1]
c3 = p[0]

x1 = (math.sqrt(4 * c2 ** 2 - 12 * c1 * c3) - 2 * c2) / (6 * c3)
para = 1/math.sqrt(x1)

print 'The optimized lattice parameter is: ' %(para)

需要注意的是，读取的晶格参数并不直接对应拟合公式中的x，而是需要通过 $x=(1/a)^2$ 的变换获得。