Difference between revisions of "Tutorial:Periodic systems"
(→Bandstructure) 
m (→Convergence in kpoints: making it work with BSD (MacOS) sed (https://unix.stackexchange.com/questions/401905/bsdsedvsgnusedandi)) 

Line 206:  Line 206:  
for nkpt in $list  for nkpt in $list  
do  do  
−  sed i "s/nk = [09]\+/nk = $nkpt/" inp  +  sed i '' "s/nk = [09]\+/nk = $nkpt/" inp 
octopus >& out$nkpt  octopus >& out$nkpt  
energy=`grep Total static/info  head 2  tail 1  cut d "=" f 2`  energy=`grep Total static/info  head 2  tail 1  cut d "=" f 2` 
Revision as of 16:01, 21 August 2018
The extension of a groundstate 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
= gsPeriodicDimensions
= 3Spacing
= 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
= yesExtraStates
= 1Output
= 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, whileLatticeParameters
defines their length.

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

KPointsGrid
: this specifies the kpoint grid to be used in the calculation. Here we employ a 2x2x2 MonkhorstPack grid with four shifts. The first line of the block defines the number of kpoints 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 kpoints. The other four lines define the shifts, one per line, expressed in reduced coordinates of the Brillouin zone. Alternatively, one can also define the reciprocalspace mesh by explicitly setting the position and weight of each kpoint using theKPoints
orKPointsReduced
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 kpoints 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 reallife 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.
******************************** Grid ******************************** Simulation Box: Type = parallelepiped Lengths [b] = ( 3.627, 3.627, 3.627) Octopus will run in 3 dimension(s). Octopus will treat the system as periodic in 3 dimension(s). Lattice Vectors [b] 0.000000 5.130000 5.130000 5.130000 0.000000 5.130000 5.130000 5.130000 0.000000 Cell volume = 270.0114 [b^3] ReciprocalLattice Vectors [b^1] 0.612396 0.612396 0.612396 0.612396 0.612396 0.612396 0.612396 0.612396 0.612396 Main mesh: Spacing [b] = ( 0.484, 0.484, 0.484) volume/point [b^3] = 0.08000 # inner mesh = 3375 # total mesh = 10639 Grid Cutoff [H] = 21.095389 Grid Cutoff [Ry] = 42.190778 **********************************************************************
Here Octopus outputs some information about the cell in real and reciprocal space.
***************************** Symmetries ***************************** Space group No. 227 International: Fd3m 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 spacegroup and the symmetries found for the specified structure.
2 kpoints 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 kpoints in reduced coordinates and their weights. Since symmetries are used, only two kpoints are generated. If we had not used symmetries, we would have 32 kpoints instead.
The rest of the output is much like its nonperiodic counterpart. After a few iterations the code should converge:
*********************** SCF CYCLE ITER # 8 ************************ etot = 7.92843852E+00 abs_ev = 4.39E06 rel_ev = 1.79E05 ediff = 1.19E06 abs_dens = 3.49E05 rel_dens = 4.37E06 Matrix vector products: 77 Converged eigenvectors: 10 # State KPoint Eigenvalue [H] Occupation Error 1 1 0.257101 2.000000 (6.5E07) 1 2 0.184777 2.000000 (9.3E07) 2 2 0.079243 2.000000 (9.4E07) 2 1 0.010873 2.000000 (8.4E07) 3 2 0.023247 2.000000 (8.2E07) 4 2 0.073689 2.000000 (8.6E07) 3 1 0.128013 2.000000 (9.0E07) 4 1 0.128013 2.000000 (9.6E07) 5 2 0.209032 0.000000 (9.3E07) 5 1 0.232110 0.000000 (9.2E07) Density of states:         % %%%%%%%%% ^ Elapsed time for SCF step 8: 0.10 ********************************************************************** Info: Writing states. 2018/08/06 at 17:15:24 Info: Finished writing states. 2018/08/06 at 17:15:24 Info: SCF converged in 8 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 static directory:
 dosXXXX.dat: the bandresolved density of states (DOS);
 totaldos.dat: the total DOS (summed over all bands);
 totaldosefermi.dat: the Fermi Energy in a format compatible with totaldos.dat.
Of course you can tune the output type and format in the same way you do in a finitesystem calculation.
Convergence in kpoints
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 kpoints 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" for nkpt in $list do sed i '' "s/nk = [09]\+/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:
#nk Total Energy 2 7.92843852 4 7.93452207 6 7.93462084 8 7.93462226
As you can see, the total energy is converged to within 0.0001 hartree for nk = 6.
Bandstructure
We now proceed with the calculation of the bandstructure of Si. In order to compute a bandstructure, we must perform a nonselfconsistent calculation, using the density of a previous groundstate calculation. So the first step is to obtain the initial groundstate. To do this, rerun the previous input file, but changing the number of kpoints to nk=6. Next, modify the input file such that it looks like this:
CalculationMode
= unoccPeriodicDimensions
= 3Spacing
= 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
= 10ExtraStatesToConverge
= 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 nonselfconsistent calculation, so we use the unoccupied 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 nonselfconsistent 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 theKPointsGrid
block. The first row describes how many kpoints will be used to sample each segment. The next rows are the coordinates of the kpoints from which each segment starts and stops. In this particular example, we chose the following path: LGamma, GammaX, XGamma using a sampling of 101015 kpoints.

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.19984750 0.10463050 0.11059562 0.11059567 0.21320524 0.27708534 0.27708537 0.43517458 0.54466879 0.54466894 0.57089044 0.57496796 0.57496801 0.63339690 0.02640129 0.45000000 0.00000000 0.00000000 0.20506492 0.09711271 0.11125581 0.11125586 0.21395854 0.27788925 0.27788927 0.43640709 0.51950978 0.51950989 0.55510464 0.60282729 0.60282738 0.64782337 0.05280258 0.40000000 0.00000000 0.00000000 0.21725636 0.07804331 0.11323409 0.11323414 0.21620689 0.28007561 0.28007563 0.43967737 0.48696242 0.48696251 0.52612903 0.64322288 0.64322300 0.66832900 ...
The first column is the coordinate of the kpoint along the path. The second, third, and fourth columns are the reduced coordinates of the kpoint. 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 groundstate 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.