Determining bulk structures:PROPERTIES:

What we typically mean by determining bulk structures includes the following:

• What is the most stable crystal structure for a material?

• What is the lattice constant of fcc Cu?

• What are the lattice parameters and internal atom parameters for TiO₂?

All of these questions can often be addressed by finding the volume, shape and atomic positions that minimize the total energy of a bulk system. This is true at 0K. At higher temperatures, one must consider minimizing the free energy, rather than the internal energy.

fcc/bcc crystal structures

The fcc and bcc structures are simple. They only have one degree of freedom: the lattice constant. In this section we show how to calculate the equilibrium volume of each structure, and determine which one is more stable. We start with the fcc crystal structure of Cu. We will manually define the crystal structure based on the definitions in Kittel cite:kittel (Chapter 1).

from vasp import Vasp
from ase import Atom, Atoms
import numpy as np

# fcc
LC = [3.5, 3.55, 3.6, 3.65, 3.7, 3.75]
fcc_energies = []
ready = True
for a in LC:
    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]]))

    calc = Vasp(label=f'bulk/Cu-{a}',
                xc='PBE',
                encut=350,
                kpts=[8, 8, 8],
                atoms=atoms)

    e = atoms.get_potential_energy()
    fcc_energies.append(e)

calc.stop_if(None in fcc_energies)

import matplotlib.pyplot as plt
plt.plot(LC, fcc_energies)
plt.xlabel('Lattice constant ($\AA$)')
plt.ylabel('Total energy (eV)')

print('#+tblname: cu-fcc-energies')
print(r'| lattice constant ($\AA$) | Total Energy (eV) |')
for lc, e in zip(LC, fcc_energies):
    print(f'| {lc} | {e} |')

Use the data in the table above to plot the total energy as a function of the lattice constant. Fit a cubic polynomial to the data, and find the volume that minimizes the total energy.

If you want to know the lattice constant that gives the lowest energy, you would fit an Equations of State to the data. Here is an example using mod:ase.utils.eos.

from vasp import Vasp
from ase.eos import EquationOfState
LC = [3.5, 3.55, 3.6, 3.65, 3.7, 3.75]
energies = []
volumes = []
for a in LC:
    calc = Vasp(label=f'bulk/Cu-{a}')
    atoms = calc.load_atoms()
    volumes.append(atoms.get_volume())
    energies.append(atoms.get_potential_energy())

calc.stop_if(None in energies)

eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()

print(f'''
v0 = {v0} A^3
E0 = {e0} eV
B  = {B} eV/A^3''')

eos.plot('images/Cu-fcc-eos.png')

Before we jump into the bcc calculations, let us consider what range of lattice constants we should choose. The fcc lattice is close-packed, and the volume of the primitive cell is \(V = 1/4 a^3\) or about 11.8 ų/atom. The volume of the equilibrium bcc primitive cell will probably be similar to that. The question is: what bcc lattice constant gives that volume? The simplest way to answer this is to compute the answer. We will make a bcc crystal at the fcc lattice constant, and then compute the scaling factor needed to make it the right volume.

from ase import Atom, Atoms
import numpy as np
a = 3.61 # lattice constant

atoms = Atoms([Atom('Cu', [0,0,0])],
              cell=0.5 * a*np.array([[ 1.0,  1.0, -1.0],
                                     [-1.0,  1.0,  1.0],
                                     [ 1.0, -1.0,  1.0]]))

print(f'BCC lattice constant = {a * (11.8 / atoms.get_volume())**(1./3.):1.3f} Ang')

Now we run the equation of state calculations.

from vasp import Vasp
from ase import Atom, Atoms
import numpy as np

LC = [2.75, 2.8, 2.85, 2.9, 2.95, 3.0]

for a in LC:
    atoms = Atoms([Atom('Cu', [0, 0, 0])],
                  cell=0.5 * a * np.array([[ 1.0,  1.0, -1.0],
                                           [-1.0,  1.0,  1.0],
                                           [ 1.0, -1.0,  1.0]]))

    calc = Vasp(label=f'bulk/Cu-bcc-{a}',
                xc='PBE',
                encut=350,
                kpts=[8, 8, 8],
                atoms=atoms)
    print(calc.potential_energy)

Finally, we will compare the two crystal structures.

from vasp import Vasp

# bcc energies and volumes
bcc_LC = [2.75, 2.8, 2.85, 2.9, 2.95, 3.0]
bcc_volumes = []
bcc_energies = []
for a in bcc_LC:
    calc = Vasp(label=f'bulk/Cu-bcc-{a}')
    atoms = calc.load_atoms()
    bcc_volumes.append(atoms.get_volume())
    bcc_energies.append(atoms.get_potential_energy())

# fcc energies and volumes
fcc_LC = [3.5, 3.55, 3.6, 3.65, 3.7, 3.75]
fcc_volumes = []
fcc_energies =[]
for a in fcc_LC:
    calc = Vasp(label=f'bulk/Cu-{a}')
    atoms = calc.load_atoms()
    fcc_volumes.append(atoms.get_volume())
    fcc_energies.append(atoms.get_potential_energy())

import matplotlib.pyplot as plt
plt.plot(fcc_volumes, fcc_energies, label='fcc')
plt.plot(bcc_volumes, bcc_energies, label='bcc')

plt.xlabel('Atomic volume ($\AA^3$/atom)')
plt.ylabel('Total energy (eV)')
plt.legend()

# print table of data
print('#+tblname: bcc-data')
print('#+caption: Total energy vs. lattice constant for BCC Cu.')
print('| Lattice constant (\AA$^3$) | Total energy (eV) |')
print('|-')
for lc, e in zip(bcc_LC, bcc_energies):
    print(f'| {lc} | {e} |')

Lattice constant (Å\(^3\))	Total energy (eV)
2.75	-3.59937543
2.8	-3.67930795
2.85	-3.71927399
2.9	-3.72637899
2.95	-3.70697046
3.0	-3.66645678

Use the data for FCC and BCC Cu to plot the total energy as a function of the lattice constant.

Note we plot the energy vs. atomic volume. That is because the lattice constants of the two crystal structures are very different. It also shows that the atomic volumes in the two structures are similar.

What can we say here? The fcc structure has a lower energy than the bcc structure, so we can conclude the fcc structure is more favorable. In fact, the fcc structure is the experimentally found structure for Cu. Some caution is in order; if you run these calculations at a \(4 \times 4 \times 4\) \(k\)-point grid, the bcc structure is more stable because the results are not converged!

Compute the energy vs. volume for fcc and bcc Cu for different $k$-point grids. Determine when each result has converged, and which structure is more stable.

What can we say about the relative stability of fcc to hcp? Nothing, until we calculate the hcp equation of state.

Optimizing the hcp lattice constant

The hcp lattice is more complicated than the fcc/bcc lattices because there are two lattice parameters: \(a\) and \(c\) or equivalently: \(a\) and \(c/a\). We will start by making a grid of values and find the set of parameters that minimizes the energy. See Figure ref:fig:ru-e-ca.

from ase.lattice.hexagonal import HexagonalClosedPacked
from vasp import Vasp
import matplotlib.pyplot as plt

atoms = HexagonalClosedPacked(symbol='Ru',
                              latticeconstant={'a': 2.7,
                                               'c/a': 1.584})

a_list = [2.5, 2.6, 2.7, 2.8, 2.9]
covera_list = [1.4, 1.5, 1.6, 1.7, 1.8]

for a in a_list:
    energies = []
    for covera in covera_list:

        atoms = HexagonalClosedPacked(symbol='Ru',
                              latticeconstant={'a': a,
                                               'c/a': covera})

        wd = f'bulk/Ru/{a:1.2f}-{covera:1.2f}'

        calc = Vasp(wd,
                    xc='PBE',
                    # the c-axis is longer than the a-axis, so we use
                    # fewer kpoints.
                    kpts=[6, 6, 4],
                    encut=350,
                    atoms=atoms)

        energies.append(atoms.get_potential_energy())
    if not None in energies:
        plt.plot(covera_list, energies, label=rf'a={a} $\AA$')

plt.xlabel('$c/a$')
plt.ylabel('Energy (eV)')
plt.legend()
plt.savefig('images/Ru-covera-scan.png')

It looks like there is a minimum in the a=2.7 Å curve, at a \(c/a\) ratio of about 1.6. We can look at the same data in a contour plot which shows more clearly there is minimum in all directions near that point (Figure ref:fig:ru-contourf).

from vasp import Vasp
import matplotlib.pyplot as plt
import numpy as np

x = [2.5, 2.6, 2.7, 2.8, 2.9]
y = [1.4, 1.5, 1.6, 1.7, 1.8]

X,Y = np.meshgrid(x, y)
Z = np.zeros(X.shape)

for i,a in enumerate(x):
    for j,covera in enumerate(y):
        wd = f'bulk/Ru/{a:1.2f}-{covera:1.2f}'
        calc = Vasp(wd)
        Z[i][j] = calc.potential_energy

calc.stop_if(None in Z)

cf = plt.contourf(X, Y, Z, 20,
                  cmap=plt.cm.jet)

cbar = plt.colorbar(cf)
cbar.ax.set_ylabel('Energy (eV)')

plt.xlabel('$a$ ($\AA$)')
plt.ylabel('$c/a$')

plt.legend()
plt.savefig('images/ru-contourf.png')
plt.show()

Complex structures with internal degrees of freedom

A unit cell has six degrees of freedom: the lengths of each unit cell vector, and the angle between each vector. There may additionally be internal degrees of freedom for the atoms. It is impractical to try the approach used for the hcp Ru on anything complicated. Instead, we rely again on algorithms to optimize the unit cell shape, volume and internal degrees of freedom. It is usually not efficient to make a wild guess of the geometry and then turn VASP loose on to optimize it. Instead, the following algorithm works pretty well.

1. Find the volume (at constant shape, with relaxed ions) that minimizes the total energy ([BROKEN LINK: incar:ISIF]=2). The goal here is to just get an idea of where the right volume is.

2. Using the results from step 1 as a starting point, perform a set of calculations at constant volume around the minimum from step 1, but the shape and internal atom positions are allowed to change ([BROKEN LINK: incar:ISIF]=4).

3. Finally, do a final calculation near the minimum energy allowing the volume to also change. ([BROKEN LINK: incar:ISIF]=3).

This multistep process is pretty reasonable to get a converged structure pretty quickly. It is not foolproof, however, and if you have materials such as graphite it may not work well. The problem with graphite is that it is a layered compound that is held together by weak van der waal type forces which are not modeled well by typical GGA functionals. Thus the change in energy due to a volume change is larger in the plane of the graphite sheet than in the direction normal to the sheet. With a typical GGA, the sheets may just move apart until they do not interact any more.

We will illustrate the process on a well-behaved system (rutile TiO₂) which has two lattice parameters and one internal degree of freedom. There are a few subtle points to mention in doing these calculations. The VASP manual recommends that you set incar:PREC to ‘high’, and that ENCUT be set to 1.3*max(ENMAX) of the pseudopotentials. This is necessary to avoid problems caused by small basis sets when the volume changes, and Pulay stress. It is important to ensure that the energies are reasonably converged with respect to k-point grids. Hence, it can be a significant amount of work to do this right! Let us start with determining the ENCUT value that is appropriate for TiO₂.

grep ENMAX $VASP_PP_PATH/potpaw_PBE/Ti/POTCAR
grep ENMAX $VASP_PP_PATH/potpaw_PBE/O/POTCAR

ENMAX  =  178.330; ENMIN  =  133.747 eV
ENMAX  =  400.000; ENMIN  =  300.000 eV

ENMAX  =  178.330; ENMIN  =  133.747 eV
ENMAX  =  400.000; ENMIN  =  300.000 eV

According to the manual, we should use ENCUT = 1.3*400 = 520 eV for good results.

Now we consider the k-point convergence. The lattice vectors of the rutile TiO₂ structure are not all the same length, which means it is not essential that we use the same number of k-points in each direction. For simplicity, however, we do that here.

# step 1 frozen atoms and shape at different volumes
from ase import Atom, Atoms
import numpy as np
from vasp import Vasp
import matplotlib.pyplot as plt

'''
create a TiO2 structure from the lattice vectors at
http://cst-www.nrl.navy.mil/lattice/struk/c4.html
This site does not exist anymore.
'''
a = 4.59  # experimental degrees of freedom.
c = 2.96
u = 0.3  # internal degree of freedom!

#primitive vectors
a1 = a * np.array([1.0, 0.0, 0.0])
a2 = a * np.array([0.0, 1.0, 0.0])
a3 = c * np.array([0.0, 0.0, 1.0])

atoms = Atoms([Atom('Ti', [0., 0., 0.]),
               Atom('Ti', 0.5 * a1 + 0.5 * a2 + 0.5 * a3),
               Atom('O', u * a1 + u * a2),
               Atom('O', -u * a1 - u * a2),
               Atom('O', (0.5 + u) * a1 + (0.5 - u) * a2 + 0.5 * a3),
               Atom('O', (0.5 - u) * a1 + (0.5 + u) * a2 + 0.5 * a3)],
              cell=[a1, a2, a3])

KPOINTS = [2, 3, 4, 5, 6, 7, 8]
energies = []

ready = True
for k in KPOINTS:
    calc = Vasp(label=f'bulk/tio2/kpts-{k}',
                encut=520,
                kpts=[k, k, k],
                xc='PBE',
                sigma=0.05,
                atoms=atoms)

    energies.append(atoms.get_potential_energy())

calc.stop_if(None in energies)

plt.plot(KPOINTS, energies)
plt.xlabel('number of k-points in each vector')
plt.ylabel('Total energy (eV)')
plt.savefig('images/tio2-kpt-convergence.png')

A k-point grid of \(5 \times 5 \times 5\) appears suitable for reasonably converged results. Now we proceed with step 1: Compute the total energy of the unit cell allowing internal degrees of freedom to relax, but keeping a constant cell shape.

# step 1 frozen atoms and shape at different volumes
from ase import Atom, Atoms
import numpy as np
from vasp import Vasp
import matplotlib.pyplot as plt

'''
create a TiO2 structure from the lattice vectors at
http://cst-www.nrl.navy.mil/lattice/struk/c4.html
'''
a = 4.59  # experimental degrees of freedom.
c = 2.96
u = 0.3  # internal degree of freedom!

#primitive vectors
a1 = a * np.array([1.0, 0.0, 0.0])
a2 = a * np.array([0.0, 1.0, 0.0])
a3 = c * np.array([0.0, 0.0, 1.0])

atoms = Atoms([Atom('Ti', [0., 0., 0.]),
               Atom('Ti', 0.5 * a1 + 0.5 * a2 + 0.5 * a3),
               Atom('O', u * a1 + u * a2),
               Atom('O', -u * a1 - u * a2),
               Atom('O', (0.5 + u) * a1 + (0.5 - u) * a2 + 0.5 * a3),
               Atom('O', (0.5 - u) * a1 + (0.5 + u) * a2 + 0.5 * a3)],
              cell=[a1, a2, a3])

v0 = atoms.get_volume()
cell0 = atoms.get_cell()

factors = [0.9, 0.95, 1.0, 1.05, 1.1] #to change volume by

energies, volumes = [], []

ready = True
for f in factors:
    v1 = f * v0
    cell_factor = (v1 / v0)**(1. / 3.)

    atoms.set_cell(cell0 * cell_factor, scale_atoms=True)

    calc = Vasp(label=f'bulk/tio2/step1-{f:1.2f}',
                encut=520,
                kpts=[5, 5, 5],
                isif=2, # relax internal degrees of freedom
                ibrion=1,
                nsw=50,
                xc='PBE',
                sigma=0.05,
                atoms=atoms)

    energies.append(atoms.get_potential_energy())
    volumes.append(atoms.get_volume())

calc.stop_if(None in energies)

plt.plot(volumes, energies)
plt.xlabel('Vol. ($\AA^3)$')
plt.ylabel('Total energy (eV)')

print('#+tblname: tio2-vol-ene')
print('#+caption: Total energy of TiO_{2} vs. volume.')
print('| Volume ($\AA^3$) | Energy (eV) |')
print('|-')
for v, e in zip(volumes, energies):
    print(f'| {v} | {e} |')

Volume (\(\AA^3\))	Energy (eV)
56.1254185488	-51.81932158
59.2434971663	-52.46118036
62.361576	-52.76017908
65.4796549456	-52.80043775
68.5977335623	-52.64628895

Now, we know the minimum energy volume is near 64 ų. You could at this point fit an equation of state to find that minimum. However, we now want to use these initial starting points for a second round of optimization where we allow the unit cell shape to change, at constant volume: ISIF=4.

from vasp import Vasp

calc = Vasp(label='bulk/tio2/step1-0.90')
calc.clone('bulk/tio2/step2-0.90')

print(calc.set(isif=4))
print(calc.calculation_required(atoms, ['potential_energy']))

{'isif': 4}
True

from vasp import Vasp

factors = [0.9, 0.95, 1.0, 1.05, 1.1]  # to change volume by

energies1, volumes1 = [], []  # from step 1
energies, volumes = [], []  # for step 2
ready = True
for f in factors:
    calc = Vasp(label=f'bulk/tio2/step1-{f:1.2f}')
    atoms = calc.load_atoms()
    energies1.append(atoms.get_potential_energy())
    volumes1.append(atoms.get_volume())

    calc.clone(f'bulk/tio2/step2-{f:1.2f}')
    calc.set(isif=4)
    # You have to get the atoms again.
    atoms = calc.load_atoms()

    energies.append(atoms.get_potential_energy())
    volumes.append(atoms.get_volume())

print(energies, volumes)
calc.stop_if(None in energies)

import matplotlib.pyplot as plt
plt.plot(volumes1, energies1, volumes, energies)
plt.xlabel('Vol. ($\AA^3)$')
plt.ylabel('Total energy (eV)')
plt.legend(['step 1', 'step 2'], loc='best')
plt.savefig('images/tio2-step2.png')

The take away message here is that the total energy slightly decreases when we allow the unit cell shape to change, especially for the larger unit cell deformation. This has little effect on the minimum volume, but would have an effect on the bulk modulus, which is related to the curvature of the equation of state. At this point, you could fit an equation of state to the step 2 data, and estimate the volume at the minimum volume, and recalculate the total energy at that volume.

An alternative is a final calculation with ISIF=3, which optimizes the unit cell volume, shape and internal coordinates. It looks like the calculation at bulk/tio2/step2-1.05 is close to the minimum, so we will use that as a starting point for the final calculation.

from vasp import Vasp
import os

calc = Vasp(label='bulk/tio2/step2-1.05')

# Check if source calculation exists before cloning
source_dir = calc.directory
if os.path.exists(os.path.join(source_dir, 'CONTCAR')) or os.path.exists(os.path.join(source_dir, 'POSCAR')):
    calc.clone('bulk/tio2/step3')
    
    # Load atoms from the cloned directory and set new parameters
    calc = Vasp(label='bulk/tio2/step3')
    atoms = calc.load_atoms()
    calc.set(isif=3)
    calc.atoms = atoms
    calc.calculate()
    print(calc)

    import spglib
    print(f'\nThe spacegroup is {spglib.get_spacegroup(atoms)}')
else:
    print('Source calculation bulk/tio2/step2-1.05 not yet complete. Run previous cells first.')

#+begin_example

################ VASP CALCULATION SUMMARY *

This is the final result. You can see that the forces on all the atoms are less than 0.01 eV/Å, and the stress is also very small. The final volume is close to where we expect it to be based on steps 1 and 2. The space group is still correct. The lattice vectors are close in length to the experimentally known values, and the angles between the vectors has not changed much. Looks good!

As a final note, the VASP manual recommends you do not use the final energy directly from the calculation, but rather run a final calculation with incar:ISMEAR set to -5. Here we examine the effect.

from vasp import Vasp
import os

calc = Vasp(label='bulk/tio2/step3')

# Check if step3 calculation exists
if not os.path.exists(os.path.join(calc.directory, 'CONTCAR')) and not os.path.exists(os.path.join(calc.directory, 'POSCAR')):
    print('Step3 calculation not yet complete. Run previous cells first.')
else:
    atoms = calc.load_atoms()
    print('default ismear: ', atoms.get_potential_energy())

    calc.clone('bulk/tio2/step4')
    calc.set(ismear=-5,
             nsw=0)
    atoms = calc.load_atoms()
    print('ismear=-5:      ', atoms.get_potential_energy())

Step3 calculation not yet complete. Run previous cells first.

The difference here is on the order of a meV. That does not seem significant here. I suspect the recommended practice stems from early days when much smaller ENCUT values were used and changes in the number of basis functions were more significant.

Effect of XC on bulk properties

The exchange correlation functional can significantly affect computed bulk properties. Here, we examine the effect on the bulk lattice constant of Pd (exp. 3.881). An excellent review of this can be found in cite:mattsson-084714. We examine several functionals. The xc keyword in Vasp is used to select the POTCARs. Let us consider the LDA functional first.

from vasp import Vasp
from ase import Atom, Atoms
from ase.utils.eos import EquationOfState
import numpy as np

LC = [3.75, 3.80, 3.85, 3.90, 3.95, 4.0, 4.05, 4.1]

volumes, energies = [], []
for a in LC:
    atoms = Atoms([Atom('Pd', (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]]))
    calc = Vasp(label=f'bulk/Pd-LDA-{a}',
                encut=350,
                kpts=[12, 12, 12],
                xc='LDA',
                atoms=atoms)

    e = atoms.get_potential_energy()
    energies.append(e)
    volumes.append(atoms.get_volume())

calc.stop_if(None in energies)

eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print(f'LDA lattice constant is {(4*v0)**(1./3.):1.3f} Ang^3')

For a GGA calculation, it is possible to specify which functional you want via the incar:GGA tag. This tag was designed to use the LDA POTCAR files, but with a GGA functional. We will consider four different functionals here.

from vasp import Vasp
from ase import Atom, Atoms
from ase.utils.eos import EquationOfState
import numpy as np

LC = [3.75, 3.80, 3.85, 3.90, 3.95, 4.0, 4.05, 4.1]

GGA = {'AM': 'AM05',
       'PE': 'PBE',
       'PS': 'PBEsol',
       'RP': 'RPBE'}

for key in GGA:
    volumes, energies = [], []
    for a in LC:
        atoms = Atoms([Atom('Pd', (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]]))
        calc = Vasp(label=f'bulk/Pd-GGA-{key}-{a}',
                    encut=350,
                    kpts=[12, 12, 12],
                    xc='LDA',
                    gga=key,
                    atoms=atoms)

        e = atoms.get_potential_energy()
        energies.append(e)
        volumes.append(atoms.get_volume())

    if None not in energies:
        eos = EquationOfState(volumes, energies)
        v0, e0, B = eos.fit()
        print(f'{GGA[key]:6s} lattice constant is {(4*v0)**(1./3.):1.3f} Ang^3')
    else:
        print(energies, LC)
        print(f'{GGA[key]} is not ready')

These results compare very favorably to those in cite:mattsson-084714. It is typical that LDA functionals underestimate the lattice constants, and that GGAs tend to overestimate the lattice constants. PBEsol and AM05 were designed specifically for solids, and for Pd, these functionals do an exceptional job of reproducing the lattice constants. RPBE is particularly bad at the lattice constant, but it has been reported to be a superior functional for reactivity cite:hammer1999:improv-pbe.