# Getting started with GPAW¶

In this exercise we will calculate structures and binding energies for simple molecules.

## 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 `EMT calculator`

to optimize the structure of H_{2}.

```
# creates: h2.emt.traj
from ase import Atoms
from ase.calculators.emt import EMT
from ase.optimize import QuasiNewton
system = Atoms('H2', positions=[[0.0, 0.0, 0.0],
[0.0, 0.0, 1.0]])
calc = EMT()
system.calc = calc
opt = QuasiNewton(system, trajectory='h2.emt.traj')
opt.run(fmax=0.05)
```

This is the first ASE script we have seen so far, so 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

`EMT`

calculator is created. A*calculator*can evaluate quantities such as energies and forces on a collection of atoms. There are different kinds of calculators, and EMT is a particularly simple one. The calculator is associated with the`Atoms`

object by calling`atoms.calc = calc`

.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.emt.traj`

in the working
directory. Use the command **ase gui** to view the
trajectory file, showing each step of the optimization.

## Structure optimization of H_{2}O with EMT and GPAW¶

Adapt the above script as needed and calculate the structure of a
H_{2}O molecule using the EMT 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 empirical EMT potential is fast, but not very accurate for
molecules in particular. We therefore want to perform this
calculation in GPAW instead. GPAW uses real-space grids to represent
density and wavefunctions, and the grids exist in a cell. For this
reason you must set a cell for the `Atoms`

object. As a
coarse value let us use a 6 Ångström cell:

```
atoms.set_cell((6.0, 6.0, 6.0))
atoms.center()
```

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.

Instead of importing and using EMT, we now use GPAW:

```
from gpaw import GPAW
...
calc = GPAW(mode='fd')
...
```

Make a copy of your script and adapt it to GPAW, then recalculate the
structure of H_{2}O (make sure to choose a new filename for the
trajectory file).

During the calculation a lot of text is printed to the terminal. 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 and compare bond lenghts and bond
angles to the EMT result. Bond lengths and angles are shown
automatically if you select two or three atoms at a time.

## Atomization energies¶

Now that we know the structure of H_{2}O, 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 H_{2} 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='fd', 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 H_{2}O.

## 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=`

, where *name*)

is a string such as
*name*`'LDA'`

, `'PBE'`

or `'RPBE'`

.

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