diff git a/doc/devel/ase_optimize/agts.py b/doc/devel/ase_optimize/agts.py
new file mode 100644
index 0000000000000000000000000000000000000000..13adc95353a0f24eeeb5484875131404173b7a98
 /dev/null
+++ b/doc/devel/ase_optimize/agts.py
@@ 0,0 +1,16 @@
+# Creates: systems.db
+
+
+def workflow():
+ from myqueue.job import Job
+ return [Job('agts.py'),
+ Job('run_tests_emt.py', deps='agts.py'),
+ Job('run_tests.py+1@8x1d', deps='agts.py'),
+ Job('run_tests.py+2@8x1d', deps='agts.py'),
+ Job('analyze.py',
+ deps=['run_tests_emt.py', 'run_tests.py+1', 'run_tests.py+2'])]
+
+
+if __name__ == '__main__':
+ from ase.optimize.test.systems import create_database
+ create_database() # creates systems.db
diff git a/doc/devel/ase_optimize/analyze.py b/doc/devel/ase_optimize/analyze.py
new file mode 100644
index 0000000000000000000000000000000000000000..3abd24737f680c0ecc043c3807f3aa3d0555cdc7
 /dev/null
+++ b/doc/devel/ase_optimize/analyze.py
@@ 0,0 +1,13 @@
+# Creates: emtiteration.csv, lcaotime.csv, systems.csv
+import ase.db
+from ase.optimize.test.analyze import analyze
+
+analyze('resultsemt.db', 'emt')
+analyze('resultslcao.db', 'lcao')
+
+db1 = ase.db.connect('systems.db')
+
+with open('systems.csv', 'w') as f:
+ print('testname,description', file=f)
+ for row in db1.select():
+ print('{},{}'.format(row.formula, row.description), file=f)
diff git a/doc/devel/ase_optimize/ase_optimize.rst b/doc/devel/ase_optimize/ase_optimize.rst
index 5844b14a672bab22bdc05302f66353c7bbc54964..b786f7fdfd14cb8263c0e4cdf9a83eea1de93b51 100644
 a/doc/devel/ase_optimize/ase_optimize.rst
+++ b/doc/devel/ase_optimize/ase_optimize.rst
@@ 3,7 +3,9 @@
===============
Optimizer tests
===============
This page shows benchmarks of optimizations done with our different optimizers.
+
+This page shows benchmarks of optimizations done with ASE's different
+:mod:`optimizers `.
Note that the iteration number (steps) is not the same as the number of force
evaluations. This is because some of the optimizers uses internal line searches
or similar.
@@ 14,107 +16,38 @@ Different optimizers may perform the same number of steps, but along a different
path, so the time spent on calculation of energy/forces may be different
due to different convergence of the selfconsistent field.
G2
==
PBE relaxation of molecules from the G2 set.

On the plots: the number of optimizer force calls (stats), the total run time,
the systems with the largest number of optimizer force calls, and the number of
systems for which optimization failed.

In the corresponding tables: the most common value of the total energy
("relaxed energy"), and the differences (optimizer  "relaxed energy").
The dot sign denotes the above difference below a threshold
(same as the printed precision of "relaxed energy" in eV),
and "N/A" denotes an optimization failure.
Only the systems that fail to converge or converge to a
total energy above the threshold are given in the tables.

Calculator used: GPAW mode='lcao' (see :git:`~doc/devel/ase_optimize/g2_dzp.py`)

Limit of optimizer steps: 25

.. csvtable::
 :file: g2_dzp.csv
 :header: optimizer, failed, time, steps


N2Cu
====
Relaxation of Cu surface.

Calculator used: EMT

.. csvtable::
 :file: N2Cusurf.csv

N2 adsorption on relaxed Ru surface

Calculator used: EMT

.. csvtable::
 :file: N2CuN2.csv

Cu_bulk
=======
Bulk relaxation of Cu where atoms has been rattled.

Calculator used: EMT

.. csvtable::
 :file: Cu_bulk.csv

CO_Au111
========
CO adsorption on Au

Calculator used: EMT

.. csvtable::
 :file: CO_Au111.csv
H2
==
Geometry optimization of gasphase molecule.
+Test systems
+============
Calculator used: EMT
+These are the test systems (:download:`systems.db`):
.. csvtable::
 :file: H2emt.csv

Calculator used: GPAW
+ :file: systems.csv
+ :headerrows: 1
.. csvtable::
 :file: H2gpaw.csv
C5H12
=====
Geometry optimization of gasphase molecule.
+EMT calculations
+================
Calculator used: GPAW (lcao)
+Calculation done with :class:`~ase.calculators.emt.EMT`. Number of steps:
.. csvtable::
 :file: C5H12gpaw.csv

nanoparticle
============
Adsorption of a NH on a Pd nanoparticle.

Calculator used: GPAW (lcao)
+ :file: emtiterations.csv
+ :headerrows: 1
.. csvtable::
 :file: nanoparticle.csv
NEB
=======
Diffusion of gold atom on Al(100) surface.
+GPAWLCAO calculations
+======================
Calculator used: EMT
+Parameters::
.. csvtable::
 :file: nebemt.csv
+ GPAW(mode='lcao',
+ basis='dzp',
+ kpts={'density': 2.0})
Calculator used: GPAW (lcao)
+Absolute time relative to fastest optimizer:
.. csvtable::
 :file: nebgpaw.csv
+ :file: lcaotime.csv
+ :headerrows: 1
diff git a/doc/devel/ase_optimize/g2_dzp.agts.py b/doc/devel/ase_optimize/g2_dzp.agts.py
deleted file mode 100644
index 0b421b099d666756e5eedde158b9981034d1e881..0000000000000000000000000000000000000000
 a/doc/devel/ase_optimize/g2_dzp.agts.py
+++ /dev/null
@@ 1,6 +0,0 @@
from myqueue.job import Job


def workflow():
 jobs = [Job('g2_dzp.py+{}@4x13h'.format(i)) for i in range(10)]
 return jobs + [Job('g2_dzp_csv.py', deps=jobs)]
diff git a/doc/devel/ase_optimize/g2_dzp.py b/doc/devel/ase_optimize/g2_dzp.py
deleted file mode 100644
index 9ec6894103d55c3948aedf9b7e6c006da2c5411f..0000000000000000000000000000000000000000
 a/doc/devel/ase_optimize/g2_dzp.py
+++ /dev/null
@@ 1,40 +0,0 @@
import time

import ase.db
import ase.optimize
from ase.collections import g2

from gpaw import GPAW, Mixer


optimizers = ['BFGS', 'BFGSLineSearch',
 'FIRE', 'GoodOldQuasiNewton',
 'LBFGS', 'LBFGSLineSearch']

con = ase.db.connect('g2_dzp.db')
for name, atoms in zip(g2.names, g2):
 if len(atoms) == 1:
 continue
 atoms.center(vacuum=3.5)
 for optimizer in optimizers:
 id = con.reserve(name=name, optimizer=optimizer)
 if id is None:
 continue
 mol = atoms.copy()
 mol.calc = GPAW(mode='lcao',
 basis='dzp',
 h=0.17,
 xc='PBE',
 mixer=Mixer(0.05, 2),
 txt='{0}{1}.txt'.format(name, optimizer))
 Optimizer = getattr(ase.optimize, optimizer)
 opt = Optimizer(mol)
 t = time.time()
 try:
 opt.run(fmax=0.05, steps=50)
 except Exception as ex:
 print(name, optimizer, ex)
 continue
 con.write(mol, name=name, optimizer=optimizer,
 steps=opt.nsteps, time=time.time()  t)
 del con[id]
diff git a/doc/devel/ase_optimize/g2_dzp_csv.py b/doc/devel/ase_optimize/g2_dzp_csv.py
deleted file mode 100644
index 672c07e7a47a8dabcea74bcc4f487084befdcd5e..0000000000000000000000000000000000000000
 a/doc/devel/ase_optimize/g2_dzp_csv.py
+++ /dev/null
@@ 1,31 +0,0 @@
# Creates: g2_dzp.csv
from __future__ import print_function
import numpy as np
import ase.db

con = ase.db.connect('g2_dzp.db')
N = 6 # number of optimizers
M = len(con) // N # number of molecules
data = []
for row in con.select():
 data.append((row.optimizer,
 row.name,
 row.get('time', 42),
 row.get('steps', 42),
 row.get('energy', 42),
 row.get('fmax', 42)))
data.sort()
optimizers = [x[0] for x in data[::M]]
data = np.array([x[2:] for x in data]).reshape((N, M, 4))
e0 = data[:, :, 2].min(0)
results = []
for opt, d in zip(optimizers, data):
 ok = (d[:, 3] < 0.05) & (d[:, 2] < e0 + 0.1)
 failed = M  sum(ok)
 time, niter = d[ok, :2].mean(0)
 results.append((failed, time, niter, opt))

with open('g2_dzp.csv', 'w') as f:
 for failed, time, niter, opt in sorted(results):
 print('{0}, {1}, {2:.1f}, {3:.1f}'.format(opt, failed, time, niter),
 file=f)
diff git a/doc/devel/ase_optimize/run_tests.py b/doc/devel/ase_optimize/run_tests.py
new file mode 100644
index 0000000000000000000000000000000000000000..45957a38c6debeacc9501eccc4d6be94bc75c4c9
 /dev/null
+++ b/doc/devel/ase_optimize/run_tests.py
@@ 0,0 +1,22 @@
+import ase.db
+from ase.optimize.test.test import (test_optimizer, all_optimizers,
+ get_optimizer)
+from gpaw import GPAW
+
+
+db1 = ase.db.connect('systems.db')
+db = ase.db.connect('resultslcao.db')
+
+
+def lcao(txt):
+ return GPAW(mode='lcao',
+ basis='dzp',
+ kpts={'density': 2.0},
+ txt=txt)
+
+
+systems = [row.toatoms() for row in db1.select()]
+
+for opt in all_optimizers:
+ optimizer = get_optimizer(opt)
+ test_optimizer(systems, optimizer, lcao, 'lcao', db)
diff git a/doc/devel/ase_optimize/run_tests_emt.py b/doc/devel/ase_optimize/run_tests_emt.py
new file mode 100644
index 0000000000000000000000000000000000000000..d8342d13dc45ed8c1620584912bdb18ab3b44cb3
 /dev/null
+++ b/doc/devel/ase_optimize/run_tests_emt.py
@@ 0,0 +1,14 @@
+import ase.db
+from ase.calculators.emt import EMT
+from ase.optimize.test.test import (test_optimizer, all_optimizers,
+ get_optimizer)
+
+
+db1 = ase.db.connect('systems.db')
+db = ase.db.connect('resultsemt.db', serial=True)
+
+systems = [row.toatoms() for row in db1.select() if row.formula != 'C5H12']
+
+for opt in all_optimizers:
+ optimizer = get_optimizer(opt)
+ test_optimizer(systems, optimizer, EMT, 'emt', db)
diff git a/doc/devel/ase_optimize/systems.py b/doc/devel/ase_optimize/systems.py
deleted file mode 100644
index e3790549b518be00b364459d600baba0f1d6e27e..0000000000000000000000000000000000000000
 a/doc/devel/ase_optimize/systems.py
+++ /dev/null
@@ 1,23 +0,0 @@
from ase.db import connect
from ase.optimize.test.systems import create_database

from gpaw import GPAW, PW


def create_database_gpaw():
 create_database()
 db = connect('systems.db')
 systems = [row.toatoms() for row in db.select()]
 db = connect('systemsgpaw.db')
 for atoms in systems:
 if atoms.number_of_lattice_vectors != 3:
 atoms.center(vacuum=3.5)
 atoms.calc = GPAW(mode=PW(500), # 'lcao',
 # basis='dzp',
 kpts={'density': 2.0},
 txt=None)
 db.write(atoms)


if __name__ == '__main__':
 create_database_gpaw()
diff git a/doc/exercises/rpa/rpa.rst b/doc/exercises/rpa/rpa.rst
index 045a48fce2cb42e48c7e27ecf4d5a39890f9afe6..79dde5165efbeea53a4659db96203788ce20fc3e 100644
 a/doc/exercises/rpa/rpa.rst
+++ b/doc/exercises/rpa/rpa.rst
@@ 4,16 +4,16 @@
RPA calculation of the cohesive energy of Si
============================================
In this exercise we will use GPAW to calculate the cohesive energy of
silicon. When calculating total energies, we shall split the
exchangecorrelation energy into an "exact" exchange contribution and
calculate the remaining correlation energy in the randomphase
approximation (RPA). A comprehesive introduction to the RPA can be
+In this exercise we will use GPAW to calculate the cohesive energy of
+silicon. When calculating total energies, we shall split the
+exchangecorrelation energy into an "exact" exchange contribution and
+calculate the remaining correlation energy in the randomphase
+approximation (RPA). A comprehesive introduction to the RPA can be
found in [Ren]_, and you are also advised to look at :ref:`this page `.
.. note::

 RPA calculations are typically rather timeconsuming, so in
+
+ RPA calculations are typically rather timeconsuming, so in
this exercise we shall cut a few corners!
We obtain the cohesive energy `E_\text{coh}` of silicon as
@@ 22,7 +22,7 @@ We obtain the cohesive energy `E_\text{coh}` of silicon as
E_\text{coh} = E_\text{at}  0.5 E_\text{bulk}
where `E_\text{at}` is the total energy of an isolated silicon atom, and
+where `E_\text{at}` is the total energy of an isolated silicon atom, and
`E_\text{bulk}` is the total energy per primitive unit cell of bulk silicon.
Do you know where the factor of 0.5 comes from?
@@ 32,27 +32,27 @@ In DFT, we partition the total energy as
E_\text{DFT} = E_\text{Kin} + E_\text{ie} + E_\text{ii} + E_\text{Hartree} + E_\text{XC}
The purpose of this exercise is to explore different approximations for
+The purpose of this exercise is to explore different approximations for
the exchangecorrelation energy, `E_\text{XC}`.
PBE cohesive energy  bulk
==========================
We start with a generalizedgradient approximation to `E_\text{XC}` using the
+We start with a generalizedgradient approximation to `E_\text{XC}` using the
PBE functional:
.. math::
E_\text{XC} = E_\text{PBE}
The script :download:`si.pbe.py` calculates the total
energy of bulk silicon. The only parameters we need to choose are the
planewave cutoff (i.e. number of plane waves in our basis set to describe
the electron wavefunctions), the kpoint sampling of the Brillouin zone,
and the lattice parameter of Si. We could calculate the
lattice parameter ourselves from first principles, but in order to compare
to previous calculations, we instead choose the experimental value of
+The script :download:`si.pbe.py` calculates the total
+energy of bulk silicon. The only parameters we need to choose are the
+planewave cutoff (i.e. number of plane waves in our basis set to describe
+the electron wavefunctions), the kpoint sampling of the Brillouin zone,
+and the lattice parameter of Si. We could calculate the
+lattice parameter ourselves from first principles, but in order to compare
+to previous calculations, we instead choose the experimental value of
5.421 Ang [Harl]_.
Make sure you understand what the script is doing, and then try running
@@ 76,8 +76,8 @@ side length `L`, and increase `L` until the replicas no longer see each other.
Unfortunately, the larger the value of `L`, the more plane waves we have,
which slows down the calculation considerably.
Therefore, for the purpose of this exercise, we shall not actually perform the
calculations on the isolated Si atom  instead just provide the numbers as
+Therefore, for the purpose of this exercise, we shall not actually perform the
+calculations on the isolated Si atom  instead just provide the numbers as
reference data. In the next section a sample script will be given to show how to
generate the following numbers:
@@ 93,17 +93,17 @@ generate the following numbers:
12.0 0.85109735188
======== ===================
The first column gives the side length (in Angstroms) of the simulation cell
+The first column gives the side length (in Angstroms) of the simulation cell
containing the isolated atom, and the second gives the total
energy in eV.
From the above data and your own calculations, calculate the cohesive energy
of silicon using the PBE functional to describe exchangecorrelation.
Compare your result to the value of 4.55 eV reported in
+From the above data and your own calculations, calculate the cohesive energy
+of silicon using the PBE functional to describe exchangecorrelation.
+Compare your result to the value of 4.55 eV reported in
[Olsen]_.
.. note::

+
The total energy delivered by GPAW is not an absolute value, but rather
given with respect to a
reference energy. It turns out that in this case, the reference
@@ 121,18 +121,18 @@ which is
.. math::
E_\text{XC} = E_\text{EXX}
An expression for the exact exchange energy `E_\text{EXX}` can be found e.g. in
+An expression for the exact exchange energy `E_\text{EXX}` can be found e.g. in
equation (9) of [Olsen]_. The main points to note are that:
* it is fully nonlocal  to get the energy we must integrate over `\mathbf{r}`
 and `\mathbf{r}'`, which is expensive.
+ and `\mathbf{r}'`, which is expensive.
* it requires knowledge of the wavefunctions, not just
 the density, which again makes it more expensive to compute.
+ the density, which again makes it more expensive to compute.
* in the formalism used here we calculate `E_\text{EXX}` nonselfconsistently;
 that is, we use one approximation for the exchangecorrelation energy
 (PBE) to obtain the wavefunctions, then use these wavefunctions to
+* in the formalism used here we calculate `E_\text{EXX}` nonselfconsistently;
+ that is, we use one approximation for the exchangecorrelation energy
+ (PBE) to obtain the wavefunctions, then use these wavefunctions to
construct the exchange energy under a different
approximation. As a result, this method is described as EXX\@PBE; had we
used LDA to obtain the wavefunctions, we would have EXX\@LDA etc.
@@ 145,9 +145,9 @@ from ``exx.py`` in our script. The ``calculate`` method performs the
calculation of the exchange energy, while the ``get_total_energy`` method
returns the total energy of our system with `E_\text{XC}=E_\text{EXX}`.
The script :download:`si.pbe+exx.py` calculates the total
energy of bulk Si in the EXX\@PBE approximation. The calculation
proceeds in two parts  first, a PBE calculation which is identical
+The script :download:`si_pbe_exx.py` calculates the total
+energy of bulk Si in the EXX\@PBE approximation. The calculation
+proceeds in two parts  first, a PBE calculation which is identical
to that of the previous section. Second, the exchange
part. This part is much slower, and it is a good idea to run on a few
processors  it takes about 5 minutes on 4 CPUs.
@@ 176,16 +176,16 @@ following output in ``pbe_and_exx_energies.txt``::
12.0 0.852243293624 9.75141580104
13.0 0.852570610869 9.75125973558
Note that :download:`atom/si.atom.pbe+exx.py` also contains
some additional tweaking not required for the bulk calculation,
+Note that :download:`atom/si.atom.pbe+exx.py` also contains
+some additional tweaking not required for the bulk calculation,
most importantly spinpolarization; by Hund's
rules, we expect a spinpolarized atom to be more stable than the
nonspinpolarized case.
You now have enough information to calculate the cohesive energy in
the EXX\@PBE approximation. Compare your value to that of 2.82 eV given in
[Olsen]_. This number is dramatically different to
the experimental value of 4.68 eV, and highlights the danger of
+[Olsen]_. This number is dramatically different to
+the experimental value of 4.68 eV, and highlights the danger of
neglecting correlation in solids!
@@ 216,21 +216,21 @@ states which turns out to work rather well (more details below).
Like for exact exchange, the first part of our RPA calculation is performed
at the PBE level to obtain the ground state density. We then use this density
to obtain the wavefunctions both of the occupied and some of the unoccupied
states. The script :download:`si.rpa_init_pbe.py` performs
this step; note it is essentially identical to
+states. The script :download:`si.rpa_init_pbe.py` performs
+this step; note it is essentially identical to
:download:`si.pbe.py` apart from the allimportant
``diagonalize_full_hamiltonian`` line. However note that we have reduced
the kpoint grid to a 4x4x4 sampling.
Having performed this step (which should take ~1 minute on 4 CPUs) we now
calculate the correlation energy using :download:`si.rpa.py`,
which imports the ``RPACorrelation`` class from ``rpa.py``. All the
computational details are read from the ``bulk.gpw`` file; the only input
we need specify is the number of plane waves used to describe `\chi_0`.
Here we give a list of values, which means that the correlation energy
will be calculated for each planewave cutoff (in eV). The reason for
this procedure is described below. Note that in principle we also need
to specify the number of unoccupied bands used in the construction of
+calculate the correlation energy using :download:`si.rpa.py`,
+which imports the ``RPACorrelation`` class from ``rpa.py``. All the
+computational details are read from the ``bulk.gpw`` file; the only input
+we need specify is the number of plane waves used to describe `\chi_0`.
+Here we give a list of values, which means that the correlation energy
+will be calculated for each planewave cutoff (in eV). The reason for
+this procedure is described below. Note that in principle we also need
+to specify the number of unoccupied bands used in the construction of
`\chi_0`  however here this choice is made by the code,
and sets the number of bands to equal the number of plane waves describing `\chi_0`.
Now, run :download:`si.rpa.py` (4 minutes, 4 CPUs).
@@ 252,14 +252,14 @@ the correlation energy by over 1 eV.
In order to converge the correlation energy, we should increase the planewave
cutoff describing `\chi_0` (and implicitly, the number of empty states).
However it is noted in [Harl]_ that for the electron
+However it is noted in [Harl]_ that for the electron
gas, one expects the correlation energy to scale as
.. math::
E_\text{RPA}(E_{cut}) = E_\text{RPA}(\infty) + A E_{cut}^{1.5}
where `E_{cut}` is the planewave cutoff describing `\chi_0`. Empirically, this
expression seems to work beyond the electron gas.
+expression seems to work beyond the electron gas.
Test this expression for silicon by plotting the correlation energy against
`E_{cut}^{1.5}`; the intercept of the straight line should give
@@ 281,12 +281,12 @@ to obtain the extrapolated correlation energy for the single atom.
Combining the correlation energies with the EXX\@PBE calculations of the
previous section, you should now be able to calculate the cohesive energy
of silicon using exact exchange and the RPA correlation energy.
+of silicon using exact exchange and the RPA correlation energy.
* Compare the result of using the extrapolated correlation energies with that
at a fixed cutoff of 164 eV.
* Compare your value to that of 4.32 eV given in [Olsen]_
+* Compare your value to that of 4.32 eV given in [Olsen]_
and the experimental value, 4.68 eV.
@@ 300,12 +300,12 @@ generalizedgradient PBE functional. Indeed, according to table VII of
range of materials. The strength of the RPA lies in its ability to describe
longrange correlation effects, e.g. in systems exhibiting van der Waals bonds.
Unfortunately, the complexity of these systems does not allow us to study them
in a quick exercise like this one. Nonetheless the procedure of calculating the
total energy employed above is exactly the same when applied to more complicated
+in a quick exercise like this one. Nonetheless the procedure of calculating the
+total energy employed above is exactly the same when applied to more complicated
systems.
In order to get a consistent, highquality description of both longrange
and shortrange correlation it is desirable to move beyond the RPA 
+In order to get a consistent, highquality description of both longrange
+and shortrange correlation it is desirable to move beyond the RPA 
but that's another story...
diff git a/doc/images.py b/doc/images.py
index 70f2614d4193da4ef6090fb5a6a1e0079ce73678..9e8105f6fe32154604f2e7eb3322977cb4d66259 100644
 a/doc/images.py
+++ b/doc/images.py
@@ 17,6 +17,8 @@ except ImportError:
from urllib.request import urlopen
from urllib.error import HTTPError
import os
+from pathlib import Path
+
srcpath = 'http://wiki.fysik.dtu.dk/gpawfiles'
agtspath = 'http://wiki.fysik.dtu.dk'
@@ 148,16 +150,15 @@ get('static', ['NOMAD_Logo_supported_by.png'])
def setup(app):
 # Get png files and other stuff from the AGTS scripts that run
+ # Get png and csv files and other stuff from the AGTS scripts that run
# every weekend:
 from gpaw.test.big.agts import AGTSQueue
 queue = AGTSQueue()
 queue.collect()
names = set()
 for job in queue.jobs:
 if not job.creates:
+ for path in Path().glob('**/*agts.py'):
+ line1 = path.read_text().split('\n', 1)[0]
+ if not line1.startswith('# Creates:'):
continue
 for name in job.creates:
+ for name in line1.split(':')[1].split(','):
+ name = name.strip()
if name in names:
raise RuntimeError(
'The name {0!r} is used in more than one place!'
@@ 165,8 +166,8 @@ def setup(app):
names.add(name)
# the files are saved by the weekly tests under agtspath/agtsfiles
# now we are copying them back to their original run directories
 path = os.path.join(job.dir, name)
 if os.path.isfile(path):
+ file = path.with_name(name)
+ if file.is_file():
continue
 print(path, 'copied from', agtspath)
 get('agtsfiles', [name], job.dir, source=agtspath)
+ print(file, 'copied from', agtspath)
+ get('agtsfiles', [name], str(file.parent), source=agtspath)
diff git a/doc/tutorials/lattice_constants/al.agts.py b/doc/tutorials/lattice_constants/al.agts.py
index 25b4c452c78ae62ade6d4cc3c3a278a8ff5fe9a8..9293431daf49a05756d0ef05e13e92652bac4d3f 100644
 a/doc/tutorials/lattice_constants/al.agts.py
+++ b/doc/tutorials/lattice_constants/al.agts.py
@@ 1,39 +1,6 @@
# Creates: Al_conv_ecut.png, Al_conv_k.png
from myqueue.job import Job
def workflow():
 return [
 Job('al.py@8x12h'),
 Job('al.agts.py', deps=['al.py'])]


if __name__ == '__main__':
 import pylab as plt
 from ase.eos import EquationOfState
 from ase.io import read

 def fit(filename):
 configs = read(filename + '@:')
 volumes = [a.get_volume() for a in configs]
 energies = [a.get_potential_energy() for a in configs]
 eos = EquationOfState(volumes, energies)
 v0, e0, B = eos.fit()
 return (4 * v0)**(1 / 3.0)

 cutoffs = range(200, 501, 50)
 a = [fit('Al%d.txt' % ecut) for ecut in cutoffs]
 plt.figure(figsize=(6, 4))
 plt.plot(cutoffs, a, 'o')
 plt.axis(ymin=4.03, ymax=4.05)
 plt.xlabel('Planewave cutoff energy [eV]')
 plt.ylabel('lattice constant [Ang]')
 plt.savefig('Al_conv_ecut.png')

 kpoints = range(4, 17)
 plt.figure(figsize=(6, 4))
 a = [fit('Al%02d.txt' % k) for k in kpoints]
 plt.plot(kpoints, a, '')
 plt.xlabel('number of kpoints')
 plt.ylabel('lattice constant [Ang]')
 plt.savefig('Al_conv_k.png')
+ return [Job('al.py@8x12h'),
+ Job('al_analysis.py', deps=['al.py'])]
diff git a/doc/tutorials/lattice_constants/al_analysis.py b/doc/tutorials/lattice_constants/al_analysis.py
new file mode 100644
index 0000000000000000000000000000000000000000..84615dad245c82b25701719b7271afa100bd06f5
 /dev/null
+++ b/doc/tutorials/lattice_constants/al_analysis.py
@@ 0,0 +1,32 @@
+# Creates: Al_conv_ecut.png, Al_conv_k.png
+import matplotlib.pyplot as plt
+
+from ase.eos import EquationOfState
+from ase.io import read
+
+
+def fit(filename):
+ configs = read(filename + '@:')
+ volumes = [a.get_volume() for a in configs]
+ energies = [a.get_potential_energy() for a in configs]
+ eos = EquationOfState(volumes, energies)
+ v0, e0, B = eos.fit()
+ return (4 * v0)**(1 / 3.0)
+
+
+cutoffs = range(200, 501, 50)
+a = [fit('Al%d.txt' % ecut) for ecut in cutoffs]
+plt.figure(figsize=(6, 4))
+plt.plot(cutoffs, a, 'o')
+plt.axis(ymin=4.03, ymax=4.05)
+plt.xlabel('Planewave cutoff energy [eV]')
+plt.ylabel('lattice constant [Ang]')
+plt.savefig('Al_conv_ecut.png')
+
+kpoints = range(4, 17)
+plt.figure(figsize=(6, 4))
+a = [fit('Al%02d.txt' % k) for k in kpoints]
+plt.plot(kpoints, a, '')
+plt.xlabel('number of kpoints')
+plt.ylabel('lattice constant [Ang]')
+plt.savefig('Al_conv_k.png')
diff git a/doc/tutorials/lattice_constants/lattice_constants.rst b/doc/tutorials/lattice_constants/lattice_constants.rst
index 913c55580cf810f71c12e7e544373a5e3399ee95..66377c3e3e5d4016488a3c1f5a6dc52a863bb339 100644
 a/doc/tutorials/lattice_constants/lattice_constants.rst
+++ b/doc/tutorials/lattice_constants/lattice_constants.rst
@@ 34,4 +34,4 @@ with respect to number of **k**points:
.. image:: Al_conv_k.png
(see also :download:`analysis script `).
+(see also :download:`analysis script `).