Vacuum reference
Finding a reference to compare between experiments or calculations can be challenging. Experimentally, there is usually a standard reference taken, e.g., the standard hydrogen electrode (SHE) under standard conditions. Computationally, it can be difficult to even compare the energies of cells if a reference is not carefully taken. The local potential cannot be used for a cell because it is relative to the individual cell. Instead, it is more useful to find a common reference, e.g., the vacuum.
We describe a procedure for finding a vacuum reference for the Fe2+/Fe3+ half-cell by taking a water slab and referencing its local potential (relative to the vacuum) relative to the O 1s states of the oxygen atoms in its water molecules. In each of the solvated ion cells, the O 1s level of the water molecules is also taken as a reference, allowing a vacuum reference to be taken for the energies of each of the solvated cells [1].
| Mind: Although we take the absolute reference for the electrochemical potential here, it is also possible to do this for other examples, such as for electron emission or catalytic processes (swap out the redox level for the valence band maximum - VBM). |
Introduction
Theory
During the Fe2+/Fe3+ half-cell reaction, one electron is transferred from the reservoir. The Helmholtz free energy difference ΔA is:
[math]\displaystyle{ \Delta A = \int_0^1 \langle U_1 - U_0 \rangle_{\lambda} d\lambda - \mu n }[/math]
where [math]\displaystyle{ \langle U_1 - U_0 \rangle }[/math] is the potential energy difference between an interacting system ([math]\displaystyle{ \lambda=1 }[/math]) and non-interacting system ([math]\displaystyle{ \lambda=0 }[/math]) taken from thermodynamic integration with coupling parameter [math]\displaystyle{ \lambda }[/math], [math]\displaystyle{ \mu }[/math] is the reference chemical potential, and [math]\displaystyle{ n }[/math] is the number of electrons involved in the reaction.
The chemical potential μ is fixed at 0 eV, the vacuum level. However, the vacuum potential is not directly accessible under bulk conditions; instead, a water slab is used as a reference to the vacuum. The O1s level is used to allow comparison between the water in the bulk ([math]\displaystyle{ \langle \epsilon_{1s,bulk} \rangle_{\lambda} }[/math]) and slab ([math]\displaystyle{ \langle \epsilon_{1s,slab} \rangle }[/math]).
| Tip: The energy of the 1s level relative to the vacuum shows analogy to calculating the work function, where a surface electron is removed to a point in the vacuum. |
[math]\displaystyle{ \mu n }[/math] can then be expressed in terms of the O1s levels, according to the following equation (cf. Fig. 1):
[math]\displaystyle{ e\Delta \bar{\phi}_{\lambda} = \mu - \langle \epsilon_{1s,slab} \rangle + \langle \epsilon_{1s,bulk} \rangle_{\lambda} }[/math]
For the Fe2+/Fe3+ half-cell reaction, [math]\displaystyle{ \langle \epsilon_{1s,bulk} \rangle_{\lambda} }[/math] becomes [math]\displaystyle{ \frac{[\langle \epsilon_{1s,Fe^{3+}} \rangle + \langle \epsilon_{1s,Fe^{2+}} \rangle]}{2} }[/math]:
[math]\displaystyle{ e \Delta \bar{\phi} = \mu - \langle \epsilon_{1s,slab} \rangle + [\langle \epsilon_{1s,Fe^{3+}} \rangle + \langle \epsilon_{1s,Fe^{2+}} \rangle]/2 }[/math]
Calculation
To obtain these three terms, we perform calculations on three different systems: Fe3+ in 64 H2O (Fe3P_64H2O), Fe2+ in 64 H2O (Fe2P_64H2O), H2O slab with 128 molecules (128H2O_slab). I.e., [math]\displaystyle{ [\mathrm{Fe}(\mathrm{H}_2\mathrm{O})_n]^{3+} }[/math], [math]\displaystyle{ [\mathrm{Fe}(\mathrm{H}_2\mathrm{O})_n]^{2+} }[/math], and the water slab.
Fe3P_64H2O), top center: Fe2+ in 64 H2O (Fe2P_64H2O), top right: bulk H2O with 64 molecules (64H2O_bulk), and bottom: H2O slab with 128 molecules (128H2O_slab).The 128H2O_slab calculations are more expensive, since so many structures (~3000) are required to achieve convergence. Only the H2O molecules' O1s levels in the center of the slab are used, shaded in grey in Fig. 2, representing the bulk. This only needs to be done once to determine the chemical potential relative to vacuum, and can then be used for any half-cell reaction within 64H2O_bulk and with the same method. The Fe2P/Fe3P_64H2O calculations, on the other hand, need to be repeated for each reaction.
Procedure
The procedure is split into two parts, the ab initio part:
- Perform an MD calculation (MLFF)
- Randomly select structures
- Calculate the energy of the O1s states and local potential (GGA)
and the post-processing:
- Calculate the average of the O1s energy levels
- Calculate the O1s level relative to the vacuum
- Calculate the difference between the local potential in the slab and the bulk
Input files
Two MD calculations need to be run in parallel between two different systems: Fe3+ in 64 H2O (Fe3P_64H2O) and Fe2+ in 64 H2O (Fe2P_64H2O). I.e., [math]\displaystyle{ [\mathrm{Fe}(\mathrm{H}_2\mathrm{O})_n]^{3+} }[/math] and [math]\displaystyle{ [\mathrm{Fe}(\mathrm{H}_2\mathrm{O})_n]^{2+} }[/math].
POSCARs
The POSCAR files can be found in the calculating the redox potential overview page.
INCAR
The INCAR files are provided in the text and discussed there.
KPOINTS
The Gamma-point only is used for the KPOINTS file:
Gamma-point only 0 Monkhorst Pack 1 1 1 0 0 0
POTCAR
Standard POTCAR files are used throughout:
PAW_PBE H 15Jun2001PAW_PBE O 08Apr2002PAW_PBE Fe_sv 23Jul2007
Ab initio
Step 1: Perform an MD calculation (MLFF)
Perform a molecular dynamics (MD) simulation for each system using a Langevin thermostat. The following INCAR file can be used:
#Molecular dynamics IBRION = 0 ISYM = 0 NSW = 100000 POTIM = 1.0 TEBEG = 298 TEEND = 298 MDALGO = 3 ISIF = 2 LANGEVIN_GAMMA = 10.0 10.0 10.0 # one for each atomic species POMASS = 2.0 16.0 55.847 RANDOM_SEED = 248489752 0 0 #Machine learning #Must monitor/change ML_LMLFF = .TRUE. # switches on machine learning ML_MODE = run
In this case, we have performed a 100000 fs MD simulation (100000 steps with a time step of 1 fs), using a Langevin thermostat.
Step 2: Randomly select structures
From each of the MD calculations in Step 1, select 100-200 random structures (for the slab, ~3000 should be used), i.e., POSCAR files, excluding the first 20000 steps to allow for time to equilibrate. Save these structures to three directories (one for each system), i.e., Fe3P_64H2O, Fe2P_64H2O, and 128H2O_slab.
Using the following Python script, you can select 200 random structures from each directory, omitting the first 20000 MD steps to allow equilibration:
from ase.io import read
from ase.io import write
import numpy as np
path1 = "$PATH_TO_MD_DIRECTORIES/"
path2 = "$PATH_TO_CORELEVEL_DIRECTORIES/"
systems = ["128H2O_slab", "Fe2P_64H2O", "Fe3P_64H2O"]
def random_md_structures(path1, path2, systems, n, cutoff):
# Read entire XDATCAR trajectory (all steps)
for system in systems:
atoms_list = read(path1 + system + "/XDATCAR", index=":")
n_struc = len(atoms_list)
n_random = np.random.randint(cutoff, n_struc, size=n)
# Write selected frames back to POSCAR format
for i, j in enumerate(n_random):
#print(i+1,j)
write(path2 + system + "/POSCAR." + str(i+1), atoms_list[j], format="vasp")
random_md_structures(path1, path2, systems, 200, 20000)
which generates 200 POSCAR.m files.
Step 3: Calculate the energy of the O1s states and local potential (GGA)
For each of the structures that you have generated, perform a first-principles (FP) calculation with ICORELEVEL = 1, LVHAR = .TRUE., and WRT_POTENTIAL = hartree ionic. The energies of the O1s levels, and the local potential will be written, which can later be used to find the vacuum reference. An example INCAR for RPBE+D3 is provided below:
ENCUT = 520.0 EDIFF = 1E-5 GGA = RP IVDW = 11 ISMEAR = 0 SIGMA = 0.10 PREC = Normal LREAL = A NELMIN = 4 ALGO = All NELM = 1000 LCHARG = .FALSE. LWAVE = .FALSE. ICHARG = 2 LVHAR = .TRUE. WRT_POTENTIAL = hartree ionic ICORELEVEL = 1 IBRION = -1 ISYM = 0 NSW = 1
In addition, for the charged systems, Fe2+:
ISPIN = 2 NELECT = 526 MAGMOM = 128*0 64*0 1*4
and Fe3+:
ISPIN = 2 NELECT = 525 MAGMOM = 128*0 64*0 1*5
An example bash script is given below for executing the calculations. Make sure to use the supporting Python scripts (O1s_level.py and vacuum_ref.py).
Important: For the slab, only the central 10 Å is used for calculating the average local potential and the average O1s levels. For the redox systems (Fe3P_64H2O and Fe2P_64H2O), only the O1s levels outside of a radius of half the cell parameter are included.
|
O1s_level.py scriptimport numpy as np
import os
import sys
def pbc_distance(r1, r2, L):
"""Compute minimum-image distance between r1 and r2 in a periodic box of size L."""
delta = (r1 - r2) % L
delta = np.where(delta > L/2, delta - L, delta)
return delta
def o1s_av(path, outcar, element_of_interest, lower_limit=0.0, upper_limit=50.0, system="bulk"):
o1s, pos, n_species, species, cell = [], [], [], [], []
n_ions, count = 0, 0
with open(path + "/" + outcar, "r") as f:
for line in f:
# Total number of ions
if "NIONS =" in line:
n_ions = float(line.split()[-1])
# Number of ions of each species
elif "ions per type =" in line:
n_species = line.split()[4:]
# Name of each species
elif "POSCAR =" in line:
species = line.split()[2:]
# Positions of all ions
elif "position of ions in cartesian coordinates (Angst)" in line:
for ion in range(0, int(n_ions)):
a = [float(x) for x in f.readline().split()]
pos.append(a)
count += 1
elif "Primitive cell" in line:
#print(line)
f.readline()
f.readline()
f.readline()
f.readline()
cell.append(float(f.readline().split()[0]))
cell.append(float(f.readline().split()[1]))
cell.append(float(f.readline().split()[2]))
elif "1s" in line:
o1s.append(float(line.split()[-1]))
'''Filters through positions so that only the O ions which are between
the upper and lower limit along the z-axis are taken, i.e., in the "bulk" of the slab.
If it contains the redox species, all O ions within a sphere of half the cell size
are excluded. By default, all positions are taken.'''
elements_dict = dict(zip(species, n_species))
n_skip = 0
index_in_limits, core_in_limits = [], []
for key in elements_dict.keys():
n_count=int(elements_dict[key])
if system == "redox":
#print(key)
pass
if key == element_of_interest:
m = 0
for ion in pos[n_skip:n_skip+n_count]:
if system == "bulk":
core_in_limits.append(o1s[m])
elif system == "slab":
if lower_limit <= ion[2] < upper_limit:
index_in_limits.append(m)
core_in_limits.append(o1s[m])
elif system == "redox":
cell_av = np.average([cell[0],cell[1],cell[2]])
cell_cubic = [cell[0],cell[1],cell[2]]
x = pbc_distance(np.array(pos[-1]), np.array(ion), np.array(cell_cubic))
magnitude = np.sqrt(x.dot(x))
if magnitude >= cell_av/2:
#print(cell_av/2, magnitude, x, ion, pos[-1], o1s[m])
index_in_limits.append(m+n_skip)
core_in_limits.append(o1s[m])
m +=1
n_skip += int(elements_dict[key])
return(np.average(core_in_limits))
def o1s_avg_all(path, element_of_interest, lower_limit=0.0, upper_limit=50.0, system="bulk"):
o1s_all = []
for outcar in (x for x in os.listdir(path + "/") if "OUTCAR." in x):
#if len(o1s_all) < 100:
o1s_all.append(o1s_av(path, outcar, element_of_interest, lower_limit, upper_limit, system))
#if bulk_slab == "slab":
# print(o1s_all[-1])
return(np.average(o1s_all))
# Expecting: python3 vacuum_ref.py z_mid z_diff system
if len(sys.argv) != 5:
print("Usage: python3 vacuum_ref.py z_mid z_diff system element")
sys.exit(1)
# Read inputs
z_mid = float(sys.argv[1])
z_diff = float(sys.argv[2])
system = sys.argv[3]
element_of_interest = sys.argv[4]
# Validate system type
if system not in ("bulk", "slab", "redox"):
print('Error: system must be either "bulk" or "slab" or "redox"')
sys.exit(1)
lower_limit, upper_limit = (z_mid - z_diff), (z_mid + z_diff)
O1s_avg = o1s_avg_all(".", element_of_interest, lower_limit, upper_limit, system)
print(O1s_avg)
vacuum_ref.py scriptimport py4vasp
import numpy as np
import sys
def local_pot_avg(path, z_lower, z_upper, z_mid, system="bulk"):
# Load calculation
calc = py4vasp.Calculation.from_path(path)
pot_dict = calc.potential.to_dict()
# Local potential (ionic + hartree)
pot_local = pot_dict["ionic"] + pot_dict["hartree"]
# Plane-average over x,y
pot_avg = np.mean(pot_local, axis=(0, 1)) # shape: (NZ,)
# Lattice length along z
a_3 = pot_dict["structure"]["lattice_vectors"][2][2]
NZ = pot_avg.shape[0]
dz = a_3 / NZ
# Original z-grid
z = np.arange(NZ) * dz
# Step 1: Extend potential with periodic images (-1, 0, 1)
extended_pot = np.zeros(3*NZ)
extended_z = np.zeros(3*NZ)
for shift in (-1, 0, 1):
base = (shift + 1) * NZ
for k in range(NZ):
extended_z[base + k] = (k + shift*NZ) * dz - (z_mid - 0.5*a_3)
extended_pot[base + k] = pot_avg[k]
# Step 2: Interpolate back onto the original z-grid
pot_interp = np.zeros(NZ)
for j in range(NZ):
# Find bracketing points
K = int((z[j] + z_mid - 0.5*a_3) / dz) + NZ # shift to center of extended array
K = max(0, min(K, 3*NZ-2)) # avoid overflow
x0, x1 = extended_z[K], extended_z[K+1]
y0, y1 = extended_pot[K], extended_pot[K+1]
pot_interp[j] = y0 + (y1 - y0) * (z[j] - x0) / (x1 - x0)
# Step 3: Find the vacuum reference and apply the vacuum shift to the potential
vac_ref = pot_interp[0] # when the slab is in the middle of the cell
#pot_interp -= vac_ref
#print(system)
if system in ("bulk", "redox"):
# Step 4: Compute average over entire cell
pot_range_avg = pot_interp.mean()
#print(f"Average potential across bulk cell is {pot_range_avg:.10f}")
elif system == "slab":
# Step 4: Compute average over z_mid ± z_diff
mask = (z >= z_lower) & (z <= z_upper)
pot_range_avg = pot_interp[mask].mean()
#print(f"Average potential from z={z_lower:.3f} to {z_upper:.3f} is {pot_range_avg:.10f}")
return z, pot_interp, pot_range_avg
# Expecting: python3 vacuum_ref.py z_mid z_diff system
if len(sys.argv) != 4:
print("Usage: python3 vacuum_ref.py z_mid z_diff system")
sys.exit(1)
# Read inputs
z_mid = float(sys.argv[1])
z_diff = float(sys.argv[2])
system = sys.argv[3]
# Validate system type
if system not in ("bulk", "slab", "redox"):
print('Error: system must be either "bulk" or "slab"')
sys.exit(1)
# Example: average potential
z_lower, z_upper = (z_mid - z_diff), (z_mid + z_diff)
path = "."
calc_temp = py4vasp.Calculation.from_path(path)
grid = np.shape(calc_temp.potential.to_dict()["ionic"])[2]
tot_avg, z, pot_avg = np.zeros(int(grid)), np.zeros(int(grid)), np.zeros(int(grid))
pot_avg_range = 0.0
z, pot_interp, pot_range_avg = local_pot_avg(path, z_lower, z_upper, z_mid, system)
#print("system, z_mid, z_diff, path, z_lower, z_upper")
#print(system, z_mid, z_diff, path, z_lower, z_upper)
print(pot_range_avg)
#!/bin/bash
VASP="mpirun -np 32 $PATH_TO_VASP/vasp_gam"
rm -rf vac_ref.dat
echo "N_struc O1s_level local_potential_avg" >> vac_ref.dat
for a in POSCAR.*; do
echo $a
cp $a POSCAR
$VASP
# Use if you don't want to run the Python scripts at the same time as the calculations
#cp OUTCAR OUTCAR.$suffix
#cp vaspout.h5 vaspout.$suffix.h5
#cp LOCPOT LOCPOT.$suffix
suffix="${a#*.}"
# Use the first for the Fe2+/3+ systems, and the second for the bare slab
echo $suffix `python3 O1s_level.py 25.0 5.0 redox O` `python3 vacuum_ref.py 25.0 5.0 redox` >> vac_ref.dat
#echo $suffix `python3 O1s_level.py 25.0 5.0 slab O` `python3 vacuum_ref.py 25.0 5.0 slab` >> vac_ref.dat
done
An alternative script exists if you are taking the local potential from LOCPOT, rather than vaspout.h5, though it has not been as extensively tested:
read_locpot.py scriptimport numpy as np
def read_locpot_single(filename, ZMID, DIFF):
"""
Process a single LOCPOT file (equivalent to NCONF=1 case of the Fortran code).
All arrays should be NumPy arrays.
"""
IREFZ, LSHIFT, MITYP = 0, 1, 0
global LOCPOTZ, RHOZ, POSION
print("---------------LOCPOT--------------------")
print("Reading single LOCPOT:", filename)
NWRITE = 5
TOTAL, NCOL, NMOD = 0, 0, 0
ZMIN = ZMID - DIFF
ZMAX = ZMID + DIFF
A = np.zeros((3, 3))
NITYP = 0
with open(filename, "r") as f:
# VASP header (8 lines)
f.readline()
f.readline()
for i in range(3):
A[i] = np.fromstring(f.readline(), sep=' ')
f.readline()
NITYP = np.fromstring(f.readline(), sep=' ', dtype=int)
print(NITYP)
NTYP = len(NITYP)
f.readline()
#for i in range(NTYP):
# for j in range(int(NITYP[i])):
# f.readline()
#print("first", f.readline())
# Allocate POSION now that NTYP and NITYP are known
POSION = np.zeros((3, int(max(NITYP)), NTYP))
# Read atomic positions
for j in range(NTYP):
for k in range(NITYP[j]):
x, y, z = map(float, f.readline().split())
POSION[0, k, j] = x
POSION[1, k, j] = y
POSION[2, k, j] = z
f.readline()
NX, NY, NZ = np.fromstring(f.readline(), sep=' ', dtype=int)
LOCPOTZ = np.zeros(NZ)
RHOZ = np.zeros(NZ)
TOTAL = NX * NY * NZ
NCOL = TOTAL // NWRITE
NMOD = TOTAL % NWRITE
MITYP = NITYP.max()
DX = A[0,0] / NX
DY = A[1,1] / NY
DZ = A[2,2] / NZ
# RHOZ accumulation (only if j == 1 in Fortran, zero-based j == 1)
for k in range(NITYP[j]):
for L in (-1, 0, 1):
XYZ = np.zeros(3)
XYZ[0] = POSION[0, k, j] * A[0,0]
XYZ[1] = POSION[1, k, j] * A[1,1]
XYZ[2] = (POSION[2, k, j] + L) * A[2,2]
if j == 1: # Fortran J==2
iz = int(XYZ[2] * NZ / A[2,2]) + 1
if 1 <= iz <= NZ:
RHOZ[iz-1] += (18.0 / 6.022045e23 / (A[0,0]*A[1,1]*DZ*1e-24))
# Skip blank lines before potential data
#f.readline()
#f.readline()
# Read LOCPOT grid (flattened)
LOCPOT_HELP = np.zeros(TOTAL)
idx = 0
for _ in range(NCOL):
vals = list(map(float, f.readline().split()))
LOCPOT_HELP[idx:idx+NWRITE] = vals
idx += NWRITE
if NMOD > 0:
vals = list(map(float, f.readline().split()))
LOCPOT_HELP[idx:idx+NMOD] = vals
# Average LOCPOT over x,y for each z-plane
LOCPOTZ_HELP1 = np.zeros(NZ)
for j in range(NZ):
plane = LOCPOT_HELP[j*NX*NY:(j+1)*NX*NY]
LOCPOTZ_HELP1[j] = plane.mean()
# Build extended arrays for interpolation
LOCPOTZ_HELP2 = np.zeros(3*NZ)
Z_HELP2 = np.zeros(3*NZ)
for shift in (-1, 0, 1):
base = (shift + 1) * NZ
for k in range(NZ):
L = base + k
Z_HELP2[L] = (L * DZ) - A[2,2] - (ZMID - 0.5*A[2,2])
LOCPOTZ_HELP2[L] = LOCPOTZ_HELP1[k]
# Interpolate back to correct Z grid
LOCPOTZ_HELP3 = np.zeros(NZ)
Z_HELP3 = np.arange(NZ) * DZ
for j in range(NZ):
K = int((Z_HELP3[j] + A[2,2] + ZMID - 0.5*A[2,2]) * NZ / A[2,2])
K = max(0, min(K, 3*NZ-2))
x0, x1 = Z_HELP2[K], Z_HELP2[K+1]
y0, y1 = LOCPOTZ_HELP2[K], LOCPOTZ_HELP2[K+1]
LOCPOTZ_HELP3[j] = y0 + (y1 - y0)*(Z_HELP3[j] - x0)/(x1 - x0)
# Apply shift if required
if LSHIFT == 1:
LOCPOTZ_VAC = LOCPOTZ_HELP3[IREFZ-1]
SHIFTZ = -LOCPOTZ_VAC
else:
SHIFTZ = 0.0
LOCPOTZ_HELP3 += SHIFTZ
# Compute average potential over ZMIN–ZMAX
mask = (Z_HELP3 >= ZMIN) & (Z_HELP3 <= ZMAX)
LOCPOTZ_AV = LOCPOTZ_HELP3[mask].mean()
print("Potential =", LOCPOTZ_AV)
# This is the final averaged LOCPOT vs z
return LOCPOTZ_HELP3, RHOZ
read_locpot_single("LOCPOT", 18.175, 3.5)
Important: Make sure to first download and save the O1s_level.py and vacuum_ref.py scripts. Calculating these values while running the jobs will save a non-trivial amount of disk space. Alternatively, you can save all of the OUTCAR and vaspout.h5 files to process later.
|
These scripts will create a vac_ref.dat file containing the number, the average 1s energy levels for O in water molecules, and the average local potential for each structure:
N_struc O1s_level local_potential_avg 1 -509.526392 -1.6969455464219334 2 -509.4371388571428 -1.5132796708524134 3 -509.30019053061227 -1.0734347143520142 4 -509.4201625058023 -1.973256948325608 ...
These will be required in the next step.
Post-processing
Step 4: Calculate the average of the O1s energy levels
Once the FP calculations have finished, extract the O1s energy levels [math]\displaystyle{ \epsilon_{1s} }[/math]. You can find them listed in OUTCAR after the core state eigenenergies. For each structure, average these energy levels over all O, except for 128H2O_slab, where only O's within the central 10 [math]\displaystyle{ \AA }[/math] of the cell should be used to mimic the bulk. With the average O1s for each individual structure, average over all structures in the ensemble, to get a single value for each system.
If you used the run script structure that we suggested above, then you will have a vac_ref.dat file. In the first column, there is the structure number, followed by the average O1s level for water molecules, then the average local potential [math]\displaystyle{ \mu }[/math].
You have now obtained the O1s level. Our values for different numbers of structures compared to the literature values are compiled in the following table:
| System | GGA (RPBE+D3) | Lit. [1], a | ||
| No. of structures | 100 | 200 | 400 | |
| [math]\displaystyle{ \langle \epsilon_{1s,slab} \rangle }[/math], 128H2O slab | -509.43 | -509.46 | -509.45 | -511.23 |
| [math]\displaystyle{ \langle \epsilon_{1s,Fe^{3+}} \rangle }[/math], Fe3+ + 64H2O (λ = 0) | -507.45 | -507.45 | -507.45 | -507.43 |
| [math]\displaystyle{ \langle \epsilon_{1s,Fe^{2+}} \rangle }[/math], Fe2+ + 64H2O (λ = 1) | -507.50 | -507.50 | -507.50 | -507.50 |
aThe number of structures used in Jinnouchi et al. is 200 for bulk water and redox systems, and 3000 for the water slab, cf. Supplementary Table 4 for the averaged 1s levels of O atoms in Jinnouchi et al. for reference [1].
There is a small difference between the literature for the bulk and redox systems using 100 structures, which is narrowed when using 200. The difference for the water slab is much larger, as the reference with respect to the vacuum has not been given, hence the missing ~1.8 eV. The average O1s level is converged for the water slab, which tells us that it is the local potential that requires many structures to converge.
Step 5: Calculate the O1s level relative to the vacuum
The value of the local potential [math]\displaystyle{ \mu }[/math] must be included to obtain an absolute redox potential; otherwise, comparing between different cells (as in this example) or between different systems cannot be done as no common reference exists. It is important to include this because the local potential means that the vacuum reference is not equal to 0 eV, it is more negative ~-1.8 eV:
128H2O_slab). Increasing from 100 to 400 structures partially smoothes out the potential, but much larger numbers of structures are required to smooth it out completely (cf. the right-hand side of Fig. 1b in [1]). You can see that the local potential is ~-1.8 eV. This is the missing term from the previous section.You can calculate the local potential using LVHAR = .TRUE.. This saves the potential to LOCPOT, from which you could perform all of the following steps to calculate the average. We instead used WRT_POTENTIAL = hartree ionic to save the local potential to vaspout.h5 and analyse using py4vasp. This gives
| System | GGA (RPBE+D3) | ||
| No. of structures | 100 | 200 | 400 |
| [math]\displaystyle{ -\langle \mu_{slab} \rangle }[/math], 128H2O slab | -1.686 | -1.735 | -1.729 |
| [math]\displaystyle{ -\langle \mu_{Fe^{3+}} \rangle }[/math], Fe3+ + 64H2O (λ = 0) | 0.239 | 0.091 | 0.089 |
| [math]\displaystyle{ -\langle \mu_{Fe^{2+}} \rangle }[/math], Fe2+ + 64H2O (λ = 1) | -0.905 | -0.680 | -0.734 |
The average local potential only needs to be considered relative to the vacuum, i.e., for the centre of the water slab. The bulk and redox systems do not need to be corrected (cf. Fig. 3 below, and Fig. 1c and Eq. 7 of [1]).
The updated O1s energy for the water slab is then much closer to the literature:
| System | GGA (RPBE+D3) | Lit. [1], a | ||
| No. of structures | 100 | 200 | 400 | Std POTCAR |
| [math]\displaystyle{ -\mu + \langle \epsilon_{1s,slab} \rangle }[/math], 128H2O slab | -511.12 | -511.20 | -511.18 | -511.23 |
aThe number of structures used in Jinnouchi et al. varies from 200 for bulk water and redox systems, to 3000 for the water slab, cf. Supplementary Table 4 for the averaged 1s levels of O atoms in Jinnouchi et al. for reference [1].
Step 6: Calculate the difference between the local potential in the slab and the bulk
The redox level [math]\displaystyle{ e \Delta \bar{\phi} }[/math] (in V) can now be obtained (note that e=1 for this half-cell reaction):
[math]\displaystyle{ e \Delta \bar{\phi} = \frac{[\langle \epsilon_{1s,Fe^{3+}} \rangle + \langle \epsilon_{1s,Fe^{2+}} \rangle]}{2} - \langle \epsilon_{1s,slab} \rangle + \mu }[/math]
[math]\displaystyle{ = \frac{[(-507.45) + (-507.50)]}{2} - (-509.45 + -1.73) = 3.71 }[/math] (this example)
[math]\displaystyle{ e \Delta \bar{\phi}_{Lit.} = \frac{[(-507.43) + (-507.50)]}{2} - (-511.23) = 3.74 }[/math] (for Std. POTCAR provided by R. Jinnouchi; GW POTCARs in lit. [1])
We agree to within 0.1 eV of the literature, having used only 400 structures for the slab, rather than the 3000 in the paper. Having obtained [math]\displaystyle{ e \Delta \bar{\phi} }[/math], you now have the reference to the vacuum, necessary for calculating the absolute redox potentials.
In the next stage, you will perform TI from Fe3+ to Fe2+ using MLFFs.
| Tip: You can find a more detailed explanation of calculating the core levels in this paper: [2] |
Recommendations and advice
- We recommend reading the local potential from vaspout.h5, rather than LOCPOT, as you can take advantage of the HDF5 tools for data extraction.
- 3000 slab calculations will create a lot of large OUTCAR, vaspout.h5, and LOCPOT files (several 100 GBs). Test your setup using the Python scripts in your run script to avoid this.
Related tags and articles
- Tags
- Files
- Theory
- How-tos
- Construction:Thermodynamic integration: redox potential
- Construction:Thermodynamic integration between MLFFs
- Construction:Thermodynamic integration between MLFF and GGAs
References
- ↑ a b c d e f g h i R. Jinnouchi, F. Karsai, and G. Kresse, Machine learning-aided first-principles calculations of redox potentials, npj Comput. Mater. 10, 107 (2024).
- ↑ R. Jinnouchi, F. Karsai, and G. Kresse, Absolute standard hydrogen electrode potential and redox potentials of atoms and molecules: machine learning aided first principles calculations, Chem. Sci. 16, 2335 (2025).