Skip to content

Fix pheasy anharmonic fitting and add fit options - #1558

Merged
JaGeo merged 12 commits into
materialsproject:mainfrom
hrushikesh-s:pheasy-anharmonic-fixes
Oct 2, 2026
Merged

JaGeo merged 12 commits into
materialsproject:mainfrom
hrushikesh-s:pheasy-anharmonic-fixes

Conversation

@hrushikesh-s

@hrushikesh-s hrushikesh-s commented Sep 24, 2026 •

Copy link
Copy Markdown
Collaborator

Summary

This PR fixes the anharmonic force constant fit in the pheasy workflow and adds options for it. The harmonic fit is unchanged, apart from a fixed LASSO seed.

Bug fixes:

  • The number of anharmonic supercells was always 20. A hard-coded line overwrote both the ALM count and num_disp_anhar.
  • The anharmonic set used the same random seed as the harmonic set, so it repeated the harmonic directions. It now uses random_seed + 1.
  • The anharmonic fit used pheasy's default estimator, which is OLS. It now uses LASSO.
  • ALM counted the anharmonic force constants with settings that did not match the pheasy fit. The cutoffs were not given per pair of elements, and the many-body terms differed from the --nbody passed to pheasy. The count now uses the same cutoffs and nbody as the fit.
  • If the harmonic fit wrote no FORCE_CONSTANTS, the anharmonic fit was skipped without an error. It now raises an error.
  • Failed pheasy commands in the anharmonic fit were ignored. They now run with check=True.
  • ALM was required even when the harmonic force constants came from finite displacements. It is now imported only where free force constants are counted.
  • Several docstrings did not match the code, for example the min_length default and the supercell kwargs.

Features:

  • anhar_max_order, 3 (default) or 4.
  • anhar_fit_methods, "cocktail" (default), "one-shot", or both. Cocktail keeps fc2 fixed to the harmonic fit. One-shot fits all orders together and writes to a one_shot folder.
  • anhar_alpha_min, the lowest LASSO penalty in the cross-validation search. A warning is raised if the chosen penalty lands on either end of the search.
  • The number of anharmonic supercells is set from the number of free force constants, with 100 force equations per free force constant. The floor is 20 and the limit is 600. An explicit num_disp_anhar is used as given.
  • The ALM counting code is moved to common/jobs/phonons.py as _get_num_irreducible_fcs and _get_num_harmonic_supercells. The pheasy and hiPhive workflows both use these helpers, so the hiPhive job no longer keeps its own copy.
  • duecredit entries for pheasy, ALM and hiPhive.
  • pheasy built its own supercell from POSCAR. That supercell can put an atom on a cell face one lattice vector away, so its atoms were not in the order of the displacement and force matrices. All pheasy commands now read the supercell from SPOSCAR with --scell SPOSCAR.
  • The short-cutoff pheasy refit, run when the harmonic fit has imaginary modes, overwrote FORCE_CONSTANTS and its result was never used. It now runs in its own folder with the displacement and force matrices, and its force constants are used.
  • With pheasy's default --tol of 1e-4, the anharmonic LASSO fits at small penalties stopped before they converged. Cross-validation then picked a penalty that could change from one machine to another. The anharmonic fit now runs with --tol 1e-8, and the same penalty is picked on every machine.

Removals:

  • renorm_phonon and renorm_temp. renorm_temp was never read. The renormalization step only reran the fit with pheasy's default OLS and did nothing with the result, so no renormalization was done. I will soon be creating a new PR with temperature-dependent effective phonons.
  • cal_ther_cond, ther_cond_mesh and ther_cond_temp. The phono3py command used flags that phono3py 4.5 does not have. The failure was ignored, and no result reached the output document.
  • Passing any of these options now raises a TypeError.

Behaviour changes:

  • The harmonic LASSO fit now passes --seed, so repeated runs give the same force constants.
  • By default only third-order force constants are fitted. Before, fourth order was always included.
  • The number of anharmonic supercells goes up from a fixed 20, often to about 100 or more.
  • The anharmonic fit uses LASSO instead of OLS.
  • The finite-displacement harmonic path no longer needs ALM.

Tests:

  • tests/common/jobs/test_pheasy.py (new) runs both jobs end to end with EMT forces on fcc Cu. It covers the dataset split, the supercell count, the penalty warnings, the cocktail and one-shot fits, and a fourth-order fit.
  • tests/vasp/flows/test_pheasy.py checks the input validation. It also checks that the maker passes the same settings to both jobs.

Additional dependencies introduced (if any)

  • pheasy is pinned to commit 9f24162 of my pheasy fork (gitlab.com/hpsahasrabuddhe/pheasy). It has two changes that are not yet in pheasy. pheasy reads the supercell given with --scell, and it builds the sensing matrix for all displaced supercells at once. The second change gives the same force constants and is about 27 times faster for a 432-atom Si supercell. I will move the pin to a pheasy release once these changes are merged there.

TODO (if any)

  • The anharmonic force constants are written to files in the job folder. They are not stored in the output document yet.
  • A follow-up PR will add a thermal expansion (CTE) workflow on top of this one.

Checklist

  • 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

JaGeo commented Sep 24, 2026

Copy link
Copy Markdown
Member

@hrushikesh-s thanks! Did you already test this on a smaller database? 😀 This would be great to know for additional context

Comment thread src/atomate2/common/flows/pheasy.py
@hrushikesh-s

Copy link
Copy Markdown
Collaborator Author

@JaGeo , yes, correct -- @ZKC19940412 is running the benchmarking tests on these, and will post the results in a few days for you to take a look 👍

@leslie-zheng

Copy link
Copy Markdown
Contributor

@JaGeo , yes, correct -- @ZKC19940412 is running the benchmarking tests on these, and will post the results in a few days for you to take a look 👍

Hi Hrushikesh, if you want, Maybe you also can add me into this branch, I can provide some good strategies and coding to make the anharmonic FCs fitting, finite-temperature phonon or phonon-related properties more efficient, accurate and useful for users, especially for high-throughput calculations and screening. Otherwise, if the phonon-related code is too slow and tedious, it could not attract too many users to try. best.

@hrushikesh-s

Copy link
Copy Markdown
Collaborator Author

@leslie-zheng , sounds good! Pls feel free to make edits/comments on this PR as you think necessary.

@hrushikesh-s

Copy link
Copy Markdown
Collaborator Author

Hi @JaGeo, this is ready for your review. I updated the description with the changes since the last round. The main one is that the anharmonic LASSO fit now runs with --tol 1e-8. The force constants do not change with the bound.

On your question about testing: we checked the fits against finite displacements with phono3py for six materials and three MACE potentials. The results are in the tutorial in #1559.

Comment thread src/atomate2/common/flows/hiphive.py
Comment thread src/atomate2/common/flows/pheasy.py
Comment thread src/atomate2/common/flows/pheasy.py Outdated
@JaGeo

JaGeo commented Oct 2, 2026

Copy link
Copy Markdown
Member

Thanks. Where exactly did you adapt the tutorial? At least, I don't see a new output there for anharmonic data.

Everything else looks good

@hrushikesh-s

Copy link
Copy Markdown
Collaborator Author

Thanks. Where exactly did you adapt the tutorial? At least, I don't see a new output there for anharmonic data.

Everything else looks good

The tutorial for anharmonic CTE calcs is in PR #1559

@JaGeo

JaGeo commented Oct 2, 2026

Copy link
Copy Markdown
Member

@hrushikesh-s one last thing: should we add a note in the documentation that the workflow might undergo further changes and has not been as widely tested as the harmonic one?

@hrushikesh-s

Copy link
Copy Markdown
Collaborator Author

@hrushikesh-s one last thing: should we add a note in the documentation that the workflow might undergo further changes and has not been as widely tested as the harmonic one?

Yes, let me add that in now!

@hrushikesh-s

Copy link
Copy Markdown
Collaborator Author

@JaGeo , done!

@JaGeo
JaGeo enabled auto-merge (squash) October 2, 2026 15:42
@JaGeo

JaGeo commented Oct 2, 2026

Copy link
Copy Markdown
Member

Will be merged as soon as all tests pass

@JaGeo
JaGeo merged commit 57adeb3 into materialsproject:main Oct 2, 2026
17 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants