Getting started with GPAW

In this exercise we will relax structures of simple molecules and calculate their binding energies.

At first let’s have a look at a simple GPAW calculation for obtaining forces on the H2 molecule:

from ase import Atoms
from gpaw import GPAW

d = 0.74
a = 6.0

atoms = Atoms('H2',
              positions=[(0, 0, 0),
                         (0, 0, d)],
              cell=(a, a, a))
atoms.center()

calc = GPAW(mode='pw', txt='h2.txt')
atoms.calc = calc
print(atoms.get_forces())

In case this is the first ASE script we have seen so far, a few comments are in order:

  • At the top is a series of import statements. These load the Python modules we are going to use.

  • An Atoms object is created, specifying an initial (possibly bad) guess for the atomic positions.

  • An GPAW calculator is created. A calculator can evaluate quantities such as energies and forces on a collection of atoms.

  • We create the GPAW calculator in plane-wave (pw) mode, with output to h2.txt.

  • The calculator is associated with the Atoms object by calling atoms.calc = calc.

  • Calling atoms.get_forces() triggers the DFT calculation for H2 and forces on the atoms are return (here: printed).

Performing a structure optimization

A structure optimization, also called a relaxation, is a series of calculations used to determine the minimum-energy structure of a given system. This involves multiple calculations of the atomic forces \(\mathbf F^a = -\tfrac{\partial E}{\partial \mathbf R^a}\) with respect to the atomic positions \(\mathbf R^a\) as the atoms are moved downhill according to an optimization algorithm.

The following script uses the GPAW calculator to optimize the structure of H2.

from ase import Atoms
from gpaw import GPAW
from ase.optimize import QuasiNewton

d = 1.0
a = 6.0

atoms = Atoms('H2',
              positions=[(0, 0, 0),
                         (0, 0, d)],
              cell=(a, a, a))
atoms.center()

calc = GPAW(mode='pw', txt='h2.txt')
atoms.calc = calc

opt = QuasiNewton(atoms, trajectory='h2.traj')
opt.run(fmax=0.05)

Apart from the commands in the static calculation above:

  • An optimizer is created and associated with the Atoms object. It is also given an optional argument, trajectory, which specifies the name of a file into which the positions will be saved for each step in the geometry optimization.

  • Finally the call opt.run(fmax=0.05) will run the optimization algorithm until all atomic forces are below 0.05 eV per Ångström.

Run the above structure optimization.

This will print the (decreasing) total energy for each iteration until it converges, leaving the file h2.traj in the working directory. Use the command ase gui to view the trajectory file, showing each step of the optimization.

Structure optimization of H2O

Adapt the above script as needed and calculate the structure of a H2O molecule using the GPAW calculator. Note that water is not a linear molecule. If you start with a linear molecule, the minimization may not be able to break the symmetry. Be sure to visualize the final configuration to check that it is reasonable.

The cell must be centered in order to prevent atoms from lying too close to the boundary, as the boundary conditions are zero by default.

During the calculation a lot of text is printed ot the output file, or if we haven’t specified any to the terminal (standard output). This includes the parameters used in the calculation: Atomic positions, grid spacing, XC functional (GPAW uses LDA by default) and many other properties. For each iteration in the self-consistency cycle one line is printed with the energy and convergence measures. After the calculation the energy contributions, band energies and forces are listed.

Use ase gui to visualize the structure and measure bond lengths and bond angles. Bond lengths and angles are shown automatically if you select two or three atoms at a time.

The following script performs the structure optimization for H2O: h2o_opt.py.

Atomization energies

Now that we know the structure of H2O, we can calculate other interesting properties like the molecule’s atomization energy.

The atomization energy of a molecule is equal to the total energy of the molecule minus the sum of the energies of each of its constituent isolated atoms. For example, the atomization energy of H2 is \(E[\mathrm{H}_2] - 2 E[\mathrm H]\).

GPAW calculations are by default spin-paired, i.e. the spin-up and spin-down densities are assumed to be equal. As this is not the case for isolated atoms, it will be necessary to instruct GPAW to do something different:

calc = GPAW(mode='pw', hund=True)

With the hund keyword, Hund’s rule is applied to initialize the atomic states, and the calculation will be made spin-polarized.

Write a script which calculates the total energy of the isolated O and H atoms, and calculate the atomization energy of H2O.

The following script calculates the atomization energy for H2O: h2o.py.

Exchange and correlation functionals

So far we have been using GPAW’s default parameters. The default exchange-correlation functional is LDA. This is not very accurate, and in particular overestimates atomization energies. You can specify different XC functionals to the calculator using GPAW(xc=name), where name is a string such as 'LDA', 'PBE' or 'RPBE'.

Calculate the atomization energy of H2O with LDA and PBE (just reuse the geometry from the LDA optimization, i.e. do not repeat the minimization).