Tutorial:Periodic systems

From OctopusWiki
Jump to: navigation, search

The extension of a ground-state calculation to a periodic system is quite straightforward in Octopus, but a few subtleties must be considered.

Needed input variables

You already know that you can start a ground state calculation in Octopus with a tiny input file basically including only the coordinates of the atoms in the system. In order to make the system periodic you basically need only two extra variables:

  1. The first variable is PeriodicDimensions which must be set equal to the number of dimensions you want to consider as periodic, namely:
    • 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.
  2. The second variable is the block that defines the k points. You can use either KPointsGrid, to generate a simple Monkhorst-Pack grid, or KPoints or KPointsReduced to customize the reciprocal-space mesh. In the former case you only input the desired number and shift of k-points along each axis in the Brillouin zone; in the latter cases you can explicitly set the position and weight of each k-point.

Warnings and comments

  • At the moment Octopus only accepts BoxShape = parallelepiped for periodic systems.
  • The number of k-points specified in input is meant to be the minimum number of points in each direction in the reciprocal space for the full Brillouin zone. The effective number of points is of course adjusted by the code when symmetries are used to shrink the Brillouin zone to its irreducible portion.
  • You must 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.

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. But 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 obtained automatically, since the value of PoissonSolver in a periodic calculation is chosen according to PeriodicDimensions. See also the variable documentation for PoissonSolver, and the section on the Sodium chain for an example.

Now we can show how some simple periodic calculations are performed.

Examples

Bulk Silicon

Let us follow this simple test. The input is

CalculationMode = gs

PeriodicDimensions = 3

Spacing = 0.6

a = 10.3

BoxShape = parallelepiped
Lsize = a/2

%Coordinates
  "Si" |   0.0       | 0.0       | 0.0 
  "Si" |   a/2       | a/2       | 0.0
  "Si" |   a/2       | 0.0       | a/2
  "Si" |   0.0       | a/2       | a/2 
  "Si" |   a/4       | a/4       | a/4 
  "Si" |   a/4 + a/2 | a/4 + a/2 | a/4
  "Si" |   a/4 + a/2 | a/4       | a/4 + a/2
  "Si" |   a/4       | a/4 + a/2 | a/4 + a/2 
%

%KPointsGrid
  2   | 2   | 2
  1/2 | 1/2 | 1/2
%
KPointsUseSymmetries = yes

ExtraStates = 1
Output = dos
OutputFormat = xyz

Octopus now outputs some info about the cell in real and reciprocal space:

******************************** Grid ********************************
Simulation Box:
  Type = parallelepiped
  Lengths [b] = (   5.150,   5.150,   5.150)
  Octopus will run in 3 dimension(s).
  Octopus will treat the system as periodic in 3 dimension(s).

  Lattice Vectors [b]
   10.300000    0.000000    0.000000
    0.000000   10.300000    0.000000
    0.000000    0.000000   10.300000
  Cell volume =          1092.7270 [b^3]
  Reciprocal-Lattice Vectors [b^-1]
    0.610018    0.000000    0.000000
    0.000000    0.610018    0.000000
    0.000000    0.000000    0.610018
Main mesh:
  Spacing [b] = ( 0.606, 0.606, 0.606)    volume/point [b^3] =      0.22242
  # inner mesh =       4913
  # total mesh =      12729
  Grid Cutoff [H] =    13.442905    Grid Cutoff [Ry] =    26.885811
**********************************************************************

and about the symmetries it finds for the specified structure:

***************************** Symmetries *****************************
Space group No.  227
International: Fd-3m
Schoenflies: Oh^7
Identity has a fractional translation     0.000000    0.500000    0.500000
Identity has a fractional translation     0.500000    0.000000    0.500000
Identity has a fractional translation     0.500000    0.500000    0.000000
Disabling fractional translations. System appears to be a supercell.
  Index                Rotation matrix                      Fractional translations
    1 :     1   0   0     0   1   0     0   0   1      0.000000    0.000000    0.000000
    2 :    -1   0   0     0   1   0     0   0  -1      0.000000    0.000000    0.000000
    3 :    -1   0   0     0  -1   0     0   0   1      0.000000    0.000000    0.000000
    4 :     1   0   0     0  -1   0     0   0  -1      0.000000    0.000000    0.000000
    5 :     0   0  -1    -1   0   0     0   1   0      0.000000    0.000000    0.000000
    6 :     0   0   1     1   0   0     0   1   0      0.000000    0.000000    0.000000
    7 :     0   0  -1     1   0   0     0  -1   0      0.000000    0.000000    0.000000
    8 :     0   0   1    -1   0   0     0  -1   0      0.000000    0.000000    0.000000
    9 :     0  -1   0     0   0   1    -1   0   0      0.000000    0.000000    0.000000
   10 :     0  -1   0     0   0  -1     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 :     0   1   0     0   0  -1    -1   0   0      0.000000    0.000000    0.000000
   13 :     0   0   1     0  -1   0    -1   0   0      0.000000    0.000000    0.000000
   14 :     0   0  -1     0  -1   0     1   0   0      0.000000    0.000000    0.000000
   15 :     0   0  -1     0   1   0    -1   0   0      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   0     0   0  -1     0  -1   0      0.000000    0.000000    0.000000
   18 :    -1   0   0     0   0   1     0  -1   0      0.000000    0.000000    0.000000
   19 :     1   0   0     0   0   1     0   1   0      0.000000    0.000000    0.000000
   20 :    -1   0   0     0   0  -1     0   1   0      0.000000    0.000000    0.000000
   21 :     0   1   0    -1   0   0     0   0  -1      0.000000    0.000000    0.000000
   22 :     0   1   0     1   0   0     0   0   1      0.000000    0.000000    0.000000
   23 :     0  -1   0    -1   0   0     0   0   1      0.000000    0.000000    0.000000
   24 :     0  -1   0     1   0   0     0   0  -1      0.000000    0.000000    0.000000
Info: The system has    24 symmetries that can be used.
**********************************************************************

We are indeed using a cubic supercell with 8 atoms rather than the fcc primitive cell with 2 atoms, which is why we see messages about fractional translations.

Finally the list of the k-points (reduced by symmetry) and their weights appear in reduced coordinates. Use of symmetries was turned on by KPointsUseSymmetries = yes, saving some time.

   1 k-points generated from parameters :
 ---------------------------------------------------
    n =    2    2    2      s =  0.50  0.50  0.50

 index |    weight    |             coordinates              |
     1 |     1.000000 |    0.250000    0.250000    0.250000  |

The rest of the output is like its non-periodic counterpart:

*********************** SCF CYCLE ITER #    9 ************************
 etot = -3.16853104E+01 abs_ev   =  6.44E-05 rel_ev   =  5.45E-05
                        abs_dens =  2.32E-04 rel_dens =  7.26E-06
Matrix vector products:    187
Converged eigenvectors:     17

#  State  KPoint  Eigenvalue [H]  Occupation    Error
      1       1       -0.259203    2.000000   (9.9E-07)
      2       1       -0.188770    2.000000   (9.8E-07)
      3       1       -0.188672    2.000000   (8.7E-07)
      4       1       -0.188672    2.000000   (9.5E-07)
      5       1       -0.086775    2.000000   (9.9E-07)
      6       1       -0.086741    2.000000   (8.4E-07)
      7       1       -0.086741    2.000000   (8.4E-07)
      8       1        0.004936    2.000000   (8.3E-07)
      9       1        0.018504    2.000000   (8.8E-07)
     10       1        0.018537    2.000000   (9.0E-07)
     11       1        0.018537    2.000000   (8.9E-07)
     12       1        0.065601    2.000000   (9.1E-07)
     13       1        0.065601    2.000000   (8.7E-07)
     14       1        0.065604    2.000000   (8.6E-07)
     15       1        0.118898    2.000000   (8.9E-07)
     16       1        0.118898    2.000000   (8.8E-07)
     17       1        0.199921    0.000000   (8.5E-07)

Density of states:

----------------------------------------------------------------------
----------------------------------------------------------------------
----------------------------------------------------------------------
----------------------------------------------------------------------
----------------------------------------------------------------------
----------------------------------------------------------------------
----------------------------------------------------------------------
-----------%--------------%---------------%------%--------------------
-----------%--------------%---------------%------%-------%------------
%----------%--------------%-------------%-%------%-------%-----------%
                                                         ^


Elapsed time for SCF step     9:          1.30
**********************************************************************

  • When the calculation stops in the static directory more files are found than in the non-periodic case:
    • 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.
  • You might also want to refine the k-point mesh, and you can do so by adding FromScratch = no (to use the previously calculated density, though not necessarily the wavefunctions, as a starting point), and increasing the number of k-points in KPointsGrid. When there is only one k-point (as in the previous run), the code does not bother writing bandstructure files. But now with more k-points, you will get these files in static:
    • bands-gp.dat: the bands in the selected format (default is gp for gnuplot). Control format via the variables OutputBandsGnuplotMode and OutputBandsGraceMode.
    • bands-efermi.dat: the Fermi energy in a format compatible with bands-gp.dat
  • The unoccupied orbitals are also calculated in the same way as for a finite system, via CalculationMode = unocc.

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

Summary

In summary, bear in mind that the periodicity of a system in Octopus is controlled by various parameters working together:

  1. using the PeriodicDimensions variable, which imposes periodic boundary conditions (and, by default, a Poisson solver) in the chosen number of dimensions;
  2. using explicitly the PoissonSolver variable, when you want to force a particular truncation scheme for the long-range part of the Coulomb interaction -- in general, you should leave this alone and use the default though;
  3. choosing an appropriate grid for the k-points via KPointsGrid, KPoints, or KPointsReduced;
  4. tuning the Lsize, which ultimately defines the size of your unit cell (or super-cell).

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); http://link.aps.org/abstract/PRB/v73/e205119

Previous Tutorial:Large systems: the Fullerene molecule - Next Tutorial:Geometry optimization

Back to Tutorial