Example 04: Automated calculation of melting temperature#

In Example 03, we calculated the melting temperature of Cu manually by running a solid and a liquid reversible-scaling sweep and looking for the free-energy crossing. Calphy can do all of this automatically with mode: melting_temperature, including iterating until the search converges.

The EAM potential we use is Mishin et al., Phys. Rev. B 63, 224106 (2001). Its melting temperature is roughly 1270–1290 K (the experimental value for Cu is 1357 K).

Input#

mode: melting_temperature needs very little structural input — no explicit temperature (an initial guess is taken from the element’s experimental melting point via mendeleev, or you can set it with melting_temperature.guess), and the search drives the lattice via reference_phase internally. The full input used in this example:

calculations:
- element: Cu
  mass: 63.546
  repeat:
  - 10
  - 10
  - 10
  md:
    timestep: 0.001
  mode: melting_temperature
  n_equilibration_steps: 10000
  n_iterations: 1
  n_print_steps: 100
  n_switching_steps: 25000
  pair_coeff: '* * ../potentials/Cu01.eam.alloy Cu'
  pair_style: eam/alloy
  equilibration_control: berendsen
  queue:
    commands:
    - conda activate calphy
    cores: 4
    scheduler: local
  phase_transition_detection:
    mode: recover
  melting_temperature:
    step: 400

Two blocks are particularly relevant:

  • phase_transition_detection: { mode: recover } lets the solid and liquid sub-sweeps safely truncate themselves if either phase loses stability mid-sweep, rather than producing contaminated free-energy curves. At 4000 atoms (repeat: [10,10,10] fcc) the default peak_threshold: 12 is appropriate; smaller cells should lower it (see Phase-transition detection).

  • melting_temperature.step: 400 sets a ±400 K initial window around the guess (1357 K for Cu via mendeleev → window [958, 1758] K). A larger half-width is helpful when MD hysteresis (superheating / supercooling) is ~100–200 K at this system size — both phases need enough room on each side of \(T_m\) to sample cleanly.

Run the calculation as usual:

calphy -i input.yaml

How the search works#

Each attempt runs one solid sweep (\(T_\mathrm{min} \to T_\mathrm{max}\)) and one liquid sweep (\(T_\mathrm{max} \to T_\mathrm{min}\)). Both sweeps produce free-energy curves that may be truncated by phase_transition_detection: recover at the boundary where the corresponding phase loses stability. The truncation points give a bracket on \(T_m\) directly:

  • The solid sweep’s last clean temperature T_sol_max is an upper bound (solid stayed stable up to here).

  • The liquid sweep’s last clean temperature T_liq_min is a lower bound (liquid stayed stable down to here).

There are three outcomes per attempt:

  1. Overlap with crossing — the two FE curves overlap in \([T_\mathrm{liq\_min}, T_\mathrm{sol\_max}]\) and \(G_\mathrm{solid}(T) - G_\mathrm{liquid}(T)\) changes sign there. The crossing is \(T_m\) → done.

  2. Overlap, no crossing — both phases coexist metastably across the overlap; one is uniformly lower-FE. Each curve is fit with the 6-term CALPHAD Gibbs polynomial \(G(T) = a + b\,T + c\,T\ln T + d\,T^2 + e\,T^3 + f/T\), and \(G_\mathrm{sol} = G_\mathrm{liq}\) is solved for the extrapolated \(T_m\). The window is recentred on that prediction and another attempt is run. Two consecutive CALPHAD predictions agreeing within 25 K stop the search.

  3. No overlap (gap) — \(T_\mathrm{sol\_max} < T_\mathrm{liq\_min}\). \(T_m\) is bracketed inside the gap. Recentre on the gap midpoint and retry.

A fourth outcome — failure of the averaging stage at the window endpoints (the streaming monitor flipped the phase before any sweep data could be collected) — raises a clear error advising you to set melting_temperature.guess closer to the expected \(T_m\).

Each attempt also reports its energy dissipation; an attempt where either phase silently sampled a hidden transition is flagged as contaminated and skipped (after two consecutive contaminated attempts the search aborts with mitigation advice).

Output#

Two files of interest end up next to the input:

  • ``<run>.log`` — full rotating log. grep STATE <run>.log distils the key milestones:

    STATE: Temperature range of 957.770000-1757.770000 K
STATE: Tm = 1285.24 K +/- 0.24 K

  • ``<run>.report.yaml`` — machine-readable result. Contains the converged \(T_m\) and its uncertainty, the experimental reference value, and the full per-attempt history with windows, brackets (\(T_\mathrm{sol\_max}\), \(T_\mathrm{liq\_min}\)), case decisions, and per-phase energy dissipations:

    status: converged_extrapolation
tm: 1285.24
tmerr: 0.24
tm_experimental: 1357.77
attempts: 2
max_attempts: 5
step: 400.0
history:
- attempt: 0
  tmin: 957.77
  tmax: 1757.77
  ediss_solid: 1.3e-04
  ediss_liquid: 2.3e-04
  T_sol_max: 1207.85
  T_liq_min: 958.00
  case: "2: overlap, no crossing → CALPHAD extrapolation"
  tpred: 1285.47
- attempt: 1
  ...

For this Cu EAM, the search converged on \(T_m = 1285.24 \pm 0.24\) K in two attempts (~6 minutes wall-clock on 4 cores), in line with the published value for this potential.

FAQs#

  • How can I tune the initial guess temperature?
    Set melting_temperature.guess (in Kelvin) inside the calculation block. If omitted, calphy uses the experimental melting point of element[0] from mendeleev.
  • How can I tune the width of the temperature window?
    Set melting_temperature.step; this is the half-width of the initial (and recentred) windows. For systems with significant MD hysteresis (small cells, slow nucleation potentials) prefer ±400 K or more; for well-converged FCC metals at ≥4000 atoms ±200 K is usually enough.
  • How many attempts will it use?
    Capped by melting_temperature.attempts (default 5). Most clean cases converge in 1–2 attempts via case 1 (overlap crossing) or case-2 extrapolation stability.
  • What if the system undergoes a solid-solid phase transition before melting?
    Calphy picks the ground-state lattice automatically. For elements with a solid-solid transition below \(T_m\) (e.g. HCP → BCC → liquid for Ti) you must specify the higher-temperature solid phase explicitly with the lattice keyword so the solid sweep starts in the correct phase.
  • How can I calculate :math:`T_m` at non-zero pressure?
    Add pressure: <bar> to the calculation block. Setting melting_temperature.guess helps too — the experimental \(T_m\) from mendeleev is for ambient pressure.
  • The run aborted with a ``contaminated`` error. What does that mean?
    Energy dissipation on one of the sub-sweeps was too high (> \(10^{-3}\) eV/atom), which means the phase-transition detector silently missed a transition during that sweep. Lower phase_transition_detection.peak_threshold (e.g. 6 for small cells), grow the system size, or raise n_switching_steps.