[WIP] Add PheasyCalc class for higher order phonon calculations

[leslie-zheng]

Summary

Major changes:

• feature 1: Create Pheasy branch to extract second, third and fourth, up to sixth order force constants.
• feature 2: Implement FourPhonon to calculate the thermal conductivity considering 3ph and 4ph scattering processes.
• fix 1: ...

Todos

If this is work in progress, what else needs to be done?

• feature 2:
• fix 2:

Checklist

• Google format doc strings added. Check with ruff.
• Type annotations included. Check with mypy.
• Tests added for new features/fixes.
• If applicable, new classes/functions/modules have duecredit @due.dcite decorators to reference relevant papers by DOI (example)

Tip: Install pre-commit hooks to auto-check types and linting before every commit:

pip install -U pre-commit
pre-commit install

[shyuep]

@rul048 You will take charge of reviewing PRs. But for now, I will be the one to merge them.

[shyuep]

@rul048 When are you going to review this? Are you actually interested in maintaining matcalc?

[rul048]

Here should be PheasyCalc not PhononCalc.

[rul048]

If the commented out codes here are not needed, we could remove them.

[rul048]

Same as above. If the commented out codes here are not needed, we could remove them.

[rul048]

I recommend not using (plenty of) command lines internally unless necessary. If unavoidable, please give the command a clear and descriptive name that reflects its purpose.

[shyuep]

Automated PR review generated by Claude (Opus 4.7) on behalf of @shyuep. Based purely on the diff and current GitHub Actions state — no code was executed locally.

Thanks for the contribution, @leslie-zheng — adding PheasyCalc, AlamodeCalc, and FourPhononCalc for higher-order force constants and 4-phonon thermal conductivity is a substantial and welcome capability. However, the PR in its current [WIP] state has several blockers and needs significant cleanup before it can be merged.

Blockers:

1. Hardcoded user-specific binary paths. Multiple subprocess call sites embed personal absolute paths:

   • _alamode.py: subprocess.run("mpirun -n 1 /home/jzheng4/alamode/_build/alm/alm alamode.in", shell=True, check=False) (and the alamode_anhar.in counterpart).
   • _fourphonon.py: subprocess.run(["mpirun", "-n", "1", "/home/jzheng4/FourPhonon-sampling_method/ShengBTE"], check=True).

   FourPhononCalc already exposes path_to_shengbte / path_to_fourphonon constructor params — please actually use them in the subprocess calls. AlamodeCalc should expose an analogous path_to_alamode parameter (or read from $ALAMODE / which alm). No code path should hardcode /home/jzheng4/....

2. shell=True + check=False. All subprocess.call(..., shell=True) and subprocess.run(..., shell=True, check=False) calls silently swallow failures. If alm, pheasy, or ShengBTE fail (missing binary, bad input file, MPI error), the calculator continues with stale/missing output files and returns a result that looks valid but is wrong. Use shell=False with argv lists, check=True, and capture stderr for diagnostics. shell=True is also an injection risk if any path or argument is ever derived from user input.

3. Hardcoded composition-specific identifiers in _alamode.py:

   • f.write(" PREFIX = EuZnAs_harmonic\n") and FC2XML = EuZnAs_harmonic.xml are hardcoded for the EuZnAs test system. Derive from the structure's reduced formula or accept a prefix parameter.
   • f.write("&cutoff\n *-* 18 10 7.5\n/") hardcodes cutoff radii. Parameterize using the existing cutoff_distance_cubic / quartic / etc. attributes (which are currently defined but unused in this code path).
   • L1_ALPHA = 3.82925e-06, CONV_TOL = 1.0e-8 — surface as __init__ parameters with documented defaults.

4. CI is fully red and the branch is dirty (CONFLICTING). Lint, test, and pre-commit.ci all fail; rebase onto main and run ruff check --fix . && ruff format . and mypy per @rul048's earlier ask.

5. Patch coverage is 14% (235 uncovered lines). External-binary calls are hard to unit-test, but the pure-Python logic absolutely should be — in particular _parse_xml, _get_forceconstants_xml, _write_fc2_phonopy, the input-file generators, and the unit conversion logic (Bohr/Rydberg). Add tests with small fixture XML/FORCE_CONSTANTS files.

Other issues:

• Many duplicated import subprocess statements at function scope (e.g. inside if self.calc_anharmonic:); move all imports to the module top.
• Large blocks of commented-out debug code (# subprocess.run(["mpirun", ..., # write_force_constants_to_hdf5(...), # print(np.shape(fc_full))). Please delete; git is the history.
• pheasy>=0.0.2 is added as a hard dependency in pyproject.toml, but pheasy is very early-stage and not on conda-forge. Make it optional via an extras group (e.g. pheasy = ["pheasy>=0.0.2"]) and import it lazily inside the calculator, raising a clear ImportError with install instructions if missing — this is how matcalc handles phono3py-only paths.
• AlamodeCalc.__init__ has 26 parameters; consider grouping anharmonic-only options into a dedicated dataclass / dict, and ensure every parameter is actually consumed (several look unused in the current code path).
• _calc_forces is duplicated semantically across _phonon.py and the new _alamode.py; extract to a shared helper.
• Variable name disp_array = np.array(disp_array) rebinds list → array; rename for clarity. Also for i, (disp, force) in enumerate(zip(...)) doesn't use i — drop the enumerate.
• The try: import f90nml except ImportError: f90nml = None pattern is fine, but I don't see f90nml used anywhere in the visible diff — confirm it's needed or drop.

Process suggestion: this is ~1650 lines across three independent calculators. I'd strongly recommend splitting:

• PR 1: PheasyCalc only (LASSO-fitted harmonic FCs).
• PR 2: AlamodeCalc (higher-order FCs, depends on PheasyCalc patterns).
• PR 3: FourPhononCalc (4-phonon thermal conductivity).

Each can then be reviewed and tested in isolation, which is much more tractable.

I am happy to do another pass once the blockers above are addressed and CI is green. The underlying capability — anharmonic FCs and 4-phonon transport via foundation potentials — is exactly the direction matcalc should be growing.

[shyuep]

🤖 Automated PR review (generated by Claude on behalf of @shyuep)

• ⚠️ pre-commit.ci - pr failing on the head commit. Needs to be cleaned up before merging.
• ⚠️ mergeable_state: dirty — branch has conflicts with main. Needs a rebase before merging.
• ⚠️ 1650 additions / 5 files / 125 commits, but the PR description still has empty Todos and all checklist items unchecked. WIP marker is honest, but this is too large to merge in its current state.
• ⚠️ Codecov reports patch coverage at 14.2% — _pheasy.py at 11.8%, _fourphonon.py at 16.5%. Two whole new public modules effectively untested. At minimum the public API (PheasyCalc, FourPhonon) needs smoke-test coverage.
• ⚠️ No new dependencies declared in pyproject.toml. If the external pheasy / fourphonon packages or binaries are required, add them (as extras at minimum) so the install path is auditable.
• 💡 @rul048's April 2025 review listed concrete next steps (ruff, mypy, tests, deps, docstrings). Confirm those are addressed before requesting re-review.
• 💡 Consider splitting: (1) PheasyCalc higher-order force-constant extraction and (2) FourPhonon 3ph+4ph thermal-conductivity workflow. Reviewable in isolation.

---

Generated by Claude Code