Add --site-mixture/--split-site-mixture options to velph init

src/phelel/velph/cli/init/cmd_init.py

@@ -180,6 +180,30 @@
         f'default="{PrimitiveCellChoice.STANDARDIZED.value}")'
     ),
 )
+@click.option(
+    "--site-mixture",
+    "site_mixture",
+    type=str,
+    default=None,
+    help=(
+        "[Experimental] Per-atom concentration weights in input-cell atom "
+        'order, e.g., "0.5 0.5 0.5 0.5". Weights of co-located atoms (same '
+        "fractional position) must sum to 1.0; an isolated atom must be 1.0. "
+        "Cannot be combined with --magmom. "
+        f"(site_mixture: str, default={VelphInitParams.site_mixture})"
+    ),
+)
+@click.option(
+    "--split-site-mixture",
+    "split_site_mixture",
+    is_flag=True,
+    default=None,
+    help=(
+        "[Experimental] Keep co-located atoms as separate weighted species "
+        "instead of merging them (currently the only supported scheme). "
+        f"(split_site_mixture: bool, default={VelphInitParams.split_site_mixture})"
+    ),
+)
 @click.option(
     "--supercell-matrix",
     "supercell_matrix",
@@ -256,6 +280,8 @@ def cmd_init(
     plusminus: bool | None,
     primitive_cell_choice: Literal["standardized", "reduced"] | None,
     project_folder: str,
+    site_mixture: str | None,
+    split_site_mixture: bool | None,
     supercell_dimension: tuple[int, int, int] | None,
     supercell_matrix: tuple[int, int, int, int, int, int, int, int, int] | None,
     symmetrize_cell: bool | None,
@@ -304,6 +330,8 @@ def cmd_init(
         phelel_nosym=phelel_nosym,
         plusminus=_plusminus,
         primitive_cell_choice=primitive_cell_choice,
+        site_mixture=site_mixture,
+        split_site_mixture=split_site_mixture,
         supercell_dimension=supercell_dimension,
         supercell_matrix=supercell_matrix,
         symmetrize_cell=symmetrize_cell,

src/phelel/velph/cli/init/init.py

@@ -25,7 +25,9 @@
 
 from phonopy.structure.atoms import PhonopyAtoms
 from phonopy.structure.cells import (
+    apply_site_mixture,
     estimate_supercell_matrix,
+    generate_standardized_cells,
     get_supercell,
     shape_supercell_matrix,
 )
@@ -43,7 +45,6 @@
 )
 from phelel.velph.templates import default_template_dict
 from phelel.velph.utils.structure import (
-    generate_standardized_cells,
     get_primitive_cell,
     get_reduced_cell,
     get_symmetry_dataset,
@@ -166,6 +167,20 @@ def _run_init(
         magmom_vals = VaspIncar().expand(vip.magmom.strip())
         input_cell.magnetic_moments = magmom_vals
 
+    #
+    # Apply site-mixture per-atom concentration weights (non-merge scheme).
+    #
+    if vip.site_mixture is not None:
+        if vip.magmom is not None:
+            raise click.ClickException(
+                "--site-mixture cannot be combined with --magmom."
+            )
+        weights = [float(x) for x in vip.site_mixture.split()]
+        try:
+            input_cell = apply_site_mixture(input_cell, weights, symprec=vip.tolerance)
+        except ValueError as e:
+            raise click.ClickException(str(e)) from e
+
     #
     # Define cells and find crystal symmetry.
     #
@@ -571,7 +586,7 @@ def _get_cells(
                 "sym_dataset must be SpglibDataset or SpglibMagneticDataset."
             )
         unitcell, _primitive, tmat = generate_standardized_cells(
-            sym_dataset, tolerance=tolerance
+            input_cell, sym_dataset, symprec=tolerance
         )
         if find_primitive:
             primitive = _primitive
@@ -1512,12 +1527,17 @@ def _get_cell_toml_lines(
         magnetic_moments = [None] * len(unitcell.symbols)
     else:
         magnetic_moments = unitcell.magnetic_moments
-    for i, (s, v, m, mag) in enumerate(
+    if unitcell.mixture_weights is None:
+        weights = [None] * len(unitcell.symbols)
+    else:
+        weights = unitcell.mixture_weights
+    for i, (s, v, m, mag, w) in enumerate(
         zip(
             unitcell.symbols,
             unitcell.scaled_positions,
             masses,
             magnetic_moments,
+            weights,
             strict=True,
         )
     ):
@@ -1532,6 +1552,8 @@ def _get_cell_toml_lines(
             else:
                 mag_str = f"[ {mag[0]:.8f}, {mag[1]:.8f}, {mag[2]:.8f} ]"
             lines.append(f"magnetic_moment = {mag_str}")
+        if w is not None:
+            lines.append(f"weight = {w:.15f}")
     return lines
 
 

src/phelel/velph/cli/utils.py

@@ -89,6 +89,8 @@ class VelphInitParams:
     kspacing: float = 0.1
     kspacing_dense: float = 0.05
     magmom: str | None = None
+    site_mixture: str | None = None
+    split_site_mixture: bool = False
     phelel_displacement_options: DisplacementOptions | None = None
     phonopy_displacement_options: DisplacementOptions | None = None
     phono3py_displacement_options: DisplacementOptions | None = None
@@ -122,6 +124,8 @@ class VelphInitOptions:
     kspacing_dense: float | None = None
     magmom: str | None = None
     max_num_atoms: int | None = None
+    site_mixture: str | None = None
+    split_site_mixture: bool | None = None
     phelel_nosym: bool | None = None
     plusminus: bool | Literal["auto"] | None = True
     primitive_cell_choice: Literal["standardized", "reduced"] | None = None

src/phelel/velph/utils/structure.py

@@ -5,8 +5,7 @@
 import numpy as np
 import spglib
 from numpy.typing import NDArray
-from phonopy.interface.vasp import sort_positions_by_symbols
-from phonopy.structure.atoms import PhonopyAtoms, get_atomic_data
+from phonopy.structure.atoms import PhonopyAtoms
 from phonopy.structure.cells import (
     get_primitive,
     get_primitive_matrix_by_centring,
@@ -15,21 +14,6 @@
 from spglib import SpglibDataset, SpglibMagneticDataset
 
 
-def generate_standardized_cells(
-    sym_dataset: SpglibDataset | SpglibMagneticDataset,
-    tolerance: float = 1e-5,
-) -> tuple[PhonopyAtoms, PhonopyAtoms, NDArray]:
-    """Return standardized unit cell and primitive cell."""
-    convcell = _get_standardized_unitcell(sym_dataset)
-    pmat = _get_primitive_matrix_from_dataset(sym_dataset)
-    if (np.abs(pmat - np.eye(3)) < 1e-8).all():
-        primitive = convcell
-    else:
-        primitive = get_primitive(convcell, primitive_matrix=pmat, symprec=tolerance)
-
-    return convcell, primitive, pmat
-
-
 def get_primitive_cell(
     cell: PhonopyAtoms,
     sym_dataset: SpglibDataset | SpglibMagneticDataset,
@@ -38,7 +22,7 @@ def get_primitive_cell(
     """Return primitive cell and transformation matrix.
 
     This primitive cell is generated from the input cell without
-    rigid rotation in contrast to `_get_standardized_unitcell`.
+    rigid rotation in contrast to the spglib-standardized cell.
 
     """
     tmat = sym_dataset.transformation_matrix
@@ -95,43 +79,3 @@ def _get_primitive_matrix_from_dataset(
         assert spg_type is not None
         centring = spg_type.international_short[0]
     return get_primitive_matrix_by_centring(centring)
-
-
-def _get_standardized_unitcell(
-    dataset: SpglibDataset | SpglibMagneticDataset,
-) -> PhonopyAtoms:
-    """Return conventional unit cell.
-
-    This conventional unit cell can include rigid rotation with respect to
-    input unit cell for which symmetry was analized.
-
-    Parameters
-    ----------
-    cell : PhonopyAtoms
-        Input cell.
-    dataset : SpgliDataset
-        Symmetry dataset of spglib.
-
-    Returns
-    -------
-    PhonopyAtoms
-        Convetional unit cell.
-
-    """
-    std_positions = dataset.std_positions
-    std_types = dataset.std_types
-    _, _, _, perm = sort_positions_by_symbols(std_types, std_positions)
-    atom_data = get_atomic_data().atom_data
-    if isinstance(dataset, SpglibDataset):
-        return PhonopyAtoms(
-            cell=dataset.std_lattice,
-            scaled_positions=std_positions[perm],
-            symbols=[atom_data[n][1] for n in std_types[perm]],
-        )
-    else:
-        return PhonopyAtoms(
-            cell=dataset.std_lattice,
-            scaled_positions=std_positions[perm],
-            symbols=[atom_data[n][1] for n in std_types[perm]],
-            magnetic_moments=dataset.std_tensors[perm],
-        )

test/velph/cli/init/test_cmd_init.py

@@ -9,6 +9,7 @@
 from collections.abc import Callable
 from typing import Literal
 
+import click
 import numpy as np
 import pytest
 import tomli
@@ -118,6 +119,78 @@ def test_run_init_read_cell_and_magmom():
     np.testing.assert_allclose(pcell.magnetic_moments, [1, -1])
 
 
+def _site_mixture_cell() -> PhonopyAtoms:
+    """Return a CsCl-like cell with a co-located Ge/Sn site and a Te site."""
+    return PhonopyAtoms(
+        symbols=["Ge", "Sn", "Te"],
+        cell=np.eye(3) * 4.0,
+        scaled_positions=[[0, 0, 0], [0, 0, 0], [0.5, 0.5, 0.5]],
+    )
+
+
+def _run_init_site_mixture(options: VelphInitOptions) -> dict:
+    toml_lines = _run_init(
+        _site_mixture_cell(), options, velph_template_fp=io.BytesIO(b"")
+    )
+    assert toml_lines is not None
+    return tomli.loads("\n".join(toml_lines))
+
+
+@pytest.mark.parametrize("symmetrize_cell", [False, True])
+def test_run_init_site_mixture(symmetrize_cell: bool):
+    """Test --site-mixture writes per-atom weights into velph.toml.
+
+    Weights survive both the default (find_primitive) path and the
+    --symmetrize-cell standardization path, and a pure site carries an
+    explicit weight of 1.0.
+
+    """
+    velph_dict = _run_init_site_mixture(
+        VelphInitOptions(
+            site_mixture="0.5 0.5 1.0",
+            split_site_mixture=True,
+            symmetrize_cell=symmetrize_cell,
+            supercell_dimension=(2, 2, 2),
+        )
+    )
+    for cell_key in ("unitcell", "primitive_cell"):
+        cell = load_phonopy_yaml(velph_dict[cell_key]).unitcell
+        assert cell is not None
+        assert cell.symbols == ["Ge", "Sn", "Te"]
+        assert cell.mixture_weights is not None
+        np.testing.assert_allclose(cell.mixture_weights, [0.5, 0.5, 1.0])
+
+
+def test_run_init_site_mixture_via_template():
+    """Test site_mixture / split_site_mixture flow through [init.options]."""
+    template = (
+        b'[init.options]\nsite_mixture = "0.5 0.5 1.0"\nsplit_site_mixture = true\n'
+    )
+    toml_lines = _run_init(
+        _site_mixture_cell(),
+        VelphInitOptions(supercell_dimension=(2, 2, 2)),
+        velph_template_fp=io.BytesIO(template),
+    )
+    assert toml_lines is not None
+    velph_dict = tomli.loads("\n".join(toml_lines))
+    unitcell = load_phonopy_yaml(velph_dict["unitcell"]).unitcell
+    assert unitcell is not None
+    assert unitcell.mixture_weights is not None
+    np.testing.assert_allclose(unitcell.mixture_weights, [0.5, 0.5, 1.0])
+
+
+def test_run_init_site_mixture_with_magmom_raises():
+    """Test --site-mixture cannot be combined with --magmom."""
+    with pytest.raises(click.ClickException):
+        _run_init_site_mixture(
+            VelphInitOptions(
+                site_mixture="0.5 0.5 1.0",
+                magmom="1 1 1",
+                supercell_dimension=(2, 2, 2),
+            )
+        )
+
+
 def test_run_init_without_max_num_atoms(
     nacl_cell: PhonopyAtoms,
 ):