Computational parameters that are important for bulk structures#
k-point convergence#
In the section on molecules, we learned that the total energy is a function of the planewave cutoff energy (ENCUT) used. In bulk systems that is true also. There is also another calculation parameter you must consider, the k-point grid. The k-point grid is a computational tool used to approximate integrals of some property, e.g. the electron density, over the entire unit cell. The integration is performed in reciprocal space (i.e. in the Brillouin zone) for convenience and efficiency, and the k-point grid is where the property is sampled for the integration. The higher the number of sampled points, the more accurately the integrals are approximated.
We will typically use a Monkhorst-Pack cite:PhysRevB.13.5188 \(k\)-point grid, which is essentially a uniformly spaced grid in the Brillouin zone. Another less commonly used scheme is the Chadi-Cohen k-point grid cite:PhysRevB.8.5747. The Monkhorst-Pack grids are specified as \(n1 \times n2 \times n3\) grids, and the total number of k-points is \(n1 \cdot n2 \cdot n3\). The computational cost is linear in the total number of k-points, so a calculation on a \(4 \times 4 \times 4\) grid will be roughly 8 times more expensive than on a \(2 \times 2 \times 2\) grid. Hence, one seeks again to balance convergence with computational tractability. Below we consider the k-point convergence of fcc Ag.
from ase.lattice.cubic import FaceCenteredCubic
from vasp import Vasp
import numpy as np
atoms = FaceCenteredCubic('Ag')
KPTS = [2, 3, 4, 5, 6, 8, 10]
TE = []
for k in KPTS:
calc = Vasp(label=f'bulk/Ag-kpts-{k}',
xc='PBE',
kpts=[k, k, k], # specifies the Monkhorst-Pack grid
encut=300,
atoms=atoms)
TE.append(atoms.get_potential_energy())
if None in TE:
calc.abort()
import matplotlib.pyplot as plt
# consider the change in energy from lowest energy state
TE = np.array(TE)
TE -= TE.min()
plt.plot(KPTS, TE)
plt.xlabel('number of k-points in each dimension')
plt.ylabel('Total Energy (eV)');
Based on this figure, we need at least a \(6 \times 6 \times 6\) k-point grid to achieve a convergence level of at least 50 meV. Note: the k-point convergence is not always monotonic like it is in this example, and sometimes very dense grids (e.g. up to \(20 \times 20 \times 20\)) are needed for highly converged properties such as the density of states in smaller unit cells. Oscillations in the total energy are typical, and it can be difficult to get high levels of convergence. The best practices are to use the same k-point sampling grid in energy differences where possible, and dense (high numbers of k-points) otherwise. It is important to check for convergence in these cases.
As unit cells get larger, the number of k-points required becomes smaller. For example, if a \(1 \times 1 \times 1\) fcc unit cell shows converged energies in a \(12 \times 12 \times 12\) k-point grid, then a \(2 \times 2 \times 2\) fcc unit cell would show the same level of convergence with a \(6 \times 6 \times 6\) k-point grid. In other words, doubling the unit cell vectors results in a halving of the number of k-points.
Sometimes you may see k-points described as k-points per reciprocal atom. For example, a \(12 \times 12 \times 12\) k-point grid for a primitive fcc unit cell would be 1728 k-points per reciprocal atom. A \(2 \times 2 \times 2\) fcc unit cell has eight atoms in it, or 0.125 reciprocal atoms, so a \(6 \times 6 \times 6\) k-point grid has 216 k-points in it, or 216/0.125 = 1728 k-points per reciprocal atom, the same as we discussed before.
In the k-point convergence example above, we used a \(6 \times 6 \times 6\) k-point grid on a unit cell with four atoms in it, leading to 864 k-points per reciprocal atom. If we had instead used the primitive unit cell, we would need either a \(9 \times 9 \times 9\) or \(10 \times 10 \times 10\) k-point grid to get a similar level of accuracy. In this case, there is no exact matching of k-point grids due to the difference in shape of the cells.
Effect of SIGMA#
In the self-consistent cycle of a DFT calculation, the total energy is minimized with respect to occupation of the Kohn-Sham orbitals. At absolute zero, a band is either occupied or empty. This discrete occupation results in discontinuous changes in energy with changes in occupation, which makes it difficult to converge. One solution is to artificially broaden the band occupancies, as if they were occupied at a higher temperature where partial occupation is possible. This results in a continuous dependence of energy on the partial occupancy, and dramatically increases the rate of convergence. [BROKEN LINK: incar:SIGMA] and [BROKEN LINK: incar:ISMEAR] affect how the partial occupancies of the bands are determined.
Some rules to keep in mind:
The smearing methods were designed for metals. For molecules, semiconductors and insulators you should use a very small SIGMA (e.g. 0.01).
Standard values for metallic systems is SIGMA=0.1, but the best SIGMA may be material specific.
The consequence of this finite temperature is that additional bands must be included in the calculation to allow for the partially occupied states above the Fermi level; the number of extra bands depends on the temperature used. An example of the maximum occupancies of the bands for an Cu bulk as a function of SIGMA is shown in Figure ref:fig:sigma-occ. Obviously, as SIGMA approaches 0, the occupancy approaches a step function. It is preferable that the occupancy of several of the highest bands be zero (or at least of order \(1\times 10^{-8}\)) to ensure enough variational freedom was available in the calculation. Consequently, it is suggested that fifteen to twenty extra bands be used for a SIGMA of 0.20. In any case, it should be determined that enough bands were used by examination of the occupancies. It is undesirable to have too many extra bands, as this will add computational time.
Below we show the effect of SIGMA on the band occupancies.
from vasp import Vasp
from ase import Atom, Atoms
import matplotlib.pyplot as plt
import numpy as np
a = 3.61
atoms = Atoms([Atom('Cu', (0, 0, 0))],
cell=0.5 * a * np.array([[1.0, 1.0, 0.0],
[0.0, 1.0, 1.0],
[1.0, 0.0, 1.0]])).repeat((2, 2, 2))
SIGMA = [0.001, 0.05, 0.1, 0.2, 0.5]
for sigma in SIGMA:
calc = Vasp(label=f'bulk/Cu-sigma-{sigma}',
xc='PBE',
encut=350,
kpts=[4, 4, 4],
ismear=-1,
sigma=sigma,
nbands=9 * 8,
atoms=atoms)
if calc.potential_energy is not None:
nbands = calc.parameters['nbands']
nkpts = len(calc.get_ibz_k_points())
occ = np.zeros((nkpts, nbands))
for i in range(nkpts):
occ[i, :] = calc.get_occupation_numbers(kpt=i)
max_occ = np.max(occ, axis=0) #axis 0 is columns
plt.plot(range(nbands), max_occ, label=f'$\\sigma = {sigma}$')
plt.xlabel('band number')
plt.ylabel('maximum occupancy (electrons)')
plt.ylim([-0.1, 2.1])
plt.legend(loc='best');
The number of bands#
In the last figure, it is evident that due to the smearing of the electronic states you need to have extra bands to accommodate the electrons above the Fermi level, and the higher the incar:SIGMA value is, the more bands you need. You need enough bands so that the highest energy bands are unoccupied, and VASP will give you a warning that looks like this:
-----------------------------------------------------------------------------
| |
| ADVICE TO THIS USER RUNNING 'VASP/VAMP' (HEAR YOUR MASTER'S VOICE ...): |
| |
| Your highest band is occupied at some k-points! Unless you are |
| performing a calculation for an insulator or semiconductor, without |
| unoccupied bands, you have included TOO FEW BANDS!! Please increase |
| the parameter NBANDS in file 'INCAR' to ensure that the highest band |
| is unoccupied at all k-points. It is always recommended to |
| include a few unoccupied bands to accelerate the convergence of |
| molecular dynamics runs (even for insulators or semiconductors). |
| Because the presence of unoccupied bands improves wavefunction |
| prediction, and helps to suppress 'band-crossings.' |
| Following all k-points will be listed (with the Fermi weights of |
| the highest band given in paranthesis) ... : |
| |
| 6 (-0.01472) |
| 8 (-0.01413) |
| 13 (-0.01733) |
| 14 (-0.01838) |
| |
| The total occupancy of band no. 49 is -0.00932 electrons ... |
| |
-----------------------------------------------------------------------------
We tell VASP the number of bands to use with the incar:NBANDS keyword. VASP will set the NBANDS automatically if you do not provide a value, but this is in general bad practice (even though it is often done in this book!). There are a few general guidelines for setting NBANDS. First we recognize that a band can only have two electrons in it (one spin up, and one spin down) in an calculation without spin-polarization, or one electron per band for a spin-polarized calculation (note that spin-polarization doubles the number of bands). There absolutely must be enough bands to accommodate all the electrons, so the minimum number of bands is int(ceil(nelectrons/2)).
Here is an example of what this equation does.
import numpy as np
print(int(np.ceil(50 / 2.)))
print(int(np.ceil(51 / 2.)))
25
26
However, due to the smearing, the minimum number of bands is almost never enough, and we always add more bands. The default behavior in VASP is:
| non-spin polarized|NELECT/2 + NIONS/2| | spin-polarized|0.6*NELECT + NMAGIONS|
These do not always work, especially for small molecular systems where NIONS/2 may be only 1, or transition metals where it may be necessary to add up to 2*NIONS extra bands.
To figure out how many bands you need, it is necessary to know how many electrons are in your calculation. The func:Vasp.getvalenceelectrons provides this for you. Alternatively, you can look in the [BROKEN LINK: *Recommended%20values%20for%20ENCUT%20and%20valence%20electrons%20for%20different%20POTCAR%20files] for a table listing the number of valence electrons for each POTCAR file. Armed with this information you can set NBANDS the way you want.
from vasp import Vasp
from ase import Atom, Atoms
atoms = Atoms([Atom('Cu', [0.000, 0.000, 0.000])],
cell= [[1.818, 0.000, 1.818],
[1.818, 1.818, 0.000],
[0.000, 1.818, 1.818]])
calc = Vasp(label='bulk/alloy/cu',
xc='PBE',
encut=350,
kpts=(13, 13, 13),
nbands=9,
ibrion=2,
isif=4,
nsw=10,
atoms=atoms)
# Cu has 11 valence electrons (3d10 4s1) per the POTCAR
print('Cu valence electrons: 11')
print(calc.potential_energy)
Cu valence electrons: 11
-3.73437228
For this calculation we need at least 6 bands (11/2=5.5 which is rounded up to 6) and we need to include some extra bands. The default rule would only add half a band, which is not enough. We add three additional bands. This system is so small it does not substantially increase the computational cost.
If you are too trifling to do that much work, you can use the Vasp.set_nbands to automatically set the number of bands. This function takes an argument N to set the number of bands to N, or an argument f to set the NBANDS according to the formula \(nbands = int(nelectrons/2 + len(atoms)*f)\). The default value of f is 1.5. If you want the default VASP behavior, set f=0.5. For transition metals, it may be required that f=2. This function does not consider whether the calculation is spin-polarized or not. Here is an example of using func:Vasp.setnbands.
from vasp import Vasp
from ase import Atom, Atoms
atoms = Atoms([Atom('Cu', [0.000, 0.000, 0.000])],
cell=[[1.818, 0.000, 1.818],
[1.818, 1.818, 0.000],
[0.000, 1.818, 1.818]])
calc = Vasp(label='bulk/alloy/cu',
xc='PBE',
encut=350,
kpts=(13, 13, 13),
ibrion=2,
isif=4,
nsw=10,
atoms=atoms)
calc.set_nbands(f=7)
calc.write_input() # you have to write out the input for it to take effect
print(calc)
Vasp('/home/jovyan/dft-book/notebooks/04-bulk-systems/bulk/alloy/cu', xc='PBE')
Note the defaults that were set.
from vasp import Vasp
from ase import Atom, Atoms
atoms = Atoms([Atom('Cu', [0.000, 0.000, 0.000])],
cell=[[1.818, 0.000, 1.818],
[1.818, 1.818, 0.000],
[0.000, 1.818, 1.818]])
calc = Vasp(label='bulk/alloy/cu-setnbands',
xc='PBE',
encut=350,
kpts=(13, 13, 13),
ibrion=2,
isif=4,
nsw=10,
atoms=atoms)
calc.set_nbands(f=3)
calc.write_input()
print(calc)
Vasp('/home/jovyan/dft-book/notebooks/04-bulk-systems/bulk/alloy/cu-setnbands', xc='PBE')
You are, of course, free to use any formula you want to set the number of bands. Some formulas I have used in the past include:
NBANDS = 0.65*NELECT + 10
NBANDS = 0.5*NELECT + 15
etc…