Crystal binding and elastic properties

Why do atoms form solids at all, and why do different solids hold together with such different strengths? Solid argon melts at 84 K; sodium chloride at 1074 K; tungsten at 3695 K. The answer lies in the balance between attraction and repulsion. At large separations, atoms attract each other through one of several mechanisms; at short separations, their electron clouds overlap and the Pauli principle produces a steep repulsion. The crystal settles at the spacing where the two balance. This lesson classifies the main types of bonding, calculates the binding energy of two important cases, a van der Waals crystal and an ionic crystal, from simple model potentials, and connects the curvature of the energy curve to the elastic stiffness of the solid. CHEM-241 introduced lattice energies through the Born–Haber cycle; here we compute them from the potentials directly.

Types of bonding

Solids are conventionally grouped into five classes by the dominant bonding mechanism, although real materials often mix them.

Van der Waals, or molecular, crystals are formed by noble gases and by molecules with closed electron shells, such as \(\text{N}_2\), \(\text{CH}_4\) and many organic compounds. The attraction comes from fluctuating dipoles: the instantaneous dipole of one atom polarizes its neighbour, and the induced dipoles attract with an energy proportional to \(-1/r^6\). The binding is weak, around 0.01 to 0.1 eV per atom, and these crystals melt at low temperatures.

Ionic crystals, such as NaCl, MgO and CaF₂, form when an electropositive atom gives up electrons to an electronegative one, leaving closed-shell ions held together by the Coulomb attraction of opposite charges. Cohesive energies are several eV per ion pair, and the crystals are hard, brittle and insulating.

Covalent crystals, such as diamond, silicon and germanium, are held together by shared electron pairs in directional bonds. The directionality produces open structures with low coordination, like the tetrahedral diamond structure of lesson 2, and cohesive energies of several eV per atom. Diamond, at 7.4 eV per atom, is among the most strongly bound of all solids.

Metallic crystals, such as sodium, copper and iron, share their valence electrons among all atoms. The attraction comes from the lowering of the electrons' kinetic energy when they spread out over the whole crystal, together with their Coulomb attraction to the ion cores. Metallic bonding is not directional, so metals tend to adopt close-packed structures and are ductile. Cohesive energies range from about 1 eV per atom for the alkali metals to 8.9 eV for tungsten, where partly filled \(d\) shells add covalent character.

Hydrogen-bonded crystals, such as ice, are held together by a hydrogen atom shared between two electronegative atoms. The bond energy is about 0.2 to 0.4 eV. Hydrogen bonds give ice its open structure and its lower density than liquid water, and they hold the two strands of DNA together.

The Lennard-Jones crystal

The simplest solids to model quantitatively are the noble-gas crystals. Two noble-gas atoms at distance \(r\) interact through the Lennard-Jones potential

where the \(r^{-6}\) term is the van der Waals attraction derived by Fritz London in 1930, and the \(r^{-12}\) term is a convenient, steep model of the repulsion from overlapping electron clouds. The parameters \(\varepsilon\) and \(\sigma\) are fitted to the properties of the gas; for argon, \(\varepsilon = 0.0104\) eV and \(\sigma = 0.340\) nm. In a crystal with closest-atom spacing \(r\), the energy per atom is half the sum of \(\phi\) over all other atoms (half, because each pair is shared by two atoms):

where \(r_{0j}\) is the distance from a reference atom to atom \(j\). The lattice sums \(A_n\) depend only on the structure. For fcc, the 12 closest neighbours contribute 12 to each, and farther shells add the rest: \(A_{12} = 12.13\) and \(A_6 = 14.45\). Minimizing \(u\) with respect to \(r\) gives

Each atom in the crystal is bound by 8.6 times the depth of the pair potential, because it has 12 close neighbours and many more distant ones; the six bonds per atom (12 neighbours, each bond shared) would give only \(6\varepsilon\).

Elastic stiffness

The curvature of the energy at its minimum determines how hard the solid is to compress. The bulk modulus is defined by \(B = -V\,\partial P/\partial V\), and since \(P = -\partial U/\partial V\) at zero temperature,

For a cubic crystal the volume per atom is proportional to \(r^3\), so \(B\) can be written in terms of \(d^2u/dr^2\). For the Lennard-Jones fcc crystal the result is \(B = 75\,\varepsilon/\sigma^3\). The program evaluates the lattice sums directly, adds a correction for the slowly converging \(r^{-6}\) tail, and applies all these formulas to argon.

import math, itertools
# Lennard-Jones solid: lattice sums for fcc, equilibrium spacing, cohesive energy, bulk modulus.
N = 12
fcc = [(i, j, k) for i, j, k in itertools.product(range(-N, N + 1), repeat=3)
       if (i + j + k) % 2 == 0 and 0 < i*i + j*j + k*k <= N * N]          # spherical cutoff of radius N
# positions in units of a/2, so the nearest-atom distance is sqrt(2) in these units
dnn = math.sqrt(2)
A6 = sum((dnn / math.sqrt(i*i + j*j + k*k)) ** 6 for i, j, k in fcc)
A12 = sum((dnn / math.sqrt(i*i + j*j + k*k)) ** 12 for i, j, k in fcc)
# the r^-6 sum converges slowly: add the smooth tail beyond the cutoff (one atom per 2 units of volume)
A6 += 2 * math.pi * dnn ** 6 / (3 * N ** 3)
A12 += 2 * math.pi * dnn ** 12 / (9 * N ** 9)
print(f"fcc lattice sums: A12 = {A12:.4f} (exact 12.1319), A6 = {A6:.4f} (exact 14.4539)")
# u(r) = 2 eps [A12 (s/r)^12 - A6 (s/r)^6] per atom
r0 = (2 * A12 / A6) ** (1 / 6)
u0 = -A6 ** 2 / (2 * A12)
# bulk modulus from B = V d2U/dV2 with V = r^3 / sqrt(2) per atom
def U(r): return 2 * (A12 / r ** 12 - A6 / r ** 6)
h = 1e-4
d2U = (U(r0 + h) - 2 * U(r0) + U(r0 - h)) / h ** 2
V0 = r0 ** 3 / math.sqrt(2)
B = r0 ** 2 * d2U / (9 * V0)
print(f"equilibrium r0 = {r0:.4f} sigma, energy per atom = {u0:.3f} eps, bulk modulus = {B:.1f} eps/sigma^3")
eps, sig = 0.0104 * 1.602176634e-19, 0.340e-9     # argon gas-phase parameters
print(f"argon: r0 = {r0 * sig * 1e9:.3f} nm (measured 0.376), cohesive energy = {-u0 * 0.0104:.3f} eV (measured 0.080),"
      f" B = {B * eps / sig ** 3 / 1e9:.2f} GPa (measured 2.7)")
print()
# Ionic crystal: Madelung constant of NaCl by Evjen's method (charge-weighted boundary)
def madelung(n):
    s = 0.0
    for i, j, k in itertools.product(range(-n, n + 1), repeat=3):
        if (i, j, k) == (0, 0, 0): continue
        w = 1.0
        for c in (i, j, k):
            if abs(c) == n: w *= 0.5
        s += w * (-1) ** (i + j + k + 1) / math.sqrt(i*i + j*j + k*k)
    return s
for n in (2, 4, 8, 12):
    print(f"Evjen cube of side {2*n:2d}: Madelung constant {madelung(n):.6f}")
# Born-Mayer energy per ion pair: U(r) = -alpha k e^2 / r + B exp(-r/rho), rho = 0.0321 nm fitted to compressibility
ke2 = 1.439964   # eV nm
alpha, r0n, rho = 1.747565, 0.2820, 0.0321
Uc = -alpha * ke2 / r0n * (1 - rho / r0n)
print(f"NaCl lattice energy (Born-Mayer) = {Uc:.2f} eV per ion pair = {-Uc * 96.485:.0f} kJ/mol (Born-Haber 787)")

fcc lattice sums: A12 = 12.1319 (exact 12.1319), A6 = 14.4538 (exact 14.4539)
equilibrium r0 = 1.0902 sigma, energy per atom = -8.610 eps, bulk modulus = 75.2 eps/sigma^3
argon: r0 = 0.371 nm (measured 0.376), cohesive energy = 0.090 eV (measured 0.080), B = 3.19 GPa (measured 2.7)

Evjen cube of side  4: Madelung constant 1.751769
Evjen cube of side  8: Madelung constant 1.747721
Evjen cube of side 16: Madelung constant 1.747574
Evjen cube of side 24: Madelung constant 1.747567
NaCl lattice energy (Born-Mayer) = -7.91 eV per ion pair = 763 kJ/mol (Born-Haber 787)

The predictions for argon are close: spacing 0.371 nm against 0.376 nm measured, cohesive energy 0.090 eV against 0.080 eV, bulk modulus 3.2 GPa against 2.7 GPa. All three errors have the same origin. Argon atoms are light, and quantum zero-point motion, ignored in the static calculation, pushes them slightly apart and reduces the binding. For neon, lighter still, the correction is larger, and for helium it is so large that the solid does not form at all at ordinary pressure: liquid helium stays liquid down to absolute zero unless compressed to about 25 atmospheres.

Ionic crystals and the Madelung constant

In an ionic crystal the dominant interaction is the Coulomb energy of the ions. For NaCl, take a sodium ion as reference. Its six closest neighbours are chlorine ions at distance \(r\), its next 12 neighbours are sodium ions at \(\sqrt2\,r\), then 8 chlorine ions at \(\sqrt3\,r\), and so on. The Coulomb energy of the reference ion is

where \(\alpha\) is the Madelung constant, after Erwin Madelung, who first evaluated it in 1918. The series as written does not converge: its terms grow, because the number of ions in each shell increases. The sum is conditionally convergent, and its value depends on the order of summation, a subtlety with physical meaning: a crystal with a net charged surface layer would have a different energy. Harald Evjen showed in 1932 how to sum it reliably: add up neutral cubes of increasing size, giving ions on the faces half weight, those on edges a quarter and those at corners an eighth, so that each cube is electrically neutral. The program does this and converges quickly to \(\alpha = 1.747565\). Other structures have other constants: 1.763 for caesium chloride, 1.638 for zinc blende.

The Coulomb attraction must be balanced by repulsion when the closed shells of neighbouring ions overlap. Max Born and Joseph Mayer modelled it in 1932 with an exponential, \(Be^{-r/\rho}\) per closest-neighbour pair, where the range \(\rho \approx 0.03\) nm is fitted to the measured compressibility. Minimizing the total energy, \(U(r) = -\alpha e^2/(4\pi\varepsilon_0 r) + zBe^{-r/\rho}\), and eliminating \(B\) with the equilibrium condition gives

The repulsion reduces the binding by only about one-ninth, because it is so short-ranged. For NaCl, with \(r_0 = 0.282\) nm, the program gives \(-7.91\) eV per ion pair, or 763 kJ per mole, compared with 787 kJ per mole from the thermochemical Born–Haber cycle. The difference of about 3 parts in 100 comes mostly from the van der Waals attraction between the ions and from zero-point energy, both omitted. That a model with one fitted parameter reproduces a lattice energy this well is a strong confirmation of the ionic picture.

Elastic constants and sound

Real crystals respond to stress differently in different directions. In a cubic crystal, Hooke's law involves three independent elastic constants, conventionally \(C_{11}\), \(C_{12}\) and \(C_{44}\), relating stresses to strains. The bulk modulus is \((C_{11} + 2C_{12})/3\); \(C_{44}\) measures resistance to shear. For copper, \(C_{11} = 168\), \(C_{12} = 121\) and \(C_{44} = 75\) GPa, giving \(B = 137\) GPa. Measuring the speeds of sound along different crystal directions determines all three: a longitudinal wave along \([100]\) travels at \(\sqrt{C_{11}/\rho}\), and a transverse wave at \(\sqrt{C_{44}/\rho}\). These long-wavelength sound waves are the low-frequency limit of the lattice vibrations of lesson 5.

A useful comparison across materials is the relation between stiffness and bond strength. Diamond, with very strong covalent bonds and light atoms, has a bulk modulus of 443 GPa and a longitudinal sound speed of about 18 km/s. Sodium, with weak metallic bonding, has a bulk modulus of only about 7 GPa. Solid argon, at 2.7 GPa, is softer than many polymers. As a rough rule, \(B\) scales as the cohesive energy divided by the atomic volume, which is why compounds of light, strongly bonded atoms such as boron, carbon and nitrogen are the stiffest known materials.