体积弛豫 | Eastsheng's Wiki

[toc]

引言

在对晶胞进行体积弛豫时，必须格外小心。如果使用的平面波截断能（ENCUT）或 k 点网格不足，会导致能量-体积曲线出现不连续
这是因为使用了未收敛的基组，产生了非物理的力，即 Pulay 应力。这种应力会使晶胞结构偏离平衡位置，通常导致体积缩小。在体积不变的弛豫中，Pulay 应力可以忽略，但在体积弛豫中则非常重要。

如何消除 Pulay 应力

有多种方法可以修正 Pulay 应力。通常，通过增大 ENCUT 来使用更大的基组，直至计算结果收敛：
1. 将 ENCUT 设为默认截断能的 1.3 倍，或在 VASP.4.4 中使用 PREC = High（注意：这是推荐的最低标准）。
2. 使用默认截断能重新运行 VASP，获得最终的弛豫位置和晶胞参数。
3. 进一步增加 ENCUT 并重复步骤 1 和 2，直到结构不再变化（即收敛）。
注意：这种弛豫会沿所有晶格矢量进行。如果只想弛豫其中一个而保持其他固定（例如表面弛豫中的真空层），应使用 LATTICE_CONSTRAINTS 标签。它可以约束正交晶胞的任意晶格矢量，或非正交晶胞晶格矩阵中的特定行列。

方法一：直接体积弛豫

该方法通过对一个结构进行一系列精确的弛豫计算来完成。

步骤 1：初始弛豫

从初始结构开始，使用高斯或 Methfessel-Paxton 展宽（ISMEAR = 0 或 1）进行弛豫（ISIF = 7；IBRION = 1, 2, 或 3）。

注意：虽然也可使用 ISIF = 3, 6, 8，但不建议同时弛豫多个自由度。示例 INCAR 设置：

1 2	ISMEAR = 0 ; SIGMA = 0.05 # 对于金属，改为 ISMEAR = 1 ; SIGMA = 0.2 ISIF = 7 ; IBRION = 2 ; NSW = 20

步骤 2：再次弛豫

将步骤 1 生成的 CONTCAR 复制为 POSCAR，再次进行弛豫，设置：

1
2
3

ISMEAR = 0 ; SIGMA = 0.05   # 对于金属，改为 ISMEAR = 1 ; SIGMA = 0.2
ISIF = 7 ; IBRION = 2 ; NSW = 20
ISTART = 1   # 保持截断能不变，但根据新晶胞体积更新平面波基组

关键：必须设置 ISTART = 1，这样基组会更新而截断能保持恒定。如果使用 ISTART = 2，基组保持不变，从而保留了人为的 Pulay 应力——这与步骤 1 完全相同的糟糕设置。

步骤 3：单点计算获得精确能量

将展宽方法改为四面体法（ISMEAR = -5），并进行单点计算（不弛豫）。这会给出非常精确的能量。
1
2
ISMEAR = -5 # 四面体法
NSW = 0

（可选）步骤 4：如果仍有问题

如果能量-体积曲线仍然呈现锯齿状，可尝试以下额外步骤：
- 步骤 4a：进一步提高 ENCUT，或通过设置 PREC = Accurate 改进 FFT 网格。
- 步骤 4b：使用 PSTRESS 修正 Pulay 应力（避免增加截断能带来的计算成本）。Pulay 应力对体积和离子构型依赖较弱，主要由成分和基组决定。可以在输出中找到一个很好的估计值，例如：
  1
  external pressure = -100.29567 kB
- 用高截断能与默认截断能计算得到的压力差来修正。方法如下：
  1. 分别用默认截断能和高截断能做两个单点计算。
  2. 从 OUTCAR 中提取两次计算的外部压力值，例如：
  1
  2
  external pressure = -1311.32 kB (默认截断)
  external pressure = -95.66 kB (高截断)
  1. 计算压力差，作为 Pulay 应力的近似值：-1215.66 kB。
  2. 在 INCAR 中设置 PSTRESS = -1215.66。
    这样可以在默认截断能下获得接近高截断能的结构。重申：仅当无法使用高截断能时才用 PSTRESS，因为它只改善结构，不能提高能量精度。

方法二：状态方程拟合

这是避免直接体积弛豫的另一种方法。

对一系列固定体积下的结构进行离子位置和晶胞形状的弛豫，然后将这些结果拟合到状态方程（如 Murnaghan 方程）。

由于 Pulay 应力几乎是各向同性的，它只会给应力张量的对角元添加一个常数值，因此固定体积下的弛豫可以给出非常准确的结构。

此方法还有一个优点：同时得到体模量，并且可以使用默认截断能安全地进行。

步骤 0：构建初始结构

使用晶胞（unit cell），如果已知，采用实验晶格参数a

步骤 1：生成不同体积的结构

从初始结构出发，创建晶格参数分别为 0.90a, 0.95a, 1.0a, 1.05a, 1.10a 的五个结构（以金刚石为例，实验晶格参数 3.567 Å，则 X 值分别为 3.210, 3.389, 3.567, 3.745, 3.924 Å）。示例 POSCAR（金刚石）：
1
2
3
4
5
6
7
8
9
cubic diamond
X
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
2
Direct
-0.125 -0.125 -0.125
0.125 0.125 0.125

步骤 2：计算每个结构的能量

对每个体积的结构进行单点能计算（无需弛豫）。使用以下基本 INCAR 设置：
1
2
ISMEAR = 0 ; SIGMA = 0.05 # 金属用 ISMEAR = 1 ; SIGMA = 0.2
IBRION = 0
如需极高的能量精度，可使用四面体法（ISMEAR = -5，不适用于金属）。

步骤 3：提取晶格参数和能量

写一个简单的 bash 脚本从每个计算的 OUTCAR 中提取数据：

#!/bin/bash
rm -f loop_lattice_constant.dat
for a in */; do
  latt=$(grep "ALAT" "$a"/OUT* | awk '{print $3}')
  E=$(grep "free  e" "$a"/OUTCAR | awk '{print $5}')
  echo "$latt $E" >> loop_lattice_constant.dat
done

步骤 4：拟合状态方程

将晶格参数与能量绘图，并用 Murnaghan 方程（或其他如 Birch-Murnaghan 方程）进行拟合。以下是一个 Python 脚本示例（点击展开）：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def murnaghan(V, E0, V0, B0, Bp):
    """Murnaghan equation of state"""
    return E0 + (B0 * V / Bp) * (((V0 / V)**Bp) / (Bp - 1) + 1) - (B0 * V0) / (Bp - 1)

# 读取数据：volume, energy (假设已转换为体积)
data = np.loadtxt('loop_lattice_constant.dat')
volumes = data[:,0]**3  # 如果数据是晶格参数，需转换为体积
energies = data[:,1]

# 拟合
popt, _ = curve_fit(murnaghan, volumes, energies, p0=[-10, 20, 100, 4])
E0, V0, B0, Bp = popt

# 绘图
V_fit = np.linspace(min(volumes), max(volumes), 100)
E_fit = murnaghan(V_fit, *popt)
plt.plot(volumes, energies, 'o', label='DFT')
plt.plot(V_fit, E_fit, '-', label='Murnaghan fit')
plt.xlabel('Volume (Å³)')
plt.ylabel('Energy (eV)')
plt.legend()
plt.show()

print(f"Equilibrium volume: {V0:.3f} Å³")
print(f"Bulk modulus: {B0:.3f} GPa")

在脚本中，我们显式定义了 Murnaghan 拟合。其他工具（如 pymatgen）也内置了状态方程。
说明：优化后的晶格参数对应 Murnaghan 拟合的最小值。示例中用了 5 个点，建议使用更多点以改进拟合。更多细节可参考硅的教程（使用一系列晶格参数）。

参考

[1] https://vasp.at/wiki/Volume_relaxation

← 声子色散关系计算（VASP） PERTURBO 教程 →