DeltaCodesDFT - Comparing Solid State DFT Codes, Basis Sets and Potentials

Note

This exercise is currently broken. Please skip it.

The webpage https://molmod.ugent.be/deltacodesdft provides a method for measuring the precision of a given calculation method against a chosen reference method (computational or experimental) for parameters of the equation of state (see Bulk aluminum) of elementary solids.

When performing any benchmark calculations, especially involving a large number of systems, it is important to be aware of the fact that we, humans tend to do mistakes.

Therefore the motto of this exercise is taken from Karl Popper’s “All life is problem solving”: the novelty in the scientific approach is that we actively seek to eliminate our attempted solutions.

In this exercise, in addition to the traditional, error prone method of writing output files generated using separate scripts on disk, we write the results into a database. All calculations are performed using one script, and therefore not only sharing the results with other researches is easy (by granting the access to the database) but also the precise method of performing the calculations should be shared (by presenting the script). Please consult the introduction to the ase.db module for details.

You can find out more about “reproducible science” with ASE in the following talks: Emacs + org-mode + python in reproducible research or How Python & the iPython notebook can revamp quantum chemical reseach.

We will compare PBE numbers from GPAW with http://www.wien2k.at/ for K, Ca and Ti. We use default PAW-datasets, a plane-wave cutoff of 340 eV, a k-point density of 3.5 Å and a Fermi-Dirac distribution with a 0.1 eV temperature. Copy this dcdft_gpaw.py to a place in your file area:

from time import time
import numpy as np
import ase.db
from ase.test.tasks.dcdft import DeltaCodesDFTCollection as Collection
from gpaw import GPAW, PW, FermiDirac

c = ase.db.connect('dcdft.db')

ecut = 340
kptdensity = 3.5
width = 0.10

collection = Collection()

for name in ['K', 'Ca', 'Ti']:
    atoms = collection[name]
    cell = atoms.get_cell()
    # Loop over volumes:
    for n, x in enumerate(np.linspace(0.98, 1.02, 5)):
        id = c.reserve(name=name, x=x)
        if id is None:
            # This calculation has been or is being done:
            continue

        atoms.set_cell(cell * x, scale_atoms=True)
        atoms.calc = GPAW(txt='%s-%d.txt' % (name, n),
                          mode=PW(ecut),
                          xc='PBE',
                          kpts={'density': kptdensity},
                          parallel={'band': 1},
                          occupations=FermiDirac(width))

        t1 = time()
        atoms.get_potential_energy()
        t2 = time()

        # Write to database:
        c.write(atoms, name=name, x=x, time=t2 - t1,
                ecut=ecut, kptdensity=kptdensity, width=width)

        del c[id]

Read the script and try to understand it. Run the script by typing:

$ python3 dcdft_gpaw.py

It should take about 15 minutes to run the script. Note that you can start several instances of the script simultaneously in order to speed things up.

The script will generate .txt files and an SQLite3 database file. Watch the progess as the calculations run:

$ ase db dcdft.db -c +x,time

Examine the equation of state (see Bulk aluminum) using ase gui:

$ ase gui dcdft.db@name=Ca

Note

The PBE reference values from https://molmod.ugent.be/deltacodesdft are:

element

\(V\) [Å\(^3\)/atom]

\(B\) [GPa]

K

73.68

3.6

Ca

42.20

17.1

Ti

17.39

112.2

Extract the results from the database in order to calculate the parameters of the equation of state:

 K  73.6852   3.6070  19.7231
Ca  42.5903  24.2170 -13.2987
Ti  17.3542 113.9593   3.5966

and use the script available from https://molmod.ugent.be/deltacodesdft to calculate the Delta factors.