Fix Grüneisen support for Phonopy 4

Co-authored-by: Cursor <cursoragent@cursor.com>

Author: Zikkying

README.md

@@ -798,7 +798,7 @@ APEX integrates parts of [dflow-phonon](https://github.com/Chengqian-Zhang/dflow
 | `approach` | String | `"linear"` | VASP phonon method: `"linear"` or `"displacement"`. |
 | `supercell_size` | Sequence[Int] | `[2, 2, 2]` | Supercell dimensions. |
 | `MESH` | Sequence[Int] | `None` | Reciprocal-space mesh (e.g., `[8, 8, 8]`). |
-| `PRIMITIVE_AXES` | String | `None` | Custom primitive axes definition (`"0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0"`). |
+| `PRIMITIVE_AXES` | String | `"P"` | Primitive axes used consistently by phonopy and phonoLAMMPS; accepts `"P"`, `"AUTO"`, or a 3×3 matrix. |
 | `BAND` | String | `None` | Band path definition (falls back to SeeK-path when omitted). |
 | `BAND_LABELS` | String | `None` | Labels for band segments. |
 | `BAND_POINTS` | Integer | `51` | Number of sampling points per segment. |
@@ -822,8 +822,8 @@ APEX supports Grüneisen workflows based on phonon calculations at multiple stra
 | `primitive` | Bool | `false` | Reduce to primitive cell before phonon calculation. |
 | `approach` | String | `"linear"` | Phonon workflow approach; VASP Grüneisen currently uses linear response. |
 | `supercell_size` | Sequence[Int] | `[2, 2, 2]` | Phonon supercell dimensions. |
-| `MESH` | Sequence[Int] | `None` | Reciprocal-space mesh for mode summation. |
-| `PRIMITIVE_AXES` | String | `None` | Custom primitive axes definition. |
+| `MESH` | Sequence[Int] | `[20, 20, 20]` | Three positive reciprocal-space mesh dimensions for mode summation. |
+| `PRIMITIVE_AXES` | String | `"P"` | Primitive axes used consistently for displacement generation and post-processing. |
 | `BAND` | String | `None` | Band path definition (falls back to SeeK-path when omitted). |
 | `BAND_LABELS` | String | `None` | Labels for band segments. |
 | `BAND_POINTS` | Integer | `51` | Number of sampling points per segment. |

apex/core/property/Gruneisen.py

@@ -32,6 +32,7 @@
 }
 
 DEFAULT_SUPERCELL = [2, 2, 2]
+DEFAULT_MESH = [20, 20, 20]
 THZ_TO_K = 47.99243073366221
 KB_EV_PER_K = 8.617333262145e-5
 
@@ -49,7 +50,8 @@ def __init__(self, parameter, inter_param=None):
         self.seekpath_from_original = parameter["seekpath_from_original"]
         parameter["seekpath_param"] = parameter.get("seekpath_param", {})
         self.seekpath_param = parameter["seekpath_param"]
-        parameter["MESH"] = parameter.get("MESH", None)
+        mesh = parameter.get("MESH")
+        parameter["MESH"] = list(DEFAULT_MESH) if mesh is None else mesh
         self.MESH = parameter["MESH"]
         parameter["PRIMITIVE_AXES"] = parameter.get("PRIMITIVE_AXES", "P")
         self.PRIMITIVE_AXES = parameter["PRIMITIVE_AXES"]
@@ -88,6 +90,17 @@ def __init__(self, parameter, inter_param=None):
         self._validate()
 
     def _validate(self):
+        if (
+            not isinstance(self.MESH, (list, tuple))
+            or len(self.MESH) != 3
+            or any(
+                not isinstance(value, int)
+                or isinstance(value, bool)
+                or value <= 0
+                for value in self.MESH
+            )
+        ):
+            raise ValueError("gruneisen MESH must contain exactly 3 positive integers")
         if len(self.volume_strains) < 3:
             raise ValueError("volume_strains must contain at least 3 points")
         if len(self.temperatures) < 1:
@@ -326,8 +339,11 @@ def _make_abacus_confs(self, path_to_work: str, path_to_equi: str) -> List[str]:
                         "ATOM_NAME =%s\n"
                         % "".join(f" {name}" for name in vasp_utils.get_poscar_types("POSCAR"))
                     )
-                    if self.PRIMITIVE_AXES:
-                        fp.write(f"PRIMITIVE_AXES = {self.PRIMITIVE_AXES}\n")
+                    primitive_axes = Phonon.primitive_axes_config_value(
+                        self.PRIMITIVE_AXES
+                    )
+                    if primitive_axes:
+                        fp.write(f"PRIMITIVE_AXES = {primitive_axes}\n")
                 subprocess.check_call(
                     Phonon.phonopy_setup_command("setting.conf --abacus -d"),
                     shell=True,
@@ -682,10 +698,10 @@ def _build_band_payload(self, poscar_path: str) -> dict:
                 self.supercell_size[2],
             )
         )
-        if self.MESH:
-            lines.append("MESH = %s %s %s" % (self.MESH[0], self.MESH[1], self.MESH[2]))
-        if self.PRIMITIVE_AXES:
-            lines.append(f"PRIMITIVE_AXES = {self.PRIMITIVE_AXES}")
+        lines.append("MESH = %s %s %s" % (self.MESH[0], self.MESH[1], self.MESH[2]))
+        primitive_axes = Phonon.primitive_axes_config_value(self.PRIMITIVE_AXES)
+        if primitive_axes:
+            lines.append(f"PRIMITIVE_AXES = {primitive_axes}")
         lines.append(f"BAND = {band_string}")
         if band_labels:
             lines.append(f"BAND_LABELS = {band_labels}")
@@ -801,6 +817,16 @@ def _ensure_mesh_yaml(self, task_dir: str) -> None:
             raise FileNotFoundError(f"FORCE_CONSTANTS not found in {task_dir}")
         if not os.path.isfile(band_conf):
             raise FileNotFoundError(f"band.conf not found in {task_dir}")
+        with open(band_conf, "r") as fp:
+            has_mesh = any(
+                line.partition("=")[0].strip().upper() == "MESH"
+                for line in fp
+            )
+        if not has_mesh:
+            raise ValueError(
+                f"band.conf in {task_dir} does not define MESH; "
+                f"Gruneisen requires MESH={self.MESH}. Regenerate the task inputs."
+            )
         if not os.path.isfile(poscar):
             raise FileNotFoundError(f"POSCAR not found in {task_dir}")
         os.chdir(task_dir)

apex/core/property/Phonon.py

@@ -5,6 +5,7 @@
 import shutil
 import re
 import subprocess
+from fractions import Fraction
 from typing import List, Dict, Any
 
 import dpdata
@@ -47,6 +48,8 @@ def phonopy_writefc_commands(arguments: str) -> List[str]:
         phonopy_command = Phonon.phonopy_command(arguments)
         if setup_command == phonopy_command:
             return [setup_command]
+        if "--dim" in arguments:
+            return [setup_command]
         return [setup_command, phonopy_command]
 
     @staticmethod
@@ -103,15 +106,17 @@ def phonopy_load_commands(
                     ),
                 )
             )
-            commands.append(
-                Phonon.phonopy_command(
-                    Phonon._join_phonopy_arguments(
-                        extra_args,
-                        '--dim=%s -c %s' % (Phonon._format_dim(supercell_size), cell_file),
-                        config_file,
+            if not shutil.which("phonopy-init"):
+                commands.append(
+                    Phonon.phonopy_command(
+                        Phonon._join_phonopy_arguments(
+                            extra_args,
+                            '--dim=%s -c %s'
+                            % (Phonon._format_dim(supercell_size), cell_file),
+                            config_file,
+                        )
                     )
                 )
-            )
         else:
             commands.append(
                 Phonon.phonopy_command(
@@ -137,10 +142,51 @@ def phonopy_writefc_load_commands(
 
     @staticmethod
     def primitive_axes_setup_argument(primitive_axes: str | None) -> str:
+        primitive_axes = Phonon.primitive_axes_config_value(primitive_axes)
         if not primitive_axes:
             return ""
         return f"--pa {primitive_axes}"
 
+    @staticmethod
+    def primitive_axes_config_value(primitive_axes: Any) -> str | None:
+        if primitive_axes is None:
+            return None
+        if isinstance(primitive_axes, str):
+            value = primitive_axes.strip()
+            if not value or value.upper() == "AUTO":
+                return None
+            if value.upper() == "P":
+                return "P"
+            tokens = value.replace(",", " ").split()
+        elif isinstance(primitive_axes, (list, tuple)):
+            tokens = []
+            for item in primitive_axes:
+                if isinstance(item, (list, tuple)):
+                    tokens.extend(item)
+                else:
+                    tokens.append(item)
+        else:
+            raise ValueError(
+                "PRIMITIVE_AXES must be 'P', 'AUTO', or a 3x3 numeric matrix"
+            )
+        if len(tokens) != 9:
+            raise ValueError(
+                "PRIMITIVE_AXES must be 'P', 'AUTO', or contain exactly 9 values"
+            )
+        try:
+            return " ".join(str(float(Fraction(str(token)))) for token in tokens)
+        except (ValueError, ZeroDivisionError) as exc:
+            raise ValueError("PRIMITIVE_AXES contains a non-numeric value") from exc
+
+    @staticmethod
+    def primitive_axes_phonolammps_argument(primitive_axes: Any) -> str:
+        config_value = Phonon.primitive_axes_config_value(primitive_axes)
+        if config_value is None:
+            return ""
+        if config_value == "P":
+            config_value = "1 0 0 0 1 0 0 0 1"
+        return f"-pa {config_value}"
+
     @staticmethod
     def run_first_success(commands: List[str], required_file: str | None = None) -> None:
         errors = []
@@ -154,6 +200,12 @@ def run_first_success(commands: List[str], required_file: str | None = None) ->
                 return
             errors.append(FileNotFoundError(f"{required_file} was not created by: {command}"))
         if errors:
+            if required_file is not None:
+                attempted = "; ".join(commands)
+                raise RuntimeError(
+                    f"{required_file} was not created after trying: {attempted}. "
+                    f"Last error: {errors[-1]}"
+                ) from errors[-1]
             raise errors[-1]
 
     @staticmethod
@@ -253,16 +305,30 @@ def _ensure_deepmd_plugin_loaded(self, input_text: str) -> str:
     def _build_phonolammps_run_command(self) -> str:
         dim_x, dim_y, dim_z = self.supercell_size
         command_template = self.phonolammps_run_command
+        primitive_axes = self.primitive_axes_phonolammps_argument(
+            self.PRIMITIVE_AXES
+        )
         if not command_template:
-            return f"phonolammps in.lammps -c POSCAR --dim {dim_x} {dim_y} {dim_z} "
-        return command_template.format(
+            return self._join_phonopy_arguments(
+                f"phonolammps in.lammps -c POSCAR --dim {dim_x} {dim_y} {dim_z}",
+                primitive_axes,
+            )
+        command = command_template.format(
             input_file="in.lammps",
             poscar="POSCAR",
             dim=f"{dim_x} {dim_y} {dim_z}",
             dim_x=dim_x,
             dim_y=dim_y,
             dim_z=dim_z,
+            primitive_axes=primitive_axes,
         )
+        if (
+            primitive_axes
+            and "{primitive_axes}" not in command_template
+            and not re.search(r"(?:^|\s)(?:-pa|--primitive_axis)(?:\s|=)", command)
+        ):
+            command = self._join_phonopy_arguments(command, primitive_axes)
+        return command
 
     def make_confs(self, path_to_work, path_to_equi, refine=False):
         path_to_work = os.path.abspath(path_to_work)
@@ -408,8 +474,9 @@ def make_confs(self, path_to_work, path_to_equi, refine=False):
                     ret += "MESH = %s %s %s\n" % (
                         self.MESH[0], self.MESH[1], self.MESH[2]
                     )
-                if self.PRIMITIVE_AXES:
-                    ret += "PRIMITIVE_AXES = %s\n" % self.PRIMITIVE_AXES
+                primitive_axes = self.primitive_axes_config_value(self.PRIMITIVE_AXES)
+                if primitive_axes:
+                    ret += "PRIMITIVE_AXES = %s\n" % primitive_axes
                 ret += "BAND = %s\n" % self.BAND
                 if self.BAND_LABELS:
                     ret += "BAND_LABELS = %s\n" % self.BAND_LABELS

apex/skills/apex-flow/reference/properties.md

@@ -233,7 +233,7 @@ These parameters appear in most property configurations:
 | `PRIMITIVE_AXES` | optional | str | `"P"` | Phonopy PRIMITIVE_AXES |
 | `BAND` | optional | str | `null` | Explicit band path |
 | `BAND_LABELS` | optional | list[str] | `null` | High-symmetry labels |
-| `phonolammps_run_command` | optional | str | `null` | Custom phonoLAMMPS command |
+| `phonolammps_run_command` | optional | str | `null` | Custom phonoLAMMPS command; supports `{primitive_axes}`, otherwise the matching `-pa` argument is appended |
 
 **cal_setting defaults**: `relax_pos=true`, `relax_shape=false`, `relax_vol=false`, `cal_type="static"`
 
@@ -494,15 +494,17 @@ VASP uses Langevin–Parrinello–Rahman NpT; ABACUS uses Nose–Hoover-style Np
 | Parameter | Required? | Type | Default | Description |
 |-----------|-----------|------|---------|-------------|
 | `supercell_size` | optional | list[int] | `[2,2,2]` | Phonon supercell |
+| `MESH` | optional | list[int] | `[20,20,20]` | Three positive reciprocal-space mesh dimensions used for mode summation |
 | `volume_strains` | **REQUIRED** | list[float] | `[-0.02,-0.01,0.0,0.01,0.02]` | Must include 0.0, ≥3 points |
 | `temperatures` | **REQUIRED** | list[float] | `[100,200,300,400,500]` | Temperature points for evaluation |
-| `alpha_mode` | optional | str | `"full"` | `"full"` or `"debye"` |
+| `alpha_mode` | optional | str | `"sign_only"` | `"sign_only"` or `"full"` |
 
 **Complete working default**:
 ```json
 {
     "type": "gruneisen",
     "supercell_size": [2, 2, 2],
+    "MESH": [20, 20, 20],
     "volume_strains": [-0.02, -0.01, 0.0, 0.01, 0.02],
     "temperatures": [100, 200, 300, 400, 500],
     "alpha_mode": "full"

apex/skills/apex-flow/scripts/generate_config.py

@@ -127,6 +127,7 @@
     "gruneisen": {
         "type": "gruneisen",
         "supercell_size": [2, 2, 2],
+        "MESH": [20, 20, 20],
         "volume_strains": [-0.02, -0.01, 0.0, 0.01, 0.02],
         "temperatures": [100, 200, 300, 400, 500],
         "alpha_mode": "full",

apex/skills/apex-flow/scripts/validate_inputs.py

@@ -210,6 +210,22 @@ def validate_properties(properties: list, interaction_type: str) -> list:
                 errors.append(f"{prefix}: decohesive requires 'min_slab_size'")
 
         elif prop_type == "gruneisen":
+            mesh = prop.get("MESH")
+            if mesh is None:
+                mesh = [20, 20, 20]
+            if (
+                not isinstance(mesh, (list, tuple))
+                or len(mesh) != 3
+                or any(
+                    not isinstance(value, int)
+                    or isinstance(value, bool)
+                    or value <= 0
+                    for value in mesh
+                )
+            ):
+                errors.append(
+                    f"{prefix}: gruneisen 'MESH' must contain 3 positive integers"
+                )
             if "volume_strains" not in prop:
                 errors.append(f"{prefix}: gruneisen requires 'volume_strains'")
             else:

tests/test_gruneisen.py

@@ -45,6 +45,10 @@ def valid_gruneisen_params(**overrides):
         ({"volume_strains": [-0.02, 0.0, 0.0, 0.02]}, "duplicates"),
         ({"temperatures": [50, 10]}, "temperatures must be strictly increasing"),
         ({"volume_strains": [-0.02, 0.0, 0.03]}, "symmetric"),
+        ({"MESH": []}, "MESH"),
+        ({"MESH": [20, 20]}, "MESH"),
+        ({"MESH": [20, 0, 20]}, "MESH"),
+        ({"MESH": [20, True, 20]}, "MESH"),
     ],
 )
 def test_gruneisen_validation_rejects_low_cost_invalid_options(overrides, message):
@@ -67,6 +71,10 @@ def test_gruneisen_validation_rejects_low_cost_invalid_options(self):
             ({"volume_strains": [-0.02, 0.0, 0.0, 0.02]}, "duplicates"),
             ({"temperatures": [50, 10]}, "temperatures must be strictly increasing"),
             ({"volume_strains": [-0.02, 0.0, 0.03]}, "symmetric"),
+            ({"MESH": []}, "MESH"),
+            ({"MESH": [20, 20]}, "MESH"),
+            ({"MESH": [20, 0, 20]}, "MESH"),
+            ({"MESH": [20, True, 20]}, "MESH"),
         ]
         for overrides, message in cases:
             with self.subTest(overrides=overrides):
@@ -110,6 +118,7 @@ def test_task_param_defaults(self):
         self.assertEqual(task_param["supercell_size"], [2, 2, 2])
         self.assertEqual(task_param["approach"], "linear")
         self.assertEqual(task_param["PRIMITIVE_AXES"], "P")
+        self.assertEqual(task_param["MESH"], [20, 20, 20])
 
     def test_validation_rejects_invalid_schema(self):
         with self.assertRaises(ValueError):
@@ -166,6 +175,10 @@ def fake_check_call(command, shell):
             self.assertGreater(volume_data["volume"], 0.0)
             self.assertGreater(volume_data["volume_per_atom"], 0.0)
             self.assertTrue((Path(task_dir) / "band.conf").is_file())
+            self.assertIn(
+                "MESH = 20 20 20",
+                (Path(task_dir) / "band.conf").read_text(),
+            )
 
         self.assertTrue((self.target_path / "band_path.json").is_file())
 
@@ -266,6 +279,23 @@ def fake_check_call(command, shell):
         self.assertIn("-c POSCAR-unitcell", calls[1])
         self.assertIn("--nomeshsym", calls[1])
 
+    def test_ensure_mesh_yaml_rejects_legacy_band_conf_without_mesh(self):
+        task_dir = self.work_root / "missing_mesh" / "task.000000"
+        task_dir.mkdir(parents=True)
+        (task_dir / "FORCE_CONSTANTS").write_text("fake force constants\n")
+        (task_dir / "band.conf").write_text("BAND = 0 0 0  0.5 0 0\n")
+        (task_dir / "POSCAR").write_text(self.source_path.read_text())
+
+        gruneisen = Gruneisen(
+            valid_gruneisen_params(),
+            inter_param={"type": "deepmd"},
+        )
+        with self.assertRaisesRegex(
+            ValueError,
+            r"does not define MESH.*MESH=\[20, 20, 20\]",
+        ):
+            gruneisen._ensure_mesh_yaml(str(task_dir))
+
     def test_sign_only_compute_lower_from_vasp_displacement_manifest(self):
         gruneisen = Gruneisen(
             {
@@ -516,6 +546,7 @@ def test_post_process_prepares_phonon_run_inputs_for_lammps(self):
         self.assertIn("/root/.dp1s/bin/lmp -in in.relax.lammps", run_script)
         self.assertIn("python3 convert_relax_dump_to_poscar.py dump.relax POSCAR.relaxed type_map.json", run_script)
         self.assertIn("phonolammps in.lammps -c POSCAR --dim 2 2 2", run_script)
+        self.assertIn("-pa 1 0 0 0 1 0 0 0 1", run_script)
         self.assertIn("cp POSCAR.relaxed POSCAR", run_script)
 
     def test_lammps_backward_files_for_gruneisen(self):