Surface adsorption study using the ASE database

In this tutorial we will adsorb C, N and O on 7 different FCC(111) surfaces with 1, 2 and 3 layers and we will use database files to store the results.

See also

The ase.db module documentation.

cu1o cu2o cu3o

Bulk

First, we calculate the equilibrium bulk FCC lattice constants for the seven elements where the EMT potential works well:

from ase.build import bulk
from ase.calculators.emt import EMT
from ase.eos import calculate_eos
from ase.db import connect

db = connect('bulk.db')
for symb in ['Al', 'Ni', 'Cu', 'Pd', 'Ag', 'Pt', 'Au']:
    atoms = bulk(symb, 'fcc')
    atoms.calc = EMT()
    eos = calculate_eos(atoms)
    v, e, B = eos.fit()  # find minimum
    # Do one more calculation at the minimu and write to database:
    atoms.cell *= (v / atoms.get_volume())**(1 / 3)
    atoms.get_potential_energy()
    db.write(atoms, bm=B)

Run the bulk.py script and look at the results:

$ python3 bulk.py
$ ase db bulk.db -c +bm  # show also the bulk-modulus column
id|age|formula|calculator|energy| fmax|pbc|volume|charge|   mass|   bm
 1|10s|Al     |emt       |-0.005|0.000|TTT|15.932| 0.000| 26.982|0.249
 2| 9s|Ni     |emt       |-0.013|0.000|TTT|10.601| 0.000| 58.693|1.105
 3| 9s|Cu     |emt       |-0.007|0.000|TTT|11.565| 0.000| 63.546|0.839
 4| 9s|Pd     |emt       |-0.000|0.000|TTT|14.588| 0.000|106.420|1.118
 5| 9s|Ag     |emt       |-0.000|0.000|TTT|16.775| 0.000|107.868|0.625
 6| 9s|Pt     |emt       |-0.000|0.000|TTT|15.080| 0.000|195.084|1.736
 7| 9s|Au     |emt       |-0.000|0.000|TTT|16.684| 0.000|196.967|1.085
Rows: 7
Keys: bm
$ ase gui bulk.db

The bulk.db is an SQLite3 database in a single file:

$ file bulk.db
bulk.db: SQLite 3.x database

If you want to see what’s inside you can convert the database file to a json file and open that in your text editor:

$ ase db bulk.db --insert-into bulk.json
Added 0 key-value pairs (0 pairs updated)
Inserted 7 rows

or, you can look at a single row like this:

$ ase db bulk.db Cu -j
{"1": {
 "calculator": "emt",
 "energy": -0.007036492048371201,
 "forces": [[0.0, 0.0, 0.0]],
 "key_value_pairs": {"bm": 0.8392875566787444},
 ...
 ...
}

The json file format is human readable, but much less efficient to work with compared to a SQLite3 file.

Adsorbates

Now we do the adsorption calculations (run the ads.py script).

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

db1 = connect('bulk.db')
db2 = connect('ads.db')


def run(symb, a, n, ads):
    atoms = fcc111(symb, (1, 1, n), a=a)
    add_adsorbate(atoms, ads, height=1.0, position='fcc')

    # Constrain all atoms except the adsorbate:
    fixed = list(range(len(atoms) - 1))
    atoms.constraints = [FixAtoms(indices=fixed)]

    atoms.calc = EMT()
    opt = BFGS(atoms, logfile=None)
    opt.run(fmax=0.01)
    return atoms


for row in db1.select():
    a = row.cell[0, 1] * 2
    symb = row.symbols[0]
    for n in [1, 2, 3]:
        for ads in 'CNO':
            atoms = run(symb, a, n, ads)
            db2.write(atoms, layers=n, surf=symb, ads=ads)

We now have a new database file with 63 rows:

$ ase db ads.db -n
63 rows

These 63 calculations only take a few seconds with EMT. Suppose you want to use DFT and send the calculations to a supercomputer. In that case you may want to run several calculations in different jobs on the computer. In addition, some of the jobs could time out and not finish. It’s a good idea to modify the script a bit for this scenario. We add a couple of lines to the inner loop:

for row in db1.select():
    a = row.cell[0, 1] * 2
    symb = row.symbols[0]
    for n in [1, 2, 3]:
        for ads in 'CNO':
            id = db2.reserve(layers=n, surf=symb, ads=ads)
            if id is not None:
                atoms = run(symb, a, n, ads)
                db2.write(atoms, layers=n, surf=symb, ads=ads)
                del db2[id]

The reserve() method will check if there is a row with the keys layers=n, surf=symb and ads=ads. If there is, then the calculation will be skipped. If there is not, then an empty row with those keys-values will be written and the calculation will start. When done, the real row will be written and the empty one will be removed. This modified script can run in several jobs all running in parallel and no calculation will be done twice.

In case a calculation crashes, you will see empty rows left in the database:

$ ase db ads.db natoms=0 -c ++
id|age|user |formula|pbc|charge| mass|ads|layers|surf
17|31s|jensj|       |FFF| 0.000|0.000|  N|     1|  Cu
Rows: 1
Keys: ads, layers, surf

Delete them, fix the problem and run the script again:

$ ase db ads.db natoms=0 --delete
Delete 1 row? (yes/No): yes
Deleted 1 row
$ python ads.py  # or sbatch ...
$ ase db ads.db natoms=0
Rows: 0

Reference energies

Let’s also calculate the energy of the clean surfaces and the isolated adsorbates (refs.py):

from ase import Atoms
from ase.calculators.emt import EMT
from ase.db import connect
from ase.build import fcc111

db1 = connect('bulk.db')
db2 = connect('ads.db')


def run(symb, a, n):
    atoms = fcc111(symb, (1, 1, n), a=a)
    atoms.calc = EMT()
    atoms.get_forces()
    return atoms


# Clean slabs:
for row in db1.select():
    a = row.cell[0, 1] * 2
    symb = row.symbols[0]
    for n in [1, 2, 3]:
        id = db2.reserve(layers=n, surf=symb, ads='clean')
        if id is not None:
            atoms = run(symb, a, n)
            db2.write(atoms, id=id, layers=n, surf=symb, ads='clean')

# Atoms:
for ads in 'CNO':
    a = Atoms(ads)
    a.calc = EMT()
    a.get_potential_energy()
    db2.write(a)

$ python refs.py
$ ase db ads.db -n
87 rows

Say we want those 24 reference energies (clean surfaces and isolated adsorbates) in a refs.db file instead of the big ads.db file. We could change the refs.py script and run the calculations again, but we can also manipulate the files using the ase db tool. First, we move over the clean surfaces:

$ ase db ads.db ads=clean --insert-into refs.db
Added 0 key-value pairs (0 pairs updated)
Inserted 21 rows
$ ase db ads.db ads=clean --delete --yes
Deleted 21 rows

and then the three atoms (pbc=FFF, no periodicity):

$ ase db ads.db pbc=FFF --insert-into refs.db
Added 0 key-value pairs (0 pairs updated)
Inserted 3 rows
$ ase db ads.db pbc=FFF --delete --yes
Deleted 3 rows
$ ase db ads.db -n
63 rows
$ ase db refs.db -n
24 rows

Analysis

Now we have what we need to calculate the adsorption energies and heights (ea.py):

from ase.db import connect

refs = connect('refs.db')
db = connect('ads.db')

for row in db.select():
    ea = (row.energy -
          refs.get(formula=row.ads).energy -
          refs.get(layers=row.layers, surf=row.surf).energy)
    h = row.positions[-1, 2] - row.positions[-2, 2]
    db.update(row.id, height=h, ea=ea)

Here are the results for three layers of Pt:

$ python3 ea.py
$ ase db ads.db Pt,layers=3 -c formula,ea,height
formula|    ea|height
Pt3C   |-3.715| 1.504
Pt3N   |-5.419| 1.534
Pt3O   |-4.724| 1.706
Rows: 3
Keys: ads, ea, height, layers, surf

Note

While the EMT description of Ni, Cu, Pd, Ag, Pt, Au and Al is OK, the parameters for C, N and O are not intended for real work!