Solvated Jellium Method (SJM)

Please familiarize yourself with the theory of the SJM approach before starting this tutorial, by reading the SJM documentation page.

A simple Au(111) slab

As a simple example, we’ll examine the calculation of a Au slab with a water overlayer calculated at a potential of -0.2 V versus SHE. To speed the computation, we’ll use a tiny slab with a single water molecule over the top.

To use the solvated jellium approach, instead of importing GPAW as our calculator object, we import SJM. All of the SJM-specific parameters are fed to the calculator as a single dictionary called sj. In our example script, the sj dictionary contains only a single parameter: the target potential; the rest of the SJM parameters are kept at their default values. Keep in mind that the absolute potential has to be provided, where the value of the SHE potential on an absolute scale is approximately 4.4 V.

Since SJM utilizes an implicit solvent to screen the field, we will need to specify some solvent parameters as well. SJM relies on GPAW’s Continuum Solvent Model (CSM) to provide this screening, and any solvation parameters should be fed directly to the calculator, outside of the sj dictionary. (The SJM calculator is a subclass of the SolvationGPAW calculator, so any solvation parameters will be passed straight through to SolvationGPAW.) Reasonable values to simulate room-temperature water are used in the script below. In practice, since the purpose of the implicit solvent is just to screen the field, the net results—such as binding energies and barrier heights—will be relatively insensitive to the choice of implicit solvent parameters, so long as they are kept consistent. Try running the script below to find the number of electrons necessary to equilibrate the work function to 4.2 eV:

from ase.build import fcc111, molecule
from ase.units import Pascal, m

from gpaw.solvation.sjm import SJM, SJMPower12Potential
from gpaw.solvation import (
    EffectivePotentialCavity,
    LinearDielectric,
    GradientSurface,
    SurfaceInteraction)

# Build a tiny gold slab with a single water molecule above.
atoms = fcc111('Au', size=(1, 1, 3))
atoms.center(axis=2, vacuum=12.)
atoms.translate([0., 0., -4.])
water = molecule('H2O')
water.rotate('y', 90.)
water.positions += atoms[2].position + (0., 0., 4.4) - water[0].position
atoms.extend(water)

# Solvated jellium parameters.
sj = {'target_potential': 4.2}  # Desired potential

# Implicit solvent parameters (to SolvationGPAW).
epsinf = 78.36  # dielectric constant of water at 298 K
gamma = 18.4 * 1e-3 * Pascal * m
cavity = EffectivePotentialCavity(
    effective_potential=SJMPower12Potential(H2O_layer=True),
    temperature=298.15,  # K
    surface_calculator=GradientSurface())
dielectric = LinearDielectric(epsinf=epsinf)
interactions = [SurfaceInteraction(surface_tension=gamma)]

# The calculator
calc = SJM(
    # General GPAW parameters.
    mode='fd',
    txt='Au111.txt',
    gpts=(16, 16, 136),
    kpts=(9, 9, 1),
    xc='PBE',
    maxiter=1000,
    # Solvated jellium parameters.
    sj=sj,
    # Implicit solvent parameters.
    cavity=cavity,
    dielectric=dielectric,
    interactions=interactions)
atoms.calc = calc

# Run the calculation.
atoms.get_potential_energy()
atoms.write('Au111.traj')
calc.write_sjm_traces(path='sjm_traces.out')  # .out for .gitignore

If you examine the output in ‘Au111.txt’, you’ll see that the code varies the number of excess electrons until the target potential is within a tolerance of 4.2 V; the default tolerance is 0.01 V. This process usually takes a few steps on the first image, but is faster on subsequent images of a trajectory, since the changes are less dramatic and the potential–vs–charge slope is retained from previous steps. You should see that in net the routine removed only about 0.007 electrons from the simulation, as compared to a charge-neutral simulation, in order to achieve the desired work function.

You’ll notice that the output in ‘Au.txt’ contains additional information, as compared to a charge-neutral simulation:

Legendre-transformed energies (Omega = E - N mu)
  (grand-potential energies)
  N (excess electrons):   +0.006663
  mu (workfunction, eV):   +4.194593
--------------------------
Free energy:    -23.630651
Extrapolated:   -23.608083

These Legendre-transformed energies, \(\Omega = E - N \mu\), are written into any ASE trajectories created, and are the quantity returned by calc.get_potential_energy() (and therefore atoms.get_potential_energy()). As discussed in Legendre-transformed energy, these grand-potential energies are consistent with the forces in the grand-canonical scheme, and are thus compatible with methods such as saddle-point searches (e.g., NEB) and energy optimizations (e.g., BFGSLineSearch). If you’d rather have it output the traditional canonical energies to the get_potential_energy methods, you can use the sj['grand_output'] keyword.

In the last line of the script above, we wrote out some additional information in the form of ‘sjm traces’. These are traces of the electrostatic potential, the background charge, and the solvent cavity across the \(z\) axis (that is, \(xy\)-averaged). It’s a good idea to take a look at these, and perhaps make a plot for each, when running a simulation with a new system so you can be sure it is behaving as you expect. The three files that are created are:

potential.txt:

Electrostatic potential averaged over \(xy\) and referenced to the system’s Fermi Level. The outer parts should correspond to the respective work functions and the potential at the top boundary should correspond to the potential set in the input.

cavity.txt:

The shape of the implicit solvent cavity averaged over \(xy\). Multiplying the cavity function with epsinf corresponds to the \(xy\)-averaged dielectric constant distribution in the unit cell.

background_charge.txt:

The shape of the jellium background charge averaged over \(xy\) and normalized.

Note

Alternatively, calc.write_sjm_traces(style='cube') can be used to write cube files on the 3-D grid, instead of traces.

Our traces showing the solvent region and the jellium region are plotted below. The black dots indiate the atom’s \(z\) positions. We see that both the solvent and the jellium are starting above the highest atom, so all is good. If you do see solvent appearing in a region where you do not want it, you can block it with ghost atoms or a boundary plane, as described in the SJMPower12Potential docstring. If you see solvent in the metal region, you could increase the atomic radii as described in the same docstring.

../../../_images/traces-Au111.png

(The script used to generate this plot is here: plot-traces.py.)

Relaxing the slab and adsorbate

Next, let’s perform a structural optimization at a constant potential. Use the same system as the previous example, but add an H atom to the surface as an adsorbate.

from ase import Atom
import ase.io
from ase.units import Pascal, m
from ase.optimize import BFGS
from ase.constraints import FixAtoms

from gpaw.solvation.sjm import SJM, SJMPower12Potential
from gpaw.solvation import (
    EffectivePotentialCavity,
    LinearDielectric,
    GradientSurface,
    SurfaceInteraction)

# Add an H adsorbate.
atoms = ase.io.read('Au111.traj')
atoms.append(Atom('H', atoms[2].position + (0., 0., 1.5)))

# Fix some atoms.
atoms.set_constraint(FixAtoms(indices=[0, 1]))

# Solvated jellium parameters.
sj = {'target_potential': 4.2}

# Implicit solvent parameters (to SolvationGPAW).
epsinf = 78.36  # dielectric constant of water at 298 K
gamma = 18.4 * 1e-3 * Pascal * m
cavity = EffectivePotentialCavity(
    effective_potential=SJMPower12Potential(H2O_layer=True),
    temperature=298.15,  # K
    surface_calculator=GradientSurface())
dielectric = LinearDielectric(epsinf=epsinf)
interactions = [SurfaceInteraction(surface_tension=gamma)]

# The calculator
calc = SJM(
    # General GPAW parameters.
    mode='fd',
    txt='Au111-H-seq.txt',
    gpts=(16, 16, 136),
    kpts=(9, 9, 1),
    xc='PBE',
    maxiter=1000,
    # Solvated jellium parameters.
    sj=sj,
    # Implicit solvent parameters.
    cavity=cavity,
    dielectric=dielectric,
    interactions=interactions)
atoms.calc = calc

# Run the calculation.
opt = BFGS(atoms, trajectory='qn-Au111-H-seq.traj',
           logfile='qn-Au111-H-seq.log')
opt.run()

If you examine the output, you should see this behaves the same as any other relaxation job, except the potential is constant, rather than the number of electrons.

A faster route: simultaneous potential and structure optimization

It is often faster to optimize the structure and the potential simultaneously. That is, instead of waiting until the potential is perfectly equilibrated before taking an ionic step, you can take an ionic step and adjust the number of electrons simultaneously. To do this, set a loose potential tolerance and turn the always_adjust keyword to True. An example follows.

from ase import Atom
import ase.io
from ase.units import Pascal, m
from ase.optimize import BFGS
from ase.constraints import FixAtoms

from gpaw.solvation.sjm import SJM, SJMPower12Potential
from gpaw.solvation import (
    EffectivePotentialCavity,
    LinearDielectric,
    GradientSurface,
    SurfaceInteraction)

# Add an H adsorbate.
atoms = ase.io.read('Au111.traj')
atoms.append(Atom('H', atoms[2].position + (0., 0., 1.5)))

# Fix some atoms.
atoms.set_constraint(FixAtoms(indices=[0, 1]))

# Solvated jellium parameters.
sj = {'target_potential': 4.2,
      'tol': 0.2,
      'always_adjust': True}

# Implicit solvent parameters (to SolvationGPAW).
epsinf = 78.36  # dielectric constant of water at 298 K
gamma = 18.4 * 1e-3 * Pascal * m
cavity = EffectivePotentialCavity(
    effective_potential=SJMPower12Potential(H2O_layer=True),
    temperature=298.15,  # K
    surface_calculator=GradientSurface())
dielectric = LinearDielectric(epsinf=epsinf)
interactions = [SurfaceInteraction(surface_tension=gamma)]

# The calculator
calc = SJM(
    # General GPAW parameters.
    mode='fd',
    txt='Au111-H-sim.txt',
    gpts=(16, 16, 136),
    kpts=(9, 9, 1),
    xc='PBE',
    maxiter=1000,
    # Solvated jellium parameters.
    sj=sj,
    # Implicit solvent parameters.
    cavity=cavity,
    dielectric=dielectric,
    interactions=interactions)
atoms.calc = calc

# Run the calculation.
opt = BFGS(atoms, trajectory='qn-Au111-H-sim.traj',
           logfile='qn-Au111-H-sim.log')
opt.run()

# Tighten the tolerances again.
sj['tol'] = 0.01
sj['always_adjust'] = False
sj['slope'] = None
calc.set(sj=sj)
opt = BFGS(atoms, trajectory='qn-Au111-H-sim-1.traj',
           logfile='qn-Au111-H-sim.log')
opt.run()

Some things to note:

  • At the end of the script, we tightened the criteria back up, in order to guarantee precise results consistent with a normal, sequential calculation (assuming it found the same local minimum).

  • It’s a good idea to use the BFGS, and not the BFGSLineSearch, algorithm for structural optimization. This is because the forces and energy are not necessarily consistent until the optimization finishes, and that might confuse an optimizer that uses the energy. BFGSLineSearch may in general have trouble with SJM, since the tolerance set on the desired potential can lead to small inconsistencies between the forces and energy. When in doubt, use the traditional BFGS optimizer.

The plot below compares the “sequential” and “simultaneuous” optimization approaches. The open circles are potential-equilibration steps, while the filled circles represent ionic steps (in which the atoms are moved). In sequential mode, it always makes sure the potential is within tolerance before taking an ionic step, and you can see that it typically takes an extra DFT calculation (or two) before moving the atoms. In simultaneous mode, it moves the atoms with much higher frequency, and only when the potential gets very far from the target does it pause to equilibrate. Generally, this results in faster total convergence. At the end of the sequential mode, we switched back to simultaneous mode to ensure that the final result was at the expected potential; this adds a few extra steps, but overall the computational savings are typically worth it.

../../../_images/simultaneous.png

(Plot script: plot-seq-sim.py.)

Finding a barrier (NEB/DyNEB)

Perhaps the most-desired use of constant-potential calculations is to find barriers at a specified potential. This avoids a well-known problem in canonical electronic structure calculations, where the work function (that is, the electrical potential) can change by 1–2 V over the course of an elementary step. Such a potential change is obviously much greater than the experimental situation, where potentiostats hold the potential to a tolerance many orders of magnitude smaller. The SJM method allows one to calculate a reaction barrier in a manner where all images in the trajectory have identical work functions, down to a user-specified tolerance.

The discrepancy between a canonical and grand-canonical simulation is most pronounced in reactions involving the creation or destruction of an ion, as the ion collides with the constant-potential surface. Unfortunately, simulating even a simple reaction like the Volmer reaction—where a proton in the solution reacts to form an adsorbed hydrogen atom—is not possible in this tutorial, as a relatively large water layer—that is, larger than the one molecule we are using!—is needed in order to stabilize the solvated proton in the initial state. Therefore, in this tutorial, we’ll do a much simpler reaction: the diffusion of our adsorbed H atom from one type of surface site to another. Note that we should expect very little difference in the number of excess electrons in order to equilibrate the initial, transition, and final state; nevertheless, it allows us to demonstrate the use of the nudged elastic band (NEB) method with SJM.

You will need to create the final state before you run the NEB; you can go through the exercise of creating this yourself, or you can use our script (run-Au111-H-hollow.py). After you have created the final state, you can run the NEB as below; note here we use the DyNEB method which tends to save a large amount of computational time for these sorts of barriers.

import ase.io
from ase.units import Pascal, m
from ase.optimize import BFGS
from ase.mep import interpolate, DyNEB

from gpaw.solvation.sjm import SJM, SJMPower12Potential
from gpaw.solvation import (
    EffectivePotentialCavity,
    LinearDielectric,
    GradientSurface,
    SurfaceInteraction)


def make_calculator(index):
    # Solvated jellium parameters.
    sj = {'target_potential': 4.2,
          'tol': 0.01,
          'excess_electrons': -0.01232,  # guess from previous
          'slope': -50.}  # guess from previous
    # Implicit solvent parameters (to SolvationGPAW).
    epsinf = 78.36  # dielectric constant of water at 298 K
    gamma = 18.4 * 1e-3 * Pascal * m
    cavity = EffectivePotentialCavity(
        effective_potential=SJMPower12Potential(H2O_layer=True),
        temperature=298.15,  # K
        surface_calculator=GradientSurface())
    dielectric = LinearDielectric(epsinf=epsinf)
    interactions = [SurfaceInteraction(surface_tension=gamma)]
    calc = SJM(
        # General GPAW parameters.
        mode='fd',
        txt=f'gpaw-{index:d}.txt',
        gpts=(16, 16, 136),
        kpts=(9, 9, 1),
        xc='PBE',
        maxiter=1000,
        # Solvated jellium parameters.
        sj=sj,
        # Implicit solvent parameters.
        cavity=cavity,
        dielectric=dielectric,
        interactions=interactions)
    return calc


# Create the band of images, attaching a calc to each.
initial = ase.io.read('qn-Au111-H-sim-1.traj')
final = ase.io.read('qn-Au111-H-hollow-1.traj')
images = [initial]
for index in range(5):
    images += [initial.copy()]
    images[-1].calc = make_calculator(index + 1)
images += [final]
interpolate(images)

# Create and relax the DyNEB.
neb = DyNEB(images)
opt = BFGS(neb, logfile='neb.log', trajectory='neb.traj')
opt.run(fmax=0.05)

After a few iterations, the band completes and we can make a plot of the barrier, using NEBTools. (Plot script: plot-neb.py.)

../../../_images/band.png

Although the reaction modeled here is very simple and doesn’t involve much need for a constant-potential method, it nonetheless allows us to demonstrate that SJM plays well with standard barrier-search methods, like the NEB.

You’ll note that we used DyNEB instead of the regular NEB. Although NEB will also work fine, DyNEB is more efficient as it skips re-calculating individual images whose maximum force is already below the cutoff (of fmax=0.05). Note also that parallelizing a NEB over images is in general not a great strategy when using a grand-canonical method, because typically only one of the images might need potential equilibration on any individual force call of the NEB—and in parallel-image mode all of the images would need to wait for that individual image. However, running in serial allows us to take advantage of the computational efficiency of the DyNEB method, which is not possible when parallelizing over images.

Constant-charge mode

The SJM code can also run in constant-charge mode, where the user specifiies the total number of electrons in the simulation. This can be a fast way to calculate a system at several potentials—that is, if one does not care about the specific potentials, but just wants to span a range of potentials. To use constant-charge mode, just specify the number of excess_electrons in your sj dict and calculate as normal:

sj = {'excess_electrons': 0.5}
calc = SJM(sj=sj, ...)

If your intent is to run in the constant-charge ensemble—like in the use of the charge-extrapolation scheme—you may want to output the canonical, rather than the grand-potential energies. The canonical energies are consistent with the forces in constant-charge mode. To accomplish this, set sj['grand_output'] = False, like:

sj = {'excess_electrons': ...,
      'grand_output': False} :
calc = SJM(sj=sj, ...)