# Tutorial: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 %`LatticeVectors`

0.0 | 0.5 | 0.5 0.5 | 0.0 | 0.5 0.5 | 0.5 | 0.0 % a = 10.18 %`LatticeParameters`

a | a | a % %`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:

: this input variable must be set equal to the number of dimensions you want to consider as periodic. Since the system is 3D (`PeriodicDimensions`

= 3is 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.`Dimensions`

= 3

and`LatticeVectors`

: these two blocks are used to define the primitive lattice vectors that determine the unit cell.`LatticeParameters`

defines the direction of the vectors, while`LatticeVectors`

`LatticeParameters`

defines their length.

: the position of the atoms inside the unit cell, in reduced coordinates.`ReducedCoordinates`

: this specifies the`KPointsGrid`

*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 variable`nk = 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.

: 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`KPointsUseSymmetries`

= yes*k*-points is adjusted.

: we ask the code to output the density of states (`Output`

`dos`).

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. The allowed box shapes available for periodic systems is currently limited.

## 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 our system is indeed being treated as periodic in 3 dimensions.

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

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.

******************************** 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 **********************************************************************

The usual information about the grid is also printed. Note that the printed lengths of the simulation box are *along* the directions defined by the lattice vectors and they are equal to their norm.

The rest of the output is much like its non-periodic counterpart. After a few iterations the code should converge:

*********************** SCF CYCLE ITER # 11 ************************ etot = -7.92846745E+00 abs_ev = 1.20E-06 rel_ev = 4.90E-06 ediff = 1.43E-06 abs_dens = 1.66E-06 rel_dens = 2.08E-07 Matrix vector products: 79 Converged eigenvectors: 5 # State KPoint Eigenvalue [H] Occupation Error 1 1 -0.257104 2.000000 (4.2E-08) 1 2 -0.184780 2.000000 (9.4E-08) 2 2 -0.079246 2.000000 (8.4E-08) 2 1 0.010870 2.000000 (8.7E-08) 3 2 0.023244 2.000000 (7.3E-08) 4 2 0.073686 2.000000 (1.3E-07) 3 1 0.128009 2.000000 (1.0E-07) 4 1 0.128009 2.000000 (1.0E-07) 5 2 0.209030 0.000000 (1.4E-07) 5 1 0.232109 0.000000 (1.2E-07) Density of states: ----------%--------------%--------------%------%------------------%--- ----------%--------------%--------------%------%------------------%--- ----------%--------------%--------------%------%------------------%--- ----------%--------------%--------------%------%------%-----------%--- ----------%--------------%--------------%------%------%-----------%--- ----------%--------------%--------------%------%------%-----------%--- ----------%--------------%--------------%------%------%-----------%--- %---------%--------------%------------%-%------%------%-----------%--% %---------%--------------%------------%-%------%------%-----------%--% %---------%--------------%------------%-%------%------%-----------%--% ^ Elapsed time for SCF step 11: 0.05 ********************************************************************** Info: Writing states. 2021/08/27 at 08:16:18 Info: Finished writing states. 2021/08/27 at 08:16:18 Info: SCF converged in 11 iterations

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

**directory:**

`static`: the band-resolved density of states (DOS);`dos-XXXX.dat`

: the total DOS (summed over all bands);`total-dos.dat`

: the Fermi Energy in a format compatible with`total-dos-efermi.dat`.`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" for nkpt in $list do sed -i "s/nk = [0-9]\+/nk = $nkpt/g" 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:

#nk Total Energy 2 -7.92846745 4 -7.93455130 6 -7.93464987 8 -7.93465246 10 -7.93465314

As you can see, the total energy is converged to within 0.0001 hartree for `nk = 6`.

## 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 %`LatticeVectors`

0.0 | 0.5 | 0.5 0.5 | 0.0 | 0.5 0.5 | 0.5 | 0.0 % a = 10.18 %`LatticeParameters`

a | a | a % %`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:

: we are now performing a non-self-consistent calculation, so we use the`CalculationMode`

= unocc`unoccupied`calculation mode;

: this is the number of unoccupied bands to calculate;`ExtraStates`

= 10

: 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.`ExtraStatesToConverge`

= 5

: 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`KPointsPath`

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

: we have turned off the use of symmetries.`KPointsUseSymmetries`

= no

After running Octopus with this input file, you should obtain a file named ** bandstructure** inside the

**directory. This is how the first few lines of the file should look like:**

`static`# coord. kx ky kz (red. coord.), bands: 14 [H] 0.00000000 0.50000000 0.00000000 0.00000000 -0.19984690 -0.10463006 0.11059627 0.11059633 0.21320569 0.27708564 0.27708566 0.43517452 0.54466924 0.54466939 0.57089088 0.57496835 0.57496840 0.63339744 0.02640129 0.45000000 -0.00000000 -0.00000000 -0.20506433 -0.09711225 0.11125646 0.11125651 0.21395899 0.27788955 0.27788958 0.43640704 0.51951019 0.51951030 0.55510510 0.60282772 0.60282781 0.64782388 0.05280258 0.40000000 -0.00000000 -0.00000000 -0.21725579 -0.07804283 0.11323475 0.11323480 0.21620735 0.28007592 0.28007595 0.43967735 0.48696280 0.48696290 0.52612948 0.64322333 0.64322345 0.66832944 ...

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

On the right 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.