Skip to content

5.1 Equation of state and bulk modulus

These are unexecuted teaching inputs and starting settings to test. Original figures are schematics, not computed results. Use licensed VASP and PAW data, replace every placeholder, record the executable version and validate convergence.

5.1.1 Model, units and provenance

Use eV for energy, Å for length and eV/Å for force; 1 kbar = 0.1 GPa. State normalization per atom, molecule, primitive cell or simulation cell. Record PAW identifiers, release, ZVAL, ENMAX and permitted hashes; never redistribute POTCAR. SCF convergence addresses the chosen electronic problem; convergence of the target property requires separate tests.

Equation of state and bulk modulus — original conceptual schematic

Original schematic. Curves explain concepts; blank data areas await verified learner results. No calculation is claimed.

\[ P=-\frac{dE}{dV},\qquad B=V\frac{d^2E}{dV^2},\qquad x=(V_0/V)^{2/3}\\ E(V)=E_0+\frac{9V_0B_0}{16}\left[B_0\prime(x-1)^3+(x-1)^2(6-4x)\right] \]

5.1.2 Worked case: procedure, interpretation and checks

Question and intuition

How difficult is it to compress a crystal uniformly? Imagine a bowl whose horizontal coordinate is volume and whose height is energy. Its minimum gives the equilibrium volume; its curvature tells you how rapidly the energetic penalty grows when the crystal is squeezed or expanded. The bulk modulus is not the depth of the bowl and not the cohesive energy. Two materials can bind strongly yet respond differently to a small volume change.

Define pressure P(V) = −dE/dV and bulk modulus B(V) = −V dP/dV = V d²E/dV². At the zero-pressure minimum V₀, B₀ is positive for a locally stable crystal. B₀′ = (dB/dP) at P = 0 is dimensionless. Static DFT yields a zero-temperature electronic equation of state; zero-point motion and thermal expansion require later lessons.

Prerequisites and workflow

Start with a relaxed single phase, preferably cubic Si for a first exercise. Decide whether the question concerns isotropic scaling at fixed shape or a fully shape-relaxed energy at each volume. They coincide for a stable cubic phase but need not coincide for a low-symmetry crystal.

  1. Obtain a well-relaxed reference volume Vref and retain the same atom count and species ordering throughout.
  2. Choose seven to eleven trial volume ratios spanning both sides of the expected minimum, for example 0.94 through 1.06. These are sampling choices, not a recommendation for every pressure range.
  3. Scale every lattice vector by s = (V/Vref)^(1/3), keeping fractional coordinates as the starting geometry. A 3% volume change is approximately a 1% lattice change, not a 3% change of every vector.
  4. At each fixed trial cell relax internal atomic coordinates using ISIF=2. For shape-relaxed EOS, explicitly implement a volume-constrained shape procedure and document it; unconstrained ISIF=3 would erase the intended volume series.
  5. Compute a consistent high-accuracy static energy for each relaxed structure. Keep cutoff fixed across volumes and initially use the same electronic mesh topology to avoid artificial jumps.
  6. Fit E(V), inspect residuals, and repeat with a narrower fitting interval or a different EOS. Report V₀, B₀, B₀′, cell normalization, fit model, and sensitivity, not just one best-fit number.

Input and analysis

An original worksheet can store ratio, volume_A3, energy_eV, max_force_eVA, pressure_kbar, scf_ok. No fitted value should be produced from unconverged rows. For a third-order Birch–Murnaghan fit, let x = (V₀/V)^(2/3):

E(V) = E₀ + (9V₀B₀/16){B₀′(x−1)³ + (x−1)²(6−4x)}.

Use B₀ in eV/ų inside this equation, then convert to GPa. Plot E−min(E), not enormous absolute energies that hide millielectronvolt structure. Inspect pressure from the derivative as an independent check against the computed hydrostatic pressure. Fitting an arbitrary high-degree polynomial can produce an attractive curve and unphysical curvature.

Checks, failure modes, exercise

Increase cutoff and k mesh until B₀ itself stabilizes. Basis-set Pulay stress can shift the apparent equilibrium. Points must bracket a minimum; do not extrapolate B₀ from only compressed structures. A spin transition or phase change creates branches and should not be hidden inside one smooth EOS fit. B₀′ is typically more sensitive than B₀, so do not report many unsupported digits. A residual plot can expose one failed SCF point or an inconsistent energy definition.

Exercise: fit the same converged series with all points and with the outermost pair removed. Which parameter changes most? Then compare B₀ from EOS with (C11+2C12)/3 for a cubic material from 5.2 Elastic tensor from finite strains. Explain why disagreement beyond numerical uncertainty can signal mismatched relaxation conditions.

5.1.3 Unexecuted inputs and analysis scaffolds

These are unexecuted teaching inputs and starting settings to test. Original figures are schematics, not computed results. Use licensed VASP and PAW data, replace every placeholder, record the executable version and validate convergence.

5.1.3.1 Input block 1

# Fixed-cell internal relaxation; merge with the common/material block
IBRION = 2
NSW    = 120
ISIF   = 2
EDIFFG = -0.003

# For the final static calculation, replace the four lines above by:
# IBRION = -1
# NSW = 0
# ISIF = 2

5.1.5 Technical sources