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

©2001 by P. Blaha, K. Schwarz and J. Luitz

Back to:

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

Eatomization = Ecrystal - x * EA - y * EB

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 RMTA=2.5 and RMTB=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 EXC 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 AxBy is the difference between the energies of this crystal and the stable phases of the elements:

Ecohesive = Ecrystal - x * EA - y * EB

If A is eg. Fe, do a bcc Fe calculation (including Vol-optimization) for EA, and so on for other elements.
If B is oxygen (hydrogen, nitrogen,...), one usually takes the O2 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 alat=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 O2 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 O2, but not for H2 or N2), but in addition optimize the X2 bond distance (the distance for O2 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/RMTO*RMTA.
• PS: Many people use the experimental formation energy of O2, 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,