Tersoff potential

How can we accurately predict the behavior of materials from the atom up? The answer lies in describing the forces between atoms, but this is far from simple. While straightforward "pair potentials" work for some simple substances, they fundamentally fail to describe the directional, shared nature of covalent bonds that form the backbone of materials like silicon and carbon. This failure is starkly revealed by their inability to predict experimental elastic constants, a discrepancy known as the Cauchy relation violation. This gap highlights the need for a more sophisticated model that understands that the strength of a bond is not fixed but depends on its local environment.

This article explores one of the most elegant and powerful solutions to this problem: the Tersoff potential. Across the following sections, you will gain a comprehensive understanding of this landmark model. First, we will dissect the ​​Principles and Mechanisms​​, revealing how the ingenious concept of "bond order" allows the potential to mimic quantum-mechanical effects in a computationally efficient form. Following that, we will journey through its diverse ​​Applications and Interdisciplinary Connections​​, seeing how the Tersoff potential serves as a computational microscope to investigate everything from crystal defects and surface reconstructions to material fracture and melting.

Let’s start our journey with a simple, beautiful idea. How do atoms hold together to form the materials we see and touch? You might imagine that atoms are like tiny, sticky balls, and the energy of the whole system is just the sum of the stickiness between every pair of them. If you have three atoms, A, B, and C, the total energy would be the energy of A-B, plus B-C, plus A-C. This is the essence of a ​​pairwise-additive potential​​. It’s wonderfully simple. The interaction between any two atoms depends only on the distance between them, $r_{ij}$, and nothing else.

This idea works remarkably well for some materials, like the noble gases. Argon atoms, for instance, are quite happy to be described this way. But what about the materials that form the backbone of our technology and our world? The silicon in a computer chip, the carbon in a diamond, the silica in a pane of glass? Here, the simple picture begins to crumble.

There is a beautiful and devastating piece of evidence against the simple pairwise idea that comes from gently squeezing a crystal. The way a crystal resists being squished or sheared is described by its ​​elastic constants​​. In a cubic crystal, like silicon or diamond, two of these constants are called $C_{12}$ and $C_{44}$. Now, it is a mathematical certainty—a direct consequence of the pairwise-additive assumption—that for any such model, these two constants must be equal: $C_{12} = C_{44}$. This is known as the ​​Cauchy relation​​.

So, we ask nature. What are the values for silicon? Experiment tells us that $C_{12}$ is about $64$ gigapascals, while $C_{44}$ is about $80$ gigapascals. They are not equal! This isn't a small experimental error; it’s a fundamental disagreement. The simple, beautiful idea is wrong. The energy of a covalent solid cannot just be a sum of pairwise distances. There must be something else at play. The forces between atoms must care not just about how far apart they are, but also about their orientation. In other words, the energy must depend on ​​bond angles​​. This is the classic signature of the covalent bond.

Why do angles matter so much? The answer lies in the quantum-mechanical nature of the covalent bond. A covalent bond is formed when atoms share electrons. But an atom doesn't have an infinite supply of electrons to give away. It has a limited "budget" of valence electrons it can use for bonding.

Think of it like a handshake. You can give one person a very firm, strong handshake. If you try to shake hands with two people at once, each handshake will be a bit weaker. If you try to shake hands with four people at once, your attention is divided, and each individual handshake is weaker still. Your total "bonding capacity" is being shared.

This is precisely the concept of ​​bond order​​. The more neighbors an atom is bonded to (its ​​coordination number​​), the weaker each individual bond becomes. A carbon atom in diamond has four neighbors, a classic $sp^3$ hybridization. In graphite, it has three neighbors ($sp^2$ hybridization). The bonds in graphite are stronger and shorter than the bonds in diamond. A model that aims to describe carbon in its many forms must capture this effect: the strength of a bond must depend on its environment.

This is the challenge. We need to build a potential where the interaction energy between two atoms, say atom $i$ and atom $j$, is not a fixed quantity but is modulated by the presence of all the other neighboring atoms. This is the definition of a ​​many-body potential​​.

How can we build such a model? One approach is to start with a simple pair potential and then explicitly add three-body terms that depend on the angles between triplets of atoms. This is the strategy of the Stillinger-Weber potential, and it works quite well for things like crystalline silicon.

But in the 1980s, Jerry Tersoff came up with a different, wonderfully subtle idea. He decided to keep the mathematical form looking like a sum over pairs, but he made the strength of each "pair" interaction dependent on the local environment.

The energy of a bond between atoms $i$ and $j$ in the Tersoff model is given by a famous expression:

Let's break this down.

• $f_c(r_{ij})$ is a ​​cutoff function​​. In a computer simulation, we can't calculate interactions between every atom and every other atom in the universe. This function ensures that the interaction smoothly goes to zero at some cutoff distance. For the simulation to be stable, the function must be "smooth" enough that the forces (which are the derivatives of the energy) don't suddenly jump, which would be like giving an atom a sudden, unphysical kick.

• $f_R(r_{ij})$ is the ​​repulsive part​​. This is simple: when atoms get too close, their electron clouds overlap and they repel each other strongly. This term, usually an exponential function like $A \exp(-\lambda_1 r_{ij})$, takes care of that.

• $f_A(r_{ij})$ is the ​​attractive part​​. This term, something like $-B \exp(-\lambda_2 r_{ij})$, represents the fundamental tendency of atoms to bond and lower their energy.

• $b_{ij}$ is the star of the show. This is the ​​bond-order parameter​​. It's a number, typically between 0 and 1, that acts like a dimmer switch on the attractive part of the potential. If $b_{ij}=1$, the bond is at full strength. If $b_{ij}$ approaches zero, the attraction is turned off. The magic of the Tersoff potential lies entirely in how this little $b_{ij}$ is calculated.

The bond-order parameter $b_{ij}$ is what makes the potential "aware" of its surroundings. It is calculated based on a quantity, $\zeta_{ij}$ (zeta), which you can think of as a measure of the "crowding" or "competition" from other neighbors around the $i-j$ bond.

The formula for the bond order itself is:

You don't have to memorize this. Just notice the crucial relationship: as the crowding term $\zeta_{ij}$ gets bigger, the denominator gets bigger, and so the bond order $b_{ij}$ gets smaller. This is our handshake analogy in mathematical form! More neighbors means more crowding, which means a weaker bond.

So how is the crowding, $\zeta_{ij}$, calculated? It's a sum over all the other neighbors, $k$, of atom $i$:

This sum has two incredibly clever features. First, it simply counts the neighbors $k$. The more terms in the sum, the larger $\zeta_{ij}$ becomes, and the weaker the $i-j$ bond gets. This correctly captures the basic coordination dependence.

Second, and most brilliantly, each neighbor's contribution is weighted by an angular function, $g(\theta_{ijk})$, which depends on the bond angle $\theta_{ijk}$ formed by atoms $j$, $i$, and $k$. This function is designed to have its minimum value (usually 1) when the angle is at the ideal, preferred angle for the material (like $109.47^\circ$ for silicon's tetrahedral bonds). For any other "wrong" angle, the value of $g(\theta_{ijk})$ is larger.

The effect is profound. A neighbor that is at a geometrically "unfavorable" angle contributes more to the crowding factor $\zeta_{ij}$ than a neighbor at an "ideal" angle. This means that forming bonds with incorrect geometry is penalized by making all the bonds involved weaker! The potential has a built-in preference for certain structures, not because it has a rigid angular rule, but because those structures maximize the overall bond strength.

This whole construction—the dimmer switch, the crowding factor, the angular function—is a beautiful piece of physical intuition built into a mathematical form. But is it just a clever trick? Or does it connect to something deeper?

It turns out it does. In a more fundamental quantum mechanical description of bonding, like the ​​tight-binding model​​, the ability of an electron to move between atom $i$ and atom $j$ is described by a "hopping integral," $H_{ij}$. The strength of the resulting bond is related to the magnitude of this integral. Crucially, the value of $H_{ij}$ is not a constant; it is also reduced by the presence of other atoms competing for the same electrons.

The Tersoff bond-order parameter $b_{ij}$ can be seen as a brilliant and computationally simple approximation for how this quantum-mechanical hopping integral is scaled by the local atomic environment. It is a classical potential that has learned a deep lesson from quantum mechanics: bonding is a collective, many-body affair.

The genius of the Tersoff potential is that it packs all of this complex physics—bond order, angular dependence, coordination effects—into a form that is still computationally manageable. Because its formulation is based on the general principles of bond formation and breakage, it is remarkably ​​transferable​​. A single set of parameters for carbon can describe the bonding not only in diamond, but also in graphite, at surfaces, around defects, and in disordered amorphous structures. This is a feat that simpler models, which are often tuned to one specific environment, cannot easily match.

Of course, no model is perfect. The Tersoff potential is designed for covalent systems and is less suitable for the delocalized bonding in metals, where models like the Embedded Atom Method (EAM) are more appropriate. Even more sophisticated (and computationally expensive) models like the Modified Embedded Atom Method (MEAM) have since been developed to try and capture an even wider range of materials with a single framework.

Nevertheless, the Tersoff potential stands as a landmark in our quest to model the material world. It shows us how a few elegant mathematical ideas, guided by profound physical intuition, can capture the complex and beautiful dance of atoms that gives rise to the world around us. It is a testament to the power of finding the right level of description—not too simple that it's wrong, and not too complex that it's unusable—to unlock the secrets of matter.

Now that we have explored the inner workings of the Tersoff potential—its elegant formulation that captures the very essence of covalent bonding—we can embark on a more exciting journey. We will see how this remarkable tool is not merely an abstract mathematical exercise, but a powerful computational microscope that allows us to witness and understand the behavior of materials across a vast landscape of science and engineering. From the subtle quantum whispers that dictate a crystal's properties to the catastrophic roar of fracture, the Tersoff potential serves as our guide.

Our first stop is the world of the crystal in its most ideal form. A central task for any interatomic potential is to correctly describe the ground state of a material. The Tersoff potential excels at this, accurately predicting the equilibrium lattice structure, bond lengths, and elastic properties of materials like silicon, carbon, and germanium. But its true genius lies in what it tells us about the nature of the forces holding the crystal together. This is a task that simpler models, like pairwise potentials, fundamentally fail at. Before a potential can be trusted for large-scale simulations, it must first prove its mettle by matching its predictions for these fundamental properties against the rigorous results of quantum mechanics, specifically Density Functional Theory (DFT). This validation process is the bedrock of modern computational materials science, ensuring our classical models are firmly anchored in physical reality.

One of the most profound insights comes from a subtle discrepancy in a crystal's elastic constants. If you imagine a crystal made of atoms connected by simple central springs (a pair potential), there is a beautiful symmetry that emerges in its elastic response, known as the Cauchy relation. For a cubic crystal, this relation predicts that two of its fundamental elastic constants, $C_{12}$ and $C_{44}$, must be equal. However, for a real covalent material like silicon, experiments definitively show that $C_{12} \neq C_{44}$. What does this violation tell us? It is a macroscopic echo of a microscopic truth: the forces are not simple springs! There is an additional stiffness against bending the bonds, a resistance to changing the angles between them. The Tersoff potential, with its explicit angular dependence, naturally breaks the Cauchy symmetry and correctly captures this essential feature of covalent bonding. The inequality $C_{12} \neq C_{44}$ is the crystal's way of telling us that its bonds have direction and character, a story the Tersoff potential was designed to translate.

Of course, no real crystal is perfect. The presence of defects, such as a missing atom—a vacancy—can dramatically alter a material's properties. These tiny imperfections are the gatekeepers of atomic motion, controlling everything from diffusion in a solid to the performance of a semiconductor device. Using the Tersoff potential, we can perform a computational experiment: we build a large, perfect crystal in our computer, then pluck one atom out and carefully measure the energy cost of this act. This quantity, the vacancy formation energy, is a cornerstone of materials science, and the ability to calculate it accurately is a crucial application of the potential. It is the first step toward understanding the rich and complex world of defects that governs the behavior of the materials we build our world with.

What happens when we slice a crystal, creating a surface? The atoms at this new frontier are suddenly in a highly unnatural state. In the bulk, each atom was content, its bonding valency satisfied by a full complement of neighbors. At the surface, these atoms are left with "dangling bonds"—unsatisfied, high-energy connections reaching out into the vacuum. The system is unhappy, and like any physical system, it will seek a state of lower energy.

This drive leads to one of the most beautiful phenomena in surface science: reconstruction. The surface atoms refuse to stay in their simple, truncated positions. They rearrange themselves, moving, twisting, and forming new bonds with each other to heal the wounds of their creation. A famous example is the (100) surface of silicon, where pairs of atoms move towards each other to form "dimers," satisfying their dangling bonds and reducing the surface energy.

The Tersoff potential captures this drama beautifully. Its core logic dictates that as an atom's coordination number decreases—as happens at a surface—its remaining bonds become stronger and shorter, because the bond order term $b_{ij}$ increases. This alone drives a relaxation. But the angular term provides the crucial choreography for the reconstruction. The atoms will not just move randomly; they will shift and rotate to achieve bond angles closer to the ideal tetrahedral angle that the potential favors. The final reconstructed surface—a complex pattern of dimers and terraces—is the equilibrium solution to a competition between bond stretching, bond bending, and the desire to eliminate dangling bonds. The potential allows us to watch this intricate atomic dance unfold and predict the final, stable configuration.

Let us now turn from creating surfaces to pulling a material apart. The Tersoff potential provides profound insights into the mechanical integrity and failure of covalent solids, connecting the atomic-scale bond-breaking process to the macroscopic field of fracture mechanics.

Before a material breaks, it often sends out warning signals. If you slowly stretch a block of silicon in a simulation, a fascinating process unfolds. Of course, the bonds stretch, which stores elastic energy. But something more subtle is happening. The strain distorts the local geometry, pushing bond angles away from their ideal tetrahedral configuration. According to the rules of the Tersoff potential, this distortion increases the local environment measure $\zeta_{ij}$, which in turn decreases the bond order $b_{ij}$. The bonds are actively getting weaker as the material is stretched.

This bond weakening leads to a phenomenon known as softening. The material's stiffness, or its resistance to further stretching, begins to decrease. We can literally "hear" this in the simulation by listening to the crystal's vibrational frequencies (its phonon modes). Like a guitar string whose pitch drops as it is overstretched, the frequencies of certain phonon modes decrease as the strain increases. The softening of a particular phonon mode to zero frequency signals a lattice instability—the point at which the crystal can no longer resist the load and is about to fail catastrophically. The Tersoff potential allows us to identify these microscopic precursors to fracture.

When the material does fail, the nature of the potential dictates the character of the fracture. The Tersoff potential describes strong, directional, and stiff covalent bonds. This results in a material that can withstand enormous stress—it has a high cohesive strength. However, once this stress is exceeded, the bonds break abruptly, and a crack propagates rapidly. This is the signature of brittle fracture. By comparing the predictions of a Tersoff potential with, for example, an Embedded Atom Method (EAM) potential designed for metals (which have non-directional, delocalized bonding), we can understand why silicon shatters like glass while copper bends and deforms. The Tersoff model, by correctly capturing the physics of covalent bonding, predicts a much smaller "process zone" of nonlinear deformation at the crack tip, leading to behavior that is very close to the ideal Griffith theory of brittle fracture.

Finally, we can use the Tersoff potential to explore what happens when we push materials to their absolute limits with temperature and high-energy impacts. By coupling the potential to the laws of Newtonian motion in a Molecular Dynamics (MD) simulation, we can effectively turn up the heat and watch a material melt.

As we increase the temperature in our simulation, the atoms vibrate more and more violently around their lattice sites. We can track the structure using the Radial Distribution Function, $g(r)$, which shows sharp peaks corresponding to the ordered shells of neighbors in the crystal. At a critical temperature, a dramatic transition occurs. The long-range order vanishes; the sharp peaks in the $g(r)$ blur out into broad humps characteristic of a liquid's short-range order. Simultaneously, we can track the atoms' motion. In the solid, their movement is caged, but in the liquid, they begin to wander freely through the simulation box—the onset of diffusion. The Tersoff potential enables us to simulate this fundamental phase transition and predict properties like the melting point and the structure of the liquid state.

Another extreme environment is found deep inside a semiconductor fabrication plant. To create modern computer chips, silicon wafers are bombarded with high-energy ions in a process called ion implantation. Understanding the damage created by these atomic-scale cannonballs is critical. What happens when a fast-moving ion slams into a silicon atom, knocking it out of its place? The energy required to do this is the threshold displacement energy, $E_d$. The Tersoff potential can be used to model this event. However, the potential was originally parameterized to describe bonding and low-energy events, not the violent, close-range collisions of high-energy particles. At very short distances, the repulsive force is better described by a different model, the Ziegler–Biersack–Littmark (ZBL) potential.

Here, we see the pragmatism and ingenuity of computational physics. Instead of searching for a single, perfect potential, scientists can create a hybrid model. They "stitch" the Tersoff and ZBL potentials together, using a smooth switching function that transitions from the ZBL description at very short distances to the Tersoff description at normal bonding distances. This hybrid approach yields much more accurate predictions for high-energy events like ion implantation, showing how different physical models can be combined to build a more complete and powerful predictive tool.

From the quiet equilibrium of a perfect crystal to the violent chaos of a particle collision, the Tersoff potential and its descendants have proven to be an indispensable tool. It is more than a set of equations; it is a lens that reveals the intricate and beautiful atomic dance that underpins the properties of the material world.