index | modules | gitlab | page source

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

## 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!