Cohesive energy

The cohesive energy is defined as the energy to separate neutral atoms in their ground electronic state from the solid at 0K at 1 atm. We will compute this for rhodium. Rh is normally an fcc metal, so we will use that structure and let VASP find the equilibrium volume for us.

from vasp import Vasp
from ase.lattice.cubic import FaceCenteredCubic
from ase import Atoms, Atom

# bulk system
atoms = FaceCenteredCubic(directions=[[0, 1, 1],
                                      [1, 0, 1],
                                      [1, 1, 0]],
                                      size=(1, 1, 1),
                                      symbol='Rh')

calc = Vasp(label='bulk/bulk-rh',
            xc='PBE',
            encut=350,
            kpts=(4, 4, 4),
            isif=3,
            ibrion=2,
            nsw=10,
            atoms=atoms)
bulk_energy = atoms.get_potential_energy()

# atomic system
atoms = Atoms([Atom('Rh',[5, 5, 5])],
              cell=(7, 8, 9))

calc = Vasp(label='bulk/atomic-rh',
            xc='PBE',
            encut=350,
            kpts=(1, 1, 1),
            atoms=atoms)
atomic_energy = atoms.get_potential_energy()


cohesive_energy = atomic_energy - bulk_energy
print('The cohesive energy is {0:1.3f} eV'.format(cohesive_energy))

The cohesive energy is 6.183 eV

According to Kittel cite:kittel, the cohesive energy of Rh is 5.75 eV. There are a few reasons we may have discrepancy here:

1. The k-point grid used in the bulk state is not very dense. However, you can see below that the total energy is pretty converged by a \(6 \times 6 \times 6\) \(k\)-point grid.

2. We did not check for convergence with the planewave cutoff.

3. We neglected spin on the atomic state. Rh in the atomic state has this electronic structure: [Kr] 4d8 5s1 and is a doublet.

First we consider the k-point convergence. index:convergence!KPOINTS

from vasp import Vasp

calc = Vasp(label='bulk/atomic-rh')
atomic_energy = calc.potential_energy

calc = Vasp(label='bulk/bulk-rh')
atoms = calc.load_atoms()

kpts = [3, 4, 6, 9, 12, 15, 18]

calcs = [Vasp(label='bulk/bulk-rh-kpts-{0}'.format(k),
                xc='PBE',
                encut=350,
                kpts=(k, k, k),
                atoms=atoms)
         for k in kpts]

energies = [calc.potential_energy for calc in calcs]

calcs[0].stop_if(None in energies)

for k, e in zip(kpts, energies):
    print('({0:2d}, {0:2d}, {0:2d}):'
          ' cohesive energy = {1} eV'.format(k,
                                             e - atomic_energy))

( 3,  3,  3): cohesive energy = -4.80100581 eV
( 4,  4,  4): cohesive energy = -6.18383533 eV
( 6,  6,  6): cohesive energy = -6.210267630000001 eV
( 9,  9,  9): cohesive energy = -6.20521195 eV
(12, 12, 12): cohesive energy = -6.21304005 eV
(15, 15, 15): cohesive energy = -6.2130027100000005 eV
(18, 18, 18): cohesive energy = -6.21410349 eV

Using only 1 k-point for the bulk energy is a terrible approximation! It takes at least a 6 × 6 × 6 grid to get the total energy converged to less than 10 meV. Note we do not need to check the k-point convergence of the atomic state because it is surrounded by vacuum on all sides, and so there should not be any dispersion in the bands.

We will examine the magnetic state next.

from vasp import Vasp
from ase.lattice.cubic import FaceCenteredCubic
from ase import Atoms, Atom
# bulk system
atoms = FaceCenteredCubic(directions=[[0, 1, 1],
                                      [1, 0, 1],
                                      [1, 1, 0]],
                                     size=(1, 1, 1),
                                     symbol='Rh')

calc = Vasp(label='bulk/bulk-rh',
            xc='PBE',
            encut=350,
            kpts=(4, 4, 4),
            isif=3,
            ibrion=2,
            nsw=10,
            atoms=atoms)
bulk_energy = atoms.get_potential_energy()

# atomic system
atoms = Atoms([Atom('Rh',[5, 5, 5], magmom=1)],
              cell=(7, 8, 9))

calc = Vasp(label='bulk/atomic-rh-sp',
            xc='PBE',
            encut=350,
            kpts=(1, 1, 1),
            ispin=2,
            atoms=atoms)
atomic_energy = atoms.get_potential_energy()

calc.stop_if(None in [atomic_energy, bulk_energy])

cohesive_energy = atomic_energy - bulk_energy
print('The cohesive energy is {0:1.3f} eV'.format(cohesive_energy))

The cohesive energy is 5.904 eV

Again, the value in Kittel cite:kittel is 5.75 eV which is very close to this value. Finally, it is also possible there is a lower energy non-spherical atom energy; we did not check that at all (see [BROKEN LINK: *Estimating%20triplet%20oxygen%20dissociation%20energy%20with%20low%20symmetry]).