Solvated Jellium Method (SJM)

Theoretical Background

The Solvated Jellium method (SJM) is a 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 unit cell charge neutral.

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 and no electric field is present.

\[\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\]

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

# Import solvation modules
from ase.data.vdw import vdw_radii
from gpaw.solvation import (
    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


def atomic_radii(atoms):
    return [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 (E_tot + E_solv - mu*ne):
Extrpol:    -9.004478833351179
Free:    -9.025229571371211
-----------------------------------------------------------

These energies are written e.g. into trajectory files if write_grandcanonical_energy = True.

Since we set verbose = 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.

Note

Alternatively, verbose = 'cube' corresponds to True plus creation of a cube file including the dielectric function (cavity) on the 3-D grid.

Structure optimization

Any kind of constant potential structure optimization can be performed by applying ase’s built-in optimizers. Two options are given for the potential equilibration during structure optimization:

potential_equilibration_mode = ‘seq’:

Sequential mode. This is the default optimization mode, which fully equilibrates the potential after each ionic step, if the current potential differs from the target by more than dpot. This mode is generally slower (up to 1.5 time CPU hours compared to constant charge GPAW), but very reliable.

potential_equilibration_mode = ‘sim’:

Simultaneous mode. In this mode ne is constantly optimized together with the geometry. It is generally reliable (with few exceptions) and does not lead to a significant performance penalty compared to constant charge relaxations. However, after convergence it is adviced to check the final potential, since there is no guarantee for it to be within dpot. During the optimization the potential can oscillate, which mostly calms down close to convergence. One can, however, control the oscillation via max_pot_deviation. This keyword automatically triggers a tight and complete potential equilibration to the target, if the current potential is outside the given threshold.

Usage Example: Running a constant potential NEB calculation

A complete automatized script for performing a NEB calculation can be downloaded here: run_SJM_NEB.py. It can, of course, be substituted by a simpler, more manual script as can be found in the NEB tutorial

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.

Note

For paralle NEB runs potential_equilibration_mode = 'sim' should be used for efficiency, since not every image triggers a potential equilibration step in sequential mode. Such a run would lead to unnecessary idle time on some cores.

Note

For CI-NEBs (with climb = True), we advice to either set max_pot_deviation to a tighter value (e.g. 0.05) in simultaneous mode or use the sequential mode.

class gpaw.solvation.sjm.SJM(ne=0, doublelayer=None, potential=None, write_grandcanonical_energy=True, dpot=0.01, potential_equilibration_mode='seq', tiny=1e-08, max_pot_deviation=0.2, 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.

potential_equilibration_mode: str

The mode for potential relaxation (only relevant for structure optimizations):

‘seq’ (Default):

Sequential mode. In a geometry optimization, the potential will be fully equilibrated after each ionic step, if outside of dpot.

‘sim’:

Simultaneous mode. After each ionic step only ne is adjusted and the geometry optimization continues. If potential deviation from target is higher than \(max_pot_deviation\) one complete potential equilibration at constant geometry is triggered.

dpot: float

Tolerance for the deviation of the target potential in sequential mode. If the potential is outside the defined range \(ne\) will be changed in order to get inside again. Default: 0.01V.

max_pot_deviation: float

Maximum tolerance for the deviation in the potential from target potential in simultaneous mode. If the tolerance is surpassed a complete potential equilibration is triggered. Default: 0.2 V.

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 Angstrom (default: z). The default is atoms.cell[2][2] - 5.

verbose: bool or ‘cube’
True:

Write final electrostatic potential, background charge and and cavity into ASCII files.

‘cube’:

In addition to ‘True’, also write the cavity on the 3D-grid into a cube file.

write_grandcanonical_energy: bool

Write the constant potential energy into output files such as trajectory files. Default: True

Constructor for SolvationGPAW class.

Additional arguments not present in GPAW class: cavity – A Cavity instance. dielectric – A Dielectric instance. interactions – A list of Interaction instances.

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).