Add catalysis project

parent 1299c811
.. _projects:
.. _dev projects:
Ongoing Projects
from myqueue.task import task
def create_tasks():
return [task('', tmax='5h'),
if __name__ == '__main__':
from pathlib import Path
Lattice="5.408 0.0 0.0 2.704 4.68346538367 0.0 0.0 0.0 12.141568" Properties=species:S:1:pos:R:3:magmoms:R:1:tags:I:1:Z:I:1:forces:R:3 pbc="T T T" energy=-81.3838458737
Ru 0.00000000 1.56115513 5.00000000 0.00000000 0 44 0.00000000 -0.06723168 0.42809835
Ru 2.70400000 1.56115513 5.00000000 0.00000000 0 44 0.00000000 -0.15085339 -0.05304611
Ru 1.35200000 3.90288782 5.00000000 0.00000000 0 44 0.02993366 0.10911374 -0.17238154
Ru 4.05600000 3.90288782 5.00000000 0.00000000 0 44 -0.02993366 0.10911374 -0.17238154
Ru 0.00000000 -0.05001999 7.01356933 0.00000000 0 44 0.00000000 0.00162123 0.00150827
Ru 2.70400000 -0.20410314 7.17709544 0.00000000 0 44 0.00000000 -0.00874450 -0.00456742
Ru 1.33105146 2.25767274 7.10060851 0.00000000 0 44 0.00200376 -0.00126831 0.00593867
Ru 4.07694854 2.25767274 7.10060851 0.00000000 0 44 -0.00200376 -0.00126831 0.00593867
N 1.78223533 1.05565589 8.46377451 0.00000000 0 7 0.00632433 0.00230531 -0.00318617
N 3.62576467 1.05565589 8.46377451 0.00000000 0 7 -0.00632433 0.00230531 -0.00318617
:download:`neb.ipynb`, :download:``
:download:`vibrations.ipynb`, :download:``
:download:`convergence.ipynb`, :download:`convergence.json`,
This exercise studies the splitting of the |N2| molecule on a Ruthenium
surface. |N2| splitting is the critical step in ammonia synthesis, which is
the main source of biologically accessible nitrogen for fertilizers.
Note that the |N2| splitting occurs most readily at the bottom of a step on
the close-packed (0001) surface. However, to keep system sizes and computer
time down at a manageable level, we shall look at a flat surface.
Tools used in this exercise:
* Structural energy minimization.
* Nudged Elastic Band (NEB) for finding transition states.
* If you have time: Extra exercise on vibrational energy
Part 1: |N2| adsorption on a flat Ru surface
The notebook N2_on_metal.ipynb shows how to set up a molecule on a flat metal surface.
* Set up a clean metal surface.
* Relax the topmost layer (ca. 10 min running time).
- While running: study the gpaw text output, to learn about e.g. number of irreducible k-points (important for parallel simulations).
* Set up and relax a single |N2| molecule (ca. 1 min running time).
* Add the molecule standing on the metal on an on-top site.
* Relax the combined system.
Part 2: Splitting |N2|: initial and final geometry
(Begin this while the last step above runs)
The |N2| molecule will not split while standing up on an on-top site. The molecule can also bind to the surface in a flat geometry - we here ignore the barrier between the two states and just use the lying-down molecule as the initial configuration.
Create scripts setting up and energy-minimizing the initial and final structures, as described in the final part of the notebook from Part 1. Submit these scripts as parallel batch jobs. When submitting make sure that the number of CPU cores matches the number of irreducible k-point in the calculation, as k-point parallelization is much more efficient than other forms of parallelization.
XXXX Write something about dry-run and job submission here.
Part 3: Learning about Nudged Elastic Band
While the calculations from the previous step runs, you can learn about the Nudged Elastic Band method for finding transition states and barriers from the notebook NEB.ipynb
Part 4: Run a parallel NEB calculation
Prepare a script running NEB using the GPAW calculator and the initial and final states from part 2 to find the barrier for |N2| dissociation.
When doing this you should parallelize over the images in the NEB calculation. A more detailed description of how to do this can be found in the
*Exercise* part of the NEB.ipynb along with some suitable input parameters for the NEB.
Extra exercise: Vibrational energy
The notebook Vibrations.ipynb will guide you through how to calculate the vibrations of the adsorbate in the inital and final state and use the Thermochamistry module to calculate the reaction free energy.
The final part of the exercise shows what happens when you calculate the vibrations of a well-converged NEB transition state.
Extra material: Convergence test
We look at the adsorption energy and height of a nitrogen atom on a Ru(0001)
surface in the hcp site. We check for convergence with respect to:
* number of layers
* number of k-points in the BZ
* plane-wave cutoff energy
.. |N2| replace:: N\ :sub:`2`
import numpy as np
from ase import Atoms
from import hcp0001
from ase.constraints import FixAtoms
from ase.optimize import BFGSLineSearch
from gpaw import GPAW, PW, Davidson
a = 2.72
c = 1.58 * a
vacuum = 4.0
def adsorb(db, height=1.2, nlayers=3, nkpts=7, ecut=400):
"""Adsorb nitrogen in hcp-site on Ru(0001) surface.
Do calculations for N/Ru(0001), Ru(0001) and a nitrogen atom
if they have not already been done.
db: Database
Database for collecting results.
height: float
Height of N-atom above top Ru-layer.
nlayers: int
Number of Ru-layers.
nkpts: int
Use a (nkpts * nkpts) Monkhorst-Pack grid that includes the
Gamma point.
ecut: float
Cutoff energy for plane waves.
Returns height.
name = f'Ru{nlayers}-{nkpts}x{nkpts}-{ecut:.0f}'
parameters = dict(mode=PW(ecut),
poissonsolver={'dipolelayer': 'xy'},
kpts={'size': (nkpts, nkpts, 1), 'gamma': True},
# N/Ru(0001):
slab = hcp0001('Ru', a=a, c=c, size=(1, 1, nlayers))
z = slab.positions[:, 2].max() + height
x, y =[2 / 3, 2 / 3], slab.cell[:2, :2])
slab.positions[-1] = [x, y, z], axis=2) # 2: z-axis
# Fix first nlayer atoms:
slab.constraints = FixAtoms(indices=list(range(nlayers)))
id = db.reserve(name=f'N/{nlayers}Ru(0001)', nkpts=nkpts, ecut=ecut)
if id is not None: # skip calculation if already done
slab.calc = GPAW(txt='N' + name + '.txt',
optimizer = BFGSLineSearch(slab, logfile='N' + name + '.opt')
height = slab.positions[-1, 2] - slab.positions[:-1, 2].max()
db.write(slab, id=id,
name=f'N/{nlayers}Ru(0001)', nkpts=nkpts, ecut=ecut,
# Clean surface (single point calculation):
id = db.reserve(name=f'{nlayers}Ru(0001)', nkpts=nkpts, ecut=ecut)
if id is not None:
del slab[-1] # remove nitrogen atom
slab.calc = GPAW(txt=name + '.txt',
db.write(slab, id=id,
name=f'{nlayers}Ru(0001)', nkpts=nkpts, ecut=ecut)
# Nitrogen atom:
id = db.reserve(name='N-atom', ecut=ecut)
if id is not None:
# Create spin-polarized nitrogen atom:
molecule = Atoms('N', magmoms=[3])
# Remove parameters that make no sense for an isolated atom:
del parameters['kpts']
del parameters['poissonsolver']
# Calculate energy:
molecule.calc = GPAW(txt=name + '.txt', **parameters)
db.write(molecule, id=id, name='N-atom', ecut=ecut)
return height
if __name__ == '__main__':
h = 1.2
from ase.db import connect
db = connect('convergence.db')
for n in range(1, 10):
h = adsorb(db, h, n, 7, 400)
for k in range(4, 18):
h = adsorb(db, h, 2, k, 400)
for ecut in range(350, 801, 50):
h = adsorb(db, h, 2, 7, ecut)
This diff is collapsed.
import numpy as np
from import read
from ase.constraints import FixAtoms
from ase.neb import NEB
from ase.optimize import BFGS
from import Trajectory
from ase.parallel import rank, size
from gpaw import GPAW, PW
initial = read('ini.POSCAR')
final = read('fin.POSCAR', -1)
images = [initial]
N = 4 #No. of images
z = initial.positions[:,2]
constraint = FixAtoms(mask=(z < z.min() + 1.0))
j = rank*N//size
n = size // N # number of cpu's per image
for i in np.r_[0:N]:
image = initial.copy()
if i ==j:
calc = GPAW(xc='PBE', mode=PW(350), communicator = range(j*2,j*2+n), txt = '%i.txt'%i, kpts={'size':(4,4,1),'gamma':True}, convergence={'eigenstates':1e-7})
neb = NEB(images, k=0.5, parallel=True)
qn = BFGS(neb)
if rank % (size // N) == 0:
traj = Trajectory('neb%d.traj' % j, 'w', images[1 + j], master=True)
# Creates: convergence.db
import check_convergence
# creates: catalysis/n2-on-metal.ipynb
# ... and other .ipynb files
import json
from pathlib import Path
def f(path):
data = json.loads(path.read_text())
lines = [f'# Converted from {path}\n']
n = 1
for cell in data['cells']:
if cell['cell_type'] == 'code':
lines.extend(['\n', f'# In [{n}]:\n'])
for line in cell['source']:
if line.startswith('%') or line.startswith('!'):
line = '# ' + line
n += 1
def convert(path):
data = json.loads(path.read_text())
cells = []
for cell in data['cells']:
if cell['cell_type'] == 'code':
lines = cell['source']
if (lines and
lines[0].replace(' ', '').lower().startswith('#teacher')):
data['cells'] = cells
new = path.with_name('.', 2)[0] + '.ipynb')
new.write_text(json.dumps(data, indent=1))
if __name__ == '__main__':
for path in Path().glob('*/*.master.ipynb'):
.. _projects:
.. toctree::
......@@ -9,11 +9,19 @@
.. highlight:: bash
.. toctree::
:maxdepth: 1
Summer school exercises in Jupyter notebooks
The Summer School includes a number of exercises, which are partly formulated as Jupyter Notebooks. In a Jupyter Notebook you are running the calculations on the DTU central computing facilities, but the output is displayed in your local browser.
The Summer School includes a number of :ref:`projects`, which are partly
formulated as Jupyter Notebooks. In a Jupyter Notebook you are running the
calculations on the DTU central computing facilities, but the output is
displayed in your local browser.
Unfortunately, this requires some setup, which is described below.
......@@ -33,7 +41,7 @@ Windows users
.. toctree::
:maxdepth: 1
......@@ -43,9 +51,9 @@ Mac and Linux users
.. toctree::
:maxdepth: 1
......@@ -226,7 +226,7 @@ Thursday (May 23)
Activities for GPAW developers (we start at 9:00):
* Coordination of code development and discussions about the future:
Quick tour of :ref:`projects` --- what's the current status?
Quick tour of :ref:`dev projects` --- what's the current status?
* Introduction to Sphinx and reStructuredText
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment