Difference between revisions of "Tutorial:Wires and slabs"

From OctopusWiki
Jump to navigation Jump to search
 
(10 intermediate revisions by 4 users not shown)
Line 3: Line 3:
 
== Introduction ==
 
== Introduction ==
  
In the [[Tutorial:Periodic systems|Periodic systems]] tutorial, we saw that the {{variable|PeriodicDimensions|System}} input variable controls the number of dimensions to be considered as periodic. In that tutorial we only considered the case <tt>{{variable|PeriodicDimensions|System}} = 3</tt>. Lets now see in detail the different cases:
+
In the [[Tutorial:Getting started with periodic systems|Getting started with periodic systems]] tutorial, we saw that the {{variable|PeriodicDimensions|System}} input variable controls the number of dimensions to be considered as periodic. In that tutorial we only considered the case <tt>{{variable|PeriodicDimensions|System}} = 3</tt>. Let's now see in detail the different cases:
  
 
* <tt>{{variable|PeriodicDimensions|System}} = 0</tt> (which is the default) gives a finite system calculation, since Dirichlet zero boundary conditions are used at all the borders of the simulation box;
 
* <tt>{{variable|PeriodicDimensions|System}} = 0</tt> (which is the default) gives a finite system calculation, since Dirichlet zero boundary conditions are used at all the borders of the simulation box;
Line 27: Line 27:
 
}}
 
}}
 
</ref>. This is done automatically, since the value of {{variable|PoissonSolver|Hamiltonian}} in a periodic calculation is chosen according to {{variable|PeriodicDimensions|System}}. See also the variable documentation for {{variable|PoissonSolver|Hamiltonian}}.
 
</ref>. This is done automatically, since the value of {{variable|PoissonSolver|Hamiltonian}} in a periodic calculation is chosen according to {{variable|PeriodicDimensions|System}}. See also the variable documentation for {{variable|PoissonSolver|Hamiltonian}}.
 +
 +
== Atomic positions ==
 +
 +
An important point when dealing with semi-periodic systems is that it is better if the coordinates of the atoms along the aperiodic directions are centered around 0. For slabs, this means that the atoms of the slab must be centered around z=0. If this is not the case, the calculation might not use all the available symmetries and it might lead to an inefficient box.
  
 
== Sodium chain ==
 
== Sodium chain ==
Line 34: Line 38:
 
=== Ground-state ===
 
=== Ground-state ===
  
First we start we the ground-state calculation using the following input file:
+
We first start with the ground-state calculation using the following input file:
  
 
  {{variable|CalculationMode|Calculation_Modes}} = gs
 
  {{variable|CalculationMode|Calculation_Modes}} = gs
  {{variable|UnitsOutput|Execution}} = ev_angstrom
+
{{variable|UnitsOutput|Execution}} = ev_angstrom
 
  {{variable|ExperimentalFeatures|Execution}} = yes
 
  {{variable|ExperimentalFeatures|Execution}} = yes
 
   
 
   
Line 44: Line 48:
 
  {{variable|Spacing|Mesh}} = 0.3*angstrom
 
  {{variable|Spacing|Mesh}} = 0.3*angstrom
 
    
 
    
  %{{variable|Lsize|Mesh}}
+
  %{{variable|LatticeParameters|Mesh}}
   1.99932905*angstrom | 5.29*angstrom | 5.29*angstrom
+
   3.9986581*angstrom | 10.58*angstrom | 10.58*angstrom
 
  %
 
  %
 
   
 
   
Line 57: Line 61:
 
  {{variable|KPointsUseSymmetries|Mesh}} = yes
 
  {{variable|KPointsUseSymmetries|Mesh}} = yes
  
Most of these input variables were already introduced in the [[Tutorial:Periodic systems|Periodic systems]] tutorial. Just note that the ''k''-points are all along the first dimension, as that is the only periodic dimension.   
+
Most of these input variables were already introduced in the [[Tutorial:Getting started with periodic systems|Getting started with periodic systems]] tutorial. Just note that the ''k''-points are all along the first dimension, as that is the only periodic dimension.   
 +
 
 +
The output should be quite familiar, but there are some noteworthy differences.
 +
 
 +
<pre>
 +
******************************* Space ********************************
 +
Octopus will run in 3 dimension(s).
 +
Octopus will treat the system as periodic in 1 dimension(s).
 +
**********************************************************************
 +
</pre>
 +
Here we see that we are indeed running a 3D system that is periodic along one dimension.
 +
 
 +
<pre>
 +
****************************** Lattice *******************************
 +
  Lattice Vectors [A]
 +
    3.998658    0.000000    0.000000
 +
  Cell length =            3.9987 [A]
 +
  Reciprocal-Lattice Vectors [A^-1]
 +
    1.571323    0.000000    0.000000
 +
  Cell angles [degree]
 +
    alpha =  90.000
 +
    beta  =  90.000
 +
    gamma =  90.000
 +
**********************************************************************
 +
</pre>
 +
Although we specified three lattice parameters in the input file corresponding to three lattice vectors, here the code tells us that it's only using one lattice vector. This is because the lattice is only periodic along one direction and therefore it is fully determined by the first vector. The other vectors are ignored when generating the Bravais lattice, but they are still used for two other purposes. First, they allow to specify atomic positions in reduced coordinates, which is useful when copying the coordinates from other codes. Second, they specify the lenght of the paralleliped box used to generate the real-space grid. This is confirmed a few lines later:
 +
 
 +
<pre>
 +
******************************** Grid ********************************
 +
Simulation Box:
 +
  Type = parallelepiped
 +
  Lengths [A] = (  3.999,  10.580,  10.580)
 +
Main mesh:
 +
  Spacing [A] = ( 0.308, 0.300, 0.300)    volume/point [A^3] =      0.02768
 +
  # inner mesh =      15925
 +
  # total mesh =      34397
 +
  Grid Cutoff [eV] =  397.448444    Grid Cutoff [Ry] =    29.211924
 +
**********************************************************************
 +
</pre>
 +
Note how the lengths along 'y' and 'z' are the same as the lenghts of the corresponding lattice parameters specified in the input file.
  
The output should be quite familiar, but the following piece of output confirms that the code is indeed using a cutoff for the calculation of the Hartree potential, as mentioned in the Introduction:
 
 
<pre>
 
<pre>
 
****************************** Hartree *******************************
 
****************************** Hartree *******************************
 +
Input: [DressedOrbitals = no]
 
The chosen Poisson solver is 'fast Fourier transform'
 
The chosen Poisson solver is 'fast Fourier transform'
 
Input: [PoissonFFTKernel = cylindrical]
 
Input: [PoissonFFTKernel = cylindrical]
 
**********************************************************************
 
**********************************************************************
  
Input: [FFTLibrary = fftw]
 
 
Info: FFT grid dimensions      = 13 x 75 x 75
 
Info: FFT grid dimensions      = 13 x 75 x 75
 
       Total grid size          = 73125 (  0.6 MiB )
 
       Total grid size          = 73125 (  0.6 MiB )
Line 72: Line 114:
 
Info: Poisson Cutoff Radius    =  11.2 A
 
Info: Poisson Cutoff Radius    =  11.2 A
 
</pre>
 
</pre>
The cutoff used is a cylindrical cutoff and, by comparing the FFT grid dimensions with the size of the simulation box, we see that the Poisson solver is using a supercell doubled in size in the ''y'' and ''z'' directions.
+
Finally, this piece of output confirms that the code is indeed using a cutoff for the calculation of the Hartree potential, as mentioned in the Introduction. The cutoff used is a cylindrical cutoff and, by comparing the FFT grid dimensions with the size of the simulation box, we see that the Poisson solver is using a supercell doubled in size in the ''y'' and ''z'' directions.
  
 
At this point you might want to play around with the number of ''k''-points until you are sure the calculation is converged.
 
At this point you might want to play around with the number of ''k''-points until you are sure the calculation is converged.
Line 83: Line 125:
  
 
  {{variable|CalculationMode|Calculation_Modes}} = unocc
 
  {{variable|CalculationMode|Calculation_Modes}} = unocc
  {{variable|UnitsOutput|Execution}} = ev_angstrom
+
{{variable|UnitsOutput|Execution}} = ev_angstrom
 
  {{variable|ExperimentalFeatures|Execution}} = yes
 
  {{variable|ExperimentalFeatures|Execution}} = yes
 
   
 
   
Line 89: Line 131:
 
   
 
   
 
  {{variable|Spacing|Mesh}} = 0.3*angstrom
 
  {{variable|Spacing|Mesh}} = 0.3*angstrom
 +
 +
%{{variable|LatticeParameters|Mesh}}
 +
  3.9986581*angstrom | 10.58*angstrom | 10.58*angstrom
 +
%
 
    
 
    
%{{variable|Lsize|Mesh}}
 
  1.99932905*angstrom | 5.29*angstrom | 5.29*angstrom
 
%
 
 
 
  %{{variable|Coordinates|System}}
 
  %{{variable|Coordinates|System}}
 
   "Na" | 0.0 | 0.0 | 0.0
 
   "Na" | 0.0 | 0.0 | 0.0
Line 108: Line 150:
 
  {{variable|KPointsUseSymmetries|Mesh}} = no
 
  {{variable|KPointsUseSymmetries|Mesh}} = no
  
and run the code. You might notice the comments on the the LCAO in the output. What's going on? Why can't a full initialization with LCAO be done?
+
and run the code.
  
You can now plot the band structure using the data from the {{file|static/bandstructure}} file, just like in the [[Tutorial:Periodic systems|Periodic systems]] tutorial.
+
You can now plot the band structure using the data from the {{file|static/bandstructure}} file, just like in the [[Tutorial:Getting started with periodic systems|Getting started with periodic systems]] tutorial.
  
 
In the introduction we mentioned that performing a calculation with <tt>{{variable|PeriodicDimensions|System}} = 1</tt> is not quite the same as performing a <tt>{{variable|PeriodicDimensions|System}}= 3</tt> calculation with a large supercell. Lets now check this. Re-run both the ground-state and the unoccupied calculations, but setting <tt>{{variable|PeriodicDimensions|System}} = 3</tt> in the above input files. Before doing so, make sure you copy the {{file|static/bandstructure}} file to a different place so that it is not overwritten (better yet, run the new calculations in a different folder). You can see the plot of the two band structures on the right. More comments on this in ref. <ref name="cutoff"/>.
 
In the introduction we mentioned that performing a calculation with <tt>{{variable|PeriodicDimensions|System}} = 1</tt> is not quite the same as performing a <tt>{{variable|PeriodicDimensions|System}}= 3</tt> calculation with a large supercell. Lets now check this. Re-run both the ground-state and the unoccupied calculations, but setting <tt>{{variable|PeriodicDimensions|System}} = 3</tt> in the above input files. Before doing so, make sure you copy the {{file|static/bandstructure}} file to a different place so that it is not overwritten (better yet, run the new calculations in a different folder). You can see the plot of the two band structures on the right. More comments on this in ref. <ref name="cutoff"/>.
Line 116: Line 158:
 
== h-BN monolayer ==
 
== h-BN monolayer ==
  
Hexagonal boron nitride (h-BN) is an insulator widely studied which has a similar structure to graphene. Here we will describe how to get the band structure of an h-BN monolayer.
+
Hexagonal boron nitride (h-BN) is an insulator widely studied which has a similar structure to graphene. Here we will describe how to get the band structure of an h-BN monolayer.  
  
=== Ground state calculation ===
+
=== Ground-state calculation ===
A layer of h-BN is periodic in the x-y directions, but not in the z. Thus we will set {{variable|PeriodicDimensions|System}} = 2 . Here we set the bond length to 1.445*angstrom. The box size in the z direction is 2*L with 'L'  large enough to describe a monolayer in the vacuum. One should always converge the box length value. Here is the inp file for the GS calculation:  
+
We will start by calculating the ground-state using the following input file:
  
 
  {{variable|CalculationMode|Calculation_Modes}} = gs
 
  {{variable|CalculationMode|Calculation_Modes}} = gs
+
  {{variable|UnitsOutput|Execution}} = ev_angstrom
  {{variable|FromScratch|Execution}} = yes
 
 
 
  {{variable|ExperimentalFeatures|Execution}} = yes
 
  {{variable|ExperimentalFeatures|Execution}} = yes
 
   
 
   
Line 130: Line 170:
 
   
 
   
 
  {{variable|Spacing|Mesh}} = 0.20*angstrom
 
  {{variable|Spacing|Mesh}} = 0.20*angstrom
 
{{variable|BoxShape|Mesh}} = parallelepiped
 
 
   
 
   
 
  BNlength = 1.445*angstrom
 
  BNlength = 1.445*angstrom
 
  a = sqrt(3)*BNlength
 
  a = sqrt(3)*BNlength
  L=40
+
  L = 40
 
 
  %{{variable|LatticeParameters|Mesh}}
 
  %{{variable|LatticeParameters|Mesh}}
 
   a | a | L
 
   a | a | L
Line 142: Line 179:
 
   
 
   
 
  %{{variable|LatticeVectors|Mesh}}
 
  %{{variable|LatticeVectors|Mesh}}
  1   | 0         | 0.
+
  1.0 | 0.0      | 0.0
   -1/2 | sqrt(3)/2 | 0.
+
   -1/2 | sqrt(3)/2 | 0.0
  0.   | 0.       | 1.
+
  0.0 | 0.0      | 1.0
 
  %
 
  %
 
    
 
    
 
  %{{variable|ReducedCoordinates|System}}
 
  %{{variable|ReducedCoordinates|System}}
   'B' | 0.0 | 0.0 | 0.00
+
   'B' | 0.0 | 0.0 | 0.0
   'N' | 1/3 | 2/3 | 0.00
+
   'N' | 1/3 | 2/3 | 0.0
 
  %  
 
  %  
 
   
 
   
Line 159: Line 196:
 
   12  | 12  | 1
 
   12  | 12  | 1
 
  %
 
  %
+
  {{variable|KPointsUseSymmetries|Mesh}} = yes
  {{variable|ExtraStates|States}} = 2
+
 
+
Most of the file should be self-explanatory, but here is a more detailed explanation for some of the choices:
{{variable|UnitsOutput|Execution}} = ev_angstrom
+
 
 +
* <tt>{{variable|PeriodicDimensions|System}} = 2</tt>: A layer of h-BN is periodic in the ''x''-''y'' directions, but not in the ''z'' direction, so there are two periodic dimensions.
  
[[Image:Tutorial_band_structure_HBN.jpg|500px|right|Band structure for a spacing of 0.20 angstrom of a monolayer h-BN.]]
+
* <tt>{{variable|LatticeParameters|Mesh}}</tt>: Here we have set the bond length to 1.445 Å. The box size in the ''z'' direction is ''2 L'' with ''L'' large enough to describe a monolayer in the vacuum. Remember that one should always check the convergence of any quantities of interest with the box length value.
  
=== Band Structure ===
+
If you now run the code, you will notice that the cutoff used for the calculation of the Hartree potential is different than for the Sodium chain, as is to be expected:
 +
<pre>
 +
****************************** Hartree *******************************
 +
The chosen Poisson solver is 'fast Fourier transform'
 +
Input: [PoissonFFTKernel = planar]
 +
**********************************************************************
  
[[Image:Tutorial_band_gap_convergence_HBN.jpg|500px|right|Convergence of the band gap with respect to the spacing for an h-BN monolayer.]]
+
Input: [FFTLibrary = fftw]
 +
Info: FFT grid dimensions      = 13 x 13 x 225
 +
      Total grid size          = 38025 (  0.3 MiB )
 +
      Inefficient FFT grid. A better grid would be: 14 14 225
 +
Info: Poisson Cutoff Radius    =  22.5 A
 +
</pre>
  
After this GS calculation, we will perform an unocc run. This non-self consistent calculation which needs the density from the previous GS calculation.
+
=== Band Structure ===
 +
After the ground-state calculation, we will now calculate the band structure. This is the input for the non-self consistent calculation:
  
 
  {{variable|CalculationMode|Calculation_Modes}} = unocc
 
  {{variable|CalculationMode|Calculation_Modes}} = unocc
 +
{{variable|UnitsOutput|Execution}} = ev_angstrom
 +
{{variable|ExperimentalFeatures|Execution}} = yes
 +
 +
{{variable|PeriodicDimensions|System}} = 2
 +
 +
{{variable|Spacing|Mesh}} = 0.20*angstrom
 +
 +
BNlength = 1.445*angstrom
 +
a = sqrt(3)*BNlength
 +
L = 40
 +
%{{variable|LatticeParameters|Mesh}}
 +
  a | a | L
 +
%
 +
 +
%{{variable|LatticeVectors|Mesh}}
 +
  1.0 | 0.0      | 0.0
 +
  -1/2 | sqrt(3)/2 | 0.0
 +
  0.0 | 0.0      | 1.0
 +
%
 +
 
 +
%{{variable|ReducedCoordinates|System}}
 +
  'B' | 0.0 | 0.0 | 0.0
 +
  'N' | 1/3 | 2/3 | 0.0
 +
%
 +
 +
{{variable|PseudopotentialSet|System}}=hgh_lda
 +
 +
{{variable|LCAOStart|SCF}}=lcao_states
 +
 +
{{variable|ExtraStatesToConverge|States}}  = 4
 +
{{variable|ExtraStates|States}}  = 8
 
   
 
   
{{variable|ExtraStatesToConverge|States}}  = 5
 
{{variable|ExtraStates|States}}  = 10
 
 
Here the number of {{variable|ExtraStates|States}} is the number of unoccupied bands in the final band structure. The highest occupied states are extremely complicated to converge. Therefore the variable {{variable|ExtraStatesToConverge|States}} specifies how many unoccupied states are considered for stopping criterion of the non-self-consistent run.
 
 
In order to calculate the band structure along a certain path along the BZ, we will use the variable {{variable|KPointsPath|Mesh}} . Instead of using the {{variable|KPointsGrid|Mesh}} block of the GS calculation, we use during this unocc calculation:
 
 
 
  %{{variable|KPointsPath|Mesh}}
 
  %{{variable|KPointsPath|Mesh}}
   12  | 7 | 12                 # Number of k point to sample each path
+
   12  | 7   | 12   # Number of k point to sample each path
   0   | 0 | 0                 # Reduced coordinate of the 'Gamma' k point
+
   0.0 | 0.0 | 0.0  # Reduced coordinate of the 'Gamma' k point
   1/3 | 1/3 | 0                 # Reduced coordinate of the 'K' k point
+
   1/3 | 1/3 | 0.0  # Reduced coordinate of the 'K' k point
   1/2 | 0 | 0                 # Reduced coordinate of the 'M' k point
+
   1/2 | 0.0 | 0.0  # Reduced coordinate of the 'M' k point
   0   | 0 | 0                 # Reduced coordinate of the 'Gamma' k point
+
   0.0 | 0.0 | 0.0  # Reduced coordinate of the 'Gamma' k point
 
  %
 
  %
 +
{{variable|KPointsUseSymmetries|Mesh}} = no
  
 +
[[Image:Tutorial_band_structure_HBN.png|thumb|420px|right|Band structure of a monolayer h-BN.]]
  
The first row describes how many k points will be used to sample each segment. The next rows are the coordinate of k points from which each segment start and stop. In this particular example, we chose the path: Gamma-K, K-M, M-Gamma using 12-7-12 k points. The output band structure is written in static/ band structure. In Figure 1 is plotted the output band structure where blue lines represent the occupied states and the red ones the unoccupied ones.
+
In this case, we chose the following path to calculate the band structure: Gamma-K, K-M, M-Gamma, with a sampling of 12-7-12 ''k''-points.
 
 
Info: The code will run in band structure mode.
 
      No restart information will be printed.
 
 
 
Moreover, by using  {{variable|KPointsPath|Mesh}}, the wave function obtained during the previous GS calculation (and stored in the restart/ directory) will not be affected by this calculation. 
 
  
One should also make sure that the calculation is converged with respect to the spacing. Figure 2 shows the band gap for several spacing values. We find that a spacing of 0.14 Angstrom is needed in order to converge the band gap up to 0.01 eV.
+
On the right is the resulting plot of the occupied bands and the first four unoccupied bands from the {{file|static/bandstructure}} file.
  
 
==References==
 
==References==
Line 203: Line 273:
 
<references/>
 
<references/>
  
{{Tutorial_foot|series=|prev=|next=}}
+
{{Tutorial_foot|series=Periodic systems|prev=Getting started with periodic systems|next=Optical spectra of solids}}
  
 
[[Category:Advanced]]
 
[[Category:Advanced]]

Latest revision as of 09:30, 9 September 2021

In this tutorial we will explain how to use the flexibility of the real-space grid to treat systems that are periodic in only one or two dimensions. As examples we will use a Na chain and a hexagonal boron nitride (h-BN) monolayer.

Introduction

In the Getting started with periodic systems tutorial, we saw that the PeriodicDimensions input variable controls the number of dimensions to be considered as periodic. In that tutorial we only considered the case PeriodicDimensions = 3. Let's now see in detail the different cases:

  • PeriodicDimensions = 0 (which is the default) gives a finite system calculation, since Dirichlet zero boundary conditions are used at all the borders of the simulation box;
  • PeriodicDimensions = 1 means that only the x axis is periodic, while in all the other directions the system is confined. This value must be used to simulate, for instance, a single infinite wire.
  • PeriodicDimensions = 2 means that both x and y axis are periodic, while zero boundary conditions are imposed at the borders crossed by the z axis. This value must be used to simulate, for instance, a single infinite slab.
  • PeriodicDimensions = 3 means that the simulation box is a primitive cell for a fully periodic infinite crystal. Periodic boundary conditions are imposed at all borders.


It is important to understand that performing, for instance, a PeriodicDimensions = 1 calculation in Octopus is not quite the same as performing a PeriodicDimensions = 3 calculation with a large supercell. In the infinite-supercell limit the two approaches reach the same ground state, but this does not hold for the excited states of the system.

Another point worth noting is how the Hartree potential is calculated for periodic systems. In fact the discrete Fourier transform that are used internally in a periodic calculation would always result in a 3D periodic lattice of identical replicas of the simulation box, even if only one or two dimensions are periodic. Fortunately Octopus includes a clever system to exactly truncate the long-range part of the Coulomb interaction, in such a way that we can effectively suppress the interactions between replicas of the system along non-periodic axes [1]. This is done automatically, since the value of PoissonSolver in a periodic calculation is chosen according to PeriodicDimensions. See also the variable documentation for PoissonSolver.

Atomic positions

An important point when dealing with semi-periodic systems is that it is better if the coordinates of the atoms along the aperiodic directions are centered around 0. For slabs, this means that the atoms of the slab must be centered around z=0. If this is not the case, the calculation might not use all the available symmetries and it might lead to an inefficient box.

Sodium chain

Let us now calculate some bands for a simple single Na chain (i.e. not a crystal of infinite parallel chains, but just a single infinite chain confined in the other two dimensions).

Ground-state

We first start with the ground-state calculation using the following input file:

CalculationMode = gs
UnitsOutput = ev_angstrom
ExperimentalFeatures = yes

PeriodicDimensions = 1

Spacing = 0.3*angstrom
 
%LatticeParameters
 3.9986581*angstrom | 10.58*angstrom | 10.58*angstrom
%

%Coordinates
 "Na" | 0.0 | 0.0 | 0.0
%

%KPointsGrid
 9 | 1 | 1
%
KPointsUseSymmetries = yes

Most of these input variables were already introduced in the Getting started with periodic systems tutorial. Just note that the k-points are all along the first dimension, as that is the only periodic dimension.

The output should be quite familiar, but there are some noteworthy differences.

******************************* Space ********************************
Octopus will run in 3 dimension(s).
Octopus will treat the system as periodic in 1 dimension(s).
**********************************************************************

Here we see that we are indeed running a 3D system that is periodic along one dimension.

****************************** Lattice *******************************
  Lattice Vectors [A]
    3.998658    0.000000    0.000000
  Cell length =             3.9987 [A]
  Reciprocal-Lattice Vectors [A^-1]
    1.571323    0.000000    0.000000
  Cell angles [degree]
    alpha =   90.000
    beta  =   90.000
    gamma =   90.000
**********************************************************************

Although we specified three lattice parameters in the input file corresponding to three lattice vectors, here the code tells us that it's only using one lattice vector. This is because the lattice is only periodic along one direction and therefore it is fully determined by the first vector. The other vectors are ignored when generating the Bravais lattice, but they are still used for two other purposes. First, they allow to specify atomic positions in reduced coordinates, which is useful when copying the coordinates from other codes. Second, they specify the lenght of the paralleliped box used to generate the real-space grid. This is confirmed a few lines later:

******************************** Grid ********************************
Simulation Box:
  Type = parallelepiped
  Lengths [A] = (   3.999,  10.580,  10.580)
Main mesh:
  Spacing [A] = ( 0.308, 0.300, 0.300)    volume/point [A^3] =      0.02768
  # inner mesh =      15925
  # total mesh =      34397
  Grid Cutoff [eV] =   397.448444    Grid Cutoff [Ry] =    29.211924
**********************************************************************

Note how the lengths along 'y' and 'z' are the same as the lenghts of the corresponding lattice parameters specified in the input file.

****************************** Hartree *******************************
Input: [DressedOrbitals = no]
The chosen Poisson solver is 'fast Fourier transform'
Input: [PoissonFFTKernel = cylindrical]
**********************************************************************

Info: FFT grid dimensions       = 13 x 75 x 75
      Total grid size           = 73125 (   0.6 MiB )
      Inefficient FFT grid. A better grid would be: 14 75 75
Info: Poisson Cutoff Radius     =  11.2 A

Finally, this piece of output confirms that the code is indeed using a cutoff for the calculation of the Hartree potential, as mentioned in the Introduction. The cutoff used is a cylindrical cutoff and, by comparing the FFT grid dimensions with the size of the simulation box, we see that the Poisson solver is using a supercell doubled in size in the y and z directions.

At this point you might want to play around with the number of k-points until you are sure the calculation is converged.

Band structure

We now modify the input file in the following way:

Band structure for a infinite chain of Sodium atoms, calculated for a single chain (purple lines), and a 3D-periodic crystal of chains in a supercell (green lines).
CalculationMode = unocc
UnitsOutput = ev_angstrom
ExperimentalFeatures = yes

PeriodicDimensions = 1

Spacing = 0.3*angstrom

%LatticeParameters
 3.9986581*angstrom | 10.58*angstrom | 10.58*angstrom
%
 
%Coordinates
 "Na" | 0.0 | 0.0 | 0.0
%

ExtraStates = 6
ExtraStatesToConverge = 4

%KPointsPath
 14
 0.0 | 0.0 | 0.0
 0.5 | 0.0 | 0.0
%
KPointsUseSymmetries = no

and run the code.

You can now plot the band structure using the data from the static/bandstructure file, just like in the Getting started with periodic systems tutorial.

In the introduction we mentioned that performing a calculation with PeriodicDimensions = 1 is not quite the same as performing a PeriodicDimensions= 3 calculation with a large supercell. Lets now check this. Re-run both the ground-state and the unoccupied calculations, but setting PeriodicDimensions = 3 in the above input files. Before doing so, make sure you copy the static/bandstructure file to a different place so that it is not overwritten (better yet, run the new calculations in a different folder). You can see the plot of the two band structures on the right. More comments on this in ref. [1].

h-BN monolayer

Hexagonal boron nitride (h-BN) is an insulator widely studied which has a similar structure to graphene. Here we will describe how to get the band structure of an h-BN monolayer.

Ground-state calculation

We will start by calculating the ground-state using the following input file:

CalculationMode = gs
UnitsOutput = ev_angstrom
ExperimentalFeatures = yes

PeriodicDimensions = 2

Spacing = 0.20*angstrom

BNlength = 1.445*angstrom
a = sqrt(3)*BNlength
L = 40
%LatticeParameters
 a | a | L
%

%LatticeVectors
  1.0 | 0.0       | 0.0
 -1/2 | sqrt(3)/2 | 0.0
  0.0 | 0.0       | 1.0
%
 
%ReducedCoordinates
 'B' | 0.0 | 0.0 | 0.0
 'N' | 1/3 | 2/3 | 0.0
% 

PseudopotentialSet=hgh_lda

LCAOStart=lcao_states 

%KPointsGrid
  12   | 12   | 1
%
KPointsUseSymmetries = yes

Most of the file should be self-explanatory, but here is a more detailed explanation for some of the choices:

  • PeriodicDimensions = 2: A layer of h-BN is periodic in the x-y directions, but not in the z direction, so there are two periodic dimensions.
  • LatticeParameters: Here we have set the bond length to 1.445 Å. The box size in the z direction is 2 L with L large enough to describe a monolayer in the vacuum. Remember that one should always check the convergence of any quantities of interest with the box length value.

If you now run the code, you will notice that the cutoff used for the calculation of the Hartree potential is different than for the Sodium chain, as is to be expected:

****************************** Hartree *******************************
The chosen Poisson solver is 'fast Fourier transform'
Input: [PoissonFFTKernel = planar]
**********************************************************************

Input: [FFTLibrary = fftw]
Info: FFT grid dimensions       = 13 x 13 x 225
      Total grid size           = 38025 (   0.3 MiB )
      Inefficient FFT grid. A better grid would be: 14 14 225
Info: Poisson Cutoff Radius     =  22.5 A

Band Structure

After the ground-state calculation, we will now calculate the band structure. This is the input for the non-self consistent calculation:

CalculationMode = unocc
UnitsOutput = ev_angstrom
ExperimentalFeatures = yes

PeriodicDimensions = 2

Spacing = 0.20*angstrom

BNlength = 1.445*angstrom
a = sqrt(3)*BNlength
L = 40
%LatticeParameters
 a | a | L
%

%LatticeVectors
  1.0 | 0.0       | 0.0
 -1/2 | sqrt(3)/2 | 0.0
  0.0 | 0.0       | 1.0
%
 
%ReducedCoordinates
 'B' | 0.0 | 0.0 | 0.0
 'N' | 1/3 | 2/3 | 0.0
% 

PseudopotentialSet=hgh_lda

LCAOStart=lcao_states 

ExtraStatesToConverge  = 4
ExtraStates  = 8

%KPointsPath
 12  | 7   | 12   # Number of k point to sample each path
 0.0 | 0.0 | 0.0  # Reduced coordinate of the 'Gamma' k point
 1/3 | 1/3 | 0.0  # Reduced coordinate of the 'K' k point
 1/2 | 0.0 | 0.0  # Reduced coordinate of the 'M' k point
 0.0 | 0.0 | 0.0  # Reduced coordinate of the 'Gamma' k point
%
KPointsUseSymmetries = no
Band structure of a monolayer h-BN.

In this case, we chose the following path to calculate the band structure: Gamma-K, K-M, M-Gamma, with a sampling of 12-7-12 k-points.

On the right is the resulting plot of the occupied bands and the first four unoccupied bands from the static/bandstructure file.

References

  1. 1.0 1.1 C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Exact Coulomb cutoff technique for supercell calculations, Phys. Rev. B 73 205119 (2006)

Previous Getting started with periodic systems - Next Optical spectra of solids

Back to Periodic systems