WIEN2k-FAQ: How to calculate atomization and cohesive energies?

Back to:

The atomization energy or cohesive energy of a compound A_xB_y is the difference between the energy of the crystal and the energies of the free atoms:

E_atomization = E_crystal - x * E_A - y * E_B

The atomic energies for elements lighter than Ne can probably be taken from case.outputst. However, for heavier atoms, due to the scalar relativistic approximation one must do an "equivalent" atomic calculation, otherwise you can get even "positive" atomization energies (compound is unstable).

Choose a FCC cell with lattice parameters of 25 - 40 bohr (depends on the size of the atom and the required accuracy, eventually make a small orthorhombic distortion to break degeneracy of the cubic p states or eg/t2g levels) and put your atom into this cell.
Use the same RMT as in your compound.
Use a spin-polarized calculation (except for group IIa (Be,Mg,..), IIb (Zn, Cd,..) and VIIIA (He, Ne,..) closed shell elements.
Use 1 k-point (Gamma) and -fermit 0.002 during initialization.
Use an "equivalent" RKmax: Suppose you have RMT_A=2.5 and RMT_B=2.0 and used RKmax=7 in the bulk calculation. Then you should do atom A with RKmax=7/2*2.5=8.75 and atom B with RKmax=7.
Initialize using: init -b -sp -numk 1 -rkmax 8.75 -fermit 0.002
Use iterative diagonalization (runsp -it) and make sure that RKmax is not reduced due to NMATMAX (see :WAR in the corresponding scf file).

PS: Of course you must use identical E_XC and other (non-default) computational parameters !!!
PPS: In case the scf stops with QTL-B error in the very first cycle, it probably comes because in an atom EF is usually a negative number. Modify case.in1 and put eg. -0.5 as EF in case.in1.
PPPS: Since for such "empty" cells EF will be most likely at negative energies (eg. -0.3 or similar), you may have to adjust the energy parameters in case.in1. Once the scf-cycle has finished, you may need to readjust these energy parameters by hand (compare eigenvalues, fermi energy and :EP labels from your scf file).

The formation energy (at T=0 K) of a compound A_xB_y is the difference between the energies of this crystal and the stable phases of the elements:

E_cohesive = E_crystal - x * E_A - y * E_B

If A is eg. Fe, do a bcc Fe calculation (including Vol-optimization) for E_A, and so on for other elements.
If B is oxygen (hydrogen, nitrogen,...), one usually takes the O₂ molecule as reference. This leads to a bit cumbersome and expensive procedure:

Choose a FCC cell with lattice parameters of 25 - 40 bohr (depends on the size of the atom and the required accuracy), make the c axis longer by about 3 bohr and put your molecule oriented along z into this cell (eg. for PBE and a_lat=27,27,30 bohr, use "1" atom (Oxygen), "MULT=2" (0,0,.5385 and 0,0,0.4615). Make sure you have inversion symmetry !!!
The short bond distance of O₂ requires to set very small RMT values (~1.1 bohr, or for PBE 1.14 bohr is possible), which makes the calculation expensive. I recommend to restrict RKmax to 6 in that case, which leads to a matrix size of ~15000. (RKmax=6.5 leads to matrices of 19000 and requires at least 4/6 GB memory without/with iterative diagonalization. Note: iterative diagonalization reduces the cpu time of lapw1 from 7:30 min to 0:28 min on my PC !!) Make sure that RKmax is not reduced due to NMATMAX (see :WAR in the corresponding scf file). If it is reduced, either increase NMATMAX in param.inc of SRC_lapw1 and recompile; or - better - use mpi parallelization on a couple of cores, which will automatically increase NMATMAX by sqrt(N-cores). Please note: for all systems with such small RMTs one should increase GMAX to about 20 (gives a significant E-difference, although the density itself is not affected at all).
Follow the procedure for free atoms as described above (spin-polarized for O₂, but not for H₂ or N₂), but in addition optimize the X₂ bond distance (the distance for O₂ with PBE is ok, but of course not for other functionals or molecules) (runsp -it -min).
To have consistent total energies with "equivalent RKmax convergence", you must at the end (I would first do more complicated structure optimizations with "normal" RMTs in order to save time) also rerun:
- your compound AB (at minimum energy structure only) with the same small O sphere and the same RKmax (6 in the example above)
- your compound A (at optimized volume only) with RKmax = 6/RMT_O*RMT_A.
- PS: Many people use the experimental formation energy of O₂, since PBE calculations are not very accurate. If you want to stay "ab initio" and still be reasonable accurate, consider the SCAN meta-GGA.

Peter Blaha,