目前文章中的电催化理论计算，大多是基于Nørskov组于2004年发展的计算氢电极(CHE)模型进行的。
Nørskov, J. et al. Origin of the overpotential for oxygen reduction at a fuel-cell cathode. J. Phys. Chem. B 2004, 108, 17886-17892.
该模型基于电中性结构进行计算，认为中间体能量不随电极电势改变，通过电子的eU修正考虑外加电势的作用；同时质子和电子耦合加入，维持体系电中性。由于模型的简单易用性，在HER、OER、ORR、CO2RR等领域获得了广泛应用。
以HER的Volmer步为例：
* + H+ + e- = H*
电中性方法(CNM)的能量计算：
ΔGcnm = G(H*) - G(*) - G(H+) - μe
= G(H*) - G(*) - 1/2G0(H2) + 0.0592pH + |e|U(SHE)
假设G(H*)、G(*) 、G(H*) - G(*)不随U(SHE)变化
相对于RHE：
ΔGcnm = G(H*) - G(*) - 1/2G0(H2) + |e|U(RHE)
与pH无关
但是显然，当模型偏离零电荷电势时，模型必然是带电的，巨势能与电极电势间通常呈现二次型关系，并且由于反应各步模型零电荷电势与电容值的不同，这种误差也无法抵消，由此引起的电催化pH效应也无法考虑。
因此，为了更准确地计算电化学反应的能量变化，我们需要采用巨正则系综的恒电势计算。通过调整体系的电子数，进而调整体系的费米能级，以达到目标电极电势。
恒电势方法(CPM)的能量计算：
*(Q1) + H+ + (1 + Q1 - Q2)e- = H*(Q2)
电势相等时Q1通常不等于Q2，实际电荷改变量不为1
ΔGcpm = G(H*(Q2)) - G(*(Q1)) - G(H+) - (1 + Q1 - Q2)μe
= G(H*(Q2)) - G(*(Q1)) - G(H+) - μe - (Q1 - Q2)μe
= G(H*(Q2)) - G(*(Q1)) - 1/2G0(H2) + 0.0592pH + |e|U(SHE) - (Q1 - Q2)μe
Q1, Q2：*，H* 所带净电荷（正）
- (Q1 - Q2)μe：损失电子能量补偿(Qμ)，JDFTx折合入G(H*(Q2))、G(*(Q1))中
相对于RHE：
ΔGcpm = G(H*(Q2)) - G(*(Q1)) - 1/2G0(H2) + |e|U(RHE) - (Q1 - Q2)μe
ΔGcpm = [G(H*(Q2)) + Q2μe] - [G(*(Q1)) + Q2μe] - 1/2G0(H2) + |e|U(RHE)
损失电子能量补偿，并校正到真空μ=0：
-N0μ+Qμ = -(N0 - Q)μ= -Nμ
H2能量亦需校正到真空μ=0：-2μ
不同pH下，同一 U(RHE) 对应着不同的 U(SHE) ，对应于不同的表面带电状态，进而影响ΔG
开源的JDFTx软件，为巨正则系综的恒电势计算提供了一个方便的平台：
R. Sundararaman, K. Letchworth-Weaver, K.A. Schwarz, D. Gunceler, Y. Ozhabes and T.A. Arias, 'JDFTx: software for joint density-functional theory', SoftwareX 2017, 6, 278.
但是由于目前系统的教程较少，为JDFTx的使用制造了阻碍，因此在这里留下一些文字，供大家参考
【基本流程】
官方网站：http://jdftx.org/
教程与练习：Tutorials部分
关键词查询：Input file documentation部分
由于JDFTx计算较慢，通常不用JDFTx进行结构优化，而是先用VASP对模型进行几何结构优化和自由能校正
利用VASP计算产生的CONTCAR文件，和cont2jdftx.py脚本，生成contcar.ionpos和contcar.lattice：
chmod +x ~/bin/cont2jdftx.py
python ~/bin/cont2jdftx.py
再利用JDFTx在溶剂化和不同带电条件下计算单点能
为辅助收敛，建议保留波函数并读取上一步波函数，依次进行：真空中、隐式溶剂中、隐式溶剂+电势1、隐式溶剂+电势2、... ... 的计算
电荷密度、静电势、Bound电荷文件等都可以用JDFTx自带的createXSF脚本转换为VESTA可以打开的文件：
createXSF XXX.out XXX.xsf XXX.nbound
【JDFTx输入文件模板！】
#结构文件
include contcar.lattice #晶胞参数，每一列代表一个晶格矢量，原子单位制Bohr
include contcar.ionpos #原子坐标
coords-type Lattice #分数坐标，默认为笛卡尔坐标Cartesian
#赝势
#GBRV超软赝势，用$ID代替任意元素，优先使用v1.2版赝势，如果没有该元素再用v1版赝势
ion-species GBRV/$ID_pbe_v1.2.uspp
ion-species GBRV/$ID_pbe_v1.uspp
#K点
kpoint-folding 1 1 1 #n×n×n K点，默认 1 1 1
kpoint 0 0 0 1. #K点位移，默认Gamma centered
#kpoint 0.5 0.5 0.5 1. #K点位移，MP方法
#电子步
elec-ex-corr gga-PBE #指定交换-相关泛函，默认GGA-PBE泛函
spintype z-spin #自旋极化，z-spin = up/dn polarization (non-relativistic spin)
#elec-initial-magnetization 2.0 yes #总磁矩初猜，Nup-Ndw yes = magnetization fixed, no = optimize it
initial-magnetic-moments Fe 3.0 #每个磁性原子总磁矩初猜，<species> <M1> <M2> ... <Mn> [<species2> ...]
add-U Fe d 0.120905 #DFT+U，<species> <orbDesc> <UminusJ(Hartree)>
#wavefunction LCAO # Initialize state from atomic orbitals (random by default)
elec-cutoff 20 100 # 平面波截断能，电荷密度截断能，Hartree
electronic-minimize \
nIterations 1000 \ #最大迭代次数
energyDiffThreshold 1e-6 #能量收敛标准
elec-smearing Gauss 0.001837 #0.05eV Gauss Smearing
#elec-n-bands 32 #NBANDS
#converge-empty-states yes #确保计算结束时空态收敛
#density-of-states Total #输出总态密度(TDOS)
#DFT-D2色散校正：
van-der-waals
#溶剂化计算：
fluid LinearPCM #Simple continuum fluid with linear response
pcm-variant CANDLE #default: GLSSA13
fluid-solvent H2O
fluid-cation K+ 0.1 #mol/L，对于LinearPCM方法，带同电荷的离子是等价的
fluid-anion F- 0.1
#外加电极电势：
target-mu -0.171252 # U = 0.00V vs SHE (4.66V)，1Hartree = 27.2114eV
#target-mu -0.189626 # U = +0.5V vs SHE
#target-mu -0.208001 # U = +1.0V vs SHE
#target-mu -0.226376 # U = +1.5V vs SHE
#target-mu -0.152877 # U = -0.5V vs SHE
#target-mu -0.134502 # U = -1.0V vs SHE
#target-mu -0.116128 # U = -1.5V vs SHE
#输出
initial-state sol.$VAR #读取之前计算的波函数
dump-name u0.$VAR #输出文件名
dump End Ecomponents State ElecDensity Dtot BoundCharge #输出文件
#Ecomponents - 能量 - .Ecomponents
#ElecDensity - 电荷密度 - .n
#EigStats - 本征值 - .eigStats
#Dtot - 静电势 - .d_tot
#RealSpaceWfns - 实空间波函数 - .wfns_0_x.rs
#State - 打印重启计算所需的所有变量: 波函数和流体状态/填充(如果有)
#BoundCharge - Bound电荷 - .nbound
#库伦势能截断：
#相当于VASP里IDIPOL，LDIPOL，DIPOL
#分子体系：
#coulomb-interaction Isolated #Switch from periodic to isolated images mode #coulomb-truncation-embed 0 0 0 #Center of molecule;
#对于表面体系：
#coulomb-interaction Slab 001
#coulomb-truncation-embed 0 0 0
#几何结构优化：
#ionic-minimize \
# nIterations 100 \ #最大离子步数
# energyDiffThreshold 1e-6 \ #能量收敛标准
# knormThreshold 1e-4 #Threshold on RMS cartesian force #RMS力收敛标准
#频率计算：
#vibrations \
# dr 0.01 \
# centralDiff yes \ #有限位移法
# translationSym yes \ #是否计算平动和转动，分子yes，固体和表面no
# rotationSym yes \
# T 298 #温度
欢迎大家补充、指正
