diff --git a/pyproject.toml b/pyproject.toml index 8a0b3db333..dcd81afd9b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -158,6 +158,20 @@ docs = [ "jupyterlab==4.6.1", ] +# The strict forcefield groups have mutually incompatible pins (e.g. numpy>=2 +# via mattersim vs numpy==1.26.4). They are meant to be installed one at a time. +# uv builds a single universal lockfile, so these must be declared as conflicting +# extras or `uv sync` fails to resolve them together. +[tool.uv] +conflicts = [ + [ + { extra = "strict-forcefields-generic" }, + { extra = "strict-forcefields-torch-limited" }, + { extra = "strict-forcefields-e3nn-limited" }, + { extra = "strict-forcefields-numpy-limited" }, + ], +] + [tool.setuptools.package-data] atomate2 = ["py.typed"] "atomate2.vasp.sets" = ["*.yaml"] @@ -204,6 +218,9 @@ exclude_lines = [ [tool.ruff] target-version = "py310" output-format = "concise" +# vampire_caller.py is vendored verbatim from pymatgen (removed there in 2026.3.23); +# exclude it from lint/format so it stays a faithful copy of the upstream source. +extend-exclude = ["src/atomate2/vampire/vampire_caller.py"] [tool.ruff.lint] select = ["ALL"] diff --git a/src/atomate2/common/flows/exchange.py b/src/atomate2/common/flows/exchange.py new file mode 100644 index 0000000000..3d08b6a26e --- /dev/null +++ b/src/atomate2/common/flows/exchange.py @@ -0,0 +1,164 @@ +"""Flow for fitting magnetic exchange parameters and estimating Tc. + +This is a post-processing workflow and runs no DFT itself. Given magnetic +structures and their energies (e.g. the output of the magnetic-orderings +workflow), it fits a classical Heisenberg Hamiltonian and optionally runs +Vampire Monte-Carlo for the critical temperature. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING + +from jobflow import Flow, Maker + +from atomate2.common.jobs.exchange import build_exchange_doc, heisenberg_mapping +from atomate2.vampire.jobs.run_vampire import run_vampire + +if TYPE_CHECKING: + from pymatgen.core.structure import Structure + + from atomate2.common.schemas.magnetism import MagneticOrderingsDocument + + +__all__ = ["ExchangeMaker"] + + +@dataclass +class ExchangeMaker(Maker): + """Maker to fit Heisenberg exchange parameters from magnetic structures. + + Given a set of magnetic structures and their energies (per atom), this fits a + classical Heisenberg Hamiltonian via pymatgen's ``HeisenbergMapper`` to extract + shell-resolved exchange parameters ``J_ij``, fit by least squares over all + supplied orderings. + + This is a post-processing workflow and runs no DFT. It is code-agnostic: the + inputs are plain pymatgen structures and energies, typically taken from the + magnetic-orderings workflow (see :meth:`make_from_ordering_doc`). + + Originally implemented in atomate (v1) for VASP as the ``ExchangeWF``. + + Parameters + ---------- + name : str + Name of the flows produced by this Maker. + heisenberg_settings : dict + Keyword arguments for pymatgen's HeisenbergMapper, i.e. the nearest-neighbour + ``cutoff`` (Angstrom) and the distance-grouping ``tol``. + run_vampire : bool + Whether to run the Vampire Monte-Carlo step to estimate the critical + temperature. Requires a VAMPIRE binary on PATH - ``vampire-serial`` unless + ``ATOMATE2_VAMPIRE_CMD`` names another (e.g. ``srun -n 24 vampire-parallel`` + for the MPI build); the step raises a clear error if it is missing. + Defaults to True (atomate1 parity). + mc_settings : dict | None + Keyword arguments for the Vampire Monte-Carlo run (e.g. ``mc_box_size``, + ``equil_timesteps``, ``mc_timesteps``). Only used if ``run_vampire``. + """ + + name: str = "exchange" + heisenberg_settings: dict = field( + default_factory=lambda: {"cutoff": 3.0, "tol": 0.04} + ) + run_vampire: bool = True + mc_settings: dict | None = None + + def make( + self, + structures: list[Structure], + energies: list[float], + parent: Structure | None = None, + ) -> Flow: + """Make a flow to fit Heisenberg exchange parameters. + + Parameters + ---------- + structures : list[Structure] + Magnetic structures, each carrying a "magmom" site property. + energies : list[float] + Energies **per atom** (eV) corresponding to each structure. + + Returns + ------- + Flow + The exchange-parameter fitting workflow. + """ + if len(structures) != len(energies): + raise ValueError( + f"Got {len(structures)} structures but {len(energies)} energies; " + "these must be equal." + ) + + # sort so the ground state (lowest energy) is index 0 + order = sorted(range(len(energies)), key=lambda i: energies[i]) + structures = [structures[i] for i in order] + energies = [energies[i] for i in order] + + # HeisenbergMapper requires a 'magmom' site property on every structure + for idx, structure in enumerate(structures): + if not structure.site_properties.get("magmom"): + raise ValueError( + f"Structure {idx} is missing a 'magmom' site property, which is " + "required to fit a Heisenberg model." + ) + + hmap = heisenberg_mapping( + structures, energies, parent, self.heisenberg_settings + ) + jobs = [hmap] + + vampire_output = None + if self.run_vampire: + vmc = run_vampire(hmap.output, self.mc_settings) + jobs.append(vmc) + vampire_output = vmc.output + + doc = build_exchange_doc( + hmap.output, + parent_structure=parent, + vampire_output=vampire_output, + vampire_settings=self.mc_settings if self.run_vampire else None, + ) + jobs.append(doc) + + formula = structures[0].composition.reduced_formula + return Flow( + jobs=jobs, + output=doc.output, + name=f"{self.name} ({formula})", + ) + + def make_from_ordering_doc(self, doc: MagneticOrderingsDocument) -> Flow: + """Make an exchange flow from a (concrete) magnetic-orderings document. + + This is a convenience constructor for chaining onto the magnetic-orderings + workflow: it pulls the structures and per-atom energies out of an + already-computed ``MagneticOrderingsDocument`` and forwards them to + :meth:`make`. Because the structures must be inspected (sorted, validated) + when the flow is built, ``doc`` must be a resolved document, not a jobflow + output reference. + + Parameters + ---------- + doc : MagneticOrderingsDocument + A computed magnetic-orderings document (e.g. from a finished + MagneticOrderingsMaker run). + + Returns + ------- + Flow + The exchange-parameter fitting workflow. + """ + structures, energies = [], [] + + for output in doc.outputs: + structure = output.structure.copy() + # HeisenbergMapper needs magmoms; the output stores them separately + if not structure.site_properties.get("magmom"): + structure.add_site_property("magmom", output.magmoms) + structures.append(structure) + energies.append(output.energy_per_atom) + + return self.make(structures, energies, doc.parent_structure) diff --git a/src/atomate2/common/jobs/exchange.py b/src/atomate2/common/jobs/exchange.py new file mode 100644 index 0000000000..c7070cd48a --- /dev/null +++ b/src/atomate2/common/jobs/exchange.py @@ -0,0 +1,95 @@ +"""Jobs for fitting Heisenberg exchange parameters and building exchange docs. + +These replace the atomate1 ``HeisenbergModelMapping``, ``HeisenbergModelToDb`` and +``VampireToDb`` firetasks. Outputs flow between jobs through jobflow references +rather than a MongoDB collection. +""" + +from __future__ import annotations + +import logging +from typing import TYPE_CHECKING + +from jobflow import job +from pymatgen.analysis.magnetism.heisenberg import HeisenbergMapper + +from atomate2.common.schemas.exchange import ExchangeDocument + +if TYPE_CHECKING: + from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel + from pymatgen.core.structure import Structure + + from atomate2.vampire.schemas.vampire_output import VampireOutput + +logger = logging.getLogger(__name__) + + +@job(name="heisenberg mapping") +def heisenberg_mapping( + structures: list[Structure], + energies: list[float], + parent: Structure | None = None, + heisenberg_settings: dict | None = None, +) -> HeisenbergModel: + """Fit a classical Heisenberg Hamiltonian to magnetic structures and energies. + + This wraps pymatgen's ``HeisenbergMapper`` to extract shell-resolved exchange + parameters ``J_ij``, fit by least squares over all supplied orderings. + + Parameters + ---------- + structures : list[Structure] + Magnetic structures, each carrying a "magmom" site property. These should be + ordered with the ground state first (index 0). + energies : list[float] + Energies **per atom** (eV) corresponding to each structure. These are + converted to total energies internally, as required by HeisenbergMapper. + heisenberg_settings : dict or None + Keyword arguments for HeisenbergMapper, e.g. ``{"cutoff": 3.0, "tol": 0.04}``. + + Returns + ------- + HeisenbergModel + The fitted, MSONable Heisenberg model. + """ + heisenberg_settings = heisenberg_settings or {} + total_energies = [e * len(s) for s, e in zip(structures, energies, strict=True)] + hmapper = HeisenbergMapper( + structures, total_energies, parent, **heisenberg_settings + ) + return hmapper.get_heisenberg_model() + + +@job(name="build exchange doc") +def build_exchange_doc( + heisenberg_model: HeisenbergModel, + parent_structure: Structure | None = None, + vampire_output: VampireOutput | None = None, + vampire_settings: dict | None = None, +) -> ExchangeDocument: + """Assemble the final ExchangeDocument from a fitted model and optional Tc run. + + Parameters + ---------- + heisenberg_model : HeisenbergModel + The fitted Heisenberg model from :func:`heisenberg_mapping`. + parent_structure : Structure or None + The full parent structure from which the magnetic structures were derived. + This is used to store the final fitted exchange parameters in the context + of the original structure. + vampire_output : VampireOutput or None + The Vampire Monte-Carlo result, if the critical-temperature step was run. + vampire_settings : dict or None + The keyword arguments :func:`run_vampire` was called with, if it was run. + + Returns + ------- + ExchangeDocument + The final summary document. + """ + return ExchangeDocument.from_model( + heisenberg_model, + parent_structure=parent_structure, + vampire_output=vampire_output, + vampire_settings=vampire_settings, + ) diff --git a/src/atomate2/common/schemas/exchange.py b/src/atomate2/common/schemas/exchange.py new file mode 100644 index 0000000000..407ae5b181 --- /dev/null +++ b/src/atomate2/common/schemas/exchange.py @@ -0,0 +1,114 @@ +"""Schemas for magnetic exchange (Heisenberg + Vampire) calculations.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from pydantic import BaseModel, Field +from pymatgen.core.structure import Structure + +if TYPE_CHECKING: + from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel + + from atomate2.vampire.schemas.vampire_output import VampireOutput + + +class ExchangeDocument(BaseModel): + """Final document with fitted Heisenberg exchange parameters and Tc. + + This is the output of the ExchangeMaker workflow. The Heisenberg fields are + always populated; the Vampire fields are populated only if the Monte-Carlo + critical-temperature step was run. + """ + + formula: str | None = Field( + None, + description="Formula taken from pymatgen.core.structure.Structure.formula.", + ) + formula_pretty: str | None = Field( + None, + description="Cleaned representation of the formula.", + ) + parent_structure: Structure | None = Field( + None, + description="The ground-state (lowest-energy) structure used for the fit.", + ) + heisenberg_settings: dict | None = Field( + None, + description="The {cutoff, tol} settings used by the HeisenbergMapper.", + ) + vampire_settings: dict | None = Field( + None, + description="Keyword arguments passed to the Vampire Monte-Carlo run (e.g. " + "mc_box_size, equil_timesteps, mc_timesteps), if run. Unset keys fall " + "back to the VampireCaller defaults.", + ) + residual: float | None = Field( + None, + description="Root-mean-square residual of the least-squares fit that " + "produced ex_params, in meV per magnetic ion. Intensive in both the cell " + "size and the number of orderings, so it is comparable between materials.", + ) + ex_params: dict | None = Field( + None, + description="Fitted exchange parameters keyed by interaction label. The " + "J_ij are in meV/muB^2 (they multiply the raw moments); the included 'E0' " + "offset is in eV per magnetic ion.", + ) + ex_mat: dict | None = Field( + None, + description="Heisenberg Hamiltonian matrix used for the Heisenberg model fit.", + ) + heisenberg_model: dict | None = Field( + None, description="Full HeisenbergModel as a serialized dict (as_dict())." + ) + critical_temp: float | None = Field( + None, + description="Critical (Curie/Neel) temperature in Kelvin from Vampire, if run.", + ) + vampire_output: dict | None = Field( + None, description="Full VampireOutput as a serialized dict (as_dict()), if run." + ) + + @classmethod + def from_model( + cls, + heisenberg_model: HeisenbergModel, + parent_structure: Structure | None = None, + vampire_output: VampireOutput | None = None, + vampire_settings: dict | None = None, + ) -> ExchangeDocument: + """Construct an ExchangeDocument from a fitted model and optional Vampire run. + + Parameters + ---------- + heisenberg_model : HeisenbergModel + The fitted Heisenberg model from pymatgen's HeisenbergMapper. + parent_structure : Structure or None + The full ground-state structure (used for the parent_structure/formula + fields). If None, falls back to ``heisenberg_model.structures[0]``, the + full cell of the ground-state ordering (all ions retained). + vampire_output : VampireOutput or None + The Vampire Monte-Carlo result, if the critical-temperature step was run. + vampire_settings : dict or None + The keyword arguments the Vampire Monte-Carlo run was called with, if run. + """ + if parent_structure is None: + parent_structure = heisenberg_model.structures[0] + + return cls( + formula=parent_structure.formula, + formula_pretty=parent_structure.composition.reduced_formula, + parent_structure=parent_structure, + heisenberg_settings={ + "cutoff": heisenberg_model.cutoff, + "tol": heisenberg_model.tol, + }, + vampire_settings=vampire_settings if vampire_output else None, + residual=heisenberg_model.residual, + ex_params=heisenberg_model.ex_params, + ex_mat=heisenberg_model.ex_mat.to_dict(), + heisenberg_model=heisenberg_model.as_dict(), + critical_temp=vampire_output.critical_temp if vampire_output else None, + vampire_output=vampire_output.as_dict() if vampire_output else None, + ) diff --git a/src/atomate2/vampire/__init__.py b/src/atomate2/vampire/__init__.py new file mode 100644 index 0000000000..4337c6e4bc --- /dev/null +++ b/src/atomate2/vampire/__init__.py @@ -0,0 +1,13 @@ +"""Interface to the external VAMPIRE atomistic spin-dynamics code. + +This subpackage vendors pymatgen's ``VampireCaller``/``VampireOutput`` (removed +from pymatgen in 2026.3.23) so atomate2's magnetic-exchange workflow can still +estimate critical temperatures via Vampire Monte-Carlo. See +:mod:`atomate2.vampire.vampire_caller` for provenance details. +""" + +from __future__ import annotations + +from atomate2.vampire.vampire_caller import VampireCaller, VampireOutput + +__all__ = ["VampireCaller", "VampireOutput"] diff --git a/src/atomate2/vampire/jobs/__init__.py b/src/atomate2/vampire/jobs/__init__.py new file mode 100644 index 0000000000..87cc72ea82 --- /dev/null +++ b/src/atomate2/vampire/jobs/__init__.py @@ -0,0 +1 @@ +"""Jobs for running the external VAMPIRE code.""" diff --git a/src/atomate2/vampire/jobs/run_vampire.py b/src/atomate2/vampire/jobs/run_vampire.py new file mode 100644 index 0000000000..08bc8dec73 --- /dev/null +++ b/src/atomate2/vampire/jobs/run_vampire.py @@ -0,0 +1,44 @@ +"""Job for running Vampire Monte-Carlo on a fitted Heisenberg model.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from jobflow import job + +from atomate2.vampire.vampire_caller import VampireCaller + +if TYPE_CHECKING: + from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel + + from atomate2.vampire.schemas.vampire_output import VampireOutput + + +@job(name="run vampire") +def run_vampire( + heisenberg_model: HeisenbergModel, + mc_settings: dict | None = None, +) -> VampireOutput: + """Run Vampire Monte-Carlo to estimate the critical temperature. + + This wraps the (vendored) ``VampireCaller``, which shells out to the external + VAMPIRE binary named by ``ATOMATE2_VAMPIRE_CMD`` (``vampire-serial`` by + default; set e.g. ``srun -n 24 vampire-parallel`` to run the MPI build). A + clear error is raised if the binary is not found on PATH. + + Parameters + ---------- + heisenberg_model : HeisenbergModel + The fitted Heisenberg model from :func:`heisenberg_mapping`. + mc_settings : dict or None + Keyword arguments for VampireCaller, e.g. ``mc_box_size``, + ``equil_timesteps``, ``mc_timesteps``. + + Returns + ------- + VampireOutput + The Vampire Monte-Carlo result, exposing ``critical_temp``. + """ + mc_settings = mc_settings or {} + vampire_caller = VampireCaller(hm=heisenberg_model, **mc_settings) + return vampire_caller.output diff --git a/src/atomate2/vampire/schemas/__init__.py b/src/atomate2/vampire/schemas/__init__.py new file mode 100644 index 0000000000..04b1d93012 --- /dev/null +++ b/src/atomate2/vampire/schemas/__init__.py @@ -0,0 +1 @@ +"""Schemas for VAMPIRE outputs.""" diff --git a/src/atomate2/vampire/schemas/vampire_output.py b/src/atomate2/vampire/schemas/vampire_output.py new file mode 100644 index 0000000000..afa9fdeede --- /dev/null +++ b/src/atomate2/vampire/schemas/vampire_output.py @@ -0,0 +1,34 @@ +"""Schema for the output of a Vampire Monte-Carlo run. + +Vendored from pymatgen's ``pymatgen.command_line.vampire_caller`` (removed there +in 2026.3.23); see :mod:`atomate2.vampire.vampire_caller` for provenance details. +""" + +from __future__ import annotations + +from monty.json import MSONable + + +class VampireOutput(MSONable): + """Process results from a Vampire Monte Carlo simulation. + + Parses the critical temperature from the simulation output. + """ + + def __init__( + self, + parsed_out: str | None = None, + nmats: int | None = None, + critical_temp: float | None = None, + ) -> None: + """Initialize the output. + + Args: + parsed_out (str): JSON rep of parsed stdout DataFrame. + nmats (int): Number of distinct materials (1 for each specie and + up/down spin). + critical_temp (float): Monte Carlo Tc result. + """ + self.parsed_out = parsed_out + self.nmats = nmats + self.critical_temp = critical_temp diff --git a/src/atomate2/vampire/vampire_caller.py b/src/atomate2/vampire/vampire_caller.py new file mode 100644 index 0000000000..9b2f5c0da2 --- /dev/null +++ b/src/atomate2/vampire/vampire_caller.py @@ -0,0 +1,397 @@ +"""This module implements an interface to the VAMPIRE code for atomistic +simulations of magnetic materials. + +.. note:: + + This module is **vendored** from pymatgen. ``VampireCaller`` and + ``VampireOutput`` lived in ``pymatgen.command_line.vampire_caller`` until + they were removed in the "Major reorganization of pymatgen repo" + (materialsproject/pymatgen#4595, merged 2026-03-02, first released in + pymatgen 2026.3.23). atomate2's + :class:`~atomate2.common.flows.exchange.ExchangeMaker` still needs them, so + this started as a copy of the file's final pre-removal state (pymatgen + commit ``8785afd0d801``, the last commit to touch it, released as + pymatgen 2025.10.7) and has since been adapted to the reworked + ``HeisenbergModel`` API: per-ordering ``sublattice_ids`` and the ``igraph`` + interaction graph replace the old ``unique_site_ids`` dict and + ``_get_j_exc`` lookup, and the ground state cell is taken from + ``magnetic_structures[0]`` rather than ``structures[0]``, which now retains + the non-magnetic ions. The original ``avg`` option is gone with the ```` + estimate it read (``HeisenbergModel.javg``): the model now fits + shell-resolved ``J_ij`` over every ordering, and ``igraph`` already carries + them per bond in VAMPIRE's normalized-spin meV convention. + The original author is Nathan C. Frey (``ncfrey``). + It is excluded from ruff (see ``[tool.ruff] extend-exclude`` in + ``pyproject.toml``). + +This module depends on a compiled vampire executable available in the PATH. +Please download at https://vampire.york.ac.uk/download/ and +follow the instructions to compile the executable. + +If you use this module, please cite: + +"Atomistic spin model simulations of magnetic nanomaterials." +R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis +and R. W. Chantrell. J. Phys.: Condens. Matter 26, 103202 (2014) +""" + +from __future__ import annotations + +import logging +import os +import shlex +import subprocess +from shutil import which + +import pandas as pd +from monty.dev import requires +from monty.json import MSONable + +from atomate2.vampire.schemas.vampire_output import VampireOutput + +__author__ = "Luguza, ncfrey" +__version__ = "0.2" +__maintainer__ = "Luca Frey, Nathan C. Frey" +__email__ = "luca.frey@student.kit.edu, ncfrey@lbl.gov" +__status__ = "Development" +__date__ = "August 2026" + +logger = logging.getLogger(__name__) + +# Command used to launch VAMPIRE, overridable per worker like ATOMATE2_VASP_CMD. +# Defaults to the serial binary; set it to the MPI build behind a launcher to run +# in parallel, e.g. ATOMATE2_VAMPIRE_CMD="srun -n 24 vampire-parallel". +VAMP_CMD = shlex.split(os.environ.get("ATOMATE2_VAMPIRE_CMD", "vampire-serial")) +VAMP_EXE = which(VAMP_CMD[0]) if VAMP_CMD else None + + +class VampireCaller: + """Run Vampire on a material with magnetic ordering and exchange parameter + information to compute the critical temperature with classical Monte Carlo. + + Attributes: + structure (Structure): Ground state structure, magnetic ions only. + Taken from ``HeisenbergModel.magnetic_structures[0]``; + ``structures[0]`` is a different cell that keeps the + non-magnetic ions. + sublattice_ids (list[int]): Parent sublattice id of each site in the + ground state structure. + igraph (StructureGraph): Ground state graph whose edge weights are the + per-bond J_ij in meV, already in the normalized-spin convention + ``E = -sum_ J_ij e_i.e_j`` that VAMPIRE uses (the fit's + meV/muB^2 parameters times the two moments). + mat_name (str): Formula unit label for input files + mat_id_dict (dict): Maps sites to material id # for vampire + indexing. + """ + + @requires( + VAMP_EXE is not None, + f"VampireCaller requires {VAMP_CMD[0] if VAMP_CMD else 'vampire-serial'} " + "(from ATOMATE2_VAMPIRE_CMD) to be in the path." + "Please follow the instructions at https://vampire.york.ac.uk/download/.", + ) + def __init__( + self, + mc_box_size=4.0, + equil_timesteps=2000, + mc_timesteps=4000, + save_inputs=False, + hm=None, + user_input_settings=None, + ): + """user_input_settings is a dictionary that can contain: + * start_t (int): Start MC sim at this temp, defaults to 0 K. + * end_t (int): End MC sim at this temp, defaults to 1500 K. + * temp_increment (int): Temp step size, defaults to 25 K. + + Args: + mc_box_size (float): x=y=z dimensions (nm) of MC simulation box + equil_timesteps (int): number of MC steps for equilibrating + mc_timesteps (int): number of MC steps for averaging + save_inputs (bool): if True, save scratch dir of vampire input files + hm (HeisenbergModel): object already fit to low energy + magnetic orderings. + user_input_settings (dict): optional commands for VAMPIRE Monte Carlo + + Todo: + * Create input files in a temp folder that gets cleaned up after run terminates + """ + self.mc_box_size = mc_box_size + self.equil_timesteps = equil_timesteps + self.mc_timesteps = mc_timesteps + self.save_inputs = save_inputs + + if not user_input_settings: # set to empty dict + self.user_input_settings = {} + else: + self.user_input_settings = user_input_settings + + # Attributes from HeisenbergModel + if hm is None: + raise ValueError("A fitted HeisenbergModel (hm=...) is required.") + self.structure = hm.magnetic_structures[0] # ground state (magnetic ions only) + self.sublattice_ids = hm.sublattice_ids[0] # site -> parent sublattice id + self.igraph = hm.igraph # ground state graph, J_ij edge weights in meV + + # The [0] above assumes igraph was built for ordering 0, which the model + # does not record (it is get_interaction_graph's default ordering_index). + # _create_mat and _create_ucf both pair graph node indices with sites of + # self.structure, so a mismatch would silently write a broken .ucf. + if not (len(self.structure) == len(self.igraph.structure) == len(self.sublattice_ids)): + raise ValueError("HeisenbergModel igraph, magnetic structure and sublattice ids are misaligned.") + + # Full structure name before reducing to only magnetic ions + self.mat_name = hm.formula + + # Switch to scratch dir which automatically cleans up vampire inputs files unless user specifies to save them + # with ScratchDir( + # "/scratch", copy_from_current_on_enter=self.save_inputs, copy_to_current_on_exit=self.save_inputs + # ): + + # Create input files + self._create_mat() + self._create_input() + self._create_ucf() + + # Call Vampire + with subprocess.Popen([VAMP_EXE, *VAMP_CMD[1:]], stdout=subprocess.PIPE, stderr=subprocess.PIPE) as process: + _stdout, stderr = process.communicate() + stdout: str = _stdout.decode() + + if stderr: + van_helsing = stderr.decode() + if len(van_helsing) > 27: # Suppress blank warning msg + logger.warning(van_helsing) + + if process.returncode != 0: + raise RuntimeError(f"Vampire exited with return code {process.returncode}.") + + self._stdout = stdout + self._stderr = stderr + + # Process output + n_mats = max(self.mat_id_dict.values()) + parsed_out, critical_temp = VampireCaller.parse_stdout("output", n_mats) + self.output = VampireOutput(parsed_out, n_mats, critical_temp) + + def _create_mat(self): + structure = self.structure + mat_name = self.mat_name + magmoms = structure.site_properties["magmom"] + + # A vampire material is a (sublattice, spin direction) group: one mat + # per sublattice, two if it hosts both spin-up and spin-down sites. + mat_ids = {} # (sublattice id, spin sign) -> material id (1-indexed) + mat_id_dict = {} # site -> material id, for vampire inputs + for site, (sub_id, magmom) in enumerate(zip(self.sublattice_ids, magmoms, strict=True)): + group = (sub_id, magmom > 0) + mat_ids.setdefault(group, len(mat_ids) + 1) + mat_id_dict[site] = mat_ids[group] + + n_mats = len(mat_ids) + mat_file = [f"material:num-materials={n_mats}"] + + # One representative site per material for the element and moment + reps = {} + for site, mat_id in mat_id_dict.items(): + reps.setdefault(mat_id, site) + + for mat_id, site in sorted(reps.items()): + atom = structure[site].species.reduced_formula + spin = 1 if magmoms[site] > 0 else -1 + + mat_file += [f"material[{mat_id}]:material-element={atom}"] + mat_file += [ + # Bind this material to its own sublattice in the unit-cell file. + # Without it every material keeps Vampire's default + # unit-cell-category of 0, so during system generation only + # category 0 matches and *all* atoms collapse into material 1 - + # the magnetic sublattices vanish and the system never orders. + # _create_ucf writes the atom material column 0-indexed + # (mat_id - 1); Vampire stores unit-cell-category internally as + # (value - 1), so passing the 1-indexed mat_id here reproduces + # the matching 0-indexed category. + f"material[{mat_id}]:unit-cell-category={mat_id}", + f"material[{mat_id}]:damping-constant=1.0", + f"material[{mat_id}]:uniaxial-anisotropy-constant=1.0e-24", + # Only positive magmoms allowed + f"material[{mat_id}]:atomic-spin-moment={abs(magmoms[site]):.2f} !muB", + f"material[{mat_id}]:initial-spin-direction=0,0,{spin}", + ] + + mat_file = "\n".join(mat_file) + mat_file_name = f"{mat_name}.mat" + + self.mat_id_dict = mat_id_dict + + with open(mat_file_name, mode="w", encoding="utf-8") as file: + file.write(mat_file) + + def _create_input(self): + structure = self.structure + mc_box_size = self.mc_box_size + equil_timesteps = self.equil_timesteps + mc_timesteps = self.mc_timesteps + mat_name = self.mat_name + + input_script = [f"material:unit-cell-file={mat_name}.ucf"] + input_script += [f"material:file={mat_name}.mat"] + + # Specify periodic boundary conditions + input_script += [ + "create:periodic-boundaries-x", + "create:periodic-boundaries-y", + "create:periodic-boundaries-z", + ] + + # Unit cell size in Angstrom + abc = structure.lattice.abc + ucx, ucy, ucz = abc[0], abc[1], abc[2] + + input_script += [f"dimensions:unit-cell-size-x = {ucx:.10f} !A"] + input_script += [f"dimensions:unit-cell-size-y = {ucy:.10f} !A"] + input_script += [f"dimensions:unit-cell-size-z = {ucz:.10f} !A"] + + # System size in nm + input_script += [ + f"dimensions:system-size-x = {mc_box_size:.1f} !nm", + f"dimensions:system-size-y = {mc_box_size:.1f} !nm", + f"dimensions:system-size-z = {mc_box_size:.1f} !nm", + ] + + # Critical temperature Monte Carlo calculation + input_script += [ + "sim:integrator = monte-carlo", + "sim:program = curie-temperature", + ] + + # Default Monte Carlo params + input_script += [ + f"sim:equilibration-time-steps = {equil_timesteps}", + f"sim:loop-time-steps = {mc_timesteps}", + "sim:time-steps-increment = 1", + ] + + # Set temperature range and step size of simulation + start_t = self.user_input_settings.get("start_t", 0) + + end_t = self.user_input_settings.get("end_t", 1500) + + temp_increment = self.user_input_settings.get("temp_increment", 25) + + input_script += [ + f"sim:minimum-temperature = {start_t}", + f"sim:maximum-temperature = {end_t}", + f"sim:temperature-increment = {temp_increment}", + ] + + # Output to save + input_script += [ + "output:temperature", + "output:mean-magnetisation-length", + "output:material-mean-magnetisation-length", + "output:mean-susceptibility", + ] + + input_script = "\n".join(input_script) + + with open("input", mode="w", encoding="utf-8") as file: + file.write(input_script) + + def _create_ucf(self): + structure = self.structure + mat_name = self.mat_name + + abc = structure.lattice.abc + ucx, ucy, ucz = abc[0], abc[1], abc[2] + + ucf = ["# Unit cell size:"] + ucf += [f"{ucx:.10f} {ucy:.10f} {ucz:.10f}"] + + ucf += ["# Unit cell lattice vectors:"] + a1 = list(structure.lattice.matrix[0]) + ucf += [f"{a1[0]:.10f} {a1[1]:.10f} {a1[2]:.10f}"] + a2 = list(structure.lattice.matrix[1]) + ucf += [f"{a2[0]:.10f} {a2[1]:.10f} {a2[2]:.10f}"] + a3 = list(structure.lattice.matrix[2]) + ucf += [f"{a3[0]:.10f} {a3[1]:.10f} {a3[2]:.10f}"] + + nmats = max(self.mat_id_dict.values()) + + ucf += ["# Atoms num_materials; id cx cy cz mat cat hcat"] + ucf += [f"{len(structure)} {nmats}"] + + # Fractional coordinates of atoms + for site, r in enumerate(structure.frac_coords): + # Back to 0 indexing for some reason... + mat_id = self.mat_id_dict[site] - 1 + ucf += [f"{site} {r[0]:.10f} {r[1]:.10f} {r[2]:.10f} {mat_id} 0 0"] + + # J_ij exchange interaction matrix; the interaction graph carries every + # bond's J_ij (meV) as an edge weight, already in the normalized-spin + # convention VAMPIRE expects: get_interaction_graph folds the ground + # state's moments into the fitted meV/muB^2 parameters, so the weights go + # straight into the ucf. + igraph = self.igraph + n_inter = 0 + for idx in range(len(igraph.graph.nodes)): + n_inter += igraph.get_coordination_of_site(idx) + + ucf += ["# Interactions"] + ucf += [f"{n_inter} isotropic"] + + iid = 0 # counts number of interaction + for idx in range(len(igraph.graph.nodes)): + for conn in igraph.get_connected_sites(idx): + dx, dy, dz = conn.jimage # relative integer coordinates of atom j + j = conn.index # index of neighbor + + j_exc = conn.weight + + # Convert J_ij from meV to Joules + j_exc *= 1.6021766e-22 + + j_exc = str(j_exc) # otherwise this rounds to 0 + + ucf += [f"{iid} {idx} {j} {dx} {dy} {dz} {j_exc}"] + iid += 1 + + ucf = "\n".join(ucf) + ucf_file_name = f"{mat_name}.ucf" + + with open(ucf_file_name, mode="w", encoding="utf-8") as file: + file.write(ucf) + + @staticmethod + def parse_stdout(vamp_stdout, n_mats: int) -> tuple: + """Parse stdout from Vampire. + + Args: + vamp_stdout (txt file): Vampire 'output' file. + n_mats (int): Number of materials in Vampire simulation. + + Returns: + parsed_out (DataFrame): MSONable vampire output. + critical_temp (float): Calculated critical temp. + """ + names = [ + "T", + "m_total", + *[f"m_{idx + 1}" for idx in range(n_mats)], + "X_x", + "X_y", + "X_z", + "X_m", + "nan", + ] + + # Parsing vampire MC output + df_stdout = pd.read_csv(vamp_stdout, sep="\t", skiprows=9, header=None, names=names).drop("nan", axis=1) + + parsed_out = df_stdout.to_json() + + # Max of susceptibility <-> critical temp + critical_temp = df_stdout.iloc[df_stdout.X_m.idxmax()]["T"] + + return parsed_out, critical_temp diff --git a/tests/vampire/test_vampire_caller.py b/tests/vampire/test_vampire_caller.py new file mode 100644 index 0000000000..0141fe8e77 --- /dev/null +++ b/tests/vampire/test_vampire_caller.py @@ -0,0 +1,168 @@ +"""Tests for the vendored VampireCaller's input-file writers. + +These guard the ``HeisenbergModel`` contract the caller relies on: ``igraph``, +``magnetic_structures[0]`` and ``sublattice_ids[0]`` all share one site indexing, +the magnetic-only cell of the ground-state ordering. ``structures[0]`` is a +different, larger cell (it keeps the non-magnetic ions), so reading the cell from +there silently writes non-magnetic ions into the ``.ucf`` as atoms and +de-synchronises the interaction block's node indices. +""" + +from __future__ import annotations + +import subprocess +from typing import TYPE_CHECKING + +import pytest +from pymatgen.analysis.graphs import StructureGraph +from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel +from pymatgen.core import Lattice, Structure + +from atomate2.vampire import vampire_caller +from atomate2.vampire.vampire_caller import VampireCaller + +if TYPE_CHECKING: + from pathlib import Path + + +@pytest.fixture +def magnetic_structure() -> Structure: + """Two antiferromagnetically aligned Mn ions on one sublattice.""" + return Structure( + Lattice.cubic(4.0), + ["Mn", "Mn"], + [[0, 0, 0], [0.5, 0.5, 0.5]], + site_properties={"magmom": [3.0, -3.0]}, + ) + + +@pytest.fixture +def full_structure() -> Structure: + """The same ordering with its non-magnetic Al ions retained.""" + return Structure( + Lattice.cubic(4.0), + ["Mn", "Mn", "Al", "Al"], + [[0, 0, 0], [0.5, 0.5, 0.5], [0.5, 0, 0], [0, 0.5, 0]], + site_properties={"magmom": [3.0, -3.0, 0.0, 0.0]}, + ) + + +@pytest.fixture +def heisenberg_model(magnetic_structure, full_structure) -> HeisenbergModel: + """A minimal model carrying only what VampireCaller reads off it.""" + igraph = StructureGraph.from_empty_graph( + magnetic_structure, + edge_weight_name="exchange_constant", + edge_weight_units="meV", + ) + igraph.add_edge(0, 1, to_jimage=(0, 0, 0), weight=5.0) + + return HeisenbergModel( + formula="Mn3Al", + structures=[full_structure], + magnetic_structures=[magnetic_structure], + sublattice_ids=[[0, 0]], + igraph=igraph, + ) + + +class _FakePopen: + """Stand in for the vampire-serial subprocess.""" + + returncode = 0 + + def __init__(self, *_args, **_kwargs): + pass + + def __enter__(self): + return self + + def __exit__(self, *_exc): + pass + + def communicate(self): + return b"", b"" + + +def _run_caller(hm, monkeypatch, tmp_path: Path) -> VampireCaller: + """Run VampireCaller.__init__ without the vampire-serial binary. + + ``__init__`` is wrapped in ``monty.dev.requires``, which raises unless + vampire-serial is on PATH; ``__wrapped__`` (set by functools.wraps) reaches + the real initialiser so the writers stay testable without the binary. + """ + monkeypatch.chdir(tmp_path) + monkeypatch.setattr(subprocess, "Popen", _FakePopen) + monkeypatch.setattr( + vampire_caller.VampireCaller, + "parse_stdout", + staticmethod(lambda *_: ("{}", 0.0)), + ) + + caller = object.__new__(VampireCaller) + VampireCaller.__init__.__wrapped__(caller, hm=hm) + return caller + + +def _parse_ucf(text: str) -> tuple[list[str], int, list[str]]: + """Return (atom lines, declared material count, interaction lines).""" + lines = text.splitlines() + atoms_at = lines.index("# Atoms num_materials; id cx cy cz mat cat hcat") + n_atoms, n_mats = (int(value) for value in lines[atoms_at + 1].split()) + inter_at = lines.index("# Interactions") + n_inter = int(lines[inter_at + 1].split()[0]) + + atom_lines = lines[atoms_at + 2 : atoms_at + 2 + n_atoms] + inter_lines = lines[inter_at + 2 : inter_at + 2 + n_inter] + assert len(atom_lines) == n_atoms, "declared atom count != atom block length" + assert len(inter_lines) == n_inter, "declared interaction count != block length" + return atom_lines, n_mats, inter_lines + + +def test_ucf_holds_only_magnetic_ions(heisenberg_model, monkeypatch, tmp_path): + """The .ucf describes the magnetic-only cell, not the full ordering cell.""" + caller = _run_caller(heisenberg_model, monkeypatch, tmp_path) + atom_lines, _, _ = _parse_ucf((tmp_path / "Mn3Al.ucf").read_text()) + + assert len(atom_lines) == len(heisenberg_model.magnetic_structures[0]) == 2 + # The regression this guards: reading the cell from structures[0] would put + # the two non-magnetic Al ions in the .ucf as atoms. + assert len(atom_lines) != len(heisenberg_model.structures[0]) + assert caller.structure == heisenberg_model.magnetic_structures[0] + + +def test_ucf_indices_stay_in_range(heisenberg_model, monkeypatch, tmp_path): + """Material ids and interaction node indices address existing atoms/materials.""" + _run_caller(heisenberg_model, monkeypatch, tmp_path) + atom_lines, n_mats, inter_lines = _parse_ucf((tmp_path / "Mn3Al.ucf").read_text()) + + mat_ids = [int(line.split()[4]) for line in atom_lines] + assert set(mat_ids) == set(range(n_mats)), "material ids not a dense 0-based range" + + for line in inter_lines: + _, i, j = line.split()[:3] + assert 0 <= int(i) < len(atom_lines) + assert 0 <= int(j) < len(atom_lines) + + +def test_mat_file_matches_ucf_materials(heisenberg_model, monkeypatch, tmp_path): + """Both spin directions of the single sublattice become their own material.""" + _run_caller(heisenberg_model, monkeypatch, tmp_path) + mat_text = (tmp_path / "Mn3Al.mat").read_text() + _, n_mats, _ = _parse_ucf((tmp_path / "Mn3Al.ucf").read_text()) + + assert int(mat_text.splitlines()[0].split("=")[1]) == n_mats == 2 + # Each material must be pinned to its own unit-cell-category, or Vampire + # collapses every atom into material 1 and the sublattices vanish. + categories = {ln for ln in mat_text.splitlines() if "unit-cell-category" in ln} + assert len(categories) == n_mats + + +def test_misaligned_model_is_rejected(heisenberg_model, monkeypatch, tmp_path): + """A model whose igraph does not match its magnetic cell fails loudly.""" + # Exactly the pre-migration bug: the full cell paired with a graph built on + # the magnetic-only cell. + heisenberg_model.magnetic_structures = heisenberg_model.structures + + with pytest.raises(ValueError, match="misaligned"): + _run_caller(heisenberg_model, monkeypatch, tmp_path)