Cell relaxation in solids
In this tutorial, we will learn how to determine the lattice parameters of a bulk material, as well as getting the relaxed atomic coordinates in bulk materials. Different flavors of cell dynamics are available in Octopus, which we will explore here.
Determination of the cell shape and cell volume
In this tutorial, we will start by determining the equilibrium lattice parameters of bulk silicon. The input file is the following:
CalculationMode = go
BoxShape = parallelepiped
PeriodicDimensions = 3
FromScratch = yes
a = 10
Spacing = 0.5
%LatticeParameters
a | a | a
%
%LatticeVectors
0.0 | 0.5 | 0.5
0.5 | 0.0 | 0.5
0.5 | 0.5 | 0.0
%
%ReducedCoordinates
"Si" | 0.0 | 0.0 | 0.0
"Si" | 1/4 | 1/4 | 1/4
%
%KPointsGrid
2 | 2 | 2
0.5 | 0.5 | 0.5
0.5 | 0.0 | 0.0
0.0 | 0.5 | 0.0
0.0 | 0.0 | 0.5
%
KPointsUseSymmetries = yes
Eigensolver = chebyshev_filter
ConvRelDens = 1e-9
ExtraStates = 4
GOType = cell_shape
SCFCalculateStress = yes
ExperimentalFeatures = yes
We are not doing here a GS calculation, but a GO calculation (that stands for geometry optimization). The type of geometry optimization is determined by the variable GOType. In the above example, we are asking for optimizing the shape of the lattice vector, meaning both their directions and lengths. We could have equally asked for optimizing only the cell volume, as we started from an already correct space group. This has the advantage to allow using symmetries, which makes the calculation faster and more accurate. Note that several options exists for controlling the behavior of the geometry optimizer, in particular GOMethod, GOMaxIter, GOTolerance, and GOStep.
Running Octopus, we will perform consecutive cycles of ground-state calculations, after which forces and stress tensor will be computed. This information is then used to generate new atomic positions and lattice vectors. Based on this, a new real space grid is regenerated. The number of grid point is however kept constant, so it is adviced to start from a reasonably-size unit cell before performing a geometry relaxation with Octopus. It is important to note that accurate forces and stress tensors are critical in order to get precise results. This is why we request here a tight convergence of the relative density.
After every iteration the lattice vectors of the newly generated cell are displayed in the output. At each step, the geometry of the system is also printed in the geom/ folder, containing for the structure in XYZ format geom/go.0001.xyz , a modified XYZ format with reduced coordinates geom/go.0001.xyz_red , and a XCrysden file geom/go.0001.xsf . The file geom/optimization.log contains the evolution of the energy, force, as well as cell parameters versus the GO iterations.
The effect of the relaxation is clear. We started with the original lattice vectors
Lattice Vectors [b]
0.000000 5.000000 5.000000
5.000000 0.000000 5.000000
5.000000 5.000000 0.000000
Lattice Vectors [b]
0.000000 5.076457 5.076457
5.076457 0.000000 5.076457
5.076457 5.076457 0.000000
What about the convergence of this value? Because of the real-space grid and in particular the so-called egg-box effect, special care is required for the convergence of forces and stress in grid spacing. To illustrate this, we can repeat the previous calculation for a grid spacing of 0.4 Bohr, yielding a lattice parameter of 10.173 Bohr, showing that the above values was not fully converged. As an exercise you might want to converge the grid spacing. Just keep in mind that calculations with smaller grid spacing might take quite some time to run.
Full geometry relaxation of bulk tellurium
We now look at the more challenging relaxation of a tellurium crystal, in which the crystal symmetry does not fix the reduced coordinates of all atoms. In particular, in this system, the atomic coordinates are defined by an internal coordinate, defined as uu in the following input file, that we want to determine alongside with the lattice parameters. This calculation is a bit more demanding than for silicon as more degrees of freedom are involved, and will require few minutes to be completed.
CalculationMode = go
ExperimentalFeatures = yes
FromScratch = yes
PeriodicDimensions = 3
%LatticeParameters
4.4*angstrom | 4.4*angstrom | 5.9*angstrom
90 | 90 | 120
%
PseudopotentialSet = hgh_lda
uu = 0.24
%ReducedCoordinates
"Te" | uu | 0 | 0
"Te" | 0.0 | uu | 1/3
"Te" |-uu |-uu | 2/3
%
%KPointsGrid
3 | 3 | 2
%
KPointsUseSymmetries = yes
ExtraStates = 5
Eigensolver = chebyshev_filter
Spacing = 0.5
ConvRelDens = 1e-9
SCFCalculateStress = yes
GOType = cell_volume + ion_positions
Here, we prepared the input file such that the space group is known. Hence, we will only optimize the length of the lattice vectors and the atomic coordiates (GOType)=ion_positions + cell_volume.
At the end of the calculation, we obtain
Lattice Vectors [b]
8.274621 0.000000 0.000000
-4.137310 7.166032 0.000000
0.000000 0.000000 11.049474
Te 0.274231 0.000000 0.000000
Te 0.000000 0.274231 0.333333
Te -0.274231 -0.274231 -0.333333
We have thus converged to a relaxed internal coordinate of u = 0.274, with lattice parameters a = 4.448 Angstrom and c = 6.00 Angstrom. The experimental lattice parameters are u=0.2633, a = 4.451 Angstrom and c = 5.926 Angstrom, see [ Keller et al., PRB 16, 4404 (1977)]. The disagreement with experimental values is larger here, which comes from several factors that are not only due to convergence of the simulation parameters: spin-orbit coupling was not included in the simulation, and LDA is not good at describing the insulating nature of the material. More advanced functionals should therefor be employed for properly relaxing bulk Te. As an exercise you might want to test the effect for DFT+U on the equilibrium property of Te, using for instance a value of U=5eV on the p orbitals for Te.
%Species
"Te" | species_pseudo | hubbard_l | 1 | hubbard_u | 5*eV
%
DFTULevel = dft_u_empirical
The details on how to run the DFT+U calculations are given in the corresponding tutorial.