*O referencing issues in the tutorial

[zulissimeta]

The demo shows a large offset in the energy of *O from John's paper to the OCP predictions:
https://open-catalyst-project.github.io/tutorial/OCP-introduction.html
John explained this during the tutorial as PBE-RPBE disagreement, but this seems very large for such a simple adsorbate.

Double checking the energies with the OCP demo shows O Pt(111) with an energy of about 1.42 eV, compared to ~0.59 eV from the calculator in the tutorial. Using 1.42 eV (with the formation reaction energies for H2O/O) yields final (O+->*) energy of -4.46 eV, which is pretty close to John's number (-4.26 eV):
https://open-catalyst.metademolab.com/results/01b2250b-f9c4-43d8-b9b7-f987988db7e8

I'm not exactly sure what happened here, but I would guess the missing tags have something to do this with. My understanding is that our models generally sum over just the surface atoms and adsorbate, so the energy in the tutorial might correspond to summing over the entire slab.

A script to automatically set tags to reasonable values might be helpful here and elsewhere?

[abhshkdz]

My understanding is that our models generally sum over just the surface atoms and adsorbate

For GemNet-OC, the summation is over all atoms in the structure (code), but the contribution for atoms tagged as surface or adsorbate is computed differently compared to subsurface atoms, and that could be causing the difference here.

Alternatively, try out EquiformerV2 where you don't have to worry about setting tags! The 31M model should be easy to fit in memory.

[abhshkdz]

Update: I re-ran John's notebook with EquiformerV2 31M. That gets ~0.70 eV, yielding a final $O + * \rightarrow O*$ energy of ~-4.90 eV — still a wide gap

[jkitchin]

That is very close to the result from the eSCN checkpoint in the tutorial (-4.88 eV). There are a few places where 0.5 eV discrepancy could show up I think, mostly related to the O2 reference energy in the paper referenced at the top. In the tutorial all the reaction energy references are from experimental data.

[jkitchin]

Something is weird, If I use the Equiformer checkpoint, I get an energy similar to what Abishek reports above, but 0.76 eV that is not close to the OCP demo energy. I thought those would be the same?

Edit: My energy is not the same as Abishek's because I used the exp lattice constant, not the coordinates from the json file. I get the same number as him when I use those.

[zulissimeta]

Two things that seem to matter:

• Relaxing the entire slab (not constraining the bottom layers) will lower the energy significantly
• eSCN returns very different results for a 2x2xN slab vs anything larger (e.g. 3x3xN like in the OCP tutorial). This is perhaps because all of our training data is for very low coverage. It would be interesting to see if total energy models worked here!

Output:

2x2x2 Pt(111) slab: 0.966696
2x2x3 Pt(111) slab: 0.754219
2x2x4 Pt(111) slab: 0.918704
2x2x5 Pt(111) slab: 0.924884
3x3x2 Pt(111) slab: 1.670233
3x3x3 Pt(111) slab: 1.634536
3x3x4 Pt(111) slab: 1.678548
3x3x5 Pt(111) slab: 1.647686
4x4x2 Pt(111) slab: 1.633825
4x4x3 Pt(111) slab: 1.928250
4x4x4 Pt(111) slab: 1.758899
4x4x5 Pt(111) slab: 1.738886

Script:

import numpy as np
from ase.build import add_adsorbate, fcc111
from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.db import connect
from ase.optimize import BFGS

for n_repeat in [2, 3, 4]:
    for n_layers in [2, 3, 4, 5]:

        # Make slab with default lattice constant
        atoms = fcc111("Pt", (n_repeat, n_repeat, n_layers))
        atoms.cell[:, 2] = [0, 0, 20]
        atoms.set_pbc(True)

        # Set constraints so that only the top layer and adsorbate is free
        atoms.set_constraint(
            FixAtoms(
                indices=np.where(
                    [
                        atom
                        for atom in atoms
                        if atom.position[2] < atoms[-2].position[2] - 1
                    ]
                )[0]
            )
        )

        # Add adsorbate
        add_adsorbate(atoms, "O", height=1.0, position="fcc")

        # Set tags, not important if not gemnet-OC
        atoms.set_tags(np.ones(len(atoms)))
        atoms.set_calculator(calc)

        # Run relaxation
        opt = BFGS(atoms, logfile=None)
        opt.run(fmax=0.05)

        # Print output
        print(
            f"{n_repeat}x{n_repeat}x{n_layers} Pt(111) slab: {atoms.get_potential_energy() :2f}"
        )

[abhshkdz]

eSCN returns very different results for a 2x2xN slab vs anything larger (e.g. 3x3xN like in the OCP tutorial).

Yeah, it might just be this causing the discrepancy.

I was able to reproduce the demo results to get ~1.42 eV, yielding a final $O+* \rightarrow *$ energy of ~-4.19 eV with the following:

from ocdata.core import Bulk, Slab, Adsorbate, AdsorbateSlabConfig

bulk = Bulk(bulk_src_id_from_db="mp-126")
slabs = Slab.from_bulk_get_specific_millers(bulk=bulk, specific_millers=(1,1,1))
adsorbate = Adsorbate(adsorbate_smiles_from_db="*O")
heuristic_adslabs = AdsorbateSlabConfig(slabs[0], adsorbate, mode="heuristic")

atoms = heuristic_adslabs.atoms_list[0]
atoms.set_calculator(calc)
opt = BFGS(atoms)

opt.run(fmax=0.05, steps=100)
re2 = atoms.get_potential_energy()
    
nO = 0
for atom in atoms:
    if atom.symbol == 'O':
        nO += 1
        re2 += re1 + re3

print(re2 / nO)

      Step     Time          Energy         fmax
BFGS:    0 13:43:03        2.238281        1.8482
BFGS:    1 13:43:03        2.150391        1.4790
BFGS:    2 13:43:03        2.011719        1.0960
BFGS:    3 13:43:04        1.972656        0.5262
BFGS:    4 13:43:04        1.957031        0.3768
BFGS:    5 13:43:04        1.931641        0.3193
BFGS:    6 13:43:04        1.929688        0.3229
BFGS:    7 13:43:04        1.925781        0.3506
BFGS:    8 13:43:04        1.919922        0.4213
BFGS:    9 13:43:04        1.896484        0.6744
BFGS:   10 13:43:04        1.869141        0.9231
BFGS:   11 13:43:05        1.796875        0.9111
BFGS:   12 13:43:05        1.769531        1.2941
BFGS:   13 13:43:05        1.728516        1.3930
BFGS:   14 13:43:05        1.660156        0.7369
BFGS:   15 13:43:05        1.583984        0.5790
BFGS:   16 13:43:05        1.519531        0.4506
BFGS:   17 13:43:05        1.478516        0.3661
BFGS:   18 13:43:05        1.455078        0.2392
BFGS:   19 13:43:05        1.451172        0.2166
BFGS:   20 13:43:06        1.449219        0.1618
BFGS:   21 13:43:06        1.447266        0.1514
BFGS:   22 13:43:06        1.441406        0.1439
BFGS:   23 13:43:06        1.435547        0.1462
BFGS:   24 13:43:06        1.433594        0.1310
BFGS:   25 13:43:06        1.425781        0.1083
BFGS:   26 13:43:06        1.421875        0.0974
BFGS:   27 13:43:06        1.423828        0.0850
BFGS:   28 13:43:06        1.423828        0.0890
BFGS:   29 13:43:07        1.421875        0.0971
BFGS:   30 13:43:07        1.421875        0.1007
BFGS:   31 13:43:07        1.419922        0.0878
BFGS:   32 13:43:07        1.419922        0.1200
BFGS:   33 13:43:07        1.421875        0.1250
BFGS:   34 13:43:07        1.419922        0.0783
BFGS:   35 13:43:07        1.421875        0.0339
-4.188124999999999

Initial structure:

Final structure:

I'm not a 100% I set everything up properly above — would be good to double check!