[WIP] CALPHAD workflow with ATAT sqs2tdb

[hrushikesh-s]

Summary

Include a summary of major changes in bullet points:

• New workflow that fits a CALPHAD database (a TDB file) for a binary system with a force field. It uses the sqs2tdb tool of ATAT (van de Walle et al., Calphad 58, 70, 2017).
  1. sqs2tdb -cp copies the special quasirandom structures (SQS) of each lattice from the ATAT database.
  2. Each solid SQS is relaxed with a force field. Its harmonic phonons give its vibrational entropy, which sqs2tdb turns into a term linear in T in the solid mixing terms.
  3. Each liquid SQS is melted and then run with isotropic NPT MD at zero pressure. The mean potential energy of the second run is used, with its standard error from five blocks of the run.
  4. sqs2tdb -fit fits each lattice, and sqs2tdb -tdb writes one TDB file. pycalphad can read it to draw the phase diagram. With short_range_order=True, sqs2tdb -fit -sro adds its low-order CVM approximation of the short range order to FCC_A1, BCC_A2, HCP_A3 and DIAMOND_A4. It is off by default.
• New files: common/jobs/calphad.py, common/flows/calphad.py, common/schemas/calphad.py and forcefields/flows/calphad.py with CalphadMaker and from_force_field_name. SQS2TDB_CMD is added to the settings.
• Tests in tests/common/jobs/test_calphad.py and tests/forcefields/flows/test_calphad.py. The energies and vibrational entropies of the tutorial runs are in tests/test_data/common/calphad, and a test refits Ni-Re from them.
• Docs in forcefields.md, with a pointer in vasp.md.
• Tutorial tutorials/calphad_workflow.ipynb for Ni-Re, Co-Ni and Cr-V with MACE-OMAT-0-medium, MACE-MATPES-PBE-0, MACE-MATPES-r2SCAN-0 and GRACE-2L-OMAT. Part A runs the workflows and the fits (cells tagged skip-execution). Part B draws all diagrams from the 26 stored TDB files in tests/test_data/common/calphad/calphad_tdbs.json.gz, and CI runs it with nbmake. It compares the diagrams with experiment and with Zhu et al. (npj Comput. Mater. 11, 340 (2025)). It also tests the size and temperature of the liquid cell and the number of liquid runs (ten 20 ps runs and one 100 ps run of each liquid SQS).
• The fit job stops with an error when sqs2tdb would carry on without one. This covers a missing energy, a link to a folder that does not exist, an empty TDB parameter and an undefined reference energy when the stable lattice of an element is not fitted. It gives a warning for a relaxation that did not converge or has a relaxation strain above 0.1.
• The relaxation strain (as ATAT checkrelax reports it), the force convergence of each relaxation, the fraction of imaginary phonon frequencies of each solid, and the mean squared displacement and energy standard error of each liquid run are stored, so a user can check them. A lattice with an SQS above max_imaginary_fraction is fitted to the energies only, with a warning.
• The phonon step reuses the supercell, displacement and force jobs of the phonon workflow. The entropy is phonopy's sum over a q-point mesh at 3000 K without the three acoustic modes at Gamma, as in PhononBSDOSDoc since Breaking: Use phonopy's sum over the q-point mesh for the phonon thermal properties #1576. With diagonal supercells of at least 20 Å, it changed by at most 0.014 k_B/atom up to 25 Å.
• Part of Tracking: atomate2 PRs from Anubhav Jain's group #1591.

Additional dependencies introduced (if any)

• List all new dependencies needed and justify why. While adding dependencies that bring
  significantly useful functionality is perfectly fine, adding ones that add trivial
  functionality, e.g., to use one single easily implementable function, is frowned upon.
  Justify why that dependency is needed. Especially frowned upon are circular dependencies.

• No new Python dependencies. The workflow needs the existing phonons extra. ATAT is not on PyPI, and its license does not allow us to ship a modified copy. The CI therefore downloads the official ATAT tarball, checks its sha256 and builds the three ATAT programs that sqs2tdb needs. This takes about 15 s. It runs in the test-non-ase jobs and in the numpy-limited force field job, and a CI step checks that ATAT is on the PATH.

• The tutorial also needs pycalphad to draw the diagrams. It is not added to the package. The notebook CI job installs it.

TODO (if any)

If this is a work-in-progress, write something about what else needs to be done.

• The liquid volume comes from NPT at zero pressure in one run. The MPMorph EquilibriumVolumeMaker finds the same volume from a pressure-volume fit over several NVT runs. I can switch to it if you prefer.
• VASP support will come in a separate PR.
• A calphy workflow with free energies from thermodynamic integration will come in a separate PR with @ZKC19940412.

Checklist

Work-in-progress pull requests are encouraged, but please put [WIP] in the pull request
title.

Before a pull request can be merged, the following items must be checked:

• Code is in the standard Python style.
  The easiest way to handle this is to run the following in the correct sequence on
  your local machine. Start with running ruff and ruff format on your new code. This will
  automatically reformat your code to PEP8 conventions and fix many linting issues.
• Doc strings have been added in the Numpy docstring format.
  Run ruff on your code.
• Type annotations are highly encouraged. Run mypy to
  type check your code.
• Tests have been added for any new functionality or bug fixes.
• All linting and tests pass.

Note that the CI system will run all the above checks. But it will be much more
efficient if you already fix most errors prior to submitting the PR. It is highly
recommended that you use the pre-commit hook provided in the repository. Simply run
pre-commit install and a check will be run prior to allowing commits.

[JaGeo]

Btw, Do you know https://github.com/ICAMS/calphy ?

[ZKC19940412]

Hi Janine. Thanks for sharing info about calphy. Based on my preliminary understanding, it seems that calphy gets full free energies from thermodynamic integration in LAMMPS, and this PR follows the ATAT sqs2tdb approach, which takes the pure elements from SGTE and fits only the mixing terms to force field energies, so it is much cheaper and works with any force field in atomate2. Maybe I can add support for calphy based phase diagram as well in another PR? What do you think?

[JaGeo]

Hi Janine. Thanks for sharing info about calphy. Based on my preliminary understanding, it seems that calphy gets full free energies from thermodynamic integration in LAMMPS, and this PR follows the ATAT sqs2tdb approach, which takes the pure elements from SGTE and fits only the mixing terms to force field energies, so it is much cheaper and works with any force field in atomate2. Maybe I can add support for calphy based phase diagram as well in another PR? What do you think?

Yes, it is a different approach indeed. I was just wondering if you are aware of this approach. And, yes, if you are interested in implementing this, I am happy to accept a PR

[hrushikesh-s]

Yes, it is a different approach indeed. I was just wondering if you are aware of this approach. And, yes, if you are interested in implementing this, I am happy to accept a PR

Sounds good.. @ZKC19940412 and I will create a new PR for calphy, and maybe also compare it with the results of this current PR... Might be a good idea to see how much they differ in their final phase diagrams for binaries.

[ZKC19940412]

The PR is in good shape. The atomate2 workflow with mace-omat-medium reproduces the phase diagram about as well as native PhaseForge.

One suggestion: the liquid-phase MD sampling needs tens of short MD runs to converge. We should either benchmark how the number of runs affects the liquid-phase region of the phase diagram, or add a way to decide how many runs are needed automatically.

For the second option, PhaseForge has something we could consider implementing: In src/ternary_search.cpp (https://github.com/dogusariturk/PhaseForge/blob/main/src/ternary_search.cpp), they run a ternary search over the POSCAR lattice scaling factor to find the equilibrium volume. Each iteration evaluates the energy (one MD run) at two interior points of the bracket, keeps the side containing the minimum, and shrinks the bracket to 2/3 of its width until it is smaller than eps. With the default bracket of 0.7 to 1.3 and eps of 0.01, that is around 20 MD runs with a deterministic stopping rule, instead of a fixed dense grid of volumes.

One caveat: their search is adaptive sampling over volume, assuming energy vs volume is unimodal in the bracket. It cuts down runs spent scanning the volume axis, but it is not a convergence criterion for the statistical averaging of the liquid itself. So it is the adaptive-search idea we would borrow, not a drop-in solution. If our bottleneck is independent runs at the same state point, we would still need a separate stopping rule, for example based on the standard error of the averaged quantity.

[hrushikesh-s]

Thanks! I ran each 864-atom liquid SQS ten times for 20 ps, plus once for 100 ps, for Ni-Re, Co-Ni and Cr-V with all four potentials. One 20 ps run gives the liquid L0 to within 0.4 to 1.2 kJ/mol and the liquidus maximum or minimum to within 11 to 25 K, the mean of several runs scatters less roughly with one over the square root of their number, and the 100 ps runs show no drift, so the error is statistical. A volume search is not needed here, since the liquid MD runs NPT at zero pressure and finds its volume in the same run. For the stopping rule, the workflow now stores the standard error of each liquid energy from five blocks of the run, and it gives the scatter between runs to within about 30%, so a single run already shows whether more are needed. Should I keep one run per SQS as the default and point users to energy_standard_error for deciding when to average more runs, as the tutorial now does? or is there something else that you suggest to do instead?

[ZKC19940412]

Nice benchmark, and I agree the error is statistical. I think MD protocol and the error metrics as shown in the tutorial works. Once you had the vibrational entropy and short range order as options whenever the users want, I think it is close to complete shape.

[hrushikesh-s]

Nice benchmark, and I agree the error is statistical. I think MD protocol and the error metrics as shown in the tutorial works. Once you had the vibrational entropy and short range order as options whenever the users want, I think it is close to complete shape.

@ZKC19940412, the vibrational entropy is now in the workflow. I have also added short_range_order, which passes the -sro option to the sqs2tdb fit and adds its CVM correction from the same SQS energies.
Does this cover what you had in mind?

[ZKC19940412]

Also the current workflow is one to one interchangeable to VASP, with small change in the parameters in case vasp MD for liquid phase gets expensive.

[hrushikesh-s]

Hi Janine. Thanks for sharing info about calphy. Based on my preliminary understanding, it seems that calphy gets full free energies from thermodynamic integration in LAMMPS, and this PR follows the ATAT sqs2tdb approach, which takes the pure elements from SGTE and fits only the mixing terms to force field energies, so it is much cheaper and works with any force field in atomate2. Maybe I can add support for calphy based phase diagram as well in another PR? What do you think?

Yes, it is a different approach indeed. I was just wondering if you are aware of this approach. And, yes, if you are interested in implementing this, I am happy to accept a PR

@JaGeo, since version 2.1.0 calphy is under the Academic Software Licence, which allows only non-commercial academic use. Would it be fine to add a calphy workflow with calphy as an optional dependency that users install themselves? The force fields would run through LAMMPS, either MACE with ML-IAP or any ASE calculator with Rootstock, which already works with calphy (calphy#256).

This would also help #1560, since the free energy and phase transition analysis in Section 5 of its tutorial could then be done with calphy. or is there something else that you suggest to do instead?

[JaGeo]

Sounds good. I will continue tomorrow.