Tutorial:Periodic systems
From OctopusWiki
The extension of a ground-state calculation to a periodic system is quite straightforward in Octopus, but a few subtleties must be considered.
Contents |
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
-
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. -
The second variable is the block that defines the k points. You can use either
KPointsMonkhorstPack, to generate a simple grid, orKPointsto customize 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 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 aPeriodicDimensions= 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 reference [1]. 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= gsPeriodicDimensions= 3Spacing=0.6 a = 10.3BoxShape= parallelepipedLsize= 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 % %KPointsMonkhorstPack2 | 2 | 2 1/4 | 1/4 | 1/4 %ExtraStates= 1
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 weights appear 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
- When the calculation stops in the static directory 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
- 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 kpoints mesh, and you can do so by adding
FromScratch= no, and increasing the number of k points inKPointsMonkhorstPackorKPoints. - The unoccupied orbitals are also calculated in the same way as for a finite system using a suitable calculation mode like
%CalculationMode
gs | 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
Units= ev_angstromDimensions= 3PeriodicDimensions= 1 %CalculationModegs | unocc "gs_" | "unocc_" % %Species"Na" | 23.0 | spec_ps_psf | 11 | 2 | 2 % %Coordinates"Na" | 0.0 | 0.0 | 0.0 | no %BoxShape= parallelepiped %Lsize1.99932905 | 5.29 | 5.29 % %Spacing0.2 | 0.2 | 0.2 % %KPointsMonkhorstPack15 | 1 | 1 %KPointsCenterOfInversion= yesLCAODimension= 9 unocc_EigenSolverMaxIter= 400NumberUnoccStates= 4
Here we have introduced the new variable KPointsCenterOfInversion to exploit the only possible symmetry for a one dimensional system. If you set it to false the k points are produced for the whole Brillouin zone, instead of half of it. The same effect can be obtained in 3D setting KPointsUseSymmetries = no.
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 ********************** Input: [PoissonSolver = fft_cyl] Info: FFT allocated with size ( 18, 105, 105) in slot 1
After Octopus has run once hardly all of the unoccupied states will be converged. Try restarting, and/or increasing the value of EigenSolverMaxIter. You might also 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 there are several ways to control the periodicity of a system in Octopus:
- using the
PeriodicDimensionsvariable, which imposes periodic boundary conditions (and, by default, a Poisson solver) in the chosen number of dimensions; - using explicitly the
PoissonSolvervariable, when you want to force a particular truncation scheme for the long-range part of the Coulomb interaction; - choosing an appropriate grid for the k points via
KPointsMonkhorstPackorKPoints; - tuning the
Lsizes, which ultimately defines the size of your unit cell (or super-cell).
References
- ↑ 1.0 1.1 1.2 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
