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
andforc_conv_thr
. These values are used to determine when the optimization finishes (listed followingcriteria:
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.