Abacus第一性原理计算初步教程——以铁单质及其简单化合物为例

空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

Abacus第一性原理计算初步教程——以铁单质及其简单化合物为例

notebook

中文

作业

ABACUS使用教程

notebook中文作业ABACUS使用教程

why

发布于 2024-05-18

推荐镜像 :ABACUS:3.6.3

推荐机型 :c2_m4_cpu

数据集

Fe_FeO(v1)

双击即可修改

代码

文本

章节一：

密度泛函理论（Density Functional Theory）基础

在许多物理学和工程领域，取得科学和技术进步的关键在于能够从原子或分子尺度，理解并调控物质的性质。对于描述原子和分子量子行为的基本关系式(薛定湾方程)，密度泛函理论 (Density Functional Theory DFT) 是一个非常好的求解方法。

1.1 DFT起源

如果我们仅考虑量子力学基本原理，当我们将要描述规则排列原子集合的性质时，需要知道原子的能量，以及移动这些原子时它们的能量将如何变化。

考虑到原子核的质量远大于单个电子质量。原子核中的每个质子或中子要比一个电子的质量大1800倍。这意味着电子对环境变化的响应远比原子核快得多。因此，我们可以把这一问题分为两部分。

第一部分，我们固定原子核的位置，求解描述电子运动的方程组。对于在给定原子核势场中运动的一系列电子，我们可以得到能量最低的电子构型，称为电子的基态(Ground State)，而将原子核和电子分为各自独立数学问题的是Born -Oppen-heimer近似。如果我们有M个在 $R_{1}$ ,…, $R_{M}$ 位置的原子核，那么我们可以把基态能量E表示为这些原子核位置的函数，即E( $R_{1}$ ,…, $R_{M}$ )。这个方程亦称为这些原子的绝热势能面(Adiabatic Potential Energy Surface)通过该势能面可以解决移动原子时，材料的能量如何变化的问题。

在多电子与多原子核交互作用体系下，定态薛定谔方程写为 $[- \frac{h ^{2}}{2 m} i = 1 \sum N \nabla_{i}^{2} + i = 1 \sum N V (r_{i}) + i = 1 \sum N j < i \sum U (r_{i}, r_{j})] ψ = E ψ$

方括号中的三项依次分别是每个电子的动能、每个电子与所有原子核之间的作用能、不同电子之间的作用能。对于我们所选定的哈密顿量， $ψ$ 是电子波函数，它是N个电子每个电子空间坐标的函数。尽管电子波函数是所有N个电子坐标的函数，但仍可以将其近似表达为单个电子波函数的乘积，即 $ψ = ψ_{1} (r) ψ_{2} (r) \dots ψ_{N} (r)$ 。该波函数表达式就是所谓的Hartree乘积。

假如我们的研究对象是一个有22个电子的 $C O_{2}$ 分子，其全波函数是一个66维函数；再假如研究对象为 100个Pt原子的纳米团簇，其全波函数则超过了23000维，这就是为什么求解实际材料的Schrödinger方程是如此困难。

1.2 基本原理与自洽求解

密度泛函理论全部建立在由Kohn和Hohenberg所证明的两个基本数学定理，以及由Kohn和Sham在1960年代中期所推演的一套方程的基础上。其中，由Kohn和Hohenberg所证明的第一个定理是从Schrödinger方程得到的基态能量是电荷密度的唯一函数。该定理表明：在基态波函数和基态电荷密度之间，存在一一对应的关系。也就是基态能量E可以表示为E[n(r)]，其中n(r)是电荷密度。

另一个表述这一结果的方式是基态电荷密度唯一决定了基态的所有性质，包括能量和波函数。也就是说可以通过找到含有三个空间变量的电荷密度函数来求Schrödinger方程，而不用求解得到含有3N个变量的波函数。

第一个Hohenberg-Kohn定理严格证明了存在一个可用来求解Schrödinger方程的电荷密度泛函，第二个Hohenberg-Kohn定理给出了这个泛函的一个重要特征：使整体泛函最小化的电荷密度就是对应于 Schrödinger方程完全解的真实电荷密度。如果已知这个"真实的"泛函形式，那就能够通过不断调整电荷密度直到由泛函所确定的能量达到最小化，并可以找到相应的电荷密度。

将Hohenberg-Kohn定理所描述的泛函写成单电子波函数 $ψ_{i} (r)$ 的形式。能量泛函可以写为 $E [ψ_{i}] = E_{kn o w n} [ψ_{i}] + E_{XC} [ψ_{i}]$

其中 $E_{kn o w n} [ψ_{i}] = - \frac{h ^{2}}{m} i \sum \int ψ_{i}^{*} \nabla^{2} ψ_{i} d^{3} r + \int V (r) n (r) d^{3} r + \frac{e ^{2}}{2} \int\int \frac{n ( r ) n ( r ^{'} )}{∣ r - r ^{'} ∣} d^{3} r d^{3} r^{'} + E_{i o n}$

在上式右侧，依次分别是:电子的动能，电子和原子核之间的库伦作用，电子之间的库仑作用，原子核之间的库仑作用。 $E_{XC} [ψ_{i}]$ 是交换关联泛函，它所定义的是没有包括在known项中的量子力学效应。

Kohn和Sham据此给出了一套方程，求解正确的电荷密度可以表示为求解这套方程，而其中每个方程都只与一个电子有关。

K-S方程的表达式为 $[- \frac{h ^{2}}{m} \nabla^{2} + V (r) + V_{H} (r) + V_{XC} (r)] ψ_{i} (r) = ϵ_{i} ψ_{i} (r)$

在方程左侧，有三个势能项，其中第一项定义的是一个电子与所有原子核之间的相互作用；第二项也称为Hatree势能，可以写为 $V_{H} (r) = e^{2} \int \frac{n ( r ^{'} )}{∣ r - r ^{'} ∣} d^{3} r^{'}$ 这个势能描述的是一个Kohn-Sham 方程所考虑的单个电子，与该问题中全部电子所产生的总电荷密度之间的排斥作用。 $V_{XC}$ 可以在形式上可以表示为交换关联能的"泛函导数"。

为了求解Kohn-Sham方程，通常用迭代算法来处理，其过程简述如下:

(1)定义一个初始的、尝试性的电荷密度 $n (r)$ (2)求解由尝试性的电荷密度所确定的Kohn-Sham方程，得到单电子波函数 $ψ_{i} (r)$ (3)计算由第(2)单粒子波函数所确定的电荷密度 (4)比较计算得到的电荷密度 $n_{K S} (r)$ 和在求解方程时使用的电荷密度。如果两个电荷密度相同，则这就是基态电荷密度，并可将其用于计算总能。如果两个电荷密度不同，则用某种方式对尝试性电荷密度进行修正，然后再从第 (2)步重新开始。这也就是自治 (Self-Consistent)求解过程。

代码

文本

章节二：

Fe单质与FeO

本教程采用的案例主要研究对象是铁单质与氧化亚铁。此处简单介绍铁元素的性质与两种物质的性质。

铁元素的原子序数为26，原子量为55.845u，属于第8族，第4区，d区元素，电子排布为 $[A r] 3 d^{6} 4 s^{2}$ ，原子半径为126pm。铁有四种已知的同素异形体，它们的名称通常表示为 $α, γ, δ, ϵ$ 。前三种铁的同素异形体可以在常压下存在。 $α$ -铁存在于912摄氏度以下，呈体心立方晶格，与碳的固溶体称作 $α$ -铁素体。 $γ$ -Fe，存在于912到1394摄氏度之间，呈面心立方晶格，与碳的固溶体称作奥氏体，由于面心比体心排列紧密，所以 $α$ -Fe转化为 $γ$ -Fe时，体积会收缩。 $δ$ -Fe，存在于1400到1535摄氏度之间，呈体心立方晶格，与碳的固溶体称作 $δ$ -铁素体。 $ϵ$ -Fe呈密排六方结构，又称六铁。在13GPa的压力以及室温下， $α$ -铁粉末会转变为 $ϵ$ -Fe。

铁是典型的过渡金属，可以形成多种氧化态，并且有许多配合物和有机金属化合物。铁可看作是所有过渡金属的原型。铁主要的氧化态为+2和+3。在有机铁化合物中，铁的氧化态可以达到+1、0、−1甚至是−2。

氧化亚铁是一种无机物，化学式FeO。是铁的氧化物之一。其外观呈黑色粉末，由氧化态为+2价的铁与氧共价结合。它的矿物形式为方铁矿。氧化亚铁属于非整比化合物，其中铁和氧元素的比例会发生变化，范围从 $F e_{0.84} O$ 到 $F e_{0.95} O$ 。

氧化亚铁属于立方晶系，类似于氯化钠晶体结构，每个铁原子周围连接着6个氧原子形成八面体型配位，每个氧原子周围也以同样的情况连接着6个铁原子。出现非整比化合物的原因在于，由于二价铁很容易被氧化为三价铁，导致FeO中少量二价铁被替换为三价铁，占据晶格中的四面体空隙。

代码

文本

章节三：

Abacus自洽场计算(LCAO基组)

ABACUS(Atomic-orbtial Based Ab-initio Computation at UStc)是一个由中国科学技术大学量子信息与超级计算中心重点实验室-CAS网络中心开发开源的计算机代码包，旨在通过第一性原理模拟大规模电子结构。如其名称所示，ABACUS主要使用数值列表中的原子中心轨道作为基函数来扩展电子波函数。

在DFT计算中，首要任务是分析体系的势能面。通过SCF计算，可以得到一个体系的基态结构和基态能量。SCF计算也是核心过程，通过优化电荷密度来寻找势能最低点。

3.1 输入文件

在进行计算之前需要准备五组输入文件。

INPUT 包含了计算过程中所需的各种参数，定义和控制计算任务
STRU 结构文件，包含了原子种类、原子位置、晶格常数以及晶格向量等信息
KPT 包含了布里渊区积分所需的k点信息
*.upf 包含了原子的赝势信息
*.orb 包含了原子轨道的数值表示

3.1.1 STRU文件

STRU 文件包含晶格几何信息，赝势和数值轨道文件的名称和/或位置，以及关于系统的结构信息。

以用作实例计算的Fe的STRU文件为例

ATOMIC_SPECIES
Fe Fe_ONCV_PBE-1.0.upf   #label; mass; pseudo_file

NUMERICAL_ORBITAL
Fe_gga_9au_100Ry_4s2p2d1f.orb

LATTICE_CONSTANT
6                              # 1.8897259886 Bohr = 1.0 Angstrom

LATTICE_VECTORS
1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00
ATOMIC_POSITIONS
Direct

Fe
0.0
2                             
0.00 0.00 0.00 m 1 1 1 mag 3.0      # 按原子设置磁矩，对该Fe原子设置初始磁矩3.0 
0.50 0.50 0.50 m 1 1 1 mag -3.0

参数解释

ATOMIC_SPECIES

从左到右依次为原子种类、相对原子质量、赝势文件

NUMERICAL_ORBITAL

原子轨道文件，数值原子轨道只需要用于LCAO计算

LATTICE_CONSTANT

系统的晶格常数，单位为玻尔，1.8897259886 Bohr=1.0 Angstrom

LATTICE_VECTORS

晶胞的晶格矢量，是一个3乘3的矩阵。这里给出的晶格矢量是按晶格常数缩放的。

ATOMIC_POSITIONS

本节规定了单个原子的位置和其他信息。第一行表示给定原子位置的方法，包含以下选项：

Direct：下面原子位置的坐标将是分数坐标。

Cartesian：笛卡尔坐标，单位为“LATTICE_CONSTANT”。

Cartesian_au：以玻尔为单位的笛卡尔坐标，与LATTICE_CONSTANT=1.0的笛卡尔坐标设置相同。

Cartesian_angstrom：笛卡尔坐标，单位为埃，与LATTICE_CONSTANT=1.8897261254577828的笛卡尔坐标设置相同。

三个数字，取值为0或1，控制原子在几何弛豫计算中的移动方式。数字0意味着这个原子不允许该方向上移动，1意味着该方向上允许移动

设置每个原子的初始磁化强度。在共线的情况下，只应给出一个数字。

代码

文本

3.1.2 INPUT文件

INPUT文件设置了所有关键参数，指导ABACUS如何进行计算和输出的内容

# INPUT_PARAMETERS
suffix 			Fe
ntype 			1
pseudo_dir 		./PP_ORB
orbital_dir 	        ./PP_ORB
ecutwfc 		100     # Rydberg
scf_thr 		1e-6    # Rydberg
basis_type 		lcao
calculation 	        scf

参数解释

suffix：系统名称，默认为ABACUS
ntype：单位晶胞中元素的种类数目
pseudo_dir：提供赝势文件的目录
orbital_dir：提供轨道文件的目录
ecutwfc：波函数展开的平面波能量截止（单位：Rydberg）
scf_thr：电荷密度收敛阈值（单位：Rydberg）
basis_type：用于展开电子波函数的基函数类型
calculation：ABACUS要执行的计算类型

INPUT 参数设置的第一原则是简洁明了。参数不是越多越好，对于大多数场景，默认值已不失为一个很好的选择。

代码

文本

3.1.3 KPT文件

最后一个必需的输入文件称为 KPT，此文件包含关于布里渊区采样的k点网格设置的信息。

K_POINTS
0           #k点总数，“0”表示自动生成
Gamma        #Monkhorst-Pack方法的类型，`Gamma'或`MP'
4 4 4 0 0 0      #前三个数字：沿着倒数向量的细分；后三个数字：网格的偏移

代码

文本

3.1.4 赝势文件与原子轨道

在ABACUS网站https://abacus.ustc.edu.cn/pseudo/list.htm 可以下载赝势文件以及与ABACUS相对应的优化原子基组。

目前，ABACUS支持UPF格式的守恒赝势。

原子基组是针对特定赝势进行优化的，这意味着当赝势发生变化时，必须重新生成原子基组。当使用原子基组时，应选择参考能量最接近计算中所用能量的基组。还要注意再次确认INPUT文件中赝势文件与原子轨道文件的路径设置正确。

代码

文本

3.2 运行计算与输出结果

进入已经准备好输入文件的文件夹，执行以下命令即可快速运行ABACUS计算

代码

文本

[ ]

cd /Fe_FeO(v1)/ABACUS

# OMP_NUM_THREADS=1 表示使用单线程，如果你的机器配置比较高，可以使用多线程，比如 4 线程，就可以写成 OMP_NUM_THREADS=4

# mpirun -n 后面的数字表示计算所使用的 CPU 核心数，这里使用 2 个核心，你可以根据你的机器配置进行修改。

OMP_NUM_THREADS=1 mpirun -n 2 abacus

代码

文本

输出文件储存在当前文件夹下的 OUT.Fe文件夹中，输出文件夹的名称为 OUT.，其中由INPUT文件中的suffix项目指定。

计算结果主要储存在OUT.Fe/running_scf.log文件中。

如果ABACUS自洽场计算成功完成，能量会输出在OUT.Fe/running_scf.log文件中：

Total Magnetism on atom:  Fe         3.032520051
Total Magnetism on atom:  Fe        -3.032224908



 Warning_Memory_Consuming allocated:  Force::dS_K 5.656379700 MB

 Warning_Memory_Consuming allocated:  Stress::dHr 5.656379700 MB

 Warning_Memory_Consuming allocated:  Stress::dSR 11.31275940 MB

 Warning_Memory_Consuming allocated:  Force::dTVNL 5.656379700 MB
 correction force for each atom along direction 1 is -4.288746584e-14
 correction force for each atom along direction 2 is -3.215925309e-15
 correction force for each atom along direction 3 is -7.696704033e-14
------------------------------------------------------------------------------------------
TOTAL-FORCE (eV/Angstrom)                                                               
------------------------------------------------------------------------------------------
Fe1                               0.0000000000         0.0000000000         0.0000000000
Fe2                               0.0000000000         0.0000000000         0.0000000000
------------------------------------------------------------------------------------------
----------------------------------------------------------------
TOTAL-STRESS (KBAR)                                           
----------------------------------------------------------------
     -238.7512893809        -0.0000000000         0.0000000000
       -0.0000000000      -238.7512893809        -0.0000000000
        0.0000000000        -0.0000000000      -238.7512893809
----------------------------------------------------------------
 TOTAL-PRESSURE: -238.751289 KBAR



 --------------------------------------------
 !FINAL_ETOT_IS -6438.6065590232401519 eV
 --------------------------------------------


TIME STATISTICS
-------------------------------------------------------------------------------------
     CLASS_NAME                  NAME            TIME(Sec)  CALLS   AVG(Sec) PER(%)
-------------------------------------------------------------------------------------
                      total                       53.04          11   4.82   100.00
Driver                reading                      0.07           1   0.07     0.13
Input                 Init                         0.06           1   0.06     0.11
Input_Conv            Convert                      0.00           1   0.00     0.00
Driver                driver_line                 52.97           1  52.97    99.87
UnitCell              check_tau                    0.00           1   0.00     0.00
PW_Basis_Sup          setuptransform               0.01           1   0.01     0.03
PW_Basis_Sup          distributeg                  0.00           1   0.00     0.00
mymath                heapsort                     0.00           7   0.00     0.00
Symmetry              analy_sys                    0.00           1   0.00     0.00
PW_Basis_K            setuptransform               0.00           1   0.00     0.00
PW_Basis_K            distributeg                  0.00           1   0.00     0.00
PW_Basis              setup_struc_factor           0.00           1   0.00     0.01
NOrbital_Lm           extra_uniform                0.08           9   0.01     0.15
Mathzone_Add1         SplineD2                     0.00           9   0.00     0.00
Mathzone_Add1         Cubic_Spline_Interpolation   0.01           9   0.00     0.01
Mathzone_Add1         Uni_Deriv_Phi                0.07           9   0.01     0.13
ppcell_vl             init_vloc                    0.00           1   0.00     0.01
Ions                  opt_ions                    52.12           1  52.12    98.27
ESolver_KS_LCAO       Run                         37.50           1  37.50    70.71
ESolver_KS_LCAO       beforescf                    0.30           1   0.30     0.57
ESolver_KS_LCAO       beforesolver                 0.25           1   0.25     0.48
ESolver_KS_LCAO       set_matrix_grid              0.22           1   0.22     0.41
atom_arrange          search                       0.00           1   0.00     0.00
Grid_Technique        init                         0.21           1   0.21     0.40
Grid_BigCell          grid_expansion_index         0.00           2   0.00     0.01
Record_adj            for_2d                       0.01           1   0.01     0.01
Grid_Driver           Find_atom                    0.00          68   0.00     0.00
LCAO_Hamilt           grid_prepare                 0.00           1   0.00     0.00
Veff                  initialize_HR                0.00           1   0.00     0.00
OverlapNew            initialize_SR                0.00           1   0.00     0.00
EkineticNew           initialize_HR                0.00           1   0.00     0.00
NonlocalNew           initialize_HR                0.01           1   0.01     0.01
Charge                set_rho_core                 0.00           1   0.00     0.00
Charge                atomic_rho                   0.01           1   0.01     0.02
PW_Basis_Sup          recip2real                   0.05         135   0.00     0.09
PW_Basis_Sup          gathers_scatterp             0.02         135   0.00     0.03
Potential             init_pot                     0.03           1   0.03     0.05
Potential             update_from_charge           0.23          10   0.02     0.44
Potential             cal_fixed_v                  0.00           1   0.00     0.00
PotLocal              cal_fixed_v                  0.00           1   0.00     0.00
Potential             cal_v_eff                    0.23          10   0.02     0.44
H_Hartree_pw          v_hartree                    0.01          10   0.00     0.02
PW_Basis_Sup          real2recip                   0.04         143   0.00     0.08
PW_Basis_Sup          gatherp_scatters             0.01         143   0.00     0.03
PotXC                 cal_v_eff                    0.22          10   0.02     0.41
XC_Functional         v_xc                         0.22          10   0.02     0.41
Potential             interpolate_vrs              0.00          10   0.00     0.00
Symmetry              rhog_symmetry                0.04          20   0.00     0.07
Symmetry              group fft grids              0.01          20   0.00     0.02
H_Ewald_pw            compute_ewald                0.00           1   0.00     0.00
Charge_Mixing         init_mixing                  0.00           1   0.00     0.00
HSolverLCAO           solve                       36.82           9   4.09    69.41
HamiltLCAO            updateHk                    22.28         180   0.12    42.01
OperatorLCAO          init                        21.16         540   0.04    39.90
Veff                  contributeHR                19.41          18   1.08    36.59
Gint_interface        cal_gint                    38.53          29   1.33    72.63
Gint_interface        cal_gint_vlocal             19.25          18   1.07    36.29
Gint_Tools            cal_psir_ylm                 3.56       10935   0.00     6.71
Gint_k                transfer_pvpR                0.16          18   0.01     0.31
OverlapNew            calculate_SR                 0.91           1   0.91     1.72
OverlapNew            contributeHk                 0.20         180   0.00     0.39
EkineticNew           contributeHR                 0.91          18   0.05     1.72
EkineticNew           calculate_HR                 0.91           1   0.91     1.72
NonlocalNew           contributeHR                 0.54          18   0.03     1.02
NonlocalNew           calculate_HR                 0.53           1   0.53     1.00
OperatorLCAO          contributeHk                 0.29         180   0.00     0.55
HSolverLCAO           hamiltSolvePsiK              0.85         180   0.00     1.61
DiagoElpa             elpa_solve                   0.77         180   0.00     1.46
ElecStateLCAO         psiToRho                    13.68           9   1.52    25.79
elecstate             cal_dm                       0.06          10   0.01     0.12
psiMulPsiMpi          pdgemm                       0.06         200   0.00     0.12
DensityMatrix         cal_DMR                      0.38          10   0.04     0.72
Local_Orbital_wfc     wfc_2d_to_grid               0.01         200   0.00     0.01
Gint                  transfer_DMR                 0.12           9   0.01     0.22
Gint_interface        cal_gint_rho                13.14           9   1.46    24.78
Charge_Mixing         get_drho                     0.00           9   0.00     0.00
Charge                mix_rho                      0.02           8   0.00     0.05
Charge                Pulay_mixing                 0.01           8   0.00     0.01
ESolver_KS_LCAO       out_deepks_labels            0.00           1   0.00     0.00
LCAO_Deepks_Interface out_deepks_labels            0.00           1   0.00     0.00
HamiltLCAO            updateSk                     0.02          20   0.00     0.04
Force_Stress_LCAO     getForceStress              14.62           1  14.62    27.56
Forces                cal_force_loc                0.00           1   0.00     0.00
Forces                cal_force_ew                 0.00           1   0.00     0.00
Forces                cal_force_cc                 0.00           1   0.00     0.00
Forces                cal_force_scc                0.00           1   0.00     0.01
Stress_Func           stress_loc                   0.01           1   0.01     0.01
Stress_Func           stress_har                   0.00           1   0.00     0.00
Stress_Func           stress_ewa                   0.00           1   0.00     0.00
Stress_Func           stress_cc                    0.00           1   0.00     0.00
Stress_Func           stress_gga                   0.01           1   0.01     0.02
Force_LCAO_k          ftable_k                    14.59           1  14.59    27.51
Force_LCAO_k          allocate_k                   4.33           1   4.33     8.17
LCAO_gen_fixedH       b_NL_mu_new                  1.80           1   1.80     3.39
Force_LCAO_k          cal_foverlap_k               0.06           1   0.06     0.10
Force_LCAO_k          cal_edm_2d                   0.05           1   0.05     0.09
DensityMatrix         sum_DMR_spin                 0.00           1   0.00     0.00
Force_LCAO_k          cal_ftvnl_dphi_k             0.01           1   0.01     0.02
Force_LCAO_k          cal_fvl_dphi_k               6.14           1   6.14    11.57
Gint_interface        cal_gint_force               6.14           2   3.07    11.57
Gint_Tools            cal_dpsir_ylm                2.45         810   0.00     4.61
Gint_Tools            cal_dpsirr_ylm               0.64         810   0.00     1.21
Force_LCAO_k          cal_fvnl_dbeta_k_new         4.05           1   4.05     7.64
ModuleIO              write_istate_info            0.01           1   0.01     0.02
-------------------------------------------------------------------------------------

 NAME---------------|MEMORY(MB)--------
               total        92.6108
         Stress::dSR        22.6255
         Force::dS_K        11.3128
         Stress::dHr        11.3128
        Force::dTVNL        11.3128
        pvpR_reduced        10.1893
           LOC::DM_R        10.1893
TwoCenterTable: Kinetic         3.2307
TwoCenterTable: Overlap         3.2307
TwoCenterTable: Nonlocal         2.2223
 -------------   < 1.0 MB has been ignored ----------------
 ----------------------------------------------------------

 Start  Time  : Fri May 17 22:26:14 2024
 Finish Time  : Fri May 17 22:27:07 2024
 Total  Time  : 0 h 0 mins 53 secs

以上是Fe元素的LCAO基组自洽计算实例，对于PW基组的计算，只需在INPUT文件中将basis_type标签从lcao 更改为pw，同时无需在STRU文件中的NUMERICAL_ORBITAL下提供轨道文件。

代码

文本

章节四：

ABACUS实现结构优化/晶格弛豫/几何优化

4.1 结构优化

在第一性原理计算中，晶格弛豫（Lattice relaxation）指的是一种过程，即通过调整晶格内原子的位置和晶格参数，使晶体结构达到最低能量状态。

晶格弛豫的主要目的是找到晶体结构的稳定构型，这对于预测和理解材料的性能至关重要。在这个过程中，计算机模拟会根据所采用的理论框架来迭代地优化原子位置和晶格参数，最终得到一个能量最低的结构。这个结构通常被认为是晶体在实

与上一章节的电子自洽迭代中类似，我们需要准备以下五个类型地输入文件。其中，我们只需要调整INPUT文件中的设置项，其它文件与上一章节中的输入文件相同。接下来我们以上一章节的Fe单质与Fe元素的简单化合物FeO为例进行晶格结构的优化。

下面是FeO对应的INPUT文件

INPUT_PARAMETERS
suffix                  FeO
ntype                   2
nelec                   0.0
pseudo_dir              ./
orbital_dir             ./
ecutwfc                 100             # Rydberg
scf_thr                 1e-4		# Rydberg
basis_type              lcao
calculation             cell-relax	# this is the key parameter telling abacus to do a optimization calculation
force_thr_ev		0.01		# the threshold of the force convergence, in unit of eV/Angstrom
stress_thr		2		# the threshold of the stress convergence, in unit of kBar
relax_nmax		100		# the maximal number of ionic iteration steps
out_stru		1

加入项目包括

calculation: 计算类型，这里设置为"cell-relax"，表示进行晶胞优化计算。
force_thr_ev: 力收敛阈值，用于判断优化过程是否收敛。这里设置为0.01 eV/Angstrom，表示当所有原子上的力小于0.01 eV/Angstrom时，认为已经收敛。
stress_thr: 应力收敛阈值，用于判断晶胞优化过程中应力是否收敛。默认值：0.5 kBar
relax_nmax: 优化过程中的最大迭代次数，这里设置为100，表示最多进行100次迭代。
out_stru: 表示是否输出优化后的结构信息。这里设置为1，表示在优化结束后将优化后的结构信息输出到文件。

其它可用的INPUT项目

smearing_method：对于金属体系或者费米面附近较复杂的体系，电子自洽迭代方法往往不容易收敛，这个时候可以选择给电子提供一个光滑的占据函数，用来调节程序计算电荷密度和总电子数时候的电子占据函数，这对于费米面附近的电子态尤为重要。我们这里称为 "smearing method"，或者称为展宽方法。具体来说，对于金属体系可以选择mp、mp2，也可以使用gauss/gaussian。
smearing_sigma：给定展宽方法的能量范围（单位：Ry），默认是 0.015 Ry，我们按照通常情况给定 "0.01"。
mixing_type：进行新旧电荷密度混合时选用的方法，可选 "plain"、"pulay"、"broyden"，默认选择 "broyden"算法。
mixing_beta：新旧电荷密度混合时，新电荷的比例，不同系统取值不同。能带大于 1 eV 的体系，设为 0.7；能带小于1eV的金属和过渡金属，设为0.2。这个参数越大，收敛得越快，但不收敛的风险也会变大，不过取值大小并不会改变基态能量的结果。
mixing_gg0：电子迭代过程中，可能出现混合电荷密度不收敛的情况，此时可以通过一种叫 Kerker Mixing 的方法来加快收敛。这个参数代表 Kerker 方法中调整电荷密度的尺度。

代码

文本

[ ]

cd /personal/ABACUS1

mpirun -n 2 abacus

代码

文本

计算完成后，输出文件储存在当前文件夹下的OUT.FeO文件夹中，输出文件夹的名称为 OUT.，其中由INPUT文件中的suffix值指定。

最终优化的结构可以在STRU_NOW.cif和OUT.FeO/running_cell-relax.log中找到。

STEP OF RELAXATION : 2
 -------------------------------------------
 DONE : LOCAL POTENTIAL Time : 138.334282 (SEC)

 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 |                                                                    |
 | Doing symmetry analysis:                                           |
 | We calculate the norm of 3 vectors and the angles between them,    |
 | the type of Bravais lattice is given. We can judge if the unticell |
 | is a primitive cell. Finally we give the point group operation for |
 | this unitcell. We use the point group operations to do symmetry |
 | analysis on given k-point mesh and the charge density.             |
 |                                                                    |
 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<




 LATTICE VECTORS: (CARTESIAN COORDINATE: IN UNIT OF A0)
            +4.026184           +0.000000           +0.000000
            +0.000000           +4.026184           -0.000000
            +0.000000           -0.000000           +4.026184
                       right hand lattice = 1

 The lattice vectors have been changed (STRU_SIMPLE.cif)

(for optimal symmetric configuration:)
                             BRAVAIS TYPE = 1
                     BRAVAIS LATTICE NAME = 01. Cubic P (simple)
                                    ibrav = 1
                                    IBRAV = 1
                                  BRAVAIS = SIMPLE CUBIC
                       LATTICE CONSTANT A = 4.026184
                       right hand lattice = 1
(for primitive cell:)
                                    IBRAV = 3
                                  BRAVAIS = FACE CENTERED CUBIC
                       LATTICE CONSTANT A = 4.026184
optimized lattice volume: 65.265077
optimized primitive cell volume: 16.316269
Original cell was built up by 4 primitive cells.
                        ROTATION MATRICES = 48
              PURE POINT GROUP OPERATIONS = 48
                   SPACE GROUP OPERATIONS = 48
                              POINT GROUP = O_h
               POINT GROUP IN SPACE GROUP = O_h
 DONE : SYMMETRY Time : 138.474755 (SEC)


 SETUP K-POINTS
                                    nspin = 1
K-POINTS DIRECT COORDINATES
 KPOINTS     DIRECT_X     DIRECT_Y     DIRECT_Z WEIGHT
       1   0.00000000   0.00000000   0.00000000 0.0312
       2   0.00000000   0.25000000   0.00000000 0.1875
       3   0.00000000   0.00000000   0.50000000 0.0938
       4   0.25000000   0.00000000   0.25000000 0.3750
       5   0.00000000   0.25000000   0.50000000 0.3750
       6   0.00000000   0.50000000   0.50000000 0.0938
       7   0.25000000   0.25000000  -0.25000000 0.2500
       8   0.25000000   0.50000000   0.25000000 0.3750
       9   0.50000000   0.25000000   0.50000000 0.1875
      10   0.50000000   0.50000000   0.50000000 0.0312

           k-point number in this process = 10
       minimum distributed K point number = 10

K-POINTS CARTESIAN COORDINATES
 KPOINTS  CARTESIAN_X  CARTESIAN_Y  CARTESIAN_Z WEIGHT
       1   0.00000000   0.00000000   0.00000000 0.0312
       2   0.00000000   0.06209354   0.00000000 0.1875
       3   0.00000000   0.00000000   0.12418707 0.0938
       4   0.06209354   0.00000000   0.06209354 0.3750
       5   0.00000000   0.06209354   0.12418707 0.3750
       6   0.00000000   0.12418707   0.12418707 0.0938
       7   0.06209354   0.06209354  -0.06209354 0.2500
       8   0.06209354   0.12418707   0.06209354 0.3750
       9   0.12418707   0.06209354   0.12418707 0.1875
      10   0.12418707   0.12418707   0.12418707 0.0312

K-POINTS DIRECT COORDINATES
 KPOINTS     DIRECT_X     DIRECT_Y     DIRECT_Z WEIGHT
       1   0.00000000   0.00000000   0.00000000 0.0312
       2   0.00000000   0.25000000   0.00000000 0.1875
       3   0.00000000   0.00000000   0.50000000 0.0938
       4   0.25000000   0.00000000   0.25000000 0.3750
       5   0.00000000   0.25000000   0.50000000 0.3750
       6   0.00000000   0.50000000   0.50000000 0.0938
       7   0.25000000   0.25000000  -0.25000000 0.2500
       8   0.25000000   0.50000000   0.25000000 0.3750
       9   0.50000000   0.25000000   0.50000000 0.1875
      10   0.50000000   0.50000000   0.50000000 0.0312
 DONE : INIT K-POINTS Time : 138.475703 (SEC)

 Find the file, try to read charge from file.
 read in fermi energy = 0.000000
 NEW-OLD atomic charge density approx. for the potential !




 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 |                                                                    |
 | Search adjacent atoms:                                             |
 | Set the adjacent atoms for each atom and set the periodic boundary |
 | condition for the atoms on real space FFT grid. For k-dependent    |
 | algorithm, we also need to set the sparse H and S matrix element   |
 | for each atom.                                                     |
 |                                                                    |
 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<





 SETUP SEARCHING RADIUS FOR PROGRAM TO SEARCH ADJACENT ATOMS
                  longest orb rcut (Bohr) = 8.000
   longest nonlocal projector rcut (Bohr) = 2.160
              searching radius is (Bohr)) = 20.330
         searching radius unit is (Bohr)) = 1.890

 SETUP EXTENDED REAL SPACE GRID FOR GRID INTEGRATION
                          real space grid = [ 54, 54, 54 ]
                 big cell numbers in grid = [ 27, 27, 27 ]
             meshcell numbers in big cell = [ 2, 2, 2 ]
                        extended fft grid = [ 29, 29, 29 ]
                dimension of extened grid = [ 86, 86, 86 ]
                            UnitCellTotal = 125
              Atom number in sub-FFT-grid = 8
    Local orbitals number in sub-FFT-grid = 160
                                ParaV.nnr = 389332
                                     nnrg = 966984

 Warning_Memory_Consuming allocated:  pvpR_reduced 7.378 MB
                                nnrg_last = 810648
                                 nnrg_now = 966984

 Warning_Memory_Consuming allocated:  LOC::DM_R 7.378 MB
 DONE : INIT SCF Time : 140.328733 (SEC)

TOTAL-PRESSURE: 173.673386 KBAR


 Largest gradient in stress is 173.673386

 Threshold is = 2.000000

 Relaxation is not converged yet!
DIRECT COORDINATES
   atom           x                y                z           mag          vx               vy               vz       
taud_Fe1       0.0000000000     0.0000000000     0.0000000000 +0.0000     0.0000000000     0.0000000000     0.0000000000
taud_Fe2       0.0000000000     0.5000000000     0.5000000000 +0.0000     0.0000000000     0.0000000000     0.0000000000
taud_Fe3       0.5000000000     0.0000000000     0.5000000000 +0.0000     0.0000000000     0.0000000000     0.0000000000
taud_Fe4       0.5000000000     0.5000000000     0.0000000000 +0.0000     0.0000000000     0.0000000000     0.0000000000
taud_O1        0.5000000000     0.0000000000     0.0000000000 +0.0000     0.0000000000     0.0000000000     0.0000000000
taud_O2        0.5000000000     0.5000000000     0.5000000000 +0.0000     0.0000000000     0.0000000000     0.0000000000
taud_O3        0.0000000000     0.0000000000     0.5000000000 +0.0000     0.0000000000     0.0000000000     0.0000000000
taud_O4        0.0000000000     0.5000000000     0.0000000000 +0.0000     0.0000000000     0.0000000000     0.0000000000


                          Volume (Bohr^3) = 471.227218
                            Volume (A^3)) = 69.828586

 Lattice vectors: (Cartesian coordinate: in unit of a_0)
            +4.117920           +0.000000           +0.000000
            +0.000000           +4.117920           -0.000000
            +0.000000           -0.000000           +4.117920
 Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)
            +0.242841           -0.000000           -0.000000
            -0.000000           +0.242841           +0.000000
            -0.000000           +0.000000           +0.242841
 DONE : SETUP UNITCELL Time : 450.109782 (SEC)

以上为输出结果中第2步迭代的运行日志的部分内容。

其中TOTAL-PRESSURE参数为晶胞内应力，对应INPUT文件中的stress_thr参数，LARGEST GRAD (eV/A) 参数为最大原子受力，对应 INPUT 文件中的 force_thr_ev 参数，当 TOTAL-PRESSURE 小于 stress_thr 参数中设置的阈值时，且 LARGEST GRAD (eV/A) 小于 force_thr_ev 参数中设置的阈值时，认为结构优化收敛，优化结束。

对Fe单质和FeO进行如上计算，可以得到优化后的cif文件，在Bohrium中预览如下图。 alt

alt

得到FeO具有类似NaCl的晶格结构，晶格常数为4.117埃，现实中FeO由于非理想配比的存在取值在0.4301nm到0.42816nm之间；得到Fe单质形成了体心立方结构，晶格常数为2.913埃，现实中，铁在912摄氏度下以体心立方结构的存在，其结构常数为0.286nm。

代码

文本

4.2 平面波基组下的收敛性测试

4.2.1 Ecut收敛度测试

平面波作为一种可以用于描述周期性边界条件下的电子波函数和电荷密度基矢量，它们都是正交的，且可以通过一个ecutwfc参数来控制基矢量的个数。ecut代表了每一个平面波所对应的动能，ecut取得越大，基矢量就可以描述震荡得越剧烈的物理量，那么计算结果就会越精确，但同时所带来的机时成本消耗也越大。因此我们采用pw基组做真正计算前，需要对 ecut 进行测试来获得一个足够准确且效率也高的取值。

设置前文用于计算的基组为"pw"，进行ecut的收敛性测试，ecut取值范围为：20~100 Ry。

计算完成后，提取不同ecut数值下OUT.Si文件夹里running_scf.log文件中收敛的系统总能量，计算单原子能量并绘制其随着ecut的变化曲线。以Fe为例如下图所示。 alt

4.2.2 k点收敛性测试

对于平面波DFT计算而言，所需要进行的大量工作可以归纳为求算下述形式的积分值，即 $\overset{g}{ˉ} = \frac{V _{ce ll}}{8 π ^{3}} \int_{BZ} g (k) d k$ 该积分的主要特点是：它是在倒易空间中所定义的，并对布里渊区中所有可能的k值进行了积分。

对于这种积分，通过将被积函数在一系列离散点上进行估值，并在每一点上使用合适权童，再对这些函数值加和，可以近似求算该积分。而且随着加和所使用离散点个数的增加，这种数值方法可以给出更精确的结果。当所使用的点非常多，并对其取极限时，这些数值方法就收敛于积分的精确解。

在实际运算过程中，因为上述积分所消耗的DFT计算工作量非常大，所以对如何尽量快速并准确地估算这些积分变得十分关键。在这里我们也可以简单地对其进行测试和估计。

类似上一小节，改变KPT文件中各个方向的k值进行计算，计算完成后，提取各K点下OUT.Fe文件夹里running_scf.log文件中收敛的系统总能量，绘制其随K点数量的变化曲线。 alt

可以通过这样的可视化图像判断如何适当选取ecut与k值来获得尽量精确高效的计算结果。

代码

文本

notebook

中文

作业

ABACUS使用教程

notebook中文作业ABACUS使用教程

已赞1