Overview

Slides: https://jensj.gitlab.io/gpaw-2021-talk/

  • Quick introduction to GPAW?

    • Wave function representations (modes)

  • Installing the code

    • from the Python package index (PyPI)

    • from source

  • Testing the code

  • Running the code

  • Running the code in parallel

  • PAW-data and pseudo-potentials

  • How to build the documentation

  • What’s in a gpw-file?

  • Parallel distribution of data

What is GPAW?

Python package for electronic structure calculations based on the projector-augmented wave (PAW) approximation (P. E. Blöchl, Phys. Rev. B 50, 17953 (1994)).

Based on: the atomic simulation environment (ASE), NumPy and SciPy.

The wave functions can be described with:

  • Plane-waves (pw)

  • Real-space uniform grids, multi-grid methods and finite-difference approximation (fd)

  • Atom-centered basis-functions (lcao)

A simple example

from ase.build import molecule
from gpaw import GPAW
atoms = molecule('H2O')
atoms.center(vacuum=3.0)
atoms.calc = GPAW(mode=dict(name='pw',
                            ecut=400),
                  xc='PBE')
energy = atoms.get_potential_energy()

Modes

In the PAW approximation, we need to describe smooth valence states:

gpaw/wavefunctions/pw.py

\[\tilde\psi(\mathbf{r}) = \sum_{G^2/2<E_c} c_\mathbf{G} e^{i\mathbf{G}\cdot\mathbf{r}}\]

gpaw/wavefunctions/lcao.py

\[\tilde\psi(\mathbf{r}) = \sum_{a\mu} c_{a\mu} \Phi_{a\mu}(\mathbf r - \mathbf R_a)\]

gpaw/wavefunctions/fd.py

\(N_1 \times N_2 \times N_3\) grid-points:

\[\tilde\psi(\mathbf{r}_{ijk}) = \tilde\psi_{ijk},\]
\[\mathbf{r}_{ijk} = \frac{i}{N_1}\mathbf{a}_1 + \frac{j}{N_2}\mathbf{a}_2 + \frac{k}{N_3}\mathbf{a}_3\]

Which mode?

Memory consumption:

LCAO < PW < FD

Speed:

PW > LCAO > FD (small systems)

LCAO > PW > FD (large systems)

Absolute convergence:

PW > FD > LCAO

Egg-box errors:

PW < LCAO < FD

Features:

FD > LCAO > PW (old features)

PW > LCAO > FD (new features)

Important tools

  • venv: Creation of virtual environments

  • PyPI: Python package index

  • pip: Package installer for Python

Installation with pip and PyPI

$ cd /tmp
$ python3 -m venv venv
$ ls -l venv/bin/
lrwxrwxrwx 1 jensj jensj   16 maj 20 18:51 python3 -> /usr/bin/python3
lrwxrwxrwx 1 jensj jensj    7 maj 20 18:51 python -> python3
-rwxrwxr-x 1 jensj jensj  226 maj 20 18:51 pip
-rw-r--r-- 1 jensj jensj 2189 maj 20 18:51 activate
...
$ . venv/bin/activate  # /tmp/venv/bin/ is now first in $PATH

Let’s install ASE from PyPI:

$ pip install ase
Using cached ase-3.21.1-py3-none-any.whl (2.2 MB)
...
$ ls venv/lib/python3.8/site-packages/ase/
atom.py          eos.py           phasediagram.py
atoms.py         formula.py       phonons.py
autoneb.py       ga               __pycache__
...
$ which ase
/tmp/venv/bin/ase

Warning

  • Don’t mess with $PYTHONPATH

  • Don’t move or rename the venv/ folder

Installing GPAW

setup.py, siteconfig_template.py -> (siteconfig.py)

We put something like this in ~/.gpaw/siteconfig.py:

scalapack = True
fftw = True
libraries = ['blas', 'xc', 'scalapack-openmpi', 'fftw3']

Which siteconfig.py is used?

  1. the file that $GPAW_CONFIG points to

  2. a siteconfig.py file next to siteconfig_template.py

  3. ~/.gpaw/siteconfig.py

Defaults:

noblas = False
nolibxc = False
fftw = False  # use numpy.fft
scalapack = False
libvdwxc = False
elpa = False
libraries = ['xc', 'blas']
library_dirs = []
include_dirs = []
extra_link_args = []
extra_compile_args = ['-Wall', '-Wno-unknown-pragmas', '-std=c99']
$ pip install gpaw
Downloading gpaw-21.1.0.tar.gz (1.5 MB)
...
Running setup.py install for gpaw ... done
Successfully installed gpaw-21.1.0
$ which gpaw
/tmp/venv/bin/gpaw
$ ls venv/lib/python3.8/site-packages/gpaw
ae.py                 grid_descriptor.py    poisson_extended.py
ah.py                 hamiltonian.py        poisson_extravacuum.py
...
$ ls venv/lib/python3.8/site-packages/_gpaw*.so
venv/lib/python3.8/site-packages/_gpaw.cpython-38-x86_64-linux-gnu.so

The first test

$ gpaw test
...
Could not find any atomic PAW-data or pseudopotentials!

You need to set the GPAW_SETUP_PATH environment variable to point to
the directories where PAW dataset and basis files are stored.  See
https://wiki.fysik.dtu.dk/gpaw/install.html#install-paw-datasets
for details.

Let’s fix that:

$ gpaw install-data --register ~/PAWDATA
...
Selected gpaw-setups-0.9.20000.tar.gz.  Downloading...
Extracting tarball into /home/jensj/PAWDATA
Setups installed into /home/jensj/PAWDATA/gpaw-setups-0.9.20000.
Setup path registered in /home/jensj/.gpaw/rc.py.
Current GPAW setup paths in order of search priority:
   1. /home/jensj/PAWDATA/gpaw-setups-0.9.20000
Installation complete.

And now …

$ gpaw test
 ----------------------------------------------------------------------------------------------
| python-3.8.5      /tmp/venv/bin/python3                                                      |
| gpaw-21.1.0       /tmp/venv/lib/python3.8/site-packages/gpaw/                                |
| ase-3.21.1        /tmp/venv/lib/python3.8/site-packages/ase/                                 |
| numpy-1.20.3      /tmp/venv/lib/python3.8/site-packages/numpy/                               |
| scipy-1.6.3       /tmp/venv/lib/python3.8/site-packages/scipy/                               |
| libxc-4.3.4       yes                                                                        |
| _gpaw             /tmp/venv/lib/python3.8/site-packages/_gpaw.cpython-38-x86_64-linux-gnu.so |
| MPI enabled       yes                                                                        |
| OpenMP enabled    no                                                                         |
| scalapack         yes                                                                        |
| Elpa              no                                                                         |
| FFTW              yes                                                                        |
| libvdwxc          no                                                                         |
| PAW-datasets (1)  /home/jensj/PAWDATA/gpaw-setups-0.9.20000                                  |
 ----------------------------------------------------------------------------------------------
Doing a test calculation (cores: 1): ... Done

Try also gpaw --parallel=2 test.

Getting the source code

Browse online here: https://gitlab.com/gpaw/gpaw

$ git clone git@gitlab.com:gpaw/gpaw
Cloning into 'gpaw'...
remote: Enumerating objects: 477, done.
remote: Counting objects: 100% (477/477), done.
remote: Compressing objects: 100% (193/193), done.
remote: Total 142819 (delta 299), reused 448 (delta 282), pack-reused 142342
Receiving objects: 100% (142819/142819), 44.04 MiB | 1.80 MiB/s, done.
Resolving deltas: 100% (112367/112367), done.
../_images/pies.png

Install from source

$ pip install gpaw/  # not "pip install gpaw" !!!
Processing ./gpaw
...
Installing collected packages: gpaw
  Attempting uninstall: gpaw
    Found existing installation: gpaw 21.1.0
    Uninstalling gpaw-21.1.0:
      Successfully uninstalled gpaw-21.1.0
    Running setup.py install for gpaw ... done
Successfully installed gpaw-21.1.1b1

Python files and _gpaw.cpython-38-x86_64-linux-gnu.so will now be in /tmp/venv/lib/python3.8/site-packages/.

Let’s do an editable install instead:

$ pip install -e gpaw/
...
Successfully installed gpaw-21.1.1b1
$ cat venv/lib/python3.8/site-packages/easy-install.pth
/tmp/gpaw
$ ls /tmp/gpaw/_gpaw*.so
/tmp/gpaw/_gpaw.cpython-38-x86_64-linux-gnu.so

Python will now find the gpaw and _gpaw modules in your git-clone folder.

Note

Remember to pip install -e gpaw/ when the C-code changes.

Command-line interface

gpaw/cli/

$ ase build -x diamond -a 5.43 Si si.xyz
$ gpaw run -p "mode=pw,kpts={density:4},xc=PBE" -w si.gpw si.xyz
$ gpaw dos si.gpw -p -w0 -a Si-sp -r -10 2
../_images/dos.png

See also

Atomic simulation recipes (ASR): Collection of recipes (currently based on GPAW):

$ python -m asr.relax

Tip

gpaw not on $PATH? Use python -m gpaw <sub-command> instead.

Sub-commands

Sub-command

Description

run

Run calculation with GPAW.

info

Show versions of GPAW and its dependencies

test

Test GPAW installation.

install-data

Install PAW data-sets, pseudo-potential or basis sets.

python

Run the Python interpreter with some special GPAW CLI options.

completion

Add tab-completion for Bash.

sbatch

Submit a GPAW Python script via sbatch.

dataset

Create PAW dataset.

atom

Solve radial equation for an atom.

dos

Calculate (projected) density of states from gpw-file.

gpw

Manipulate/show content of GPAW-restart file.

diag

Calculate all or some eigenvectors/values for fixed H.

symmetry

Analyze symmetry.

rpa

Run RPA-correlation calculation.

Running python scripts

Serial

$ python script.py
$ gpaw python [--dry-run=N] script.py

Using MPI

$ gpaw -P8 python script.py
$ mpiexec -n 8 gpaw python script.py

Warning

Almost the same as mpiexec -n 8 python script.py but not quite!

Use $GPAW_MPI_OPTIONS to pass options to mpiexec (like --oversubscribe).

Submit to a queue

Write a script.sh file with some magic Slurm/PBS/… stuff that you can sbatch/qsub/…:

#!/bin/sh
mpiexec gpaw python script.py

Alternatives:

  1. gpaw sbatch -- [sbatch options] script.py [script options]

  2. mq submit script.py -R 8:1h (see MyQueue)

Testing

gpaw/test/, pytest.ini

Run the test-suite (no MPI):

$ pip install pytest
$ cd /tmp/gpaw
$ pytest
.................^C^C^C^C^C^C^C^C^C
$ pip install pytest-xdist
$ pytest -n4 -q
s.................s..............................................s.. [ 11%]
....................x..........................................s.... [ 22%]
...
...........sss.....sss...............s.................s.........s.. [100%]
642 passed, 34 skipped, 9 xfailed in 927.35s (0:15:27)

With MPI (2, 4 and 8 cores):

$ mpiexec -n 2 pytest -q  # don't mix with pytest-xdist!
...

All these tests run nightly on a server somewhere.

Example: gpaw/test/something/test_thing.py:

from gpaw import GPAW

def test_stuff(in_tmp_dir, gpw_files):
    """Test some feature that needs a GPAW object.

    Also write some files and make a mess.  Pytest will clean up.
    """
    calc = GPAW(gpw_files['bcc_li_lcao'])
    ...

Special GPAW-fixtures: gpaw/test/conftest.py

Doctests:

def func(x: int) -> int:
    """Silly function.

    >>> func(1)
    2
    """
    return x + 1

Bigger tests

gpaw/test/big/, doc/tutorials/

In these folders there are *agts.py scripts (advanced gpaw test system) that define workflows that we run when I remember to run them every weekend on the Niflheim supercomputer using MyQueue:

$ mq workflow -p agts.py .
Scanning 99 scripts:  |--------------------| 100.0%
Submitting 350 tasks: |--------------------| 100.0%

doc/tutorials/jellium/jellium.agts.py:

from myqueue.workflow import run

def workflow():
    """Run jellium tutorial and create fig2.png."""
    run(script='bulk.py', cores=4)
    with run(script='surface.py', cores=4, tmax='1h'):
        with run(script='sigma.py'):
            run(script='fig2.py')

Gitlab CI

.gitlab-ci.yml, pytest.ini, mypy.ini, .flake8

CI = continuous integration: runs at every push on one of GitLab’s servers.

Test summary

  • gpaw test

  • pytest on 1, 2, 4 and 8 cores (every night)

  • AGTS (every weekend)

  • GitLab-CI (every push)

C-code

c/, c/blas.c, c/fftw.c, c/mpi.c, c/elpa.c, c/blacs.c, c/xc/libxc.c, …

  • c/transformers.c: multi-grid transformations

  • c/operators.c: finite-difference stencils

  • c/lcao.c: \(\langle\phi_\mu|\tilde v|\phi_\nu\rangle\)

  • c/lfc.c: \(f(|\mathbf{r}_{ijk}-\mathbf{R}^a|)Y_{lm}(\mathbf{r}_{ijk}-\mathbf{R}^a)\) (localized functions collection)

Warning

The C-code typically gets a pointer to the first element of an np.ndarray. Run in debug-mode (python -d ...) to check memory layout and data type!

time-profile

from ase.build import molecule
from gpaw import GPAW
atoms = molecule('H2O')
atoms.center(vacuum=3.0)
atoms.calc = GPAW(mode=dict(name='pw',
                            ecut=400),
                  xc='PBE')
energy = atoms.get_potential_energy()
$ python -m cProfile -s time h2o.py > prof.txt
$ less prof.txt
...
 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   6622    0.441    0.000    0.441    0.000 {method 'calculate' of 'XCFunctional' objects}
   1299    0.295    0.000    0.295    0.000 {built-in method _gpaw.FFTWExecute}
      6    0.289    0.048    0.289    0.048 {built-in method _gpaw.FFTWPlan}
      1    0.280    0.280    0.280    0.280 broadcast_imports.py:1(<module>)
    665    0.208    0.000    0.208    0.000 {built-in method _gpaw.pwlfc_expand}
   1484    0.150    0.000    0.150    0.000 {built-in method _gpaw.mmm}
     84    0.124    0.001    0.124    0.001 {built-in method _gpaw.symmetrize}
29k/26k    0.094    0.000    0.148    0.000 {built-in method numpy.core._multiarray_umath.implement_array_function}
 449300    0.084    0.000    0.084    0.000 spline.py:46(__call__)
    132    0.081    0.001    0.311    0.002 lda.py:14(__call__)
     47    0.062    0.001    0.159    0.003 function_base.py:2184(_vectorize_call)
    588    0.051    0.000    0.051    0.000 {built-in method marshal.loads}
     88    0.049    0.001    0.278    0.003 pw.py:824(apply_pseudo_hamiltonian)
...

LibXC, FFTW, BLAS, NumPy, …, Python, …

PAW-data, PP’s and basis sets

Make your own:

$ gpaw dataset H -f PBE -w -t abc123
...
$ gzip H.abc123.PBE
$ GPAW_SETUP_PATH=. python script.py

Specify PP/PAW-data like this: calc = GPAW(..., setups={'H': 'abc123'}).

The web-page

doc/, doc/tutorials/, doc/exercises/, …

Built once a day from reStructuredText files and docstrings using Sphinx.

$ pip install sphinx-rtd-theme
$ cd doc
$ make  # sphinx-build -b html -d build/doctrees  -n . build/html
...
$ firefox build/html/index.html

We don’t want .png and .csv files in our git repository. Instead, we generate those files by running Python scripts that have a special comment in line 1:

# creates: figure.png, table.csv
...

The .rst files contain reStructuredText directives like shown here:

.. literalinclude:: script.py

.. image:: figure.png

.. csv-table::
   :file: table.csv

Inside a gpw-file

calc.write('abc.gpw')  # , mode='all')
del calc
calc = GPAW('abc.gpw')  # ... and we're back!

A .gpw file uses the well known ase.io.ulm format: Efficient and NumPy-friendly (also used for ASE trajectory files).

$ ase ulm abc.gpw
abc.gpw  (tag: "GPAW", 1 item)
item #0:
{
    atoms: {
        cell: [[0.0, 2.715, 2.715], [2.715, 0.0, 2.715], [2.715, 2.715, 0.0]],
        numbers: <ndarray shape=(2,) dtype=int64>,
        pbc: [True, True, True],
        positions: <ndarray shape=(2, 3) dtype=float64>},
    bohr: 0.5291772105638411,
    density: {
        atomic_density_matrices: <ndarray shape=(1, 182) dtype=float64>,
        density: <ndarray shape=(1, 18, 18, 18) dtype=float64>},
    gpaw_version: 21.1.1b1,
    ha: 27.211386024367243,
    hamiltonian: {
        atomic_hamiltonian_matrices: <ndarray shape=(1, 182) dtype=float64>,
        e_coulomb: -13.668977343791727,
        e_entropy: -0.0010673898921926962,
        e_external: 0.0,
        e_kinetic: 15.569058505949739,
        e_total_extrapolated: -10.788983937133224,
        e_total_free: -10.78951763207932,
        e_xc: -12.661684474376418,
        e_zero: -0.02684692996872059,
        potential: <ndarray shape=(1, 18, 18, 18) dtype=float64>,
        xc: {}},
    parameters: {
        kpts: {'density': 4},
        mode: pw,
        xc: PBE},
    results: {
        dipole: <ndarray shape=(3,) dtype=float64>,
        energy: -10.788983937133224,
        forces: <ndarray shape=(2, 3) dtype=float64>,
        free_energy: -10.78951763207932,
        magmom: 0.0,
        magmoms: <ndarray shape=(2,) dtype=float64>,
        stress: <ndarray shape=(6,) dtype=float64>},
    scf: {
        converged: True},
    version: 3,
    wave_functions: {
        coefficients: <ndarray shape=(1, 35, 8, 577) dtype=complex128>,
        eigenvalues: <ndarray shape=(1, 35, 8) dtype=float64>,
        fermi_levels: <ndarray shape=(1,) dtype=float64>,
        ha: 27.211386024367243,
        indices: <ndarray shape=(35, 577) dtype=int32>,
        kpts: {
            atommap: <ndarray shape=(24, 2) dtype=int64>,
            bz2ibz: <ndarray shape=(729,) dtype=int64>,
            bzkpts: <ndarray shape=(729, 3) dtype=float64>,
            ibzkpts: <ndarray shape=(35, 3) dtype=float64>,
            rotations: <ndarray shape=(24, 3, 3) dtype=int64>,
            translations: <ndarray shape=(24, 3) dtype=int64>,
            weights: <ndarray shape=(35,) dtype=float64>},
        occupations: <ndarray shape=(1, 35, 8) dtype=float64>,
        projections: <ndarray shape=(1, 35, 8, 26) dtype=complex128>,
        version: 2}}
>>> from ase.io.ulm import ulmopen
>>> u = ulmopen('abc.gpw')
>>> eps = u.wave_functions.eigenvalues
>>> eps.shape  # (spins, k-points, bands)
(1, 35, 8)
>>> eps[0, 0]
array([-6.6520288 ,  5.31633309,  5.3163331 ,  5.3163331 ,  7.87544082,
        7.87544083,  7.87544083,  8.68213876])

Same thing using gpaw.GPAW.get_eigenvalues():

>>> from gpaw import GPAW
>>> GPAW('abc.gpw').get_eigenvalues(spin=0, kpt=0)
array([-6.6520288 ,  5.31633309,  5.3163331 ,  5.3163331 ,  7.87544082,
        7.87544083,  7.87544083,  8.68213876])

Parallelization distribution of data

Parallel distributions (see the parallel keyword):

  • k-points

  • plane-waves, real-space grid-points, basis functions

  • electronic bands

>>> from gpaw import GPAW
>>> calc = GPAW('abc.gpw')
>>> calc.wfs.kpt_qs
[[KPoint(..., s=0, k=0, q=0)], ..., [KPoint(..., s=0, k=34, q=34)]]
>>> len(calc.wfs.kpt_qs)
35
>>> calc.wfs.kpt_qs[0][0].psit.array.shape
(8, 577)

Tip

Naming of numpy.ndarray variables: name_abc. See naming convention for indices (a, c, v, k, q, s, G, …).

The gpaw.GPAW.get_pseudo_wave_function() method will collect all the distributed pieces and return the whole thing.

>>> wfs = calc.get_pseudo_wave_function(spin=0, kpt=0, band=0)
>>> wfs.shape
(18, 18, 18)

Descriptors:

Last slide

Slides: https://jensj.gitlab.io/gpaw-2021-talk/

Source: https://gitlab.com/jensj/gpaw-2021-talk

Links to important tools worth learning:

  • Python: The language

  • venv: Creation of virtual environments

  • PyPI: Python package index

  • pip: Package instaler for Python

  • git: Distributed version control system

  • bash: Bourne Again SHell

  • Sphinx: Python documentation generator

  • pytest: Test framework

Feedback:

Extra stuff

Python -m trick

...

def main():
    parser = argparse.ArgumentParse()
    ...

if __name__ == '__main__':
    main()

Example:

python -m gpaw.hyperfine --help

Closing words

Thank you!

Feature requests room:

  • more docstrings

  • better developer docs

  • better handling of charged systems in PW-mode

  • augment_grid by default

  • better handling of large basis sets

  • make MGGA’s robust

  • real-time TDDFT with external potentials

We need your help:

code, code review, tutorials, benchmarks, bug reporting

The “big merge request” problem:

LCAO-TDDFT-k-omega code, GPU, DO, SIC -> gpaw-master

Feedback?

Please fill in this survey

Contact:

See you in GatherTown