新建
快速开始ABACUS | 计算硫（S）元素单质及其化合物的晶体结构
donglikun@dp.tech
推荐镜像 :ABACUS:3.6.0-dpdata-dpgen
推荐机型 :c2_m4_cpu
赞
目录
数据集
Abacus_DFT_S(v2)
致谢：
原notebook链接：https://bohrium.dp.tech/notebooks/16846533825
1. 元素背景知识
硫（S）是一种非金属元素，化学符号为S，原子序数为16，其电子排布式为。
硫在自然界中常以硫化物矿和硫酸盐矿的形式存在，早在3000年前就被古埃及文明发现并用于制造染料和药物。在隋朝，诞生了硝石、硫磺和木炭三元体系火药。1746年，英国化学家 J.Roebuck 发明了铅室法制造硫酸。1776年，法国化学家拉瓦锡首先确定了硫的不可分割性，认为硫是一种元素。
物理化学性质：硫是一种黄色固体，在常温下呈现为黄色结晶，具有特殊的气味。单质硫不溶于水，微溶于乙醇，易溶于二硫化碳CS2。硫有多种同素异形体，其中正交硫是室温下唯一稳定的硫的存在形式；而当硫被加热到100°C左右时会转变为。在硫晶体中，硫原子形成化学式为S8的环状八原子分子。自然界中通常以单质硫、硫化物、硫酸盐的形式存在，重要的硫化物如闪锌矿（ZnS），常见的价态包括-2、0、+4、+6等。
常用化合物：硫形成了许多重要的化合物。其中最常见的是二氧化硫（SO2），一种重要的工业气体，也是大气中的污染物之一。硫酸（H2SO4）是一种强酸，被广泛用于工业生产和化学实验中。硫化氢（H2S）是一种有毒气体，也是一种重要的化学原料。
用途：硫化物可以用于金属冶炼，硫酸是许多工业过程的重要原料，硫化氢在化学合成和环境监测中也有重要用途。此外，硫还被用于制造橡胶、烟火和染料，还可以用于治疗皮肤病。
代码
文本
2. DFT算法基本原理
密度泛函理论(Density functional theory ，DFT)是一种研究多电子体系电子结构的方法。密度泛函理论在物理和化学上都有广泛的应用，特别是用来研究分子和凝聚态的性质，是凝聚态物理计算材料学和计算化学领域最常用的方法之一。
根据Hohenberg-Kohn第一定理，对于任意相互作用的电子系统，其受到的外势与基态的电子密度一一对应。在材料计算中，这个外势主要由正电荷背景提供，与原子的不同排布方式有关。而Hohenberg-Kohn第二定理则证明了以基态密度为变量，将体系能量最小化之后就得到了基态能量。
密度泛函理论最普遍的应用是通过Kohn-Sham方法实现的。 在Kohn-Sham DFT的框架中，最难处理的多体问题（由于处在一个外部静电势中的电子相互作用而产生的）被简化成了一个没有相互作用的电子在有效势场中运动的问题。这个有效势场包括了电子动能、电子-离子吸引势（外势）、电子-电子排斥势、交换关联势能以及离子-离子排斥势。其中最难处理的是交换关联能，一般是用一些近似手法（如局域密度近似LDA）获得近似的势能分布，即赝势。获得了赝势分布就可以通过ABACUS求解系统的基态能量。
有关更多的 DFT 理论背景介绍，推荐阅读：了解第一性原理计算与密度泛函理论
在接下来的算例中，我们将利用ABACUS的平面波pw模块对硫元素的BCC结构以及闪锌矿（ZnS）的FCC结构进行计算。
代码
文本
3. ABACUS的平面波计算与收敛性测试
3.1 输入文件
在正式开始 ABACUS SCF 计算之前，请确保这些文件已经准备好并存储在工作目录中：
INPUT：包含了计算过程中所需的各种参数，定义和控制计算任务；STRU：结构文件，包含了原子种类、原子位置、晶格常数以及晶格向量等信息；KPT：包含了布里渊区积分所需的k点信息；*.upf：包含了原子的赝势信息；*.orb：包含了原子轨道的数值表示（pw计算用不到）；
在示例数据集中，我们为你准备好了格式转换的示例文件。我们可以直接访问：
代码
文本
[1]
! tree /bohr/
/bohr/ └── ABACUS-DFT-96ro └── v2 └── ABACUS_DFT ├── PP_ORB │ ├── S_ONCV_PBE-1.0.upf │ └── Zn_ONCV_PBE-1.0.upf ├── S_CELL-RELAX │ ├── INPUT │ ├── KPT │ └── STRU ├── S_SCF │ ├── INPUT │ ├── KPT │ └── STRU └── ZnS_CELL-RELAX ├── INPUT ├── KPT └── STRU 7 directories, 11 files
代码
文本
出于安全考虑，我们没有数据集所在文件夹的写入权限，因此我们将其复制到 /data/ 目录下:
代码
文本
[2]
! cp -nr /bohr/ /data/
代码
文本
进入S_SCF文件夹准备SCF计算：
代码
文本
[3]
import os
bohr_dataset_url = "/ABACUS-DFT-96ro/v2/ABACUS_DFT/S_SCF/" # url 可从左侧数据集复制
work_path = os.path.join("/data", bohr_dataset_url[1:])
os.chdir(work_path)
print(f"当前路径为：{os.getcwd()}")
当前路径为：/data/ABACUS-DFT-96ro/v2/ABACUS_DFT/S_SCF
代码
文本
3.1.1 STRU 文件
文件示例
STRU 文件包含晶格几何信息，赝势和数值轨道文件的名称和/或位置，以及关于系统的结构信息。我们可以通过 cat 命令进行查看：
代码
文本
[4]
! cat STRU
ATOMIC_SPECIES S 32.059 S_ONCV_PBE-1.0.upf NUMERICAL_ORBITAL S_gga_10au_100Ry_2s2p1d.orb LATTICE_CONSTANT 1.8897261258369282 LATTICE_VECTORS 3.16 0.00 0.00 0.00 3.16 0.00 0.00 0.00 3.16 ATOMIC_POSITIONS Direct S 0.0000000000 2 0.0000000000 0.0000000000 0.0000000000 1 1 1 mag 0.0 0.5000000000 0.5000000000 0.5000000000 1 1 1 mag 0.0
代码
文本
下面显示了 SCF 计算中的 BCC S 的 STRU 文件：
#This is the atom file containing all the information
#about the lattice structure.
ATOMIC_SPECIES
S 32.059 S_ONCV_PBE-1.0.upf # element name, atomic mass, pseudopotential file
NUMERICAL_ORBITAL
S_gga_10au_100Ry_2s2p1d.orb
LATTICE_CONSTANT
1.8897261258369282 # 1.8897259886 Bohr = 1.0 Angstrom
LATTICE_VECTORS
3.16 0.00 0.00
0.00 3.16 0.00
0.00 0.00 3.16
ATOMIC_POSITIONS
Direct #Cartesian(Unit is LATTICE_CONSTANT)
S #Name of element
0.0 #Magnetic for this element.
2 #Number of atoms
0.0 0.0 0.0 1 1 1 mag 0.0 #x,y,z, move_x, move_y, move_z
0.5 0.5 0.5 1 1 1 mag 0.0 #x,y,z, move_x, move_y, move_z
参数解释
- ATOMIC_SPECIES
从左到右依次为原子种类、相对原子质量、赝势文件
- NUMERICAL_ORBITAL
原子轨道文件
- 赝势文件
S_ONCV_PBE-1.0.upf应放在 INPUT 文件的pseudo_dir目录下（这里是..\PP_ORB目录）。
赝势文件可以从 ABACUS-Pseudopot-Nao-Square 网站 下载。
- LATTICE_CONSTANT
晶格常量，晶胞长度的缩放系数。
- LATTICE_VECTORS
晶格向量，即晶胞的 x,y,z 轴长度（实际长度需乘以晶格常量）。
- ATOMIC_POSITIONS
原子位置信息。BCC 的每个晶胞内都有 1 + 8*1/8 = 2个原子，因此只需要输入两个位置信息，剩下的位置可以通过对称性自动算出。
有关更多 STRU 文件的参数信息，请参阅：The STRU file
代码
文本
3.1.2 INPUT 文件
接下来，需要 INPUT 文件，该文件设置了所有关键参数，以指导 ABACUS 如何进行计算以及输出什么内容：
示例文件
INPUT_PARAMETERS
suffix S
calculation scf
symmetry 1
pseudo_dir ../PP_ORB
orbital_dir ../PP_ORB
basis_type pw
ecutwfc 80
scf_thr 1e-8
nbands 22
参数解释
- suffix：系统名称，默认为ABACUS
- calculation：ABACUS要执行的计算类型
- symmentry：可以利用对称性减少计算量
- pseudo_dir：提供赝势文件的目录
- orbital_dir：提供轨道文件的目录（pw用不到）
- basis_type：用于展开电子波函数的基函数类型
- ecutwfc：波函数展开的平面波能量截止（单位：Rydberg）
- scf_thr：电荷密度收敛阈值（单位：Rydberg）
- nbands：计算的Kohn-Sham电子轨道数，ABACUS目前采用
nbands = max(1.2*occupied_bands, occupied_bands + 10)。occupied_bands可以从赝势文件中的z_valence计算得到，这里z_valence = 6，每个晶胞中有2个原子，所以nbands = max(1.2*6*2, 6*2+10) = 22。如果不知道是多少，可以设置为0，让程序自己算。
参数列表始终以关键字 INPUT_PARAMETERS 开头。在 INPUT_PARAMETERS 之前的任何内容都将被忽略。以 # 或 / 开头的任何行也将被忽略。
注意：INPUT 参数设置的第一原则是：简洁明了。 参数不是越多越好，对于大多数场景，默认值已不失为一个很好的选择。请确保你设置的参数你都明确清楚其意义。
再注意：用户无法将文件名“INPUT”更改为其他名称。 布尔参数（如 out_chg）可以通过使用 True 和 False、1 和 0 或 T 和 F 来设置。大小写不敏感，因此也支持使用 true 和 false、TRUE 和 FALSE 以及 t 和 f 设置布尔值的其他偏好。
有关 INPUT 文件的所有可能参数信息，请参阅：Full List of INPUT Keywords
代码
文本
3.1.3 KPT 文件
最后一个必需的输入文件称为 KPT，此文件包含关于布里渊区采样的k点网格设置的信息。以下是一个示例：
K_POINTS
0 # k点的总数，`0'表示自动生成
Gamma # Monkhorst-Pack方法的类型，`Gamma'或`MP'
4 4 4 0 0 0 # 前三个数字：沿着倒数向量的细分
# 后三个数字：网格的偏移
代码
文本
3.2 运行计算
确认进入S_SCF文件夹后，执行以下命令可以快速运行ABACUS计算：
代码
文本
[5]
! mpirun -n 2 abacus
ABACUS v3.6.0
Atomic-orbital Based Ab-initio Computation at UStc
Website: http://abacus.ustc.edu.cn/
Documentation: https://abacus.deepmodeling.com/
Repository: https://github.com/abacusmodeling/abacus-develop
https://github.com/deepmodeling/abacus-develop
Commit: 5b263cbff (Wed Mar 27 18:42:03 2024 +0800)
Sat May 18 23:11:13 2024
MAKE THE DIR : OUT.S/
RUNNING WITH DEVICE : CPU / Intel(R) Xeon(R) Platinum
UNIFORM GRID DIM : 36 * 36 * 36
UNIFORM GRID DIM(BIG) : 36 * 36 * 36
DONE(0.0228798 SEC) : SETUP UNITCELL
DONE(0.181749 SEC) : SYMMETRY
DONE(0.419332 SEC) : INIT K-POINTS
---------------------------------------------------------
Self-consistent calculations for electrons
---------------------------------------------------------
SPIN KPOINTS PROCESSORS
1 10 2
---------------------------------------------------------
Use plane wave basis
---------------------------------------------------------
ELEMENT NATOM XC
S 2
---------------------------------------------------------
Initial plane wave basis and FFT box
---------------------------------------------------------
DONE(0.42266 SEC) : INIT PLANEWAVE
MEMORY FOR PSI (MB) : 4.33716
DONE(0.429407 SEC) : LOCAL POTENTIAL
DONE(0.469441 SEC) : NON-LOCAL POTENTIAL
DONE(0.469507 SEC) : INIT BASIS
-------------------------------------------
SELF-CONSISTENT :
-------------------------------------------
START CHARGE : atomic
DONE(0.516526 SEC) : INIT SCF
ITER ETOT(eV) EDIFF(eV) DRHO TIME(s)
CG1 -5.553750e+02 0.000000e+00 2.075e-03 3.907e+00
CG2 -5.553791e+02 -4.187408e-03 1.101e-04 6.784e-01
CG3 -5.553796e+02 -4.279199e-04 7.833e-06 7.832e-01
CG4 -5.553797e+02 -1.365589e-04 2.896e-06 8.772e-01
CG5 -5.553797e+02 2.968616e-06 9.486e-07 6.820e-01
CG6 -5.553798e+02 -6.684664e-05 2.459e-08 9.196e-01
CG7 -5.553798e+02 -1.874727e-07 4.404e-10 8.909e-01
TIME STATISTICS
-------------------------------------------------------------------------------------
CLASS_NAME NAME TIME(Sec) CALLS AVG(Sec) PER(%)
-------------------------------------------------------------------------------------
total 9.26 17 0.54 100.00
Driver reading 0.01 1 0.01 0.09
Input Init 0.00 1 0.00 0.05
Input_Conv Convert 0.00 1 0.00 0.01
Driver driver_line 9.25 1 9.25 99.91
UnitCell check_tau 0.00 1 0.00 0.00
PW_Basis_Sup setuptransform 0.00 1 0.00 0.01
PW_Basis_Sup distributeg 0.00 1 0.00 0.01
mymath heapsort 0.00 151 0.00 0.02
Symmetry analy_sys 0.00 1 0.00 0.00
PW_Basis_K setuptransform 0.00 1 0.00 0.02
PW_Basis_K distributeg 0.00 1 0.00 0.01
PW_Basis setup_struc_factor 0.00 1 0.00 0.01
ppcell_vnl init 0.00 1 0.00 0.01
ppcell_vl init_vloc 0.00 1 0.00 0.03
ppcell_vnl init_vnl 0.04 1 0.04 0.43
WF_atomic init_at_1 0.00 1 0.00 0.00
wavefunc wfcinit 0.00 1 0.00 0.00
Ions opt_ions 8.79 1 8.79 94.91
ESolver_KS_PW Run 8.79 1 8.79 94.91
H_Ewald_pw compute_ewald 0.00 1 0.00 0.02
Charge set_rho_core 0.00 1 0.00 0.00
Charge atomic_rho 0.01 1 0.01 0.10
PW_Basis_Sup recip2real 0.04 57 0.00 0.44
PW_Basis_Sup gathers_scatterp 0.02 57 0.00 0.20
Potential init_pot 0.03 1 0.03 0.32
Potential update_from_charge 0.17 8 0.02 1.85
Potential cal_fixed_v 0.00 1 0.00 0.01
PotLocal cal_fixed_v 0.00 1 0.00 0.01
Potential cal_v_eff 0.17 8 0.02 1.84
H_Hartree_pw v_hartree 0.01 8 0.00 0.15
PW_Basis_Sup real2recip 0.04 73 0.00 0.42
PW_Basis_Sup gatherp_scatters 0.01 73 0.00 0.16
PotXC cal_v_eff 0.15 8 0.02 1.67
XC_Functional v_xc 0.15 8 0.02 1.66
Potential interpolate_vrs 0.00 8 0.00 0.01
Symmetry rhog_symmetry 0.04 9 0.00 0.40
Symmetry group fft grids 0.01 9 0.00 0.09
Charge_Mixing init_mixing 0.00 1 0.00 0.00
ESolver_KS_PW hamilt2density 8.57 8 1.07 92.57
HSolverPW solve 8.52 8 1.06 92.00
Nonlocal getvnl 0.07 80 0.00 0.73
pp_cell_vnl getvnl 0.07 80 0.00 0.73
Structure_Factor get_sk 0.01 80 0.00 0.11
DiagoIterAssist diagH_subspace 1.18 70 0.02 12.72
Operator hPsi 6.49 7667 0.00 70.07
Operator EkineticPW 0.13 7667 0.00 1.37
Operator VeffPW 5.75 7667 0.00 62.08
PW_Basis_K recip2real 3.27 10897 0.00 35.29
PW_Basis_K gathers_scatterp 1.28 10897 0.00 13.78
PW_Basis_K real2recip 2.45 9137 0.00 26.41
PW_Basis_K gatherp_scatters 0.71 9137 0.00 7.68
Operator NonlocalPW 0.60 7667 0.00 6.43
Nonlocal add_nonlocal_pp 0.23 7667 0.00 2.45
DiagoIterAssist diagH_LAPACK 0.02 70 0.00 0.24
DiagoCG diag_once 6.53 80 0.08 70.50
DiagoCG_New spsi_func 0.05 15194 0.00 0.58
DiagoCG_New hpsi_func 5.48 7597 0.00 59.15
ElecStatePW psiToRho 0.58 8 0.07 6.26
Charge_Mixing get_drho 0.01 8 0.00 0.11
Charge_Mixing inner_product_recip_rho 0.00 8 0.00 0.01
Charge mix_rho 0.01 6 0.00 0.13
Charge Broyden_mixing 0.00 6 0.00 0.04
Charge_Mixing inner_product_recip_hartree 0.00 30 0.00 0.03
ModuleIO write_istate_info 0.00 1 0.00 0.01
-------------------------------------------------------------------------------------
START Time : Sat May 18 23:11:13 2024
FINISH Time : Sat May 18 23:11:23 2024
TOTAL Time : 10
SEE INFORMATION IN : OUT.S/
代码
文本
输出文件存储在当前文件夹下的OUT.S文件夹下，文件夹名字可以从INPUT文件的suffix值修改。
计算结果主要储存在OUT.S/running_scf.log中，我们这里关心的只有系统基态能量的部分，可以利用grep操作提取相关信息：
代码
文本
[6]
! grep ETOT OUT.S/running_scf.log
!FINAL_ETOT_IS -555.3797772318255 eV
代码
文本
可以看到系统基态能量为-555.3797772318255 eV，平均每个原子的能量为-277.6899 eV/atom（精确到0.1 meV/atom）。
代码
文本
3.3 收敛性测试
3.3.1 Ecut 收敛
平面波（pw）作为一种可以用于描述周期性边界条件下的电子波函数和电荷密度基矢量，它们都是正交的，且可以通过一个 ecutwfc 参数来控制基矢量的个数，ecut 其实代表了每一个平面波所对应的动能，如果 ecut 取得越大，则基矢量可以描述震荡得越剧烈的物理量（例如氧原子的 2p 轨道），那么计算结果就会越精确，但同时所带来的机时成本消耗也越大。因此我们采用 pw 基组做真正计算前，需要对 ecut 进行测试来获得一个足够准确且效率也高的取值。
保持前文用于计算的基组为"pw"，进行 ecut 的收敛性测试，修改INPUT文件中的 ecut 值，范围为：40~100 Ry，重复计算每个原子的基态能量并作Energy-ecut图如下：
可以看见随 ecut 增大，基态能量趋于收敛。在 ecut=80 Ry 时，认为基态能量收敛(与 ecut=70 Ry 的能量差小于 1 meV/atom)。
代码
文本
3.3.2 K点收敛
对具有周期性的体系进行计算时，DFT 计算实际会对第一布里渊区中不同离散化k点的单独进行计算（例如每个k点都会单独求解一次依赖于该k点的 Kohn-Sham 方程），并将计算结果进行积分。k点越多，则需要求的离散积分就越多，也将越近似连续积分的结果，但计算资源也会增加。因此，pw 和 LCAO 都有必要对 "需要多少k点" 进行测试并得出合适的k点选取方案。
修改KPT文件的k值，范围从2~12，此时除了基态能量ETOT外，我们还关注计算需要的总时间Total Time，可以用grep操作提取信息：
代码
文本
[7]
! grep secs OUT.S/running_scf.log
Total Time : 0 h 0 mins 10 secs
代码
文本
即测试算例计算所需总时间为10 s。
重复计算得到基态能量关于k值的图Energy-Nk，与计算总时间关于k值的图Time-Nk如下：
可见随着Nk的增大，能量在振荡着收敛，而计算时长在逐步爬升。最终能量在Nk=12时收敛（与Nk=11时的能量差在 1 meV/atom量级），此时的计算时间在100s左右。
因此选择12 x 12 x 12的k点可以同时获得较准确的计算结果和较高的计算效率，此时计算出的能量为Energy = -277.668 eV/atom。
另外，可以看到时间随 Nk 的变化呈现阶梯式上升的趋势，相邻的奇偶数 Nk 计算所需的时间相近。这是因为在INPUT文件中我们设置了参数symmetry = 1，所以在计算K点个数时采用了一些对称性简化，得到更少的k点；而这种简化就是阶梯性的，如下表：
代码
文本
4. 结构优化（弛豫）
结构优化是从一个一般的结构出发，利用优化算法得到基态结构的过程。对于周期性体系，优化计算可以根据需要针对体系的不同自由度来进行优化。
在ABACUS中，如果只优化原子位置，而晶胞固定不变，则需在INPUT中将calculation设置为relax;
如果需要同时做原子位置优化和晶胞优化，则需要将calculation设置为cell-relax；
此外，需要注意的参数还有原子受力的收敛标准force_thr_ev和应力收敛标准stress_thr的设置，具体会在下面展开说明；
更多的优化设置，请阅读ABACUS Documentation。关于结构优化，可参考文档中对这部分的详细说明。
4.1 输入文件
主要的输入文件与SCF计算差不多，主要是INPUT文件有所更改。先进入S_CELL-RELAX文件夹：
代码
文本
[8]
bohr_dataset_url = "/ABACUS-DFT-96ro/v2/ABACUS_DFT/S_CELL-RELAX/" # url 可从左侧数据集复制
work_path = os.path.join("/data", bohr_dataset_url[1:])
os.chdir(work_path)
print(f"当前路径为：{os.getcwd()}")
当前路径为：/data/ABACUS-DFT-96ro/v2/ABACUS_DFT/S_CELL-RELAX
代码
文本
4.1.1 INPUT 文件
利用cat操作查看INPUT文件：
代码
文本
[9]
!cat INPUT
INPUT_PARAMETERS suffix S calculation cell-relax symmetry 1 pseudo_dir ../PP_ORB orbital_dir ../PP_ORB basis_type pw ecutwfc 80 scf_thr 1e-8 nbands 22 force_thr_ev 0.01 stress_thr 0.5 out_stru 1
代码
文本
可以看到cell-relax与scf的区别主要在以下几处：
calculation：这里将计算类型设置为cell-relax（注意是中划线），即结构弛豫；force_thr_ev：力收敛阈值，用于判断优化过程是否收敛。这里设置为0.01 eV/Angstrom，表示当所有原子上的力小于0.01 eV/Angstrom时，认为已经收敛。默认值：0.0257112 eV/Angstrom(你也可以使用force_thr参数（默认值：0.001 Ry/Bohr，两者只是单位不同。)；stress_thr: 应力收敛阈值，用于判断晶胞优化过程中应力是否收敛。这里设置为2 kBar，表示当晶胞内的应力小于2 kBar时，认为已经收敛。默认值：0.5 kBar；out_stru: 表示是否输出优化后的结构信息。这里设置为1，表示在优化结束后将优化后的结构信息输出到文件。
代码
文本
4.1.2 STRU 文件
通过cat操作查看文件：
代码
文本
[10]
!cat STRU
ATOMIC_SPECIES S 32.059 S_ONCV_PBE-1.0.upf NUMERICAL_ORBITAL S_gga_10au_100Ry_2s2p1d.orb LATTICE_CONSTANT 1.8897261258369282 LATTICE_VECTORS 3.16 0.00 0.00 0.00 3.16 0.00 0.00 0.00 3.16 ATOMIC_POSITIONS Direct S 0.0000000000 2 0.0000000000 0.0000000000 0.0000000000 1 1 1 mag 0.0 0.5000000000 0.5000000000 0.5000000000 1 1 1 mag 0.0
代码
文本
可以发现与scf计算的文件没有区别。但是这里需要注意几点：
LATTICE_VECTOR中的值会随着优化的进行而改变，因此这里的值均为预设值。预设值与最终值的差别越大，优化所需要的时间越长；（预设值的选取是一门艺术）ATOMIC_POSITIONS中的值也会改变。如果不想让它的某个值改变，就可以把对应的move值从1改成0。举个例子，如果我想固定0 0 0处的原子的x y坐标不变，只变化z坐标的值，那就把后面的1 1 1改成0 0 1；（这对于单分子体系很有用）
代码
文本
4.1.3 KPT 文件
与scf一样，没有更改，可以通过cat操作查看：
代码
文本
[11]
!cat KPT
K_POINTS 0 Gamma 4 4 4 0 0 0
代码
文本
4.2 运行计算
设置好输入文件后直接运行ABACUS程序：
代码
文本
[12]
!mpirun -n 2 abacus
ABACUS v3.6.0
Atomic-orbital Based Ab-initio Computation at UStc
Website: http://abacus.ustc.edu.cn/
Documentation: https://abacus.deepmodeling.com/
Repository: https://github.com/abacusmodeling/abacus-develop
https://github.com/deepmodeling/abacus-develop
Commit: 5b263cbff (Wed Mar 27 18:42:03 2024 +0800)
Sat May 18 23:11:41 2024
MAKE THE DIR : OUT.S/
RUNNING WITH DEVICE : CPU / Intel(R) Xeon(R) Platinum
UNIFORM GRID DIM : 36 * 36 * 36
UNIFORM GRID DIM(BIG) : 36 * 36 * 36
DONE(0.0235207 SEC) : SETUP UNITCELL
DONE(0.186284 SEC) : SYMMETRY
DONE(0.41169 SEC) : INIT K-POINTS
---------------------------------------------------------
Cell relaxation calculations
---------------------------------------------------------
SPIN KPOINTS PROCESSORS
1 10 2
---------------------------------------------------------
Use plane wave basis
---------------------------------------------------------
ELEMENT NATOM XC
S 2
---------------------------------------------------------
Initial plane wave basis and FFT box
---------------------------------------------------------
DONE(0.415025 SEC) : INIT PLANEWAVE
MEMORY FOR PSI (MB) : 4.33716
DONE(0.422005 SEC) : LOCAL POTENTIAL
DONE(0.484586 SEC) : NON-LOCAL POTENTIAL
DONE(0.484645 SEC) : INIT BASIS
-------------------------------------------
STEP OF RELAXATION : 1
-------------------------------------------
START CHARGE : atomic
DONE(0.523503 SEC) : INIT SCF
ITER ETOT(eV) EDIFF(eV) DRHO TIME(s)
CG1 -5.553750e+02 0.000000e+00 2.075e-03 3.895e+00
CG2 -5.553791e+02 -4.187408e-03 1.101e-04 6.789e-01
CG3 -5.553796e+02 -4.279199e-04 7.833e-06 8.193e-01
CG4 -5.553797e+02 -1.365589e-04 2.896e-06 8.613e-01
CG5 -5.553797e+02 2.968616e-06 9.486e-07 6.911e-01
CG6 -5.553798e+02 -6.684664e-05 2.459e-08 9.102e-01
CG7 -5.553798e+02 -1.874727e-07 4.404e-10 8.910e-01
----------------------------------------------------------------
TOTAL-STRESS (KBAR)
----------------------------------------------------------------
1.1469814052 0.0000000000 -0.0000000000
0.0000000000 1.1469814052 -0.0000000000
-0.0000000000 -0.0000000000 1.1469814052
----------------------------------------------------------------
TOTAL-PRESSURE: 1.146981 KBAR
ETOT DIFF (eV) : 0.000000
LARGEST GRAD (eV/A) : 0.000000
DONE(9.507173 SEC) : SETUP UNITCELL
-------------------------------------------
STEP OF RELAXATION : 2
-------------------------------------------
DONE(9.513637 SEC) : LOCAL POTENTIAL
DONE(9.669451 SEC) : SYMMETRY
DONE(9.670139 SEC) : INIT K-POINTS
DONE(9.733629 SEC) : NON-LOCAL POTENTIAL
DONE(9.806991 SEC) : INIT SCF
ITER ETOT(eV) EDIFF(eV) DRHO TIME(s)
CG1 -5.553798e+02 0.000000e+00 5.812e-08 1.459e+00
CG2 -5.553798e+02 -8.784881e-08 8.766e-10 6.776e-01
----------------------------------------------------------------
TOTAL-STRESS (KBAR)
----------------------------------------------------------------
-0.3990839414 0.0000000000 0.0000000000
0.0000000000 -0.3990839414 -0.0000000000
0.0000000000 -0.0000000000 -0.3990839414
----------------------------------------------------------------
TOTAL-PRESSURE: -0.399084 KBAR
ETOT DIFF (eV) : -0.000013
LARGEST GRAD (eV/A) : 0.000000
TIME STATISTICS
-------------------------------------------------------------------------------------
CLASS_NAME NAME TIME(Sec) CALLS AVG(Sec) PER(%)
-------------------------------------------------------------------------------------
total 12.18 29 0.42 100.00
Driver reading 0.01 1 0.01 0.08
Input Init 0.00 1 0.00 0.04
Input_Conv Convert 0.00 1 0.00 0.00
Driver driver_line 12.17 1 12.17 99.92
UnitCell check_tau 0.00 1 0.00 0.00
PW_Basis_Sup setuptransform 0.00 1 0.00 0.01
PW_Basis_Sup distributeg 0.00 1 0.00 0.01
mymath heapsort 0.00 300 0.00 0.02
Symmetry analy_sys 9.49 1 9.49 77.90
PW_Basis_K setuptransform 0.00 1 0.00 0.01
PW_Basis_K distributeg 0.00 1 0.00 0.00
PW_Basis setup_struc_factor 0.00 2 0.00 0.02
ppcell_vnl init 0.00 1 0.00 0.01
ppcell_vl init_vloc 0.01 2 0.00 0.05
ppcell_vnl init_vnl 0.13 2 0.06 1.03
WF_atomic init_at_1 0.00 2 0.00 0.00
wavefunc wfcinit 0.00 1 0.00 0.00
Ions opt_ions 11.70 1 11.70 96.01
ESolver_KS_PW Run 11.32 2 5.66 92.94
H_Ewald_pw compute_ewald 0.00 2 0.00 0.01
Charge set_rho_core 0.00 2 0.00 0.00
Charge atomic_rho 0.02 4 0.01 0.19
PW_Basis_Sup recip2real 0.06 87 0.00 0.48
PW_Basis_Sup gathers_scatterp 0.02 87 0.00 0.20
Potential init_pot 0.05 2 0.02 0.40
Potential update_from_charge 0.23 11 0.02 1.91
Potential cal_fixed_v 0.00 2 0.00 0.01
PotLocal cal_fixed_v 0.00 2 0.00 0.01
Potential cal_v_eff 0.23 11 0.02 1.89
H_Hartree_pw v_hartree 0.02 11 0.00 0.20
PW_Basis_Sup real2recip 0.07 113 0.00 0.55
PW_Basis_Sup gatherp_scatters 0.02 113 0.00 0.21
PotXC cal_v_eff 0.20 11 0.02 1.68
XC_Functional v_xc 0.20 11 0.02 1.67
Potential interpolate_vrs 0.00 11 0.00 0.01
Symmetry rhog_symmetry 0.05 13 0.00 0.43
Symmetry group fft grids 0.01 13 0.00 0.10
Charge_Mixing init_mixing 0.00 2 0.00 0.00
ESolver_KS_PW hamilt2density 10.67 11 0.97 87.56
HSolverPW solve 10.59 11 0.96 86.97
Nonlocal getvnl 0.09 110 0.00 0.77
pp_cell_vnl getvnl 0.13 150 0.00 1.05
Structure_Factor get_sk 0.03 310 0.00 0.23
DiagoIterAssist diagH_subspace 1.33 80 0.02 10.93
Operator hPsi 8.00 9601 0.00 65.66
Operator EkineticPW 0.17 9601 0.00 1.36
Operator VeffPW 7.07 9601 0.00 58.01
PW_Basis_K recip2real 4.12 13701 0.00 33.78
PW_Basis_K gathers_scatterp 1.61 13701 0.00 13.20
PW_Basis_K real2recip 3.01 11281 0.00 24.67
PW_Basis_K gatherp_scatters 0.86 11281 0.00 7.05
Operator NonlocalPW 0.74 9601 0.00 6.10
Nonlocal add_nonlocal_pp 0.29 9601 0.00 2.37
DiagoIterAssist diagH_LAPACK 0.03 80 0.00 0.21
DiagoCG diag_once 8.17 110 0.07 67.09
DiagoCG_New spsi_func 0.07 19042 0.00 0.59
DiagoCG_New hpsi_func 6.86 9521 0.00 56.35
ElecStatePW psiToRho 0.82 11 0.07 6.71
Charge_Mixing get_drho 0.01 11 0.00 0.11
Charge_Mixing inner_product_recip_rho 0.00 11 0.00 0.01
Charge mix_rho 0.01 7 0.00 0.11
Charge Broyden_mixing 0.00 7 0.00 0.03
Charge_Mixing inner_product_recip_hartree 0.00 30 0.00 0.02
Forces cal_force_loc 0.00 2 0.00 0.03
Forces cal_force_ew 0.00 2 0.00 0.02
Forces cal_force_nl 0.05 2 0.02 0.41
Forces cal_force_cc 0.00 2 0.00 0.00
Forces cal_force_scc 0.01 2 0.00 0.08
Stress_PW cal_stress 0.31 2 0.15 2.51
Stress_Func stress_kin 0.01 2 0.00 0.07
Stress_Func stress_har 0.00 2 0.00 0.02
Stress_Func stress_ewa 0.00 2 0.00 0.03
Stress_Func stress_gga 0.02 2 0.01 0.19
Stress_Func stress_loc 0.02 2 0.01 0.19
Stress_Func stress_cc 0.00 2 0.00 0.00
Stress_Func stress_nl 0.24 2 0.12 2.00
ESolver_KS_PW init_after_vc 0.23 1 0.23 1.86
ModuleIO write_istate_info 0.00 1 0.00 0.01
-------------------------------------------------------------------------------------
START Time : Sat May 18 23:11:41 2024
FINISH Time : Sat May 18 23:11:54 2024
TOTAL Time : 13
SEE INFORMATION IN : OUT.S/
代码
文本
4.3 结构优化过程
从上面的运行输出日志中把优化步骤截取下来：
---------------------------------------------------------
DONE(0.404773 SEC) : INIT PLANEWAVE
MEMORY FOR PSI (MB) : 4.33716
DONE(0.411825 SEC) : LOCAL POTENTIAL
DONE(0.474569 SEC) : NON-LOCAL POTENTIAL
DONE(0.474766 SEC) : INIT BASIS
-------------------------------------------
STEP OF RELAXATION : 1
-------------------------------------------
START CHARGE : atomic
DONE(0.512434 SEC) : INIT SCF
ITER ETOT(eV) EDIFF(eV) DRHO TIME(s)
CG1 -5.553750e+02 0.000000e+00 2.075e-03 3.917e+00
CG2 -5.553791e+02 -4.187408e-03 1.101e-04 6.806e-01
CG3 -5.553796e+02 -4.279199e-04 7.833e-06 7.728e-01
CG4 -5.553797e+02 -1.365589e-04 2.896e-06 8.492e-01
CG5 -5.553797e+02 2.968616e-06 9.486e-07 6.705e-01
CG6 -5.553798e+02 -6.684664e-05 2.459e-08 9.107e-01
CG7 -5.553798e+02 -1.874727e-07 4.404e-10 8.876e-01
----------------------------------------------------------------
TOTAL-STRESS (KBAR)
----------------------------------------------------------------
1.1469814052 0.0000000000 -0.0000000000
0.0000000000 1.1469814052 -0.0000000000
-0.0000000000 -0.0000000000 1.1469814052
----------------------------------------------------------------
TOTAL-PRESSURE: 1.146981 KBAR
ETOT DIFF (eV) : 0.000000
LARGEST GRAD (eV/A) : 0.000000
DONE(9.455294 SEC) : SETUP UNITCELL
-------------------------------------------
STEP OF RELAXATION : 2
-------------------------------------------
DONE(9.463790 SEC) : LOCAL POTENTIAL
DONE(9.617383 SEC) : SYMMETRY
DONE(9.618123 SEC) : INIT K-POINTS
DONE(9.681088 SEC) : NON-LOCAL POTENTIAL
DONE(9.750329 SEC) : INIT SCF
ITER ETOT(eV) EDIFF(eV) DRHO TIME(s)
----------------------------------------------------------------
TOTAL-STRESS (KBAR)
----------------------------------------------------------------
-0.3990839414 0.0000000000 0.0000000000
0.0000000000 -0.3990839414 -0.0000000000
0.0000000000 -0.0000000000 -0.3990839414
----------------------------------------------------------------
TOTAL-PRESSURE: -0.399084 KBAR
ETOT DIFF (eV) : -0.000013
LARGEST GRAD (eV/A) : 0.000000
可以看到结构优化的主要评判标准是TOTAL - PRESSURE的大小，以及LARGEST GRAD的大小。前者是晶胞的内应力，对应INPUT文件的stress_thr参数；后者是最大原子受力，对应INPUT文件的force_thr_ev参数。当二者都小于设定的阈值时，认为结构优化收敛，优化结束。
在本例中，我们设定stress_thr = 0.5 kbar,force_thr_ev = 0.01 eV。在第1步优化中，TOTAL - PRESSURE = 1.147 kbar超过设定阈值，因此优化未结束；而在第2步优化中，TOTAL - PRESSURE = -0.399 kbar小于设定阈值，而且LARGEST GRAD也小于设定阈值，因此结构收敛，优化结束。
代码
文本
4.4 优化结果
优化后的结构可以在OUT.S/STRU_NOW.cif中找到：
代码
文本
[13]
!cat ./OUT.S/STRU_NOW.cif
data_none _audit_creation_method generated by ABACUS _cell_length_a 3.16178 _cell_length_b 3.16178 _cell_length_c 3.16178 _cell_angle_alpha 90 _cell_angle_beta 90 _cell_angle_gamma 90 loop_ _atom_site_label _atom_site_fract_x _atom_site_fract_y _atom_site_fract_z S 0 0 0 S 0.5 0.5 0.5
代码
文本
优化后的晶体结构长这样：
代码
文本
可以看到优化后的晶格结构仍然是BCC，但是晶格常数有所变化，由原来的3.16变成了3.16178。本例取的k值为4，取更高的k值可以算出更精确的数据，但是耗费的时间也就更多。
另外，如何验证数据的准确性是个问题。因为自然界中的硫都是以正交硫的形式存在的，但是那种结构过于复杂，我们的双核CPU计算不了。本例计算的BCC硫结构没有实验数据支撑。在这里笔者提供一个定性的验证方案：将算得的结构数据用scf计算一遍基态能量，如果算得的能量比原先的要低，可以定性认为此次结构优化的方向是正确的。
代码
文本
5. 多组分材料的结构优化
对于多组分材料的结构优化，如闪锌矿（ZnS），只需要在STRU文件中多加一个元素的赝势文件和结构位置即可。
由于闪锌矿的晶体结构比较复杂，可以在Materials Project上找到对应晶体结构的CIF或者POSCAR文件，然后通过ase接口转换成STRU文件。具体可以参考这篇Notebook：ASE-ABACUS接口使用方法介绍
在本例中，我们已经帮你转换好了STRU文件，直接用就行。
5.1 输入文件
进入ZnS_CELL-RELAX文件夹：
代码
文本
[14]
bohr_dataset_url = "/ABACUS-DFT-96ro/v2/ABACUS_DFT/ZnS_CELL-RELAX/" # url 可从左侧数据集复制
work_path = os.path.join("/data", bohr_dataset_url[1:])
os.chdir(work_path)
print(f"当前路径为：{os.getcwd()}")
当前路径为：/data/ABACUS-DFT-96ro/v2/ABACUS_DFT/ZnS_CELL-RELAX
代码
文本
可以通过cat操作查看三个输入文件INPUT,STRU和KPT：
代码
文本
[15]
!cat INPUT
INPUT_PARAMETERS suffix ZnS calculation cell-relax symmetry 1 pseudo_dir ../PP_ORB orbital_dir ../PP_ORB basis_type pw ecutwfc 80 scf_thr 1e-8 nbands 0 force_thr_ev 0.01 stress_thr 0.5 out_stru 1
代码
文本
[16]
!cat STRU
ATOMIC_SPECIES Zn 65.38 Zn_ONCV_PBE-1.0.upf S 32.06 S_ONCV_PBE-1.0.upf LATTICE_CONSTANT 1.889726 LATTICE_VECTORS 5.4432200000 0.0000000000 0.0000000000 0.0000000000 5.4432200000 0.0000000000 0.0000000000 0.0000000000 5.4432200000 ATOMIC_POSITIONS Direct Zn 0.0000000000 4 0.0000000000 0.0000000000 0.0000000000 1 1 1 mag 0.0 0.0000000000 0.5000000000 0.5000000000 1 1 1 mag 0.0 0.5000000000 0.0000000000 0.5000000000 1 1 1 mag 0.0 0.5000000000 0.5000000000 0.0000000000 1 1 1 mag 0.0 S 0.0000000000 4 0.2500000000 0.7500000000 0.7500000000 1 1 1 mag 0.0 0.2500000000 0.2500000000 0.2500000000 1 1 1 mag 0.0 0.7500000000 0.7500000000 0.2500000000 1 1 1 mag 0.0 0.7500000000 0.2500000000 0.7500000000 1 1 1 mag 0.0
代码
文本
[17]
!cat KPT
K_POINTS 0 Gamma 2 2 2 0 0 0
代码
文本
可以看到，与BCC硫输入文件的最主要区别就是加入了Zn元素的位置数据。相应的，我们还加入了Zn元素对应的赝势文件，放在../PP_ORB中。
代码
文本
5.2 运行计算
由于这个晶格结构比较复杂，如果预设值不准的话算得就会很慢。为了提高计算效率，我们可以先取一个低k值，粗略优化预设的晶格结构，然后逐步提升k值，利用上一次优化后的结构作为下一次的预设值。这样计算效率就会提高不少。
接下来我们先把k值从6改成2，可以用sed操作替换KPT文件中的内容：
代码
文本
[18]
!sed -i 's/6/2/g' KPT
!cat KPT
K_POINTS 0 Gamma 2 2 2 0 0 0
代码
文本
改完之后直接运行（预计用时250s，请耐心等待）：
代码
文本
[19]
!mpirun -n 2 abacus
ABACUS v3.6.0
Atomic-orbital Based Ab-initio Computation at UStc
Website: http://abacus.ustc.edu.cn/
Documentation: https://abacus.deepmodeling.com/
Repository: https://github.com/abacusmodeling/abacus-develop
https://github.com/deepmodeling/abacus-develop
Commit: 5b263cbff (Wed Mar 27 18:42:03 2024 +0800)
Sat May 18 23:12:14 2024
MAKE THE DIR : OUT.ZnS/
RUNNING WITH DEVICE : CPU / Intel(R) Xeon(R) Platinum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Warning: the number of valence electrons in pseudopotential > 12 for Zn: [Ar] 3d10 4s2
Pseudopotentials with additional electrons can yield (more) accurate outcomes, but may be less efficient.
If you're confident that your chosen pseudopotential is appropriate, you can safely ignore this warning.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
UNIFORM GRID DIM : 60 * 60 * 60
UNIFORM GRID DIM(BIG) : 60 * 60 * 60
DONE(0.0612782 SEC) : SETUP UNITCELL
DONE(0.226467 SEC) : SYMMETRY
DONE(0.44347 SEC) : INIT K-POINTS
---------------------------------------------------------
Cell relaxation calculations
---------------------------------------------------------
SPIN KPOINTS PROCESSORS
1 4 2
---------------------------------------------------------
Use plane wave basis
---------------------------------------------------------
ELEMENT NATOM XC
Zn 4
S 4
---------------------------------------------------------
Initial plane wave basis and FFT box
---------------------------------------------------------
DONE(0.451874 SEC) : INIT PLANEWAVE
MEMORY FOR PSI (MB) : 25.3477
DONE(0.520903 SEC) : LOCAL POTENTIAL
DONE(0.653123 SEC) : NON-LOCAL POTENTIAL
DONE(0.653214 SEC) : INIT BASIS
-------------------------------------------
STEP OF RELAXATION : 1
-------------------------------------------
START CHARGE : atomic
DONE(0.849534 SEC) : INIT SCF
ITER ETOT(eV) EDIFF(eV) DRHO TIME(s)
CG1 -2.290573e+04 0.000000e+00 4.854e-01 5.725e+01
CG2 -2.290809e+04 -2.352325e+00 2.051e-01 1.066e+01
CG3 -2.290864e+04 -5.551700e-01 1.986e-03 9.637e+00
CG4 -2.290867e+04 -2.480952e-02 8.151e-04 1.541e+01
CG5 -2.290867e+04 -3.724842e-03 4.763e-05 1.021e+01
CG6 -2.290867e+04 -3.263089e-04 2.658e-06 1.307e+01
CG7 -2.290867e+04 -3.042760e-05 1.044e-06 1.448e+01
CG8 -2.290867e+04 -2.574908e-06 2.698e-07 1.152e+01
CG9 -2.290867e+04 -1.356487e-06 9.202e-09 1.174e+01
----------------------------------------------------------------
TOTAL-STRESS (KBAR)
----------------------------------------------------------------
2.5142861550 0.0000000000 0.0000000000
0.0000000000 2.5142861550 0.0000000000
0.0000000000 0.0000000000 2.5142861550
----------------------------------------------------------------
TOTAL-PRESSURE: 2.514286 KBAR
ETOT DIFF (eV) : 0.000000
LARGEST GRAD (eV/A) : 0.000069
DONE(157.437811 SEC) : SETUP UNITCELL
-------------------------------------------
STEP OF RELAXATION : 2
-------------------------------------------
DONE(157.475393 SEC) : LOCAL POTENTIAL
DONE(157.555420 SEC) : SYMMETRY
DONE(157.555967 SEC) : INIT K-POINTS
DONE(157.690589 SEC) : NON-LOCAL POTENTIAL
DONE(158.052464 SEC) : INIT SCF
ITER ETOT(eV) EDIFF(eV) DRHO TIME(s)
CG1 -2.290821e+04 0.000000e+00 1.321e-06 1.873e+01
CG2 -2.290867e+04 -4.586716e-01 1.108e-06 1.426e+01
CG3 -2.290867e+04 -5.053294e-06 6.745e-07 1.151e+01
CG4 -2.290867e+04 -1.924784e-06 2.537e-09 9.423e+00
----------------------------------------------------------------
TOTAL-STRESS (KBAR)
----------------------------------------------------------------
-1.0248738162 -0.0000000000 0.0000000000
-0.0000000000 -1.0248738162 -0.0000000000
0.0000000000 -0.0000000000 -1.0248738162
----------------------------------------------------------------
TOTAL-PRESSURE: -1.024874 KBAR
ETOT DIFF (eV) : -0.000455
LARGEST GRAD (eV/A) : 0.000046
DONE(214.553462 SEC) : SETUP UNITCELL
-------------------------------------------
STEP OF RELAXATION : 3
-------------------------------------------
DONE(214.590282 SEC) : LOCAL POTENTIAL
DONE(214.674549 SEC) : SYMMETRY
DONE(214.675134 SEC) : INIT K-POINTS
DONE(214.809087 SEC) : NON-LOCAL POTENTIAL
DONE(215.334314 SEC) : INIT SCF
ITER ETOT(eV) EDIFF(eV) DRHO TIME(s)
CG1 -2.290936e+04 0.000000e+00 2.817e-06 1.944e+01
CG2 -2.290867e+04 6.940648e-01 3.133e-06 1.166e+01
CG3 -2.290867e+04 -1.491092e-05 1.431e-06 1.064e+01
CG4 -2.290867e+04 -3.878533e-06 3.863e-09 9.479e+00
----------------------------------------------------------------
TOTAL-STRESS (KBAR)
----------------------------------------------------------------
0.2348576997 -0.0000000000 0.0000000000
-0.0000000000 0.2348576997 0.0000000000
0.0000000000 0.0000000000 0.2348576997
----------------------------------------------------------------
TOTAL-PRESSURE: 0.234858 KBAR
ETOT DIFF (eV) : -0.000038
LARGEST GRAD (eV/A) : 0.000005
TIME STATISTICS
-------------------------------------------------------------------------------------
CLASS_NAME NAME TIME(Sec) CALLS AVG(Sec) PER(%)
-------------------------------------------------------------------------------------
total 269.14 41 6.56 100.00
Driver reading 0.01 1 0.01 0.00
Input Init 0.00 1 0.00 0.00
Input_Conv Convert 0.00 1 0.00 0.00
Driver driver_line 269.13 1 269.13 100.00
UnitCell check_tau 0.00 1 0.00 0.00
PW_Basis_Sup setuptransform 0.00 1 0.00 0.00
PW_Basis_Sup distributeg 0.00 1 0.00 0.00
mymath heapsort 0.03 3986 0.00 0.01
Symmetry analy_sys 157.41 2 78.71 58.49
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.06 3 0.02 0.02
ppcell_vnl init 0.02 1 0.02 0.01
ppcell_vl init_vloc 0.06 3 0.02 0.02
ppcell_vnl init_vnl 0.40 3 0.13 0.15
WF_atomic init_at_1 0.00 3 0.00 0.00
wavefunc wfcinit 0.00 1 0.00 0.00
Ions opt_ions 268.48 1 268.48 99.76
ESolver_KS_PW Run 261.38 3 87.13 97.12
H_Ewald_pw compute_ewald 0.02 3 0.01 0.01
Charge set_rho_core 0.00 3 0.00 0.00
Charge atomic_rho 0.16 6 0.03 0.06
PW_Basis_Sup recip2real 0.58 153 0.00 0.21
PW_Basis_Sup gathers_scatterp 0.24 153 0.00 0.09
Potential init_pot 0.34 3 0.11 0.13
Potential update_from_charge 2.13 20 0.11 0.79
Potential cal_fixed_v 0.02 3 0.01 0.01
PotLocal cal_fixed_v 0.02 3 0.01 0.01
Potential cal_v_eff 2.11 20 0.11 0.78
H_Hartree_pw v_hartree 0.21 20 0.01 0.08
PW_Basis_Sup real2recip 0.72 201 0.00 0.27
PW_Basis_Sup gatherp_scatters 0.31 201 0.00 0.12
PotXC cal_v_eff 1.88 20 0.09 0.70
XC_Functional v_xc 1.88 20 0.09 0.70
Potential interpolate_vrs 0.01 20 0.00 0.00
Symmetry rhog_symmetry 0.55 23 0.02 0.20
Symmetry group fft grids 0.13 23 0.01 0.05
Charge_Mixing init_mixing 0.00 3 0.00 0.00
ESolver_KS_PW hamilt2density 257.01 20 12.85 95.49
HSolverPW solve 256.26 20 12.81 95.22
Nonlocal getvnl 1.12 80 0.01 0.42
pp_cell_vnl getvnl 1.46 104 0.01 0.54
Structure_Factor get_sk 0.34 488 0.00 0.13
DiagoIterAssist diagH_subspace 19.72 68 0.29 7.33
Operator hPsi 194.76 19557 0.01 72.36
Operator EkineticPW 1.70 19557 0.00 0.63
Operator VeffPW 99.51 19557 0.01 36.97
PW_Basis_K recip2real 63.29 28813 0.00 23.51
PW_Basis_K gathers_scatterp 19.80 28813 0.00 7.36
PW_Basis_K real2recip 33.51 23773 0.00 12.45
PW_Basis_K gatherp_scatters 9.85 23773 0.00 3.66
Operator NonlocalPW 93.43 19557 0.00 34.71
Nonlocal add_nonlocal_pp 45.76 19557 0.00 17.00
DiagoIterAssist diagH_LAPACK 0.12 68 0.00 0.04
DiagoCG diag_once 224.72 80 2.81 83.49
DiagoCG_New spsi_func 1.31 38978 0.00 0.49
DiagoCG_New hpsi_func 179.23 19489 0.01 66.59
ElecStatePW psiToRho 9.38 20 0.47 3.48
Charge_Mixing get_drho 0.15 20 0.01 0.06
Charge_Mixing inner_product_recip_rho 0.01 20 0.00 0.00
Charge mix_rho 0.13 12 0.01 0.05
Charge Broyden_mixing 0.03 12 0.00 0.01
Charge_Mixing inner_product_recip_hartree 0.02 60 0.00 0.01
Forces cal_force_loc 0.07 3 0.02 0.03
Forces cal_force_ew 0.07 3 0.02 0.02
Forces cal_force_nl 1.15 3 0.38 0.43
Forces cal_force_cc 0.00 3 0.00 0.00
Forces cal_force_scc 0.12 3 0.04 0.05
Stress_PW cal_stress 5.68 3 1.89 2.11
Stress_Func stress_kin 0.17 3 0.06 0.06
Stress_Func stress_har 0.02 3 0.01 0.01
Stress_Func stress_ewa 0.06 3 0.02 0.02
Stress_Func stress_gga 0.17 3 0.06 0.06
Stress_Func stress_loc 0.21 3 0.07 0.08
Stress_Func stress_cc 0.00 3 0.00 0.00
Stress_Func stress_nl 5.04 3 1.68 1.87
ESolver_KS_PW init_after_vc 0.51 2 0.26 0.19
ModuleIO write_istate_info 0.00 1 0.00 0.00
-------------------------------------------------------------------------------------
START Time : Sat May 18 23:12:14 2024
FINISH Time : Sat May 18 23:16:43 2024
TOTAL Time : 269
SEE INFORMATION IN : OUT.ZnS/
代码
文本
5.3 优化结果
优化后的结构可以在OUT.ZnS/STRU_NOW.cif中找到：
代码
文本
[20]
!cat ./OUT.ZnS/STRU_NOW.cif
data_none _audit_creation_method generated by ABACUS _cell_length_a 5.44987 _cell_length_b 5.44987 _cell_length_c 5.44987 _cell_angle_alpha 90 _cell_angle_beta 90 _cell_angle_gamma 90 loop_ _atom_site_label _atom_site_fract_x _atom_site_fract_y _atom_site_fract_z Zn 0 0 0 Zn 0 0.5 0.5 Zn 0.5 0 0.5 Zn 0.5 0.5 0 S 0.25 0.75 0.75 S 0.25 0.25 0.25 S 0.75 0.75 0.25 S 0.75 0.25 0.75
代码
文本
优化后的结构长这样：
代码
文本
可以看到晶格结构没有变化，而晶格常数从5.44322变成了5.44987。继续提高k的取值，晶格常数会逐渐趋于稳定。
实验中测量到的晶格常量为5.41，比计算值低。原因可能是系统收敛程度不高，或者是赝势文件的选取不当，又或者是温度对晶体有影响。具体问题有待进一步分析。
代码
文本
点个赞吧
本文被以下合集收录
快速开始ABACUS | 计算化学元素单质/及其化合物的晶体结构
donglikun@dp.tech
更新于 2024-08-02
26 篇4 人关注
推荐阅读
公开
利用ABACUS完成DFT计算——以硫元素（S）为例陈贝宁
发布于 2024-05-13
1 赞
公开
快速开始ABACUS | 计算碲（Te）元素单质的晶体结构donglikun@dp.tech
更新于 2024-07-16