Getting started with periodic systems
The extension of a ground-state calculation to a periodic system is quite straightforward in Octopus. In this tutorial we will explain how to perform some basic calculation using bulk silicon as an example.
Input
As always, we will start with a simple input file. In this case we will use a primitive cell of Si, composed of two atoms.
CalculationMode = gs
PeriodicDimensions = 3
Spacing = 0.5
a = 10.18
%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
%
nk = 2
%KPointsGrid
nk | nk | nk
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
ExtraStates = 1
%Output
dos
%
Lets see more in detail some of the input variables:
-
PeriodicDimensions = 3
: this input variable must be set equal to the number of dimensions you want to consider as periodic. Since the system is 3D (Dimensions = 3
is the default), by setting this variable to 3 we impose periodic boundary conditions at all borders. This means that we have a fully periodic infinite crystal. -
LatticeVectors
andLatticeParameters
: these two blocks are used to define the primitive lattice vectors that determine the unit cell.LatticeVectors
defines the direction of the vectors, while LatticeParameters defines their length. -
ReducedCoordinates
: the position of the atoms inside the unit cell, in reduced coordinates. -
KPointsGrid
: this specifies the ‘‘k’'-point grid to be used in the calculation. Here we employ a 2x2x2 Monkhorst-Pack grid with four shifts. The first line of the block defines the number of ‘‘k’'-points along each axis in the Brillouin zone. Since we want the same number of points along each direction, we have defined the auxiliary variablenk = 2
This will be useful later on to study the convergence with respect to the number of ‘‘k’'-points. The other four lines define the shifts, one per line, expressed in reduced coordinates of the Brillouin zone. Alternatively, one can also define the reciprocal-space mesh by explicitly setting the position and weight of each ‘‘k’'-point using the KPoints or KPointsReduced variables. -
KPointsUseSymmetries = yes
: this variable controls if symmetries are used or not. When symmetries are used, the code shrinks the Brillouin zone to its irreducible portion and the effective number of ‘‘k’'-points is adjusted. -
Output = dos
: we ask the code to output the density of states.
Here we have taken the value of the grid spacing to be 0.5 bohr. Although we will use this value throughout this tutorial, remember that in a real-life calculation the convergence with respect to the grid spacing must be performed for all quantities of interest.
Note that for periodic systems the default value for the BoxShape
variable is parallelepiped
, although in this case the name can be misleading, as the actual shape also depends on the lattice vectors. This is the only box shape currently available for periodic systems.
Output
Now run Octopus using the above input file. Here are some important things to note from the output.
******************************* Space ********************************
Octopus will run in 3 dimension(s).
Octopus will treat the system as periodic in 3 dimension(s).
**********************************************************************
This tells us that out system is indeed being treated as periodic in 3 dimensions.
******************************** Grid ********************************
Simulation Box:
Type = parallelepiped
Lengths [b] = ( 7.198, 7.198, 7.198)
Main mesh:
Spacing [b] = ( 0.514, 0.514, 0.514) volume/point [b^3] = 0.09612
# inner mesh = 2744
# total mesh = 9192
Grid Cutoff [H] = 18.666387 Grid Cutoff [Ry] = 37.332774
**********************************************************************
Here Octopus outputs some information about the cell in real and reciprocal space.
***************************** Symmetries *****************************
Space group No. 227
International: Fd-3m
Schoenflies: Oh^7
Index Rotation matrix Fractional translations
1 : 1 0 0 0 1 0 0 0 1 0.000000 0.000000 0.000000
2 : 0 -1 1 0 -1 0 1 -1 0 0.000000 0.000000 0.000000
3 : -1 0 0 -1 0 1 -1 1 0 0.000000 0.000000 0.000000
4 : 0 1 -1 1 0 -1 0 0 -1 0.000000 0.000000 0.000000
5 : -1 1 0 -1 0 0 -1 0 1 0.000000 0.000000 0.000000
6 : 0 0 1 1 0 0 0 1 0 0.000000 0.000000 0.000000
7 : 1 -1 0 0 -1 1 0 -1 0 0.000000 0.000000 0.000000
8 : 0 0 -1 0 1 -1 1 0 -1 0.000000 0.000000 0.000000
9 : 0 -1 0 1 -1 0 0 -1 1 0.000000 0.000000 0.000000
10 : -1 0 1 -1 1 0 -1 0 0 0.000000 0.000000 0.000000
11 : 0 1 0 0 0 1 1 0 0 0.000000 0.000000 0.000000
12 : 1 0 -1 0 0 -1 0 1 -1 0.000000 0.000000 0.000000
13 : 0 0 -1 1 0 -1 0 1 -1 0.000000 0.000000 0.000000
14 : -1 1 0 -1 0 1 -1 0 0 0.000000 0.000000 0.000000
15 : 1 -1 0 0 -1 0 0 -1 1 0.000000 0.000000 0.000000
16 : 0 0 1 0 1 0 1 0 0 0.000000 0.000000 0.000000
17 : 1 0 -1 0 1 -1 0 0 -1 0.000000 0.000000 0.000000
18 : 0 -1 0 0 -1 1 1 -1 0 0.000000 0.000000 0.000000
19 : 0 1 0 1 0 0 0 0 1 0.000000 0.000000 0.000000
20 : -1 0 1 -1 0 0 -1 1 0 0.000000 0.000000 0.000000
21 : 0 1 -1 0 0 -1 1 0 -1 0.000000 0.000000 0.000000
22 : 1 0 0 0 0 1 0 1 0 0.000000 0.000000 0.000000
23 : -1 0 0 -1 1 0 -1 0 1 0.000000 0.000000 0.000000
24 : 0 -1 1 1 -1 0 0 -1 0 0.000000 0.000000 0.000000
Info: The system has 24 symmetries that can be used.
**********************************************************************
This block tells us about the space-group and the symmetries found for the specified structure.
****************************** Lattice *******************************
Lattice Vectors [b]
0.000000 5.090000 5.090000
5.090000 0.000000 5.090000
5.090000 5.090000 0.000000
Cell volume = 263.7445 [b^3]
Reciprocal-Lattice Vectors [b^-1]
-0.617209 0.617209 0.617209
0.617209 -0.617209 0.617209
0.617209 0.617209 -0.617209
Cell angles [degree]
alpha = 60.000
beta = 60.000
gamma = 60.000
**********************************************************************
Here Octopus outputs some information about the unit cell in real and reciprocal space.
Checking if the generated full k-point grid is symmetric
2 k-points generated from parameters :
---------------------------------------------------
n = 2 2 2
s1 = 0.50 0.50 0.50
s2 = 0.50 0.00 0.00
s3 = 0.00 0.50 0.00
s4 = 0.00 0.00 0.50
index | weight | coordinates |
1 | 0.250000 | 0.250000 0.000000 0.000000 |
2 | 0.750000 | -0.250000 0.250000 0.250000 |
Next we get the list of the ‘‘k’'-points in reduced coordinates and their weights. Since symmetries are used, only two ‘‘k’'-points are generated. If we had not used symmetries, we would have 32 ‘‘k’'-points instead.
The rest of the output is much like its non-periodic counterpart. After a few iterations the code should converge:
*********************** SCF CYCLE ITER # 12 ************************
etot = -7.92933294E+00 abs_ev = 2.51E-07 rel_ev = 1.02E-06
ediff = 3.50E-07 abs_dens = 1.09E-06 rel_dens = 1.36E-07
Matrix vector products: 69
Converged eigenvectors: 10
# State KPoint Eigenvalue [H] Occupation Error
1 1 -0.257170 2.000000 ( 4.6E-08)
1 2 -0.184858 2.000000 ( 8.4E-08)
2 2 -0.079342 2.000000 ( 6.6E-08)
2 1 0.010753 2.000000 ( 6.1E-08)
3 2 0.023125 2.000000 ( 6.3E-08)
4 2 0.073567 2.000000 ( 8.7E-08)
3 1 0.127870 2.000000 ( 7.1E-08)
4 1 0.127870 2.000000 ( 7.1E-08)
5 2 0.208835 0.000000 ( 8.3E-08)
5 1 0.231963 0.000000 ( 6.8E-08)
Density of states:
----------%--------------%--------------%------%------------------%---
----------%--------------%--------------%------%------------------%---
----------%--------------%--------------%------%------------------%---
----------%--------------%--------------%------%------%-----------%---
----------%--------------%--------------%------%------%-----------%---
----------%--------------%--------------%------%------%-----------%---
----------%--------------%--------------%------%------%-----------%---
%---------%--------------%------------%-%------%------%-----------%--%
%---------%--------------%------------%-%------%------%-----------%--%
%---------%--------------%------------%-%------%------%-----------%--%
^
Elapsed time for SCF step 12: 0.03
**********************************************************************
Info: Writing states. 2023/06/07 at 14:52:45
Info: Finished writing states. 2023/06/07 at 14:52:45
Info: SCF converged in 12 iterations
Info: Number of matrix-vector products: 1066
Info: Writing output to static/
Info: Finished writing information to 'restart/gs'.
Calculation ended on 2023/06/07 at 14:52:45
Walltime: 01.168s
Octopus emitted 15 warnings.
As usual, the static/info file contains the most relevant information concerning the calculation. Since we asked the code to output the density of states, we also have a few new files in the static directory:
- dos-XXXX.dat : the band-resolved density of states (DOS);
- total-dos.dat : the total DOS (summed over all bands);
- total-dos-efermi.dat : the Fermi Energy in a format compatible with total-dos.dat .
Of course you can tune the output type and format in the same way you do in a finite-system calculation.
Convergence in k-points
Similar to the convergence in spacing, a convergence must be performed for the sampling of the Brillouin zone. To do this one must try different numbers of ‘‘k’'-points along each axis in the Brillouin zone. This can easily be done by changing the value of the nk auxiliary variable in the previous input file. You can obviously do this by hand, but this is something that can also be done with a script. Here is such a script, which we will name kpts.sh :
#!/bin/bash
echo "#nk Total Energy" > kpts.log
list="2 4 6 8 10 12"
for nkpt in $list
do
echo "$nkpt"
sed -i "s/nk = [0-9]\+/nk = $nkpt/" inp
octopus >& out-$nkpt
energy=`grep Total static/info | head -2 | tail -1 | cut -d "=" -f 2`
echo $nkpt $energy >> kpts.log
rm -rf restart
done
After running the script (source kpts.sh
), you should obtain a file called kpts.log
that should look like this:
As you can see, the total energy is converged to within 0.0001 hartree for nk = 6
.
You can now play with an extended range, e.g. from 2 to 12. You should then see that the total energy is converged to less than a micro hartree.
Band-structure
We now proceed with the calculation of the band-structure of Si. In order to compute a band-structure, we must perform a non-self-consistent calculation, using the density of a previous ground-state calculation. So the first step is to obtain the initial ground-state. To do this, rerun the previous input file, but changing the number of k-points to nk=6
. Next, modify the input file such that it looks like this:
CalculationMode = unocc
PeriodicDimensions = 3
Spacing = 0.5
a = 10.18
%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
%
ExtraStates = 10
ExtraStatesToConverge = 5
%KPointsPath
10 | 10 | 15
0.5 | 0.0 | 0.0 # L point
0.0 | 0.0 | 0.0 # Gamma point
0.0 | 0.5 | 0.5 # X point
1.0 | 1.0 | 1.0 # Another Gamma point
%
KPointsUseSymmetries = no
Here are the things we changed:
-
CalculationMode = unocc
: we are now performing a non-self-consistent calculation, so we use theunoccupied
calculation mode; -
ExtraStates = 10
: this is the number of unoccupied bands to calculate; -
ExtraStatesToConverge = 5
: the highest unoccupied states are very hard to converge, so we use this variable to specify how many unoccupied states are considered for the stopping criterion of the non-self-consistent run. -
KPointsPath
: this block is used to specify that we want to calculate the band structure along a certain path in the Brillouin zone. This replaces the KPointsGrid block. The first row describes how many ‘‘k’'-points will be used to sample each segment. The next rows are the coordinates of the ‘‘k’'-points from which each segment starts and stops. In this particular example, we chose the following path: L-Gamma, Gamma-X, X-Gamma using a sampling of 10-10-15 ‘‘k’'-points. -
KPointsUseSymmetries = no
: we have turned off the use of symmetries.
After running Octopus with this input file, you should obtain a file named bandstructure inside the static directory. This is how the first few lines of the file should look like:
# coord. kx ky kz (red. coord.), bands: 14 [H]
0.00000000 0.50000000 0.00000000 0.00000000 -0.19992301 -0.10472513 0.11047004 0.11047009 0.21305293 0.27686191 0.27686193 0.43506938 0.54432395 0.54432410 0.57080931 0.57480646 0.57480651 0.63337123
0.02640129 0.45000000 0.00000000 0.00000000 -0.20513963 -0.09720871 0.11112960 0.11112965 0.21380657 0.27766704 0.27766706 0.43630273 0.51921757 0.51921768 0.55502603 0.60261231 0.60261241 0.64779943
0.05280258 0.40000000 0.00000000 0.00000000 -0.21732917 -0.07814289 0.11310604 0.11310609 0.21605590 0.27985723 0.27985725 0.43957622 0.48667851 0.48667860 0.52604769 0.64299609 0.64299621 0.66831748
0.07920387 0.35000000 0.00000000 0.00000000 -0.23165974 -0.05243307 0.11638717 0.11638722 0.21975932 0.28269147 0.28269149 0.44256402 0.45780330 0.45780338 0.49751716 0.66341701 0.68308558 0.68308573
...
The first column is the coordinate of the ‘‘k’'-point along the path. The second, third, and fourth columns are the reduced coordinates of the ‘‘k’'-point. The following columns are the eigenvalues for the different bands. In this case there are 14 bands (4 occupied and 10 unoccupied).
Below you can see the plot of the band structure.
This plot shows the occupied bands (purple) and the first 5 unoccupied bands (green). If you are using gnuplot, you can obtain a similar plot with the following command:
plot for [col=5:5+9] 'static/bandstructure' u 1:(column(col)) w l notitle ls 1
Note that when using the KPointsPath input variable, Octopus will run in a special mode, and the restart information of the previous ground-state calculation will not be altered in any way. The code informs us about this just before starting the unoccupied states iterations:
Info: The code will run in band structure mode.
No restart information will be printed.