空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

ABACUS计算模拟实例 | III. 材料平衡晶格常数计算

notebook

ABACUS

计算材料学

ABACUS使用教程

notebookABACUS计算材料学ABACUS使用教程

FermiNoodles

weiqingzhou@whu.edu.cn

更新于 2024-06-21

推荐镜像 :abacus:3.6.3-user-guide

推荐机型 :c32_m64_cpu

数据集

ABACUS-cases(v10)

ABACUS计算模拟实例 | III. 材料平衡晶格常数计算

代码

文本

©️ Copyright 2023 @ Authors
作者：Jinghang Wang (AISI电子结构团队实习生) 📨
日期：2024-6-21
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的开始连接按钮，选择 ABACUS:3.6.3-user-guide 镜像及 c32_m64_cpu 节点配置，并在连接后点击右上方的切换内核，选择Bash内核稍等片刻即可运行。本教程运行时间较长，因此不要求完成具体计算，若要完成可按该配置进行。

代码

文本

🎯 本教程旨在帮助初学者学习使用ABACUS自动计算材料的晶格常数，然后使用Brich-Murnaghan方程拟合，最终得到晶体材料的平衡晶格常数。

代码

文本

本节教程目录如下：
1. 自动计算材料晶格常数
2. Birch-Murnaghan方程拟合

代码

文本

在使用ABACUS进行计算时，所有操作都需要在指定的文件夹目录下进行。我们已经准备好了数据集ABACUS-cases，其中文件夹“3”对应本节教程。

代码

文本

首先我们将数据集中的文件复制到根目录下:

代码

文本

[2]

cp -r /bohr/ABACUS-cases-nu6t/v10/ABACUS-cases/3 .

代码

文本

[3]

cd 3 && ls

Birch_fitting.jpg  INPUT  Ni_fcc_energy.csv  STRU                 result
Birch_fitting.py   KPT    OUT.Ni             lattice_constant.sh  time.json

代码

文本

1.自动计算材料晶格常数

代码

文本

在文件夹中，Birch_fitting.py, INPUT, KPT, lattice_constant.sh是本次教程的输入文件，其余均为输出文件。

代码

文本

INPUT文件中的关键参数：
nspin 2: Ni为磁性体系，因此要打开自旋极化；
smearing_method mp: 一阶Methfessel-Paxton方法，对应VASP中的ISMEAR=1，适用于金属体系。

代码

文本

[4]

cat INPUT

INPUT_PARAMETERS RUNNING ABACUS-DFT

#Parameters (1.General)
suffix                  Ni  # suffix of OUTPUT DIR
nspin                   2   #  1/2/4 4 for SOC
symmetry                1   #  0/1  1 for open, default
# symmetry_autoclose      1   # if symmetry error: set symmetry to 0
# symmetry_prec           1e-5  # default
esolver_type            ksdft  # ksdft, ofdft, sdft, tddft, lj, dp
# dft_functional          pbe  # same as upf file, can be lda/pbe/scan/hf/pbe0/hse
ks_solver             genelpa  # default for ksdft-lcao
pseudo_dir              /abacus-develop/tests/PP_ORB
orbital_dir             /abacus-develop/tests/PP_ORB
#init_wfc                file
#init_chg                file


#Parameters (2.Iteration)
calculation             scf # scf relax cell-relax md
ecutwfc                 50 # Ry, Default: 50
scf_thr                 1e-7
scf_nmax                100
#relax_nmax              200
#relax_method            bfgs  # cg, bfgs, cg_bfgs, sd, "fire"
#force_thr_ev            0.05  # ev
# stress_thr            5

#Parameters (3.Basis)
basis_type              lcao  # lcao or pw
# kspacing              0.25 # replace KPT
# gamma_only            1  # 0/1, replace KPT

#Parameters (4.Smearing)
smearing_method         mp    # mp/gaussian/fd/fixed, mp for metal gau for semicon
smearing_sigma          0.008  # Rydberg, 0.008 for mp 0.001 for gau

#Parameters (5.Mixing)
mixing_type             broyden  # pulay/broyden, Default: broyden
#mixing_beta             0.2  # for metal: 0.05-0.4, 0 - 0.1 for difficult
#mixing_gg0              1.5  # only for transition metal

#Parameters (6.Calculation)
cal_force          1
cal_stress         1
out_stru           1  # print STRU in OUT
out_chg            1  # print CHG or not
out_wfc_lcao       1  # print wfc or not, 0: not 1: text format 2: binary format
out_bandgap        0
out_mul            0  # print Mulliken charge and mag of atom in mulliken.txt
# out_interval       1
# restart_save     auto  # false, auto, other
# restart_load     false

代码

文本

[5]

cat KPT

K_POINTS
0
Gamma
7 7 7 0 0 0

代码

文本

这里为读者提供自动计算晶格常数的lattice_constant.sh脚本：
运行该脚本即可获得不同晶格常数下的能量，晶格常数从3埃变化到4埃，用于下一步拟合。
注意：不同于VASP，abacus中的初始磁矩在STRU文件中指定，并且一共有三种指定方式，详见：ABACUS 使用教程｜磁性材料计算，这里采用按元素设定，对所有Ni原子设置初始磁矩为2.

代码

文本

[6]

cat lattice_constant.sh

astart=3.00 #Note: IN bash script, there should be no whitespace in variable assignment. e.g. astart = 30 will fail
astep=0.04
n=25

echo "lat,ene">Ni_fcc_energy.csv
for ((b=0; b<=$n; b++))
do
i=$(echo "$astart + $astep * $b" | bc) # use bash calculator
cat > STRU << !
ATOMIC_SPECIES
Ni 58.693 Ni_ONCV_PBE-1.0.upf

NUMERICAL_ORBITAL
Ni_gga_9au_100Ry_4s2p2d1f.orb

LATTICE_CONSTANT
1.889726

LATTICE_VECTORS
$i  0.000000000000  0.000000000000
0.000000000000  $i  0.000000000000
0.000000000000  0.000000000000  $i

ATOMIC_POSITIONS
Direct

Ni
2.0
4
0.0 0.0 0.0
0.5 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.5
!

mpirun -n 32 abacus >> result
energy_line=$(grep "etot" OUT.*/running_scf.log)
E_value=$(echo $energy_line | awk '{for(i=1;i<=NF;i++) if ($i=="is") print $(i+1)}')
E=$(printf "%.8f" $E_value) # Keep eight decimal places
echo "$i,$E" >> Ni_fcc_energy.csv
done

代码

文本

该过程使用32核需要运行30分钟，感兴趣的用户可以自行计算。去除下面的#号注释即可运行。

代码

文本

[7]

#. lattice_constant.sh

代码

文本

运行完毕后，晶格常数和能量都输出在Ni_fcc_energy.csv文件中：

代码

文本

[9]

cat Ni_fcc_energy.csv

lat,ene
3.00,-16768.13483400
3.04,-16770.05924100
3.08,-16771.71528400
3.12,-16773.13166800
3.16,-16774.33095700
3.20,-16775.33304400
3.24,-16776.15902800
3.28,-16776.82763200
3.32,-16777.35683300
3.36,-16777.76238600
3.40,-16778.05503200
3.44,-16778.25080900
3.48,-16778.36217100
3.52,-16778.39511200
3.56,-16778.36122700
3.60,-16778.26955900
3.64,-16778.12907600
3.68,-16777.94620800
3.72,-16777.72523800
3.76,-16777.47636900
3.80,-16777.20101800
3.84,-16776.90549900
3.88,-16776.59369800
3.92,-16776.26749900
3.96,-16775.92957500
4.00,-16775.58235100

代码

文本

2. Birch-Murnaghan方程拟合

代码

文本

很多实际的应用，如测定相变、膨胀系数、弹性、热电效应、组分变化、应力等，都要求精确地测定晶格常数。理论计算平衡晶格常数多采用固体的状态方程（equation of state, EOS），该方程表示晶体体系的总能量E随体系体积V的变化，即具有E(V)的形式。晶体的状态方程对于基础科学和应用科学具有重要的意义，比如平衡体积（V0）、体弹性模量（B₀）以及它的一阶导数（B₀'），这些可测的物理量与晶体的状态方程直接相关。在高压下状态方程可用好几种不同的函数形式来描述，比如Murnaghan方程、Birch-Murnaghan(BM)方程和普适方程等。
三阶Birch-Murnaghan方程如下：

代码

文本

$E (V) = E_{0} + \frac{9 V _{0} B _{0}}{16} ⎩ ⎨ ⎧ [(\frac{V _{0}}{V})^{2/3} - 1]^{3} B_{0}^{'} + [(\frac{V _{0}}{V})^{2/3} - 1]^{2} [6 - 4 (\frac{V _{0}}{V})^{2/3}] ⎭ ⎬ ⎫$

代码

文本

由该方程可以看出，当V=V₀时，E=E₀，这就分别是平衡晶格常数及其对应的能量。在实际操作中，我们就可以在平衡晶格常数a₀附近计算得到若干E-V数据点，经曲线拟合得到体系最低能量对应的晶格常数。采用最小二乘法拟合可以得到V₀, E₀, B₀, B₀'。拟合脚本Birch_fitting.py如下：

代码

文本

[10]

cat Birch_fitting.py

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.optimize import leastsq

# Birch-Murnaghan EOS
def Birch_Murnaghan(p,x):
        E0,V0,B0,B1=p
        return E0+(9*V0*B0/16)*((((V0/(x**3))**(2/3)-1)**3)*B1+(((V0/(x**3))**(2/3)-1)**2)*(6-4*(V0/(x**3))**(2/3)))
# error
def error(p,x,y):
        return Birch_Murnaghan(p,x)-y

def main():
        # read data
        filename="Ni_fcc_energy.csv"
        DataFrame=pd.read_csv(filename)
        x=np.array(DataFrame['lat'])
        y=np.array(DataFrame['ene'])
        # fitting by least square method
        p0=[ min(y), min(x)**3, 10, 1 ] # initial args
        para=leastsq(error, p0, args=(x,y)) # optimized args
        y_fitted=Birch_Murnaghan(para[0],x) # fitting curves

        print ("Emin=", para[0][0], "opt_lat", para[0][1]**(1/3))
        print ("B0=", para[0][2], "B1=", para[0][3])

        # plot
        plt.figure(figsize=(8,6))
        plt.scatter(x,y,20, color='r', label='ABACUS computed data')
        plt.plot(x,y_fitted, '-b', label='Birch-Murnaghan Fitting')
        plt.xlabel('Lattice constant/Å', fontsize=20)
        plt.ylabel('Energy/eV',fontsize=20)
        plt.legend(fontsize=10)
        plt.savefig("Birch_fitting.jpg", dpi=300)
        print (para[0])

if __name__=='__main__':
        main()

代码

文本

[11]

python3 Birch_fitting.py

Emin= -16778.37920200902 opt_lat 3.523123194799099
B0= 1.2606500039348976 B1= 4.8107190112115115
[-1.67783792e+04  4.37304039e+01  1.26065000e+00  4.81071901e+00]

代码

文本

拟合结果如下图所示。得到平衡晶格常数为3.523 Å，与实验值3.52 Å非常接近。通过拟合还可以得到Ni的FCC结构的体模量B₀=1.261 eV/Å，根据1.0 eV/Å = 160.2 GPa，转化后为202.0122 GPa（实验值约为180 GPa，误差在12.3%左右）。

代码

文本

代码

文本

恭喜你完成了第三个ABACUS计算模拟实例🎉

代码

文本

如果你还想尝试更多的计算实例，可以点击下方了解ABACUS计算模拟合集链接： ABACUS计算模拟实例 | 概述

代码

文本

notebook

ABACUS

计算材料学

ABACUS使用教程

notebookABACUS计算材料学ABACUS使用教程

点个赞吧

本文被以下合集收录

ABACUS@计算材料学 | 计算模拟实例

weiqingzhou@whu.edu.cn

更新于 2024-08-09

13 篇8 人关注

ABACUS计算模拟实例

FermiNoodles

更新于 2024-07-15

13 篇1 人关注