2 Classical potentials

2.1 Why start here?

Classical potentials are not the most accurate method for calculating E(r), but they are the most transparent. You can write down the function, plot it, and see the potential energy surface directly. This builds intuition before moving to methods where \(E(r)\) comes out of a more complex calculation.

2.2 The simplest model: pair potentials

Classical potentials are parameterised analytical functions that give the energy of a configuration of atoms. A mathematical form is chosen, parameters are fitted to reproduce known properties, and the resulting function is used to calculate energies for any configuration.

The simplest approach is to decompose the total energy into contributions from pairs of atoms:

\[E = \sum_{ij} V(r_{ij})\]

where \(r_{ij}\) is the distance between atoms \(i\) and \(j\), and \(V(r)\) is a pair potential. For each pair of atoms, we look up their separation, evaluate the pair potential, and sum over all pairs. This is conceptually simple and computationally cheap.

2.3 Building a potential from physics

What should \(V(r)\) look like for an ionic crystal? The dominant contribution is electrostatic. Ions carry charges, so opposite charges attract and like charges repel. For two ions with charges \(q_i\) and \(q_j\) separated by distance \(r\), the Coulomb energy is \(q_iq_j/r\). This is long-range — it falls off slowly with distance.

Electrostatics alone would cause the crystal to collapse — opposite charges attracting without limit. What prevents this is short-range repulsion. When ions get close, their electron clouds overlap, and this costs energy. The repulsion is modelled with a steeply rising function at short range. A common choice is the Born-Mayer form: \(A\exp(−r/\rho)\), an exponential that increases rapidly as \(r\) decreases.

There is also a weak, long-range attraction from van der Waals forces, typically modelled as \(−C/r^6\). For hard ions this is often a small correction, but it matters more for polarisable species.

Combining Coulomb, exponential repulsion, and dispersion gives a Buckingham-type potential — the workhorse for simulating ionic solids like oxides and halides.

2.4 Periodic systems

Most simulations aim to model bulk materials, but a computer simulation is necessarily finite — typically thousands to millions of atoms. This creates a problem: a finite cell has edges, and atoms near those edges behave differently from those in the interior. The solution is to apply periodic boundary conditions: the simulation cell tiles infinitely in all directions, so an atom near one boundary sees atoms from the opposite side as nearest neighbours. The result is a system with no surfaces — only bulk.

Periodicity introduces its own complication. Each atom now interacts not just with other atoms in its cell, but with all their periodic images — infinitely many copies stretching in every direction. The total energy becomes an infinite sum.

For short-range interactions, this infinite sum is manageable. The exponential repulsion and dispersion terms decay rapidly with distance, so contributions beyond a nanometre or so are negligible. But which atoms count as “nearby” in a periodic system? Consider atom A interacting with atom B. There is a copy of B in every periodic image — infinitely many, at different distances from A. For a short-range interaction, only the nearest copy of B is close enough to matter. This is the minimum image convention: for each pair of atoms, consider only the nearest periodic image. Combined with a cutoff distance beyond which interactions are ignored entirely, and provided this cutoff is less than half the cell length, the infinite sum reduces to a finite and well-defined calculation.

The Coulomb interaction poses a harder problem. It decays as \(1/r\) — too slowly to truncate. Cutting off the sum at any finite distance gives results that depend on where you cut, which is unphysical; the sum over periodic images is only conditionally convergent. Ewald summation resolves this by splitting the interaction into two parts. Each point charge is screened with a compensating Gaussian distribution. The screened charges interact only at short range and sum quickly in real space. The smooth Gaussian distributions sum quickly in reciprocal space. A correction removes the spurious self-interaction of each charge with its own screening cloud. The result is two absolutely convergent sums that together give the correct electrostatic energy. The reciprocal-space treatment will be familiar from diffraction.

2.5 When pair potentials aren’t enough

Pair potentials assume the energy contribution from two atoms depends only on their separation. This works well for ionic materials where bonding is non-directional. But for covalent or partially covalent materials, where bond angles matter, pair potentials are insufficient.

In a silicate, for instance, the O–Si–O angle isn’t arbitrary — there’s a preferred tetrahedral geometry. Capturing this requires three-body terms that penalise deviations from the preferred angle. For molecular systems or structures with well-defined bonded units, the “force field” approach uses harmonic springs for bond stretches, angle terms for bond angles, and torsional terms for dihedrals. More complex functional forms exist for metals (embedded atom method), semiconductors (Tersoff potentials), and reactive systems (ReaxFF).

The basic ionic model also treats ions as fixed point charges. Real ions are polarisable — their electron clouds respond to local electric fields — and may have non-spherical charge distributions. Extensions address this through explicit induced dipoles, higher-order multipoles, variable ion size and shape, or charge equilibration schemes where charges vary self-consistently subject to a total charge constraint.

Classical potentials can be much more sophisticated than the basic Buckingham form. But there is always a tradeoff: more realism means more parameters, more complexity, and more fitting required.

2.6 The parameterisation problem

A classical potential has parameters — charges, repulsion coefficients, cutoffs. Where do they come from?

Historically, parameters were fitted to reproduce experimental properties: lattice constants, elastic moduli, phonon frequencies, defect energies. The potential is only as good as the data it was fitted to, and only as transferable as that data allows. There is also a risk of overfitting: fitting five parameters to three experimental observables leaves the problem underdetermined. Parameter sets may reproduce those observables perfectly but for the wrong reasons, and fail elsewhere.

Increasingly, potentials are fitted to DFT calculations. A set of configurations is generated, their energies and forces are calculated with DFT, and the potential parameters are fitted to reproduce them. This covers regions of configuration space that experiment doesn’t probe directly, and provides enough data points to constrain all parameters.

Either way, the quality of a classical potential depends entirely on how well the parameterisation captures the relevant physics. And there is a more fundamental issue: the functional form itself bounds the accuracy. The physics that can be captured is limited by the mathematical form chosen. If the real interactions don’t match that form, no amount of parameter fitting will fix it.

2.7 Strengths and limitations

Classical potentials are fast — orders of magnitude cheaper than quantum-mechanical methods. Simulations of millions of atoms for nanoseconds are routine. With methods like particle-mesh Ewald or fast multipole for electrostatics, computational cost scales as O(N log N) with system size. Classical potentials are also interpretable: each term’s contribution is visible, and parameter–property relationships can be understood. When something goes wrong, diagnosis is possible. For some material classes — simple ionic solids, for instance — good potentials exist and have been extensively validated.

The limitations are significant. The treatment of electrostatics is typically simple: fixed point charges throughout the simulation, with no polarisation response or charge transfer. Electronic structure is implicit — electrons aren’t explicitly treated, so processes involving electronic rearrangement cannot be described.

Developing a good classical potential is hard work — designing a functional form, generating training data, fitting, validating, finding problems, and iterating. This requires substantial effort with no guarantee of success. Classical potentials also have limited transferability. A potential fitted to one material, or one region of configuration space, may fail for another. Bulk doesn’t guarantee surfaces; ambient conditions don’t guarantee high pressure. Extrapolation fails silently — the potential gives a number even outside its valid domain.

Classical potentials are the right choice for large system sizes or long timescales, for well-characterised material classes where good potentials exist, or for quick exploratory calculations. They are the wrong choice when the required accuracy exceeds what the potential can provide, or when moving outside its validated domain.