To find the minimum energy position of an atom, we could can manually move it, calculating the total energy each time, effectively finding the position where the force on it is zero. Some examples of how you can do this are given here which you can take a look through in your own time if you’re interested. Instead, in this lab we’ll be looking at how forces and stresses can be calculated within DFT at essentially no extra cost via the Hellman-Feynman theorem, and how these can be used.

In Quantum Espresso you can enable the calculation of forces and stresses by
setting `tprnfor = .true.`

and `tstress = .true.`

respectively in the
`CONTROL`

section of the input file.

## Forces in Methane

As a first example, let’s look at methane and calculate how the forces
converge with planewave energy cut-off. We have set up a template input file
as before in
`01_forces/01_methane/CH4_base.in`

. Take a
look at this now. The only new settings here are the two additional variables
mentioned above.

We also have a simple script to run this template file with energy cut-off
values from 20 to 60 Ry as `auto_run.sh`

.
Run this now and take a look at one of the output files. You’ll see before the
final timing information a section that looks like the following:

```
Forces acting on atoms (cartesian axes, Ry/au):
atom 1 type 1 force = 0.00000620 0.00000000 0.00002841
atom 2 type 2 force = 0.00000137 0.00000000 0.01218625
atom 3 type 2 force = 0.01151561 0.00000000 -0.00407450
atom 4 type 2 force = -0.00576159 -0.00997884 -0.00407008
atom 5 type 2 force = -0.00576159 0.00997884 -0.00407008
Total force = 0.024421 Total SCF correction = 0.000247
Computing stress (Cartesian axis) and pressure
total stress (Ry/bohr**3) (kbar) P= -0.06
-0.00000044 0.00000000 0.00000000 -0.06 0.00 0.00
0.00000000 -0.00000044 0.00000000 0.00 -0.06 0.00
0.00000000 0.00000000 -0.00000044 0.00 0.00 -0.06
```

So we have a list of the components of the force acting on each atom along
Cartesian axes, Then a total force and total scf correction.
**Note** the `Total force`

listed here is the square root of the sum of *all*
of the force components squared rather than the sum of the magnitudes of the
individual forces on the atoms. I’m not sure why, but it’s likely because the
number is intended more as a guide to check overall accuracy. If the
`Total SCF correction`

is comparable to the `Total force`

it usually means
you need to try to better converge the SCF cycle (via `conv_thr`

).

Following this we can see the calculated stress and corresponding pressure on the unit cell. While this number doesn’t mean much in principle for a calculation on a molecule in a box, if these numbers are not all close to zero it indicates your box size is likely not big enough.

### Convergence

Now let’s look at the convergence of the forces. As we’ve been doing in
earlier labs, we can extract the total force as a function of the energy
cut-off using `awk`

:

```
awk '/kinetic-energy/{ecut=$4} /Total force/{print ecut, $4}' *out
```

This works by saving (space delimited) field number 4 as a variable `ecut`

on
a line containing the string `kinetic-energy`

, and then on a line containing
the text `Total force`

it outputs the value of this variable along with the
text in the fourth field.

We could modify this to also output the total energy as

```
awk '/kinetic-energy/{ecut=$4}
/!.*total/{etot=$5}
/Total force/{print ecut, etot, $4}' *out
```

(You can add line breaks for clarity within an awk command if it gets long).

And add in the total pressure also with

```
awk '/kinetic-energy/{ecut=$4}
/!.*total/{etot=$5}
/Total force/{totfor=$4}
/total.*stress/{print ecut, etot, totfor, $6}' *out
```

Try running this and saving the convergence to a file.

*Task*

- Plot the fractional difference with respect to the most well converged result for each of total energy, total force and pressure as a function of energy cut-off.
- What can you see about the rate of convergence of each of these parameters?

*Task - Forces in Carbon*

- Make a directory
`01_forces/02_diamond`

. - In this directory calculate the forces and stresses on a carbon diamond system, and how they converge with energy cut-off.
- What do you notice about the forces? Why do you think this is?
- Now make a directory
`01_forces/03_graphite`

and repeat this process. - How does it compare to the diamond structure?
- Make a similar analysis of k-point convergence for the diamond case.

## Optimizing Ionic Positions

Given that we can readily calculate forces, wouldn’t it be nice if the code could automatically use these forces and find the atomic positions where these forces are zero (or at least within some tolerance of zero).

In Quantum Espresso this type of calculation, where the atomic positions are
relaxed to their minimum energy positions can be performed by selecting
`calculation = 'relax'`

in the `CONTROL`

section. There are a number of
additional variables that you can specify to control this process:

`tprnfor`

is automatically set to`.true.`

so that forces are output.- You’ll need to add an
`IONS`

section to the input. This is the only mandatory addition. There are several variables that can be specified within this section (it can be empty if you’re happy with all defaults), which control the algorithm used to find the optimal ionic positions. Consult the PW documentation for details.

### CH4

The directory `02_structure/01_methane`

contains
an input file for methane. Run this and take
a look at the output file once the calculation finishes. The interesting
bit is right near the end, just before the timing output. It should look
something like the following:

```
Forces acting on atoms (cartesian axes, Ry/au):
atom 1 type 1 force = -0.00008477 0.00000000 -0.00021366
atom 2 type 2 force = 0.00001647 0.00000000 0.00008941
atom 3 type 2 force = -0.00008533 0.00000000 -0.00000690
atom 4 type 2 force = 0.00007681 0.00007855 0.00006557
atom 5 type 2 force = 0.00007681 -0.00007855 0.00006557
Total force = 0.000318 Total SCF correction = 0.000117
SCF correction compared to forces is large: reduce conv_thr to get better values
bfgs converged in 4 scf cycles and 3 bfgs steps
(criteria: energy < 1.0E-04 Ry, force < 1.0E-03 Ry/Bohr)
End of BFGS Geometry Optimization
Final energy = -15.8757498591 Ry
Begin final coordinates
ATOMIC_POSITIONS (angstrom)
C -0.000016545 0.000000000 0.000113593
H 0.000017994 0.000000000 1.113983413
H 1.050254342 0.000000000 -0.371443564
H -0.525128396 -0.909710844 -0.371326721
H -0.525128396 0.909710844 -0.371326721
End final coordinates
```

- We have the force output. They should all be pretty close to zero now.
- You’ll likely notice a warning that the SCF correction is large and to
reduce
`conv_thr`

. We’ll try this shortly. - Then it should say that our relaxation converged (
`bfgs`

is the name of the default algorithm used) and how many steps it took. - We have the final total energy output.
- Finally we have the optimized atomic positions. Since it has moved all the atoms a bit, it’s not immediately clear if the C-H bond is longer or shorter than it was before.

Let’s try again, and address some of these things. We have a second input
file in the directory
`02_structure/02_methane`

where we’ve made a few
changes. Use `diff`

to see what’s changed.

- We’ve added two new variables to the control section:
`etot_conv_thr`

and`forc_conv_thr`

. These values are used to determine when the optimization finishes (listed following`criteria:`

in the output above). We have lowered both by a factor of 10 from their defaults. - We’ve increased the plane-wave energy cut-off to 40 Ry. This should really be tested systematically, but this is just for illustrative purposes.
- We’ve also lowered the scf convergence criteria from the default by a factor of 1000.
- We’ve added some 1s and 0s following the atomic positions. These define
multiplying factors for the various calculated force components on an atom.
- By adding three 0s following the carbon position we ensure it will be fixed at 0,0,0.
- We only allow the first H atom to move along the z-axis.
- And we only allow the second H atom to move within the x-z plane.

*Task*

- Run
`pw.x`

with this input file and check the result. - Has the C-H bond been lengthened or shortened?
- How much lower is the total energy than the starting configuration.
- Perform this calculation for energy cut-offs of 10 and 50 Ry, and see how this affects predicted C-H bond length.
- What do you think would happen if you did this calculation for carbon diamond? What about graphite?

## Optimizing Unit Cells

In a similar way to using the calculated forces to optimize the ionic positions we can also uses the calculated stresses to optimize the unit cell. In doing this you should keep in mind that it can take a higher energy-cut off to converge the stress as we saw at the start of this lab. Also, if you recall the number of planewaves depends on the cell volume, so if during the course of the calculation we change the volume we may also change the number of plane waves, or if we fix the number of plane waves we are changing the effective energy cut-off (the latter for Quantum Espresso, but different codes will handle this differently). So you need to be quite careful about this when optimizing lattice vectors.

In Quantum Espresso we can do this variable-cell relaxation by setting
`calculation = 'vc-relax'`

in the `CONTROL`

section. We must additionally
specify both an `IONS`

section as previously, along with a `CELL`

section.

### Silicon

An example input file is given in the
directory `02_forces/03_silicon`

. Take a look at this
now. You’ll notice in addition to the inputs mentioned, there’s also a fairly
high energy cut-off, and we’ve lowered the SCF convergence threshold from
the default. Try running this now and let’s look at the output.

The output is a little different in this case, since at the end of the
optimization, an `scf`

calculation is automatically performed starting from the
optimized structure. This is because Quantum Espresso fixes the basis set as
that for the original input structure during the calculation, so if the
structure has changed a lot, you may calculate a different stress when
starting from the relaxed structure. You will be able to see from the final
stress or pressure whether you should rerun your calculation.

Additionally, the cell output will likely be in Quantum Espresso’s
representation where the cell vectors are given explicitly in Bohr along with
a scaling factor `alat`

which is fixed. (In our case here `alat`

will be `A`

converted from Angstrom to Bohr). If you want to rerun your calculation you
could either input the cell using these directly, or calculate appropriate
values for the input. You may need to do the latter if you want to find the
new lattice length anyway.

*Task*

So for the silicon case, you’ll see the components of the lattice vectors have changed from 0.5 to 0.497277 or something close to that. What is our predicted silicon lattice constant?

## Summary

In this lab we have seen

- How to output forces and stresses in our calculation, and how to check these quantities have converged.
- How to use these forces in a calculation to optimize the atomic positions, where they are moved until the forces on the atoms are less then some threshold value.
- How the stresses can be used to optimize the structure, where the lattice constant is changed until the stress on the system is less than some threshold value.

### Extra: Graphene

The folder `02_structure/04_graphene`

contains an
input file for graphene. The thing
to note in this case is the additional use of `cell_dofree = '2Dxy'`

in the
`CELL`

section. We have used this to say that we only want to optimize the
cell in the xy plane. The spacing between periodic images in the z-direction
should be large enough such that the interaction is small, but we otherwise
don’t care about the stresses in this direction.

*Task*

- Starting from the provided graphene input file, find the length of the C-C
bond in graphene to 3 significant figures.
- Be sure to test your result is converged with respect to plane wave energy cut-off and k-point sampling, along with the internal tolerances being sufficiently small.