Difference between revisions of "Tutorial:Periodic systems"

From OctopusWiki
Jump to navigation Jump to search
Line 34: Line 34:
 
* You must understand that performing for instance, a {{variable|PeriodicDimensions|Mesh}}<tt> = 1</tt> calculation in Octopus is not quite the same as performing a {{variable|PeriodicDimensions|Mesh}}<tt> = 3</tt> 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.
 
* You must understand that performing for instance, a {{variable|PeriodicDimensions|Mesh}}<tt> = 1</tt> calculation in Octopus is not quite the same as performing a {{variable|PeriodicDimensions|Mesh}}<tt> = 3</tt> 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<ref>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</ref>. This is obtained automatically, since the value of {{variable|PoissonSolver|Hamiltonian}} in a periodic calculation is defaulted to the same value as the {{variable|PeriodicDimensions|Mesh}}. The truncation scheme is described in details in the reference. See also the variable documentation for {{variable|PoissonSolver|Hamiltonian}}, and the section on the [[Tutorial:Periodic systems#Na chain|Na chain]] for an example.
+
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<ref>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</ref>. This is obtained automatically, since the value of {{variable|PoissonSolver|Hamiltonian}} in a periodic calculation is defaulted to the same value as {{variable|PeriodicDimensions|Mesh}}. The truncation scheme is described in details in the reference. See also the variable documentation for {{variable|PoissonSolver|Hamiltonian}}, and the section on the [[Tutorial:Periodic systems#Na chain|Na chain]] for an example.
  
 
Now we can show how some simple periodic calculations are performed.
 
Now we can show how some simple periodic calculations are performed.

Revision as of 09:44, 15 February 2008

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, that 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 KPointsMonkhorstPack, to generate a simple grid, or KPoints to custom the reciprocal space mesh. In the former case you only input the desired number of k points along each axis in the reciprocal space, in the latter 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 effective number of kpoints in the Monkhorst-Pack grid can be slightly adjusted by the code depending on the actual box sizes.
  • 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 defaulted to the same value as PeriodicDimensions. The truncation scheme is described in details in the reference. See also the variable documentation for PoissonSolver, and the section on the Na 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.7

a = 10.2

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 
%

%KPointsMonkhorstPack
  2   | 2   | 2
  1/4 | 1/4 | 1/4
%

ExtraStates = 1
ConvAbsDens = 1e-4

Octopus now outputs some info about the size of the reciprocal unit cell

CRYSTAL: cell v    1.0000[b^3]
CRYSTAL: lattice vectors [b]
a1 =    1.00000    0.00000    0.00000
a2 =    0.00000    1.00000    0.00000
a3 =    0.00000    0.00000    1.00000
CRYSTAL: reciprocal basis [b]
b1 =    6.28319    0.00000    0.00000
b2 =    0.00000    6.28319    0.00000
b3 =    0.00000    0.00000    6.28319

 direct space matrix (some units)

    1.00000   0.00000   0.00000
    0.00000   1.00000   0.00000
    0.00000   0.00000   1.00000

 reciprocal space metrix (some units)

   39.47842   0.00000   0.00000
    0.00000  39.47842   0.00000
    0.00000   0.00000  39.47842

and about the symmetries it finds after having parsed the input file

 The crystal system is cubic        with operations: 
       1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
      25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


 the space group of the crystal is symmorphic

 Operation number    1  5  9 37 41 45



    Rotation matrices (r-lattice) and fractional translations (r-lattice)


    1     1  0  0     0  1  0     0  0  1      0.00000  0.00000  0.00000 E    
    2     0  0  1     1  0  0     0  1  0      0.00000  0.00000  0.00000 C3   
    3     0  1  0     0  0  1     1  0  0      0.00000  0.00000  0.00000 C3   
    4     0  1  0     1  0  0     0  0  1      0.00000  0.00000  0.00000 IC2  
    5     1  0  0     0  0  1     0  1  0      0.00000  0.00000  0.00000 IC2  
    6     0  0  1     0  1  0     1  0  0      0.00000  0.00000  0.00000 IC2  
  skipping symmetry check!

Finally the list of the k points and their wheights sappears both in absolute units, and relative to the primitive vectors

 ikp       weight             kpoint (relative)                kpoint (absolute)
   1     0.125000     0.125000  0.125000  0.125000     0.785398  0.785398  0.785398
   2     0.375000     0.125000  0.125000  0.625000     0.785398  0.785398 -2.356194
   3     0.375000     0.125000  0.625000  0.625000     0.785398 -2.356194  3.926991
   4     0.125000     0.625000  0.625000  0.625000     2.356194  2.356194  2.356194


   4 k points generated by program from parameters :
 ---------------------------------------------------

    n =    2    2    2      s =  0.25  0.25  0.25

The rest of the output is like its non-periodic counterpart, but the list of eigenvalues is now ordered first by the k points

Eigenvalues [H]
Kpoints [b^-1]
#k =   1, k = (    0.077000,    0.077000,    0.077000)
 #st  Spin   Eigenvalue     Occupation
   1   --    -0.200409       2.000000
   2   --    -0.086346       2.000000
   3   --    -0.085696       2.000000
   4   --    -0.059482       2.000000
   5   --     0.095024       2.000000
   6   --     0.124028       2.000000
   7   --     0.126384       2.000000
   8   --     0.196539       2.000000
   9   --     0.230058       2.000000
  10   --     0.232980       2.000000
  11   --     0.246907       2.000000
  12   --     0.289648       2.000000
  13   --     0.316069       2.000000
  14   --     0.316711       2.000000
  15   --     0.399462       2.000000
  16   --     0.410461       2.000000
  17   --     0.460077       0.000000

Finally

Info: SCF converged in   11 iterations

In the static directory now more files are found than in the non-periodic case:

  • bands-gp.dat: the bands in the selected format (default is gp for gnuplot)
  • bands-efermi.dat: the Fermi energy in a format compatible with bands-gp.dat
  • 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

TO DO: to be continued

References

  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