diff --git a/src/phelel/velph/cli/init/cmd_init.py b/src/phelel/velph/cli/init/cmd_init.py index 75412cd..0e3ba43 100644 --- a/src/phelel/velph/cli/init/cmd_init.py +++ b/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, diff --git a/src/phelel/velph/cli/init/init.py b/src/phelel/velph/cli/init/init.py index 6af8102..e453e9a 100644 --- a/src/phelel/velph/cli/init/init.py +++ b/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 diff --git a/src/phelel/velph/cli/utils.py b/src/phelel/velph/cli/utils.py index d824b93..0c1f49b 100644 --- a/src/phelel/velph/cli/utils.py +++ b/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 diff --git a/src/phelel/velph/utils/structure.py b/src/phelel/velph/utils/structure.py index 50ad20b..60dd64f 100644 --- a/src/phelel/velph/utils/structure.py +++ b/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], - ) diff --git a/test/velph/cli/init/test_cmd_init.py b/test/velph/cli/init/test_cmd_init.py index e5186c1..bc8b91d 100644 --- a/test/velph/cli/init/test_cmd_init.py +++ b/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, ):