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.db import connect
from ase.eos import calculate_eos
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.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
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.build import fcc111
from ase.calculators.emt import EMT
from ase.db import connect
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!