Solvated Jellium Method (SJM)

Theoretical Background

The Solvated Jellium method (SJM) is simple method for the simulation of electrochemical interfaces in DFT. A full description of the model can be found in [1]. It can be used like the standard GPAW calculator, meaning stable intermediates and reaction barriers can be calculated at defined electrode potential via e.g. the Nudged Elastic Band method (NEB) [2].

The basis of the model is keeping control of the electrode potential by charging the electrodes interface, while keeping the periodic unit cell charge neutral. This is done by adding a JelliumSlab in the region above the electrode surface. Doing so both electrons/holes in the SCF cycle and spatially constant counter charge are introduced, therefore keeping the net charge at 0.

Additionally, an implicit solvent [3] is introduced above the slab, which screens the electric field created by dipole consisting of electrode and counter charge.

The electrode potential is then defined as the Fermi Level (\(\mu\)) referenced to the electrostatic potential deep in the solvent, where the whole charge on the electrode has been screened.

\[\Phi_e = \Phi_w - \mu.\]

The energy used in the analysis of electrode reactions is the Grand Potential Energy

\[\Omega = E_{tot} + \Phi_e N_e\]
class gpaw.solvation.sjm.SJM(ne=0, doublelayer=None, potential=None, dpot=0.01, tiny=1e-08, verbose=False, **gpaw_kwargs)[source]

Subclass of the SolvationGPAW class which includes the Solvated Jellium method.

The method allows the simulation of an electrochemical environment by calculating constant potential quantities on the basis of constant charge DFT runs. For this purpose, it allows the usage of non-neutral periodic slab systems. Cell neutrality is achieved by adding a background charge in the solvent region above the slab

Further detail are given in http://dx.doi.org/10.1021/acs.jpcc.8b02465

Parameters:
ne: float
Number of electrons added in the atomic system and (with opposite sign) in the background charge region. At the start it can be an initial guess for the needed number of electrons and will be changed to the current number in the course of the calculation
potential: float
The potential that should be reached or kept in the course of the calculation. If set to “None” (default) a constant charge charge calculation based on the value of \(ne\) is performed.
dpot: float
Tolerance for the deviation of the input \(potential\). If the potential is outside the defined range \(ne\) will be changed in order to get inside again.
doublelayer: dict

Parameters regarding the shape of the counter charge region Implemented keys:

‘start’: float or ‘cavity_like’
If a float is given it corresponds to the lower boundary coordinate (default: z), where the counter charge starts. If ‘cavity_like’ is given the counter charge will take the form of the cavity up to the ‘upper_limit’.
‘thickness’: float
Thickness of the counter charge region in Angstrom. Can only be used if start is not ‘cavity_like’ and will be overwritten by ‘upper_limit’
‘upper_limit’: float
Upper boundary of the counter charge region in terms of coordinate in Anfstrom (default: z). The default is atoms.cell[2][2] - 5.
verbose: bool
Write final electrostatic potential, background charge and and cavity into ASCII file.

Usage Example: A simple Au(111) slab

As a usage example, given here is the calculation of a simple Au slab at a potential of -1 V versus SHE. Keep in mind that the absolute potential has to be provided, where the value of the SHE potential on an absolute scale is approx. 4.4V.

from gpaw.solvation.sjm import SJM, SJMPower12Potential

from ase.build import fcc111
from gpaw import FermiDirac
from ase.io import read

# Import solvation modules
from ase.data.vdw import vdw_radii
from gpaw.solvation import (
    SolvationGPAW,
    EffectivePotentialCavity,
    LinearDielectric,
    GradientSurface,
    SurfaceInteraction)

# Solvent parameters
u0 = 0.180  # eV
epsinf = 78.36  # Dielectric constant of water at 298 K
gamma = 0.00114843767916  # 18.4*1e-3 * Pascal* m
T = 298.15   # K
atomic_radii = lambda atoms: [vdw_radii[n] for n in atoms.numbers]

# Structure is created
atoms = fcc111('Au', size=(1, 1, 3))
atoms.center(axis=2, vacuum=10)
atoms.translate([0, 0, -2])

# SJM parameters
potential = 3.4
ne = 0.1
dpot = 0.01

# The calculator
calc = SJM(doublelayer={'upper_limit': 23},
           potential=potential,
           dpot=dpot,
           ne=ne,
           verbose=True,

           gpts=(16, 16, 136),
           poissonsolver={'dipolelayer': 'xy'},
           kpts=(9, 9, 1),
           xc='PBE',
           txt='Au_pot_%1.2f.txt' % (potential),
           occupations=FermiDirac(0.1),
           cavity=EffectivePotentialCavity(
               effective_potential=SJMPower12Potential(atomic_radii, u0),
               temperature=T,
               surface_calculator=GradientSurface()),
           dielectric=LinearDielectric(epsinf=epsinf),
           interactions=[SurfaceInteraction(surface_tension=gamma)])

# Run the calculation
atoms.set_calculator(calc)
atoms.get_potential_energy()

The output in ‘Au_pot_3.4.txt’ is extended by the grand canonical energy and contains the new part:

----------------------------------------------------------
Grand Potential Energy (Composed of E_tot + E_solv - mu*ne):
Extrpol:    -8.9735918012
Free:    -8.99437808372
-----------------------------------------------------------

These energies are written e.g. into trajectory files.

Since we set the ‘verbose’ keyword to True, the code produced three files:

elstat_potential.out:
Electrostatic potential averaged over xy and referenced to the systems Fermi Level. The outer parts should correspond to the respective work functions.
cavity.out:
The shape of the implicit solvent cavity averaged over xy.
background_charge.out:
The shape of the jellium background charge averaged over x and y.

Usage Example: Running a constant potential NEB calculation

A complete script for performing an NEB calculation can be downloaded here:

import sys
import numpy as np
from ase.parallel import world
from ase.optimize import BFGS
from ase.visualize import view
from gpaw import FermiDirac
from ase.io import write, read
from ase.units import Pascal, m
from ase.neb import NEB

# Import solvation modules
from ase.data.vdw import vdw_radii
from gpaw.solvation import (
    LinearDielectric,
    EffectivePotentialCavity,
    GradientSurface,
    SurfaceInteraction)

# SJM import
from gpaw.solvation.sjm import SJM, SJMPower12Potential

# Solvent parameters from JCP 141, 174108 (2014):
u0 = 0.180  # eV
epsinf = 78.36  # dielectric constant of water at 298 K
gamma = 18.4*1e-3 * Pascal * m  # Surface tension
T = 298.15  # Temperature
atomic_radii = lambda atoms: [vdw_radii[n] for n in atoms.numbers]

# NEB parameters
nimg = 4         # Number of neb images
potential = 4.4  # Desired potential

interpolate = True # True if it is a fresh start with interpolating IS and FS
relax_end_states = True
climb = False

if interpolate:
    initial_file = 'IS_start.traj'
    final_file = 'FS_start.traj'
    ne_IS = -0.4  # First guess of charge on IS
    ne_FS = -0.3  # First guess of charge on FS
    restart_file = None
else:
    restart_file = 'neb_start.traj'


# The calculator
def calculator():
    # Obviously this calculator should be adapted
    return SJM(poissonsolver={'dipolelayer': 'xy'},
               gpts=(48, 32, 168),
               kpts=(4, 6, 1),
               xc='PBE',
               spinpol=False,
               potential=potential,
               occupations=FermiDirac(0.1),
               maxiter=400,
               cavity=EffectivePotentialCavity(
                   effective_potential=SJMPower12Potential(
                       atomic_radii, u0, H2O_layer=True),
                   temperature=T,
                   surface_calculator=GradientSurface()),
               dielectric=LinearDielectric(epsinf=epsinf),
               interactions=[SurfaceInteraction(surface_tension=gamma)],
               )


# Setting up the images
if restart_file:
    images = read(restart_file, index='-%i:' % (nimg+2))
    try:
        # This needs a slight adaptation in ase
        ne_img = [float(image.calc.results['ne']) for image in images]
    except (AttributeError, KeyError):
        # Very bad initial guesses! Should be exchanged by actual values
        ne_img = [i/10. for i in list(range(nimg + 2))]

else:
    initial = read(initial_file)
    final = read(final_file)

    # Shifts atoms in z direction so the lowest layer is equal in all images
    initial.translate([0, 0, -initial.positions[0][2] +
                      final.positions[0][2]])

    images = [initial]
    ne_img = [ne_IS]

    for i in range(nimg):
            images.append(initial.copy())
            ne_img.append(ne_IS + (ne_FS - ne_IS) * (i+1) / float(nimg + 1))

    images.append(final)
    ne_img.append(ne_FS)

# If the endstates should be relaxed in the same run
if relax_end_states:
    if relax_end_states == 'IS':
        endstates = [0]
    elif relax_end_states == 'FS':
        endstates = [-1]
    else:
        endstates = [0, -1]

    system = ['IS', 'FS']
    for i in endstates:
        images[i].set_calculator(calculator())
        images[i].calc.set(txt=system[i]+'.txt')
        images[i].calc.ne = ne_img[i]

        qn = BFGS(images[i], logfile=system[i]+'.log',
                  trajectory=system[i]+'.traj')
        qn.run(fmax=0.03)

        write(system[i]+'_relaxed_pot_%1.2f_ne_%1.5f.traj'
              % (potential, images[i].calc.ne), images[i])
else:
    for i in [0, -1]:
        images[i].calc.ne = ne_img[i]


# Combine NEB images with their respective calculators
for i in range(1, nimg+1):
    images[i].set_calculator(calculator())
    images[i].calc.set(txt='image_%i.txt' % (i))
    images[i].calc.ne = ne_img[i]

# Run the NEB
neb = NEB(images, climb=climb)
if interpolate:
    neb.interpolate()

if world.size == 1:
    view(images)
    sys.exit()
else:
    qn = BFGS(neb, logfile='neb.log', trajectory='neb.traj')
    qn.run(fmax=0.05)


write('neb_final.traj', images)

Note

In this example the keyword ‘H2O_layer = True’ in the ‘SJM_Power12Potential’ class has been used. This keyword frees the interface between the electrode and a water layer from the implicit solvent. It is needed since the rather high distance between the two subsystems would lead to partial solvation of the interface region, therefore screening the electric field in the most interesting area.

References

[1]G. Kastlunger, P. Lindgren, A. A. Peterson, Controlled-Potential Simulation of Elementary Electrochemical Reactions: Proton Discharge on Metal Surfaces, J. Phys. Chem. C 122 (24), 12771 (2018)
[2]G. Henkelman and H. Jonsson, Improved Tangent Estimate in the NEB method for Finding Minimum Energy Paths and Saddle Points, J. Chem. Phys. 113, 9978 (2000)
[3]A. Held and M. Walter, Simplified continuum solvent model with a smooth cavity based on volumetric data, J. Chem. Phys. 141, 174108 (2014).