lmf in reduced dimensions

Summary

lmf can treat a system that is periodic in a plane and finite in the perpendicular direction, such as a slab, a monolayer or a surface, without periodic images along the third lattice vector. The cell contains the slab and a layer of vacuum, as in the usual supercell approach; in 2D mode the Coulomb interaction does not wrap around the cell along the third lattice vector, so a net dipole of the slab does not produce a spurious periodic-image electric field in the vacuum.

2D mode is switched on with the tag HAM_DIM_TWOD=1, which also selects the truncated Coulomb kernel for the electrostatic potential (token VHAR, which follows TWOD by default). Self-consistency requires some care with the density mixing: see here for the reasoning; for a comparison of 2D and 3D results in a Ni/graphene system, see here.

Table of Contents

• The HAM_DIM tag
• What changes in 2D mode
  • 2D Ewald sums
  • Truncated Coulomb kernel
  • Planar average of the potential in the vacuum
• Density mixing in a slab: special considerations
  • Why screening is important
  • etf is an alternative to elind
  • Why the Thomas-Fermi function is preferred in a slab
• Comparison with 3D
• Dependence on the vacuum thickness
• Convergence of the density mixing in a slab
  • Recommended settings
  • Origin of the floor
  • The –floorxfer switch
  • Dependence on etf and beta
    • 3D
    • 2D
    • Interpretation
  • Convergence to high precision, 2D and 3D
• References

The HAM_DIM tag

2D mode is controlled by the compound tag HAM_DIM in category HAM:

HAM_DIM[twod=1 vhar=1 margin=#]

For example, the following fragment of a ctrl file turns on the 2D mode and keeps the default margin:

HAM     ...
        dim[twod=1]

Geometry requirements:

• the first two lattice vectors PLAT(1) and PLAT(2) must lie in the xy plane and PLAT(3) must be along z.
• the slab must have a thickness no greater than $z_c = L/2$, with $L=\mathrm{ALAT}\times|\mathrm{PLAT}(3)|$, with the charge confined to the slab (the vacuum must be at least $L/4$ on each side). See below. The density tails count: leave at least about 5 a.u. beyond the outermost atoms within $z_c$; see here. The program warns if the room is less than twice the largest augmentation radius.

Restrictions: a homogeneous background charge and the dipole correction are not available when VHAR=1.

TWOD and VHAR act independently and normally you should not set VHAR: its default follows TWOD. The Bloch sums of the compensating densities that correct the plane-wave potential must use the same kernel as the plane-wave sums, so they follow VHAR, not TWOD. The geometry requirements apply when either token is 1.

The mixed settings are not normally needed. In the five-atom slab of this section (vacuum of about 32 a.u.), TWOD=0 and TWOD=1 at the same VHAR gave the same total energies, RMS DQ and forces to at least seven digits over 10 iterations, because the basis functions decay within the vacuum; TWOD=1 was 10–18% faster in that test, since the 2D lattice lists are much shorter. The results depend on VHAR, not on TWOD.

Check that the tag is read with lmf --input (prints the help for HAM_DIM_TWOD, HAM_DIM_VHAR and HAM_DIM_MARGIN) and with --show=2 (prints what was read). The program announces the two switches in the line special: 2D, 2D-ves (2D for TWOD=1 and 2D-ves for VHAR=1).

What changes in 2D mode

Three parts of the program are modified. The Bloch sums of the basis functions are controlled by TWOD; the truncated kernel, the 2D Bloch sums of the compensating densities and the planar-average replacement are controlled by VHAR.

2D Ewald sums

The envelope functions are smoothed Hankel functions. Their Bloch sums, and the Bloch sums of the Gaussian functions, are lattice sums that are evaluated by the Ewald method. In 3D the sums run over the three lattice vectors. In 2D mode the translations run only over the two in-plane lattice vectors, PLAT(1) and PLAT(2). This applies to the envelope functions when TWOD=1, and to the Gaussian and Hankel functions that enter the electrostatic potential when VHAR=1. The Ewald parameter and the real- and reciprocal-space lattice vector lists are generated in 2D mode, from the cell area, in a separate pass of the lattice generator; the 2D lists are not subsets of the 3D ones.

Truncated Coulomb kernel

The electrostatic potential of the smooth density is made in reciprocal space. In 3D, $$v(\mathbf{G}) = \frac{8\pi}{G^2} n(\mathbf{G}) \qquad \mathrm{(Ry\ units)}.$$

In 2D mode the Coulomb interaction is truncated at the length $z_c = L/2$ along z (the slab-truncated kernel of Ismail-Beigi and of Rozzi et al., Refs. 1, 2): $$v(\mathbf{G}) = \frac{8\pi}{G^2}\Bigl[1-e^{-G_\parallel z_c}\cos(G_z z_c)\Bigr] n(\mathbf{G}).$$

The $\mathbf{G}=0$ term is dropped as in 3D. This is the exact potential of an isolated slab for charge confined to a region of thickness $\le z_c$, so the slab and the vacuum must fit the cell accordingly (vacuum of at least $L/4$ on each side). The in-plane components $G_\parallel\ne0$ decay as $e^{-G_\parallel z_c}$ away from the slab, as for the true potential. The same truncation is applied to the Gaussian (compensating) charges.

The kernel is exact inside $z_c$ of all the charge. Farther away the planar average ($G_\parallel=0$) of the truncated potential is the sign-flipped image of the slab region, that is, it is not the correct vacuum potential. The next section fixes this.

Planar average of the potential in the vacuum

After the smooth electrostatic potential has been transformed back to the real-space mesh, the routine smvpln2d (in fp/smves.f) replaces the planar average of the potential ($G_\parallel=0$ component) outside the slab by its vacuum continuation. The components $G_\parallel\ne0$ are left alone.

1. The slab is bounded by the outermost atoms: the largest periodic gap in z between atoms is taken to be vacuum. The slab is extended on each side by the margin (HAM_DIM_MARGIN).
2. Inside this region the planar average is unchanged.
3. In the vacuum the planar average is replaced by the constant value it has at the edge on each side of the slab, and the two values are joined by a quintic step over the central half of the vacuum, so that the potential remains periodic and smooth.

The default margin is $\min(2r_\mathrm{MT}^\mathrm{max},(z_c-t)/2)$. The program prints the atom extent, the margin, $z_c$ and the step of the potential across the slab, for example

2D planar ves: atoms span z = 0.000 to 13.350,  margin = 4.562 (zc = 22.738),  step across slab (top-bottom) = -0.096190

The margin keeps the replacement away from the region outside the outermost atoms where the density is still appreciable. The program prints a warning if the slab plus margin do not fit within $z_c$. With a negative value of HAM_DIM_MARGIN the replacement is skipped. In tests, switching the replacement on or off changed the force correction by at most 0.4 mRy/au.

Density mixing in a slab: special considerations

The density of a slab is mixed with the same algorithm as any other lmf calculation, controlled by the string ITER_MIX, but the choice of the model dielectric function that screens the density difference matters more. For slabs we recommend the option etf (Thomas-Fermi function) rather than elind (Lindhard function), for example

mix=a2,etf=0.83,b=0.2,n=10,k=10

together with the switch –floorxfer. Tests and tables are below; this section explains the reasoning.

Why screening is important

Let $\Delta=n^\mathrm{out}-n^\mathrm{in}$ be the difference between the output and input densities of an iteration. The change in the electrostatic potential caused by $\Delta$ is $\delta V(\mathbf G)=(8\pi/G^2)\,\Delta(\mathbf G)$ (Ry units), so the small-$G$ components of $\Delta$ have the largest effect. If the new input density is simply $n^\mathrm{in}+\beta\Delta$, these long-wavelength components make the potential, and then the output density, overshoot from one iteration to the next: charge sloshing. The cell of a slab is long in the direction perpendicular to the slab, so its smallest $G_z$ is small, and charge can slosh across it. The mixing step therefore replaces $\Delta(\mathbf G)$ by $\varepsilon^{-1}(G)\,\Delta(\mathbf G)$, with a model dielectric function $\varepsilon$ that is large at small $G$ (the “Kerker” idea; see here). In the language of the code, a screening charge $S=(\varepsilon^{-1}-1)\Delta$ is added to the output density.

Without screening the iteration for the test slab is unstable even for $\beta=0.1$: with simple linear mixing the RMS DQ grew from 0.031 to 0.16 in six iterations. With screening it converges. At the smallest $G_z$ of the cell of this test, $\varepsilon\approx50$, so these modes are damped by that factor.

etf is an alternative to elind

Both options of ITER_MIX define the model dielectric function from a single parameter, a Fermi energy $k_F^2$, with the same default: the value of HAM_ELIND, which is $-0.7$ unless you set it, meaning 0.7 times the free-electron gas value for the number of valence electrons (see here). For the test slab this default is 0.83 Ry, which is why the tests with etf=0.83 use the same $k_F^2$ as the default Lindhard function. They differ only in the function: $$\varepsilon_\mathrm{L}(G)=1+\frac{k_\mathrm{TF}^2}{G^2}\,F\!\left(\frac{G}{2k_F}\right),\qquad \varepsilon_\mathrm{TF}(G)=1+\frac{k_\mathrm{TF}^2}{G^2},\qquad k_\mathrm{TF}^2=\frac{4k_F}{\pi},$$

where $F(x)=\tfrac12+\frac{1-x^2}{4x}\ln\left|\frac{1+x}{1-x}\right|$ is the Lindhard function (elind; exact for free electrons). The Thomas-Fermi function (etf) is the case $F=1$. Since $F(0)=1$, the two functions are the same at small $G$, which is the region that controls charge sloshing. They differ at $G\sim2k_F$, where $F$ has a logarithmic singularity in its derivative.

Why the Thomas-Fermi function is preferred in a slab

The screening charge $S(\mathbf G)=(\varepsilon^{-1}-1)\Delta(\mathbf G)$ is nonlocal: in real space it is the density difference convolved with the Fourier transform of $\varepsilon^{-1}-1$. Inside the slab $\Delta$ is large (up to 3×10⁻² e/a.u.³ in the first iteration), while in the vacuum it is below 8×10⁻⁵. What the screening does with the large differences in the slab therefore decides what appears in the vacuum.

• For the Thomas-Fermi function $\varepsilon^{-1}-1=-k_\mathrm{TF}^2/(G^2+k_\mathrm{TF}^2)$ is smooth, and the corresponding real-space kernel decays exponentially with the screening length $1/k_\mathrm{TF}$. The screening charge generated in the slab does not reach far into the vacuum.
• The Lindhard function is non-analytic at $G=2k_F$. The screening charge found in the tests oscillates in sign and does not decay exponentially with the distance from the slab, which is consistent with that. It deposits a small spurious density in the vacuum in every iteration.

In the first iteration of the test slab the mean screening charge in the vacuum interior is 2.2×10⁻⁶ e/a.u.³ for the Lindhard function and 8.8×10⁻⁸ for the Thomas-Fermi function with the same parameter, 25 times less. Over six iterations of simple mixing the vacuum density of the input density grows to 1.8×10⁻⁶ with the Lindhard function and to 1.5×10⁻⁷ with the Thomas-Fermi function. This growing background (the “floor”) is hardly visible in the RMS DQ, because that quantity is dominated by the large-$G$ components.

Thus the Thomas-Fermi function keeps what is wanted from the screening, the damping of the small-$G$ components, and removes most of an artifact of the Lindhard function. It does not remove the floor completely; the switch –floorxfer does the rest. The Lindhard function remains the default and is the exact free-electron response; the problem described here is specific to slabs and other cells with a large vacuum.

The value of etf acts as a strength of the damping. The tests used 0.4, 0.83 and 1.6; 0.83 is a good compromise. A large value stabilizes the iteration at the larger $\beta=0.2$ but slows the convergence of the energy: with etf=1.6 the difference ehk−ehf was 7×10⁻⁴–9×10⁻⁴ Ry after 50 iterations in 2D, against at most 5.3×10⁻⁵ Ry for the smaller values. See here.

Comparison with 3D

In 3D the periodic images of a slab with a net dipole produce a spurious electric field in the vacuum. With the truncated kernel (VHAR=1) there is no such field, so the total energy and forces differ slightly from the 3D values for such a slab, even for a thick vacuum. The following compares the two modes for a five-atom slab (two species, not mirror symmetric, with a net dipole) and a cell with $L=45.5$ a.u., in which the atoms span 13.35 a.u. along z. Both runs start from free atoms, use Anderson mixing with Thomas-Fermi screening (mix=a2,etf=0.83,b=0.2) and --floorxfer (see below); they are converged to RMS DQ of about 5×10⁻⁶ (3D, 70 iterations) and 2×10⁻⁶ (2D, 90 iterations).

Total force along z (mRy/au) including the eigenvalue correction:

• The forces agree to 0.09 mRy/au or better. The electrostatic and eigenvalue parts of the forces separately differ by up to 0.6 mRy/au; these compensate in the total.
• The total energies differ by about 0.2 mRy (ehf) and 0.16 mRy (ehk). Part of this, about 0.115 mRy, is due to the finite thickness of the vacuum in the 2D calculation; see the dependence on the vacuum thickness below. The rest is expected from the different treatments of electrostatics.
• The energy from the 2D run is still drifting slowly at iteration 90 (by about −3×10⁻⁷ Ry per iteration), so the 2D energy is not completely converged. Extrapolating this decay suggests that ehf will settle at about 1.2–2.0×10⁻⁵ Ry lower; this extrapolation is rough.
• ¹ The Fermi levels are not directly comparable, because the zero of the potential is not the same in the two modes.

Dependence on the vacuum thickness

The same calculation was done for three cells that differ only in the length of the third lattice vector: 0.8, 1.0 and 1.2 times that of the cell above. The atoms (spanning 13.35 a.u. along z) and the in-plane lattice are the same. All calculations use mix=a2,etf=0.83,b=0.2 and –floorxfer, start from free atoms and were iterated 90 times (70 for the 3D run of the middle cell, which is the one used in the comparison above). Here $t$ is the span of the atoms.

Total force along z (mRy/au):

• 3D is independent of the length of the cell in this range: the energies of the 36.4 and 45.5 a.u. cells differ by 0.03 mRy and the forces by at most 0.03 mRy/au. The meshes of these cells are different (28×28×210 and 30×30×264 divisions), so a mesh effect is smaller than that.
• 2D is not. In the shortest cell the 2D calculation is wrong: the energy is 31 mRy too high, the forces differ from 3D by up to 17 mRy/au, the Harris force is 46.7 mRy/au, and the iteration does not converge as it does in the other cells (the RMS DQ is 10⁻⁵ and the energy is still drifting by 4×10⁻⁶ Ry per iteration at iteration 90).
• Cause. The truncated Coulomb kernel is exact only for charge that is confined to a region of thickness $z_c$. In the shortest cell only 2.4 a.u. are available beyond the outermost atoms, about the radius of the augmentation sphere (2.28 a.u.), so part of the density tail lies outside that region. In the middle cell there are 4.7 a.u. and in the longest 7.0 a.u. This is consistent with the results, but we did not test it independently.
• Remaining dependence. The 2D energies of the two longer cells still differ by 0.115 mRy and the forces by up to 0.10 mRy/au; as above, the 3D results indicate that this is not a mesh effect. It is probably the same cause, so the 0.196 mRy difference between 2D and 3D for the middle cell may contain about 0.115 mRy that is due to the finite thickness of the vacuum. For the longest cell the difference from the 3D energy of the middle cell is 0.081 mRy (and the ehf of both 2D runs is still drifting by about −3×10⁻⁷ Ry per iteration).

Rule of thumb. Leave at least about 5 a.u. of vacuum beyond the outermost atoms, inside the truncation length, that is $(z_c-t)/2\gtrsim5$ a.u., with $z_c=L/2$ and $t$ the span of the atoms. In this system 2.4 a.u. is not enough, and with 4.7 a.u. the energy is within 0.12 mRy and the forces within 0.10 mRy/au of those for 7.0 a.u. We did not test longer cells, so the energy at 7.0 a.u. may itself not be fully converged with respect to the length of the cell. The program prints a warning if $(z_c-t)/2$ is smaller than twice the largest augmentation radius (4.56 a.u. here). Note that this is a condition on the density and not only on the atoms.

Convergence of the density mixing in a slab

In a calculation for a slab (a cell with a large region of vacuum, in 3D or 2D mode), the density mixing algorithm can build up a small spurious density in the vacuum, a “floor”, of order 10⁻⁵ e/a.u.³. The floor is not present in the output density of the band calculation. It is made by the mixing step, and it persists because the screened RMS difference between input and output densities (printed as RMS DQ) is almost insensitive to it, so the iteration can look converged while the vacuum density is still wrong.

This section explains where the floor comes from, how to avoid it (Anderson mixing with Thomas-Fermi screening, option etf of ITER_MIX, and the switch –floorxfer) and shows how convergence depends on the screening parameter etf and the mixing parameter beta, in 3D and in 2D. The reasoning behind etf is given above, and the mathematics of –floorxfer in the self-consistency page.

Recommended settings

• Anderson mixing with Thomas-Fermi screening, with etf of about 0.83 and beta between 0.1 and 0.2, for example in category ITER:

  mix=a2,etf=0.83,b=0.2,n=10,k=10
  

• the switch –floorxfer on the command line (see below).

Do not rely on RMS DQ alone: also check that ehk and ehf agree. Do not judge forces before the iteration has converged; in the first iterations the Harris force can be 40–80 times its converged value (see below).

Origin of the floor

The mixing step adds to the linearly mixed density the screening charge $S(\mathbf G)=(\varepsilon^{-1}(G)-1)\Delta(\mathbf G)$ (see above). It is made from the large differences $\Delta$ inside the slab and extends into the vacuum.

Tests starting from free atoms, with simple linear mixing (a0, beta=0.1), for the five-atom slab described below, in 3D:

The table shows the mean density of the input density in the vacuum (e/a.u.³); the output density of the band calculation has about 10⁻⁹ in the vacuum. With the default mixing (Anderson, Lindhard function) from free atoms the floor was 1.2×10⁻⁵, still rising, after 15 iterations.

• No screening: no floor, but not stable. Without screening (elind=0) the vacuum density stays at about 3–8×10⁻⁹, but the iteration is unstable even at beta=0.1: the RMS DQ is 0.031, 0.032, 0.038, 0.068, 0.12, 0.16 in the six iterations.
• The floor is made by the screening of the large density differences in the slab. About 99% of the screening charge in the vacuum interior at iteration 1 comes from $\Delta$ inside the slab; the difference in the vacuum contributes 100 times less. The tail of the Lindhard screening charge is 25 times larger than that of the Thomas-Fermi function (see above).
• Why it persists. At the lowest $G_z$ of the cell, $\varepsilon\approx53$, so a Fourier component of the floor is removed by the mixing at a rate of only $\beta/\varepsilon\approx0.4\%$ per iteration, a time constant of the order of 260 iterations. The RMS DQ is dominated by the larger-$G$ components, and does not see it.
• Thomas-Fermi helps but does not cure. With etf=0.83, which reproduces the Lindhard function at the relevant small $G$, the floor is generated about 12 times more slowly (table above), but an existing floor is not removed. A uniform shift cannot be used to remove it: the $\mathbf G=0$ term of the screening charge is zero, and the neutrality correction of the mixing step removes any uniform shift within one iteration.

The –floorxfer switch

The switch –floorxfer[~theta=#][~alpha=#] of lmf (see command line switches) moves the charge that the mixing step leaves in the vacuum into the dense region, conserving the total charge. Its motivation, the derivation and its properties are given here. The defaults are $\theta=10^{-4}$ and $\alpha=1$.

In the tests below the moved charge is 1.4×10⁻² to 1.9×10⁻² e in the first iteration and at most about 2×10⁻⁵ e (3D) and 7×10⁻⁵ e (2D) after 50 iterations.

Dependence on etf and beta

Every test was started from free atoms (no restart file, no mixing file) with mix=a2, –floorxfer and the parameters shown, for 50 iterations. The cell is a five-atom slab; the 3D and 2D cells are identical, so 2D differs only in the treatment of the Coulomb interaction (2D mode). In the tables DQ is the RMS DQ and ΔE = ehk − ehf (Ry).

3D

2D

Interpretation

• Smooth versus fast. beta=0.1 is smooth and monotonic for every etf but slow: the RMS DQ is 2–4×10⁻⁴ after 50 iterations in both 3D and 2D. beta=0.2 reaches a lower RMS DQ (1.3–2.3×10⁻⁵ in 3D; 2.5×10⁻⁵ to 2.3×10⁻⁴ in 2D) but not monotonically: after iteration 10 the RMS DQ rises by a factor of 2 or more in 5, 4 and 2 iterations in 3D (etf=0.4, 0.83, 1.6) and 5, 6 and 3 iterations in 2D.
• Larger etf is more stable at the higher beta but slower in the energy. A large etf damps the small-$G$ components more strongly, which suppresses the jumps at beta=0.2, but with etf=1.6 the ΔE at iteration 50 is 7×10⁻⁴–9×10⁻⁴ Ry in 2D, against at most 5.3×10⁻⁵ Ry for the smaller values. The value 0.83, which reproduces the Lindhard function at small $G$, is a good compromise.
• Energies converge faster than the RMS DQ suggests. With the floor removed, ΔE is about 10⁻⁵ Ry or less at RMS DQ of 1–5×10⁻⁴; the RMS DQ of 2×10⁻⁷ that is the default convc in these tests was not reached in any of the runs within 50 iterations.
• Forces are not meaningful in the first iterations. The maximum Harris force has a spike of up to about 1200 (3D) and 2200 (2D) mRy/au in iterations 4–9 and settles at 29–29.4 mRy/au.
• 2D versus 3D. The convergence is similar in 2D and 3D. The RMS DQ at iteration 50 is usually somewhat larger in 2D (1.2 to 2.5 times), with the exception of etf=1.6, beta=0.2, where it is 18 times larger after a late rise.

Convergence to high precision, 2D and 3D

The case etf=0.83, beta=0.2 was continued, with the mixing file removed every 10 iterations (kill=10), to 70 iterations in 3D and 90 iterations in 2D. The comparison of the final results is in Comparison with 3D.

The energy of the 2D run drifts steadily and slowly, while the RMS DQ stays at a few times 10⁻⁶. The change in ehf per iteration falls from −4.9×10⁻⁷ (iterations 51–60) to −3.4×10⁻⁷ (iterations 81–90). This is again a slow, long-wavelength mode that the RMS DQ does not see. It is not understood.

References

[1] S. Ismail-Beigi, Truncation of periodic image interactions for confined systems, Phys. Rev. B73, 233103 (2006).

[2] C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, Exact Coulomb cutoff technique for supercell calculations, Phys. Rev. B73, 205119 (2006).