Stillinger-Weber Potential for Silicon

Paper Information

Citation: Stillinger, F. H., & Weber, T. A. (1985). Computer simulation of local order in condensed phases of silicon. Physical Review B, 31(8), 5262-5271. https://doi.org/10.1103/PhysRevB.31.5262

Publication: Physical Review B, 1985

What kind of paper is this?

This is a Method paper.

Its primary contribution is the formulation of the Stillinger-Weber potential, a non-additive potential energy function designed to model tetrahedral semiconductors. While it uses simulation to explore physical properties (Discovery), the enduring legacy is the potential function itself that allowed computational chemistry to move beyond simple pair potentials for covalent materials.

The paper can be thought of as approximately 70% Method and 30% Discovery - the methodological contribution (the potential architecture) is what enabled subsequent research.

What is the motivation?

The authors aimed to simulate the melting and liquid properties of tetrahedral semiconductors (Silicon and Germanium).

• The Problem: Standard pair potentials (like Lennard-Jones) favor close-packed structures (12 nearest neighbors) and cannot stabilize the open diamond structure (4 nearest neighbors) of Silicon.
• The Gap: No existing classical potential could capture the “profound structural change” where Silicon shrinks upon melting (coordination number increases from 4 to >6) while remaining conductive.
• The Goal: To construct a potential that spans the entire configuration space - describing both the rigid crystal and the diffusive liquid - without requiring quantum mechanical calculations.

What is the novelty here?

The core novelty is the introduction of a stabilizing three-body interaction term ($v_3$) to the potential energy function.

• 3-Body Term: Explicitly penalizes deviations from the ideal tetrahedral angle ($\cos \theta_t = -1/3$).
• Unified Model: Unlike “Keating” potentials which only handle small elastic deformations, this potential handles bond breaking and reforming, allowing for the simulation of melting and liquid diffusion.
• Mapping Technique: The application of “steepest-descent mapping” to quench dynamical configurations into their underlying “inherent structures” (local minima), revealing the fundamental topology of the liquid energy landscape.

What experiments were performed?

The authors performed Molecular Dynamics (MD) simulations using the proposed potential.

• System: 216 Silicon atoms in a cubic cell with periodic boundary conditions.
• State Points: Fixed density $\rho = 2.53 \text{ g/cm}^3$ (matching experimental liquid density at melting).
• Process:
  1. Start with diamond crystal at low temperature.
  2. Systematically heat to induce spontaneous nucleation and melting.
  3. Equilibrate the liquid.
  4. Periodically map configurations to potential minima (inherent structures) using steepest descent.

What outcomes/conclusions?

• Validation: The potential successfully stabilizes the diamond structure as the global minimum at zero pressure.
• Liquid Structure: The simulated liquid pair-correlation function $g(r)$ and structure factor $S(k)$ qualitatively match experimental diffraction data, including the characteristic shoulder on the structure factor peak.
• Inherent Structure: The liquid possesses a temperature-independent inherent structure (amorphous network) hidden beneath thermal vibrations.
• Melting/Freezing Criteria: The study proposes an “Inverse Lindemann Criterion”: while crystals melt when vibration amplitude exceeds ~0.19 lattice spacings, liquids freeze when atom displacements from their inherent minima drop below ~0.30 neighbor spacings.

Reproducibility Details

Algorithms

• Integration: Equations of motion integrated using a fifth-order Gear algorithm.
• Time Step: $\Delta t = 5 \times 10^{-3} \tau$ (approx $3.83 \times 10^{-16}$ s).
• Minimization: Steepest-descent mapping utilized Newton’s method to find limiting solutions ($\nabla \Phi = 0$).

Models

To reproduce this work, one must implement the potential $\Phi = \sum v_2 + \sum v_3$ with the exact functional forms and parameters provided.

Reduced Units

• $\sigma = 0.20951 \text{ nm}$
• $\epsilon = 50 \text{ kcal/mol} = 3.4723 \times 10^{-12} \text{ erg}$

Two-Body Term ($v_2$)

$$ v_2(r_{ij}) = \epsilon A (B r_{ij}^{-p} - r_{ij}^{-q}) \exp[(r_{ij} - a)^{-1}] \quad \text{for } r_{ij} < a $$

(Vanishes for $r \geq a$)

Three-Body Term ($v_3$)

$$ v_3(r_i, r_j, r_k) = \epsilon [h(r_{ij}, r_{ik}, \theta_{jik}) + h(r_{ji}, r_{jk}, \theta_{ijk}) + h(r_{ki}, r_{kj}, \theta_{ikj})] $$

where:

$$ h(r_{ij}, r_{ik}, \theta_{jik}) = \lambda \exp[\gamma(r_{ij}-a)^{-1} + \gamma(r_{ik}-a)^{-1}] (\cos\theta_{jik} + \frac{1}{3})^2 $$

(Vanishes if distances $\geq a$)

Parameters

Evaluation

The paper evaluates the model against experimental diffraction data.

Citation

@article{stillingerComputerSimulationLocal1985,
  title = {Computer Simulation of Local Order in Condensed Phases of Silicon},
  author = {Stillinger, Frank H. and Weber, Thomas A.},
  year = 1985,
  month = apr,
  journal = {Physical Review B},
  volume = {31},
  number = {8},
  pages = {5262--5271},
  publisher = {American Physical Society},
  doi = {10.1103/PhysRevB.31.5262}
}