固定电势方法在ABACUS中的实现
©️ Copyright 2024 @ Authors
作者：孙梦琳 📨 刘裕 📨
日期：2024-10-20
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的 开始连接 按钮，选择 fcp-abacus-ase:202410 镜像及 c2_m4_cpu 节点配置，稍等片刻即可运行。由于具体计算的案例都需要比较长的计算时间，本案例不会要求完成具体计算，故使用c2_m4_cpu即可。
1. 背景
电催化体系微观尺度计算模拟是探究相关机理不可或缺的研究手段，在实验中，电化学表面反应的工作电极电位由外部电路控制，电势在反应过程中保持恒定。但电极或催化剂材料表面和溶液共同构成的电化学开放系统对计算模拟具有很高的挑战性，在传统的电子结构理论方法中，基于密度泛函理论的第一性原理计算通常处理的是正则系综，即计算模拟电化学反应通常是基于恒定的电子数，是在恒定电荷条件下进行，模拟过程中，材料表面的电荷密度以及相应的电势会发生剧烈变化。计算反应初、末态的电势往往会有很大差别，从而忽略了电势对实验中反应性的影响，很难准确描述一个带电体系基本性质。 因此，电化学界面的准确建模，对于理解电极上的电化学和电催化至关重要。这要求电极电位（费米能级）是固定的，而电子总数可以在原子水平上变化。
清华大学肖海课题组基于牛顿方法开发了一种高效的全收敛固定电势（fully converged constant-potential, FCP）算法,其能够高效地将模拟体系收敛到预设的电化学势。我们基于ASE的Python接口，使用ABACUS作为其方法的第一性原理计算器实现了恒电势计算。
首先，我们进入镜像中提供的测试算例文件，tests文件夹中提供了几个简单的例子的输入与输出。其中 Pt333-fcp-scf为静态自洽计算，Pt333-fcp-relax为结构优化计算。neb-Pt333为NEB计算。
/FCP-abacus-ase/tests
Pt333-fcp-relax/ Pt333-fcp-scf/ neb-Pt333/
设置运行脚本
from ase.calculators.FCPelectrochem import FCP
from ase.calculators.abacus import Abacus, AbacusProfile
from ase.io import read
from ase.optimize import LBFGS
# abacus parameters
pseudo_dir = r'/ABACUS/PP'
basis_dir = r'/ABACUS/ORB'
pp = {'Pt':'Pt_ONCV_PBE-1.0.upf','H':'H_ONCV_PBE-1.0.upf','O':'O_ONCV_PBE-1.0.upf'}
basis = {'Pt':'Pt_gga_7au_100Ry_4s2p2d1f.orb','H':'H_gga_6au_100Ry_2s1p.orb','O':'O_gga_7au_100Ry_2s2p1d.orb'}
kpts = [3,3,1]
profile = AbacusProfile(argv=['mpirun', '-n', '32', 'abacus'])
cal_abacus = Abacus(profile=profile, pp=pp, pseudo_dir=pseudo_dir, basis=basis, basis_dir=basis_dir, kpts=kpts,ntype=3, ecutwfc=60, scf_nmax=200, basis_type='lcao', smearing_method='gaussian', smearing_sigma=0.001,
mixing_gg0=1.5, mixing_beta=0.4, mixing_ndim=20, mixing_type='pulay',scf_thr=5.000000e-05, calculation='scf', force_thr=0.04,
cal_force=1, cal_stress=1, out_stru=1, imp_sol=0, efield_flag=1, dip_cor_flag=1, efield_amp=0, efield_dir=2, efield_pos_max=0.7, gate_flag=1,
nelec=543, zgate=0.5, vdw_method='d3_0')
# FCP parameters
cal_FCP=FCP(innercalc = cal_abacus,
fcptxt = 'log-fcp.txt',
U = -1.5,
NELECT = 543,
C = 1/1000, #1/k capacitance per A^2
FCPmethod = 'Newton-fitting',
FCPconv = 0.01,
NELECT0 = 543,
adaptive_lr = False,
work_ref = 4.4,
max_FCP_iter = 10000
)
#________________________________________________________________
atoms=read(filename='STRU', format='abacus')
atoms.calc = cal_FCP
atoms.get_potential_energy()算例中提供了一个运行脚本示例，脚本的写法需要利用ASE-ABACUS接口，即ase.calculators.abacus模块来调用ABACUS进行计算时保持一致。注意有两个位置需要修改：一是需添加 # FCP parameters 部分参数，二是设置calculator为 cal_FCP。
# FCP parameters 涉及的主要参数有：
U: 电极电势 (V vs. reference electrode)NELECT: 体系总电子数的初猜NELECT0: 零电荷电位的电子数C: 单位表面积电容的初猜（e/V/(Å)）
其中，NELECT和C的初猜值会影响计算的收敛效率。具体参数设置规范请参考 https://github.com/YuLiu98/FCP-abacus-ase，并且其README中有详细的环境配置和使用方法介绍。
这些例子在具体进行FCP计算时会耗费较长的运算时间,感兴趣的用户可以自行计算。去除下面的#注释即可运行.
我们以经典的氢还原反应中的Volmer步（H + e + * H*）为例，这是一个经典的涉及电子转移的电化学反应过程，往往需要考虑固定电势条件。提供的算例体系均为Pt(333)表面，表面上方放置一层HO分子模拟溶液体系，其中包含一个额外的溶剂化质子 H。
在电势收敛过程中需要不断调整体系电子数，我们已在 ABACUS 中实现了带有补偿电荷平板的功能，能够实现在调整电子总数的情况下保持体系中性，其位置也是可调整的。具体可参看 ABACUS 手册中 INPUT 参数介绍的Gate field (compensating charge)部分。https://abacus.deepmodeling.com/en/latest/advanced/input_files/input-main.html#gate-field-compensating-charge
以relax计算为例，执行运行脚本后，我们直接来看计算输出，FCP的计算输出相比原用ASE调用ABACUS多了几个输出文件：log-fcp.txt, tmp-log-fcp.txt, K.txt 。K.txt 输出每一步电子步体系的电容倒数，tmp-log-fcp.txt是在给定离子构型的情况下不断改变电子数做scf的相关输出。log-fcp.txt输出了做relax时每个离子步的收敛结果。
/FCP-abacus-ase/tests/Pt333-fcp-relax
loop NELECT Fermi(eV) Fermishift(eV) mu(eV) U(V) conv(V) Ewithoutentropy(eV) Ewithoutentropy_grand(eV) Etoten(eV) Etoten_grand(eV) Cpersurf(e/V/A^2) 3 543.130310 3.767617 -6.659328 -2.900 -1.508289 0.008289 -92064.163376 -92064.653252 -92064.177687 -92064.667563 0.001006 1 543.129796 3.766371 -6.660627 -2.900 -1.505744 0.005744 -92064.255032 -92064.743145 -92064.269354 -92064.757467 0.001006 2 543.128482 3.761937 -6.662747 -2.900 -1.499189 -0.000811 -92064.712250 -92065.195695 -92064.724498 -92065.207944 0.000956 1 543.128530 3.753041 -6.661112 -2.900 -1.491929 -0.008071 -92064.882292 -92065.365707 -92064.891261 -92065.374676 0.000956 1 543.129005 3.752971 -6.662390 -2.900 -1.490581 -0.009419 -92064.940663 -92065.426031 -92064.950630 -92065.435998 0.000956 2 543.132067 3.754758 -6.655688 -2.900 -1.499071 -0.000929 -92065.069203 -92065.565207 -92065.082058 -92065.578062 0.000977 1 543.132123 3.753457 -6.658992 -2.900 -1.494465 -0.005535 -92065.080510 -92065.577160 -92065.093596 -92065.590246 0.000977
从输出可以看出，在结构优化的最后一步，实现了在误差范围内Fermi + Fermishift = mu，即有效控制了体系电势，其中mu表示我们希望设置的电极电势值（vs vacuum），由输入参数中的U和 work_ref共同控制。
此外，还可以利用ATST-Tools进行固定电势的NEB计算。在安装好ATST-Tools后，https://bohrium.dp.tech/notebooks/39369325971，在source/abacus_neb.py中添加from ase.calculators.FCPelectrochem import FCP，并在脚本中加入如下代码：
class FcpNEB(AbacusNEB):
"""Customize Nudged Elastic Band calculation workflow by using FCP-abacus-ase"""
def __init__(self, init_chain, fcp_parameters, parameters, abacus='abacus',
dyneb=False, k=0.1, algorism="improvedtangent",
directory='NEB', mpi=1, omp=1, parallel=True, ) -> None:
AbacusNEB.__init__(self, init_chain, parameters, abacus, dyneb, k, algorism, directory, mpi, omp, parallel)
"""parameters for FCP"""
self.fcp_parameters = fcp_parameters
def set_calculator(self):
"""Set FCP calculators"""
os.environ['OMP_NUM_THREADS'] = f'{self.omp}'
profile = AbacusProfile(command=f"mpirun -np {self.mpi} {self.abacus}")
if self.parallel:
out_directory = f"{self.directory}-rank{world.rank}"
else:
out_directory = self.directory
calc = Abacus(profile=profile, directory=out_directory, **self.parameters)
calc_fcp = FCP(innercalc = calc, **self.fcp_parameters)
return calc_fcp
重新编译后，即可利用ATST-Tools工具实现基于ASE的固定电势的过渡态计算。运行脚本可参考/FCP-abacus-ase/tests/neb-Pt333/中的runeb.py，与SCF计算和relax计算同理，脚本中需要添加的修改为# FCP parameters 部分参数，并设置calculator为 cal_FCP。
from ase.calculators.FCPelectrochem import FCP
from ase.io import read
from ase.constraints import FixAtoms
from ase.calculators.emt import EMT
from ase.mep import NEB
from ase.optimize import BFGS
from ase.calculators.abacus import Abacus, AbacusProfile
pseudo_dir = r'/ABACUS/PP'
basis_dir = r'/ABACUS/ORB'
pp = {'Pt':'Pt_ONCV_PBE-1.0.upf','H':'H_ONCV_PBE-1.0.upf','O':'O_ONCV_PBE-1.0.upf'}
basis = {'Pt':'Pt_gga_7au_100Ry_4s2p2d1f.orb','H':'H_gga_6au_100Ry_2s1p.orb','O':'O_gga_7au_100Ry_2s2p1d.orb'}
kpts = [3,3,1]
profile = AbacusProfile(command='mpirun -n 32 abacus')
cal = Abacus(profile=profile, pp=pp, pseudo_dir=pseudo_dir, basis=basis, basis_dir=basis_dir, kpts=kpts,ntype=3, ecutwfc=60, scf_nmax=200, basis_type='lcao', smearing_method='gaussian', smearing_sigma=0.001,
mixing_gg0=1.5, mixing_beta=0.4, mixing_ndim=20, mixing_type='pulay',scf_thr=5.000000e-05, calculation='scf', force_thr=0.04,
cal_force=1, cal_stress=1, out_stru=1, imp_sol=0, efield_flag=1, dip_cor_flag=1, efield_amp=0, efield_dir=2, efield_pos_max=0.7, gate_flag=1,
nelec=543, zgate=0.5, vdw_method='d3_0')
initial = read('initial.traj')
final = read('final.traj')
images = [initial]
for i in range(3):
image = initial.copy()
#image.calc = EMT()
image.calc = Abacus(profile=profile, pp=pp, pseudo_dir=pseudo_dir, basis=basis, basis_dir=basis_dir, kpts=kpts,ntype=3, ecutwfc=60, scf_nmax=200, basis_type='lcao', smearing_method='gaussian', smearing_sigma=0.001,
mixing_gg0=1.5, mixing_beta=0.4, mixing_ndim=20, mixing_type='pulay',scf_thr=5.000000e-05, calculation='scf', force_thr=0.04,
cal_force=1, cal_stress=1, out_stru=1, imp_sol=0, efield_flag=1, dip_cor_flag=1, efield_amp=0, efield_dir=2, efield_pos_max=0.7, gate_flag=1,
nelec=543, zgate=0.5, vdw_method='d3_0')
#image.set_constraint(constraint)
images.append(image)
images.append(final)
# FCP parameters
cal_FCP=FCP(innercalc = cal,
fcptxt = 'log-fcp.txt',
U = -1.5,
NELECT = 543,
C = 1/1000, #1/k capacitance per A^2
FCPmethod = 'Newton-fitting',
FCPconv = 0.01,
NELECT0 = 543,
adaptive_lr = False,
work_ref = 4.4,
max_FCP_iter = 10000
)
for image in images:
image.calc = cal_FCP
neb = NEB(images,allow_shared_calculator=True)
neb.interpolate(apply_constraint=False)
qn = BFGS(neb, trajectory='neb.traj')
qn.run(fmax=0.3)
我们基于ASE ，构建了完全收敛的恒电势 （FCP）算法的ABACUS接口，实现了使用ABACUS作为恒电势运算的DFT计算器，预期能在ABACUS中实现更广泛的电极表面系统模拟，处理更多涉及电子转移的电化学体系。
[1] Xia, Z. & Xiao, H. Grand Canonical Ensemble Modeling of Electrochemical Interfaces Made Simple. Journal of Chemical Theory and Computation (2023) doi:10.1021/acs.jctc.3c00237.
