# Lab 5: Forces, Stresses and Structures

Back to Course Overview

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.

• 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.

• 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.

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.