Atomistic models, assumptions and validation

Atoms in a calculation follow a representation

First-principles and empirical atomistic methods use different approximations to electronic interactions and atomic motion. Cell size, boundary conditions, temperature, loading rate and numerical convergence affect what is being predicted. A detailed atomic image does not remove those assumptions.

Two calculations place an isolated pair of model atoms at the same equilibrium separation and give exactly the same stiffness for very small displacements. One predicts twice the other’s largest tensile holding force. There is no contradiction: the equilibrium fit samples only the bottom of an energy well. We will derive the difference, then ask whether smaller numerical errors or more atoms can resolve it.

All parameters and benchmark records below are invented teaching inputs. The curves are actual evaluations of the stated equations; they are not new electronic-structure calculations, molecular-dynamics measurements or discoveries of a candidate material. Energy and the meaning of stability prepares the energy comparison; Elastic response and elastic modulus prepares the distinction between a local slope and a failure response.

Choose the represented physics before solving it

An atomistic configuration specifies species, positions, cell shape and relevant charge or spin state. An interaction model assigns an energy U to that configuration. In a classical trajectory the nuclear coordinates move under forces F_i = −∂U/∂R_i and Newton’s equation m_iR̈_i = F_i. A static relaxation instead searches for a stationary configuration under specified constraints. Neither the word atomistic nor an attractive picture says which of those tasks was done.

In a common first-principles route, an electronic-structure calculation supplies the energy for each nuclear configuration. A ground-state treatment holds the nuclei fixed while solving the electronic problem, then uses the resulting energy surface for relaxation or classical nuclear motion. That separation does not describe explicitly excited-electron dynamics. Kohn–Sham density functional theory represents the electron density through orbitals and uses an exchange-correlation approximation. Its functional, treatment of core electrons, charge and spin state remain model choices. A finite plane-wave basis, reciprocal-space sampling and iterative solution tolerances introduce additional numerical choices. Increasing a basis cutoff does not change the chosen exchange-correlation approximation.

An empirical potential instead evaluates a prescribed energy function with fitted parameters, without solving the electronic problem afresh for each configuration. Modern fitted models can depend on many neighbours or use learned functions; they are not all pair potentials. The configurations and responses used to fit them matter. A model fitted to bulk elastic data may encounter different coordination, charge redistribution or bonding at a free surface, defect or separating interface.

A static potential-energy minimum is also not automatically a finite-temperature free-energy minimum. Classical nuclear motion omits nuclear quantum effects unless they are added by another treatment. Temperature, external pressure and constraints must match the claim. The following pair model deliberately removes those complications so that its assumptions can be inspected rather than hidden.

Differentiate a specified energy and keep the force sign

Let two sites lie on an x axis with separation r. Write δ = r − r₀ and use the untruncated Morse form φ(r) = D[e^(−2aδ) − 2e^(−aδ)], with D > 0 and a > 0. The parameters D, a and r₀ have units of energy, inverse length and length. At δ = 0, φ = −D; as the sites separate without bound, φ tends to zero. Therefore D is the separation energy of this isolated pair within this model.

Define f = dφ/dr = 2Da[e^(−aδ) − e^(−2aδ)]. This f is the positive external holding force needed on the right site during quasistatic tensile separation. The interaction force on that site is −f, directed back towards its partner. On compression f is negative. One electronvolt is 1.602176634 × 10⁻¹⁹ J, one ångström is 10⁻¹⁰ m, and 1 eV/Å = 1.602176634 nN.

The local tangent is d²φ/dr² = 2Da²[2e^(−2aδ) − e^(−aδ)]. At equilibrium it is k = 2Da², with units of force per length. For small δ, φ + D ≈ ½kδ² and f ≈ kδ. This pair stiffness is not Young’s modulus: no three-dimensional specimen, cross-sectional area or lateral response has been defined.

No interaction cut-off or smoothing is applied in this analytic example. A software implementation may truncate the interaction at a finite separation. Shifting energy by a constant can make the energy continuous there without removing a jump in its derivative. A force-smoothing modification changes the interaction and must be declared, especially when the predicted event explores those separations.

The holding force stops increasing when its derivative is zero. Set e^(−aδ) = 1/2 to obtain δ_peak = ln(2)/a and f_peak = Da/2. This is the maximum on one prescribed pair-separation coordinate. It is not an experimentally measured fracture stress, toughness or finite-temperature rupture load. In a larger body, relaxation, defects or another mode can intervene before that coordinate reaches its maximum.

An identical equilibrium fit leaves the separating response undetermined

Give both models r₀ = 2.50 Å. Model A has D = 1.00 eV and a = 2.00 Å⁻¹; model B has D = 4.00 eV and a = 1.00 Å⁻¹. Both have zero force at r₀ and k = 8.00 eV/Ų = 128.174 N/m. A fit that records only r₀ and k cannot distinguish them.

Their wells nevertheless differ away from equilibrium. A needs 1.00 eV to separate its pair and reaches f_peak = 1.00 eV/Å = 1.60218 nN at δ = 0.346574 Å. B needs 4.00 eV and reaches 2.00 eV/Å = 3.20435 nN at δ = 0.693147 Å. Equal local curvature does not determine the depth or width of the well.

At the same withheld displacement δ = 0.250 Å, A predicts 1.52945 nN and B predicts 2.20806 nN. Suppose an independent teaching benchmark specifies 2.20 nN for exactly this constrained pair coordinate and declares agreement to mean an absolute residual no greater than 0.10 nN. The benchmark is hypothetical, and the tolerance is a decision convention, not a confidence interval. A fails this comparison; B passes it. The parameters were fixed before using this benchmark.

This is a useful falsifiable comparison because the models predict different observables under the same condition. One passing point does not establish the entire B curve, its separation energy or any real material’s fracture behaviour. A further displacement or separation-energy comparison can still challenge B. Refitting A to the withheld point would use it for calibration and require a new independent check.

A converged numerical derivative can still converge to the wrong physical answer

To verify a force implementation, compare its reported interaction force with an independent derivative of its own energy. For the positive holding-force convention use f_h = [φ(r + h) − φ(r − h)]/(2h). Taylor expansion gives f_h = φ′(r) + h²φ‴(r)/6 + terms of order h⁴ for a smooth energy. Halving h should approximately quarter the leading truncation error while that term dominates. Too small a difference can instead amplify round-off or noise from incompletely converged electronic energies.

For A at δ = 0.250 Å, the table approaches the analytic 1.52945 nN. The final two estimates differ by only 0.001734 nN, yet the analytic result is about 0.671 nN below the supplied 2.20 nN benchmark. Refining this derivative cannot make A converge to B or to the benchmark. It verifies the solution of A, not the adequacy of A’s interaction.

A change between two numerical settings is a sensitivity observation, not by itself a rigorous bound on all remaining errors. Check the actual predicted observable, not just total energy digits. In an electronic calculation, independently examine basis cutoff, reciprocal-space sampling, iterative electronic convergence and force/structural relaxation tolerances as relevant. A smooth energy curve alone does not demonstrate converged stress or energy differences.

A published 2016 inter-code study compared energy–volume predictions within a fixed density-functional approximation for elemental crystals. Close agreement between implementations established precision for that benchmark. It did not remove the shared approximation or validate every property. The same distinction appears here in a fully transparent algebraic model.

More atoms can change the represented specimen rather than merely improve a calculation

Consider an original one-dimensional, nearest-neighbour chain made from A. Ignore all other interactions, temperature and transverse motion. An open chain of N atoms has N − 1 bonds. At separation r₀ its energy is −(N − 1)D, so energy per atom is −D(1 − 1/N). A one-dimensional periodically repeated cell with N atoms and one bond across its boundary has N bonds per cell and energy per atom −D. Use N ≥ 3 here, so the distinct boundary bond is unambiguous.

For open chains with N = 4, 8 and 16, the energies per atom are −0.7500, −0.8750 and −0.9375 eV. The periodic value is −1.0000 eV. These values are exact for the stipulated model. The difference is the contribution of the free ends per atom, not an integration error. Increasing N reduces the end fraction; replacing free ends with a periodic bond changes the boundary condition immediately.

The same example warns against using a bulk-looking energy as proof of a converged failure load. If every open-chain bond is constrained to the same extension δ, its end-to-end length is L = (N − 1)(r₀ + δ) and U = (N − 1)φ. Differentiation gives dU/dL = dφ/dr = f: the tensile end force is independent of N under that constraint. The energy per atom converges with N while this force does not change. Unequal bond extensions and localisation are excluded from this constrained path.

In a real atomistic study, periodic boundaries repeat the structure and can remove surfaces or make defects interact with their images. Cell size limits which spatial perturbations fit. Holding one defect while increasing the cell also changes its concentration. To interpret a size comparison, preserve the intended material state, imposed strain or stress, defect density and observable, or explicitly identify which of them changed. Neither more atoms nor periodicity is universally the more physical choice.

A stable trajectory step is not an accurate material prediction

Near equilibrium the pair’s small relative displacement x obeys μẍ = −kx, where μ is the reduced mass and ω² = k/μ. A central Verlet update gives x_(n+1) = [2 − (ωΔt)²]x_n − x_(n−1). Let q = ωΔt. Substitution of x_n proportional to e^(inθ) gives cos θ = 1 − q²/2. For 0 < q < 2 the ideal harmonic recurrence has bounded oscillatory solutions; at q = 2 the repeated root can produce growing solutions, and beyond it generic solutions are unstable.

Within that bounded interval, the discrete frequency is ω_num = θ/Δt = 2arcsin(q/2)/Δt. At q = 1, it is 4.7198% above ω, despite the absence of runaway growth. At q = 0.2, the error is 0.1674%. This is an exact small-oscillation numerical example, not a time-step recommendation for any material. Anharmonic motion and other modes require their own checks.

Temperature introduces a distribution of configurations and velocities. A thermostat exchanges energy to control a specified ensemble; it does not validate the force law or prove that rare events have been sampled. Report how equilibration, independent initial states and averaging were checked. In driven deformation, loading rate and heat removal affect the available relaxation and activated events. Simply running a loading trajectory with a smaller time step holds its physical loading rate fixed; changing that rate is a different sensitivity test.

The static, constrained pair maximum derived earlier contains no thermal activation, crack population or rate process. It cannot be relabelled a room-temperature strength or service lifetime. A short trajectory with no rearrangement is finite observation under that protocol, not proof that rearrangement is impossible.

A stationary structure must be tested against the perturbations in the claim

Zero force means a stationary point of the represented energy, not necessarily a minimum. Near one such point, write ΔU ≈ ½uᵀHu, where u collects allowed small displacements and H is the matrix of second derivatives, the Hessian. A negative curvature along any allowed direction lowers energy. Rigid translation produces zero modes in an isolated system; it should not be mistaken for a restoring vibration.

For an invented two-coordinate example, let ΔU = ½(8)x² − ½(2)y² eV when x and y are measured in Å. Both forces vanish at the origin. If a constraint fixes y = 0, the x test finds positive curvature 8 eV/Ų. Releasing y exposes curvature −2 eV/Ų: at x = 0 and y = 0.050 Å the energy falls by 0.00250 eV. The restricted minimum is a saddle in the expanded configuration space.

Atomic coordinates, cell shape, long-wavelength modes, spin state or a chemical rearrangement can provide missing directions in an actual candidate study. Choose the relevant thermodynamic potential and constraints when pressure, temperature or loading is imposed. A positive test over sampled modes remains conditional on the tested cell and model; it is not a proof of the global lowest state, a synthesis route or unlimited persistence. Metastability and persistence develops the distinction between an energy landscape and an observed duration.

Validation must target the relevant response

Compare with appropriate structural, energetic or mechanical observations and examine sensitivity to model choices. Agreement for equilibrium geometry does not automatically validate bond-breaking dynamics. Stability against one perturbation does not establish synthesizability or lifetime under every environment.

State a prediction as an observable for a defined configuration and protocol. For the pair example, retain r₀, the separation coordinate, boundary constraints, zero-temperature energy interpretation and force sign. Keep the independently withheld benchmark separate from the quantities used for fitting. Compare both the response and its uncertainty or declared tolerance, and state which conclusion would change if another deformation were supplied.

For a material claim, choose measurements that interrogate the proposed mechanism: equilibrium structure for geometry, elastic response for local curvature, defect or surface energies for changed coordination, and appropriate separation or fracture observations for failure. Ensure that composition, defects, temperature, direction, geometry and loading conditions correspond, or supply and test a model connecting them. A first-principles result used to fit a potential is not simultaneously an independent experimental validation of that potential.

A reproducible calculation record retains the exact model and parameters, configuration, boundaries, numerical settings, sampling and observable definition. Predictions, targets and levels of evidence explains how such a conditional record differs from phase identification and a measured specimen. Here A can be rejected by a specified pair-force comparison; B remains consistent with that one comparison. Neither record establishes a synthesized, recoverable or usable bulk candidate.

Worked examples

Try the reasoning

1. Why vary simulation-cell size?

   Show the answer and reasoning

   To examine finite-size effects.

   Artificial boundaries can influence results. Preserve or explicitly change defect concentration, surfaces, loading and the observable. The open-chain energy example has exact free-end physics, whereas a periodic cell represents a different boundary problem.

2. Does converged computation guarantee a complete physical model?

   Show the answer and reasoning

   No.

   Numerical solution and representation validity are different. A refined energy derivative can converge to a force that fails an independent comparison. A larger electronic basis also leaves the chosen exchange-correlation approximation in place.

3. At r = r₀ + 0.250 Å in A, the positive holding force is 1.529446 nN. With the left site at x = 0 and the right site on +x, what is the interaction force on the right site? Does φ′ > 0 mean the atoms repel?

   Show the answer and reasoning

   The interaction force on the right site is −1.529446 nN along x. The interaction is attractive at this stretched separation.

   The external holding force is +φ′ while the interaction force is −φ′. Increasing separation raises energy here, so the interaction opposes that increase. Confusing the holding-force convention with the interaction force reverses the interpretation.

4. Within this Morse family keep k = 8 eV/Ų and r₀ fixed, but double a from 1 to 2 Å⁻¹. Derive the changes in D, f_peak and δ_peak without reusing a table entry.

   Show the answer and reasoning

   D falls by a factor four; f_peak and δ_peak each fall by a factor two.

   D = k/(2a²), f_peak = Da/2 = k/(4a), and δ_peak = ln(2)/a. A sharper well with the same local curvature has a smaller separation energy and lower peak holding force in this family. It is not a universal material relation.

5. For constrained open chains of A with N = 4 and 8, find the energy per atom at equilibrium. Then stretch every bond by 0.250 Å. Should their end forces be multiplied by N − 1?

   Show the answer and reasoning

   At equilibrium the energies per atom are −0.7500 and −0.8750 eV. Both stretched chains have end force 1.529446 nN under the stated equal-extension constraint.

   There are N − 1 serial bonds, so U = (N − 1)φ but L = (N − 1)r. Their factors cancel in dU/dL. Adding bond forces as if they were parallel would describe another arrangement. Unequal extension/localisation and transverse modes remain excluded.

6. For the harmonic recurrence at fixed Δt, quadruple the reduced mass while keeping the interaction stiffness fixed. An original q = 1 becomes what q, and what happens to the static pair-force maximum?

   Show the answer and reasoning

   ω falls by a factor two, so q becomes 0.5. The static maximum is unchanged.

   ω = √(k/μ) depends on mass, whereas the static derivative of φ does not. The new frequency ratio is 2arcsin(0.25)/0.5 ≈ 1.010721, a 1.0721% numerical error. This diagnoses the integrator for the stipulated oscillator, not a universal isotope effect on real finite-temperature strength.

7. In the stated two-coordinate energy ΔU = 4x² − y² eV, a relaxation begins at x = y = 0 and reports zero force. If y was fixed, what new conclusion follows on allowing y = 0.050 Å? Does longer observation with y fixed resolve it?

   Show the answer and reasoning

   The energy can decrease by 0.00250 eV, so the origin is unstable in the newly allowed y direction. More time with y fixed does not test that direction.

   The Hessian eigenvalues in these coordinates are +8 and −2 eV/Ų. Zero gradient is compatible with a saddle, and a constraint can suppress a lowering direction. This example supplies no barrier, synthesis route or lifetime.

8. A potential was fitted to a quantum calculation’s equilibrium cell and elastic constants. The report repeats those fits and claims experimental fracture validation. Name a meaningful next comparison and the conditions it must retain.

   Show the answer and reasoning

   Use a withheld response relevant to separation or fracture, with specified structure, defects, temperature, direction and loading conditions, and compare the observable with its appropriate uncertainty or declared tolerance.

   Repeating calibration data checks neither an independent observable nor an experiment. A different electronic calculation may test model transfer, but it remains computational evidence; an experimental comparison must also establish a corresponding specimen and measurement model. Neither result alone demonstrates usable bulk material.

What to carry forward

• Inspect physical and numerical assumptions.
• Validate the particular predicted mechanism.
• Energy derivatives determine forces and local curvature, but the interaction representation and allowed configurations determine what those derivatives mean.
• Matching an equilibrium location and stiffness does not uniquely determine separation energy or peak holding force.
• Refinement can reduce numerical error while leaving physical discrepancy unchanged.
• Cell size, surfaces, defects, temperature, loading rate and sampling define different aspects of the represented problem.
• A conditional pair maximum or restricted minimum does not establish bulk fracture, synthesis or persistence.