Tutorial:Wires and slabs

From OctopusWiki
Revision as of 09:35, 7 August 2018 by Micael (talk | contribs)
Jump to navigation Jump to search

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

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

Our input is

UnitsOutput = ev_angstrom
ExperimentalFeatures = yes

Dimensions = 3
PeriodicDimensions = 1

FromScratch = yes

%Species
   'Na' | species_pseudo | db_file | 'PSF/Na.psf' | lmax | 2 | lloc | 2
%

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

BoxShape = parallelepiped

%Lsize
 1.99932905*angstrom | 5.29*angstrom | 5.29*angstrom
%

Spacing = 0.2*angstrom

%KPointsGrid
 15 | 1 | 1
%
KPointsUseSymmetries = yes

Later re-run for unoccupied states by adding

CalculationMode = unocc
EigensolverMaxIter = 400
ExtraStates = 4

Look at the comments on the LCAO in the output file. What's going on? Why can't a full initialization with LCAO be done?

Bands for a infinite chain of Sodium atoms, calculated for a single chain (blue line), and a 3D-periodic crystal of chains in a supercell (green line).

The following piece of output confirms that the code is actually using a cylindrical cutoff, which, in turn, requires the supercell to be doubled in size in the y and z directions:

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

Input: [FFTLibrary = fftw]
Info: FFT grid dimensions       = 20 x 105 x 105
      Total grid size           = 220500 (   1.7 MiB )
Info: Poisson Cutoff Radius     =  10.5 A

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

Just to stress the differences between a PeriodicDimensions = 1 and a PeriodicDimensions = 3 calculation let us change in the input above just the line

PeriodicDimensions = 3

and we can now compare the bands obtained in the two cases. 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

A layer of h-BN is periodic in the x-y directions, but not in the z. Thus we will set PeriodicDimensions = 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:

CalculationMode = gs

FromScratch = yes

ExperimentalFeatures = yes

PeriodicDimensions = 2

Spacing = 0.20*angstrom

BoxShape = parallelepiped

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

%LatticeParameters
 a | a | L
%

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

PseudopotentialSet=hgh_lda

LCAOStart=lcao_states 

%KPointsGrid
  12   | 12   | 1
%

ExtraStates = 2

UnitsOutput = ev_angstrom

Band Structure

Convergence of the band gap with respect to the spacing for an h-BN monolayer.

After this GS calculation, we will perform an unocc run. This non-self consistent calculation which needs the density from the previous GS calculation.

CalculationMode = unocc

ExtraStatesToConverge  = 5
ExtraStates  = 10

Here the number of ExtraStates is the number of unoccupied bands in the final band structure. The highest occupied states are extremely complicated to converge. Therefore the variable ExtraStatesToConverge 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 KPointsPath . Instead of using the KPointsGrid block of the GS calculation, we use during this unocc calculation:

%KPointsPath
 12  |  7  | 12                 # Number of k point to sample each path
 0   |  0  | 0                  # Reduced coordinate of the 'Gamma' k point
 1/3 | 1/3 | 0                  # Reduced coordinate of the 'K' k point
 1/2 |  0  | 0                  # Reduced coordinate of the 'M' k point
 0   |  0  | 0                  # Reduced coordinate of the 'Gamma' k point
%


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.

Info: The code will run in band structure mode.
     No restart information will be printed.

Moreover, by using KPointsPath, 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.

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)


Back to Tutorials