From 8142c33f5925ff84a9540626cce4afc1febd2be3 Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Mon, 1 Jun 2026 17:06:32 +0200 Subject: [PATCH 01/11] preliminary trial implementation of porting the exchange workflow from atomate1 --- pyproject.toml | 17 + src/atomate2/common/flows/exchange.py | 161 +++++++++ src/atomate2/common/jobs/exchange.py | 119 +++++++ src/atomate2/common/schemas/exchange.py | 101 ++++++ src/atomate2/vampire/__init__.py | 13 + src/atomate2/vampire/vampire_caller.py | 434 ++++++++++++++++++++++++ 6 files changed, 845 insertions(+) create mode 100644 src/atomate2/common/flows/exchange.py create mode 100644 src/atomate2/common/jobs/exchange.py create mode 100644 src/atomate2/common/schemas/exchange.py create mode 100644 src/atomate2/vampire/__init__.py create mode 100644 src/atomate2/vampire/vampire_caller.py diff --git a/pyproject.toml b/pyproject.toml index 1ebc3419b6..61e772342c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -163,6 +163,20 @@ docs = [ "jupyterlab==4.5.7", ] +# 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"] @@ -209,6 +223,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..05718afad1 --- /dev/null +++ b/src/atomate2/common/flows/exchange.py @@ -0,0 +1,161 @@ +"""Flow for fitting magnetic exchange parameters and estimating Tc. + +Ports atomate1's ``ExchangeWF`` to atomate2. This is a post-processing workflow: +it 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, + 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 + exchange parameters ``J_ij`` and the average exchange ```` (``javg``). + + 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 the external ``vampire-serial`` binary on PATH; 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``, ``avg``). 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], + ) -> 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, 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 + + # structures[0] is the full ground-state structure (Heisenberg strips + # non-magnetic atoms internally, so pass the original through for provenance) + doc = build_exchange_doc( + hmap.output, + parent_structure=structures[0], + vampire_output=vampire_output, + ) + 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) diff --git a/src/atomate2/common/jobs/exchange.py b/src/atomate2/common/jobs/exchange.py new file mode 100644 index 0000000000..c4aea5304a --- /dev/null +++ b/src/atomate2/common/jobs/exchange.py @@ -0,0 +1,119 @@ +"""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 +from atomate2.vampire import VampireCaller + +if TYPE_CHECKING: + from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel + from pymatgen.core.structure import Structure + + from atomate2.vampire import VampireOutput + + +logger = logging.getLogger(__name__) + + +@job(name="heisenberg mapping") +def heisenberg_mapping( + structures: list[Structure], + energies: list[float], + heisenberg_settings: dict | None = None, +) -> HeisenbergModel: + """Fit a classical Heisenberg Hamiltonian to magnetic structures and energies. + + This wraps pymatgen's ``HeisenbergMapper`` to extract exchange parameters + ``J_ij`` and the average exchange ```` (``javg``). + + 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, **heisenberg_settings) + return hmapper.get_heisenberg_model() + + +@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-serial`` binary. 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``, ``avg``. + + 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 + + +@job(name="build exchange doc") +def build_exchange_doc( + heisenberg_model: HeisenbergModel, + parent_structure: Structure | None = None, + vampire_output: VampireOutput | 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 ground-state structure, used to populate the document's + parent_structure/formula fields. If None, the model's (magnetic-only) + structure is used instead. + vampire_output : VampireOutput or None + The Vampire Monte-Carlo result, if the critical-temperature step was run. + + Returns + ------- + ExchangeDocument + The final summary document. + """ + return ExchangeDocument.from_model( + heisenberg_model, + parent_structure=parent_structure, + vampire_output=vampire_output, + ) diff --git a/src/atomate2/common/schemas/exchange.py b/src/atomate2/common/schemas/exchange.py new file mode 100644 index 0000000000..c0f4dbc885 --- /dev/null +++ b/src/atomate2/common/schemas/exchange.py @@ -0,0 +1,101 @@ +"""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 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.", + ) + nn_cutoff: float | None = Field( + None, description="Nearest-neighbour cutoff radius (Angstrom) used in the fit." + ) + nn_tol: float | None = Field( + None, description="Tolerance for grouping near-equal bond distances." + ) + javg: float | None = Field( + None, description="Average exchange parameter in Kelvin." + ) + ex_params: dict | None = Field( + None, description="Fitted exchange parameters J_ij keyed by interaction label." + ) + 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, + ) -> 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]``, which + is the magnetic sublattice only (non-magnetic atoms stripped). + vampire_output : VampireOutput or None + The Vampire Monte-Carlo result, if the critical-temperature step was 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, + }, + nn_cutoff=heisenberg_model.cutoff, + nn_tol=heisenberg_model.tol, + javg=heisenberg_model.javg, + ex_params=heisenberg_model.ex_params, + 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/vampire_caller.py b/src/atomate2/vampire/vampire_caller.py new file mode 100644 index 0000000000..c1891515ed --- /dev/null +++ b/src/atomate2/vampire/vampire_caller.py @@ -0,0 +1,434 @@ +"""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 is a verbatim copy of the file's final pre-removal state (pymatgen + commit ``8785afd0d801``, the last commit to touch it, released as + pymatgen 2025.10.7). The original author is Nathan C. Frey (``ncfrey``); + ``HeisenbergMapper``, on which this depends, is still shipped by pymatgen. + It is excluded from ruff (see ``[tool.ruff] extend-exclude`` in + ``pyproject.toml``) so it stays a faithful copy of the upstream source. + +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 subprocess +from shutil import which + +import pandas as pd +from monty.dev import requires +from monty.json import MSONable + +from pymatgen.analysis.magnetism.heisenberg import HeisenbergMapper + +__author__ = "ncfrey" +__version__ = "0.1" +__maintainer__ = "Nathan C. Frey" +__email__ = "ncfrey@lbl.gov" +__status__ = "Development" +__date__ = "June 2019" + +logger = logging.getLogger(__name__) + +VAMP_EXE = which("vampire-serial") + + +class VampireCaller: + """Run Vampire on a material with magnetic ordering and exchange parameter + information to compute the critical temperature with classical Monte Carlo. + + Attributes: + sgraph (StructureGraph): Ground state graph. + unique_site_ids (dict): Maps each site to its unique identifier + nn_interactions (dict): {i: j} pairs of NN interactions + between unique sites. + ex_params (dict): Exchange parameter values (meV/atom) + mft_t (float): Mean field theory estimate of critical T + 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, + "VampireCaller requires vampire-serial to be in the path." + "Please follow the instructions at https://vampire.york.ac.uk/download/.", + ) + def __init__( + self, + ordered_structures=None, + energies=None, + mc_box_size=4.0, + equil_timesteps=2000, + mc_timesteps=4000, + save_inputs=False, + hm=None, + avg=True, + 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: + ordered_structures (list): Structure objects with magmoms. + energies (list): Energies of each relaxed magnetic structure. + 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. + avg (bool): If True, simply use exchange parameter estimate. + If False, attempt to use NN, NNN, etc. interactions. + 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 + self.avg = avg + + if not user_input_settings: # set to empty dict + self.user_input_settings = {} + else: + self.user_input_settings = user_input_settings + + # Get exchange parameters and set instance variables + if not hm: + hmapper = HeisenbergMapper(ordered_structures, energies, cutoff=3.0, tol=0.02) + + hm = hmapper.get_heisenberg_model() + + # Attributes from HeisenbergModel + self.hm = hm + self.structure = hm.structures[0] # ground state + self.sgraph = hm.sgraphs[0] # ground state graph + self.unique_site_ids = hm.unique_site_ids + self.nn_interactions = hm.nn_interactions + self.dists = hm.dists + self.tol = hm.tol + self.ex_params = hm.ex_params + self.javg = hm.javg + + # 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], 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"] + + # Maps sites to material id for vampire inputs + mat_id_dict = {} + + n_mats = 0 + for key in self.unique_site_ids: + spin_up, spin_down = False, False + n_mats += 1 # at least 1 mat for each unique site + + # Check which spin sublattices exist for this site id + for site in key: + if magmoms[site] > 0: + spin_up = True + if magmoms[site] < 0: + spin_down = True + + # Assign material id for each site + for site in key: + if spin_up and not spin_down: + mat_id_dict[site] = n_mats + if spin_down and not spin_up: + mat_id_dict[site] = n_mats + if spin_up and spin_down: + # Check if spin up or down shows up first + m0 = magmoms[key[0]] + if magmoms[site] > 0 and m0 > 0: + mat_id_dict[site] = n_mats + if magmoms[site] < 0 and m0 < 0: + mat_id_dict[site] = n_mats + if magmoms[site] > 0 > m0: + mat_id_dict[site] = n_mats + 1 + if magmoms[site] < 0 < m0: + mat_id_dict[site] = n_mats + 1 + + # Increment index if two sublattices + if spin_up and spin_down: + n_mats += 1 + + mat_file = [f"material:num-materials={n_mats}"] + + for key in self.unique_site_ids: + i = self.unique_site_ids[key] # unique site id + + for site in key: + mat_id = mat_id_dict[site] + + # Only positive magmoms allowed + m_magnitude = abs(magmoms[site]) + + if magmoms[site] > 0: + spin = 1 + elif magmoms[site] < 0: + spin = -1 + else: + spin = 0 + + atom = structure[i].species.reduced_formula + + mat_file += [f"material[{mat_id}]:material-element={atom}"] + mat_file += [ + f"material[{mat_id}]:damping-constant=1.0", + f"material[{mat_id}]:uniaxial-anisotropy-constant=1.0e-24", # xx - do we need this? + f"material[{mat_id}]:atomic-spin-moment={m_magnitude:.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 + sgraph = self.sgraph + n_inter = 0 + for idx in range(len(sgraph.graph.nodes)): + n_inter += sgraph.get_coordination_of_site(idx) + + ucf += ["# Interactions"] + ucf += [f"{n_inter} isotropic"] + + iid = 0 # counts number of interaction + for idx in range(len(sgraph.graph.nodes)): + connections = sgraph.get_connected_sites(idx) + for c in connections: + jimage = c[1] # relative integer coordinates of atom j + dx = jimage[0] + dy = jimage[1] + dz = jimage[2] + j = c[2] # index of neighbor + dist = round(c[-1], 2) + + # Look up J_ij between the sites + # if case: Just use estimate + j_exc = self.hm.javg if self.avg is True else self.hm._get_j_exc(idx, j, dist) + + # 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 + + +class VampireOutput(MSONable): + """This class processes results from a Vampire Monte Carlo simulation + and parses the critical temperature. + """ + + def __init__(self, parsed_out=None, nmats=None, critical_temp=None): + """ + 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 From c731bc7e27860cae2438a9351d9be8d9f4f59472 Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Mon, 8 Jun 2026 19:21:48 +0200 Subject: [PATCH 02/11] refactored vampire specific code to fit into the atomate2 folder scheme added a ex_mat field to the ExchangeDocument --- src/atomate2/common/flows/exchange.py | 5 ++- src/atomate2/common/jobs/exchange.py | 33 ------------------- src/atomate2/common/schemas/exchange.py | 10 ++++-- src/atomate2/vampire/jobs/run_vampire.py | 32 ++++++++++++++++++ .../vampire/schemas/vampire_output.py | 17 ++++++++++ src/atomate2/vampire/vampire_caller.py | 18 +--------- 6 files changed, 59 insertions(+), 56 deletions(-) create mode 100644 src/atomate2/vampire/jobs/run_vampire.py create mode 100644 src/atomate2/vampire/schemas/vampire_output.py diff --git a/src/atomate2/common/flows/exchange.py b/src/atomate2/common/flows/exchange.py index 05718afad1..7c87aa751c 100644 --- a/src/atomate2/common/flows/exchange.py +++ b/src/atomate2/common/flows/exchange.py @@ -1,7 +1,6 @@ """Flow for fitting magnetic exchange parameters and estimating Tc. -Ports atomate1's ``ExchangeWF`` to atomate2. This is a post-processing workflow: -it runs no DFT itself. Given magnetic structures and their energies (e.g. the +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. """ @@ -16,8 +15,8 @@ from atomate2.common.jobs.exchange import ( build_exchange_doc, heisenberg_mapping, - run_vampire, ) +from atomate2.vampire.jobs.run_vampire import run_vampire if TYPE_CHECKING: from pymatgen.core.structure import Structure diff --git a/src/atomate2/common/jobs/exchange.py b/src/atomate2/common/jobs/exchange.py index c4aea5304a..feae3cc739 100644 --- a/src/atomate2/common/jobs/exchange.py +++ b/src/atomate2/common/jobs/exchange.py @@ -14,15 +14,11 @@ from pymatgen.analysis.magnetism.heisenberg import HeisenbergMapper from atomate2.common.schemas.exchange import ExchangeDocument -from atomate2.vampire import VampireCaller if TYPE_CHECKING: from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel from pymatgen.core.structure import Structure - from atomate2.vampire import VampireOutput - - logger = logging.getLogger(__name__) @@ -59,35 +55,6 @@ def heisenberg_mapping( return hmapper.get_heisenberg_model() -@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-serial`` binary. 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``, ``avg``. - - 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 - - @job(name="build exchange doc") def build_exchange_doc( heisenberg_model: HeisenbergModel, diff --git a/src/atomate2/common/schemas/exchange.py b/src/atomate2/common/schemas/exchange.py index c0f4dbc885..86a1bf6996 100644 --- a/src/atomate2/common/schemas/exchange.py +++ b/src/atomate2/common/schemas/exchange.py @@ -10,7 +10,7 @@ if TYPE_CHECKING: from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel - from atomate2.vampire import VampireOutput + from atomate2.vampire.schemas.vampire_output import VampireOutput class ExchangeDocument(BaseModel): @@ -44,10 +44,13 @@ class ExchangeDocument(BaseModel): None, description="Tolerance for grouping near-equal bond distances." ) javg: float | None = Field( - None, description="Average exchange parameter in Kelvin." + None, description="Estimated average exchange parameter in meV/atom (atom = magnetic ion) from the energy difference between the lowest energy FM and AFM orderings." ) ex_params: dict | None = Field( - None, description="Fitted exchange parameters J_ij keyed by interaction label." + None, description="Fitted exchange parameters J_ij keyed by interaction label in meV/atom (atom = 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())." @@ -95,6 +98,7 @@ def from_model( nn_tol=heisenberg_model.tol, javg=heisenberg_model.javg, 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/jobs/run_vampire.py b/src/atomate2/vampire/jobs/run_vampire.py new file mode 100644 index 0000000000..bff0028a8f --- /dev/null +++ b/src/atomate2/vampire/jobs/run_vampire.py @@ -0,0 +1,32 @@ +from atomate2.vampire.vampire_caller import VampireCaller +from atomate2.vampire.schemas.vampire_output import VampireOutput +from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel +from jobflow import job + +@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-serial`` binary. 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``, ``avg``. + + 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 \ No newline at end of file diff --git a/src/atomate2/vampire/schemas/vampire_output.py b/src/atomate2/vampire/schemas/vampire_output.py new file mode 100644 index 0000000000..fd0172eda2 --- /dev/null +++ b/src/atomate2/vampire/schemas/vampire_output.py @@ -0,0 +1,17 @@ +from monty.json import MSONable + +class VampireOutput(MSONable): + """This class processes results from a Vampire Monte Carlo simulation + and parses the critical temperature. + """ + + def __init__(self, parsed_out=None, nmats=None, critical_temp=None): + """ + 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 index c1891515ed..efa9d307bc 100644 --- a/src/atomate2/vampire/vampire_caller.py +++ b/src/atomate2/vampire/vampire_caller.py @@ -38,6 +38,7 @@ from monty.json import MSONable from pymatgen.analysis.magnetism.heisenberg import HeisenbergMapper +from atomate2.vampire.schemas.vampire_output import VampireOutput __author__ = "ncfrey" __version__ = "0.1" @@ -415,20 +416,3 @@ def parse_stdout(vamp_stdout, n_mats: int) -> tuple: critical_temp = df_stdout.iloc[df_stdout.X_m.idxmax()]["T"] return parsed_out, critical_temp - - -class VampireOutput(MSONable): - """This class processes results from a Vampire Monte Carlo simulation - and parses the critical temperature. - """ - - def __init__(self, parsed_out=None, nmats=None, critical_temp=None): - """ - 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 From df449fe03812ea1ad0125ea812adecfd699dbce2 Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Mon, 8 Jun 2026 19:23:07 +0200 Subject: [PATCH 03/11] ammend --- src/atomate2/common/jobs/exchange.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/atomate2/common/jobs/exchange.py b/src/atomate2/common/jobs/exchange.py index feae3cc739..1a2f0d74a9 100644 --- a/src/atomate2/common/jobs/exchange.py +++ b/src/atomate2/common/jobs/exchange.py @@ -18,6 +18,7 @@ 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__) From ed51af10c3b6cb5081030f56c1093db50d97dab2 Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Thu, 25 Jun 2026 16:58:21 +0200 Subject: [PATCH 04/11] Wire parent through --- src/atomate2/common/flows/exchange.py | 8 +++++--- src/atomate2/common/jobs/exchange.py | 8 ++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/atomate2/common/flows/exchange.py b/src/atomate2/common/flows/exchange.py index 7c87aa751c..d92cf3899a 100644 --- a/src/atomate2/common/flows/exchange.py +++ b/src/atomate2/common/flows/exchange.py @@ -68,6 +68,7 @@ def make( self, structures: list[Structure], energies: list[float], + parent: Structure | None = None, ) -> Flow: """Make a flow to fit Heisenberg exchange parameters. @@ -102,7 +103,7 @@ def make( "required to fit a Heisenberg model." ) - hmap = heisenberg_mapping(structures, energies, self.heisenberg_settings) + hmap = heisenberg_mapping(structures, energies, parent, self.heisenberg_settings) jobs = [hmap] vampire_output = None @@ -115,7 +116,7 @@ def make( # non-magnetic atoms internally, so pass the original through for provenance) doc = build_exchange_doc( hmap.output, - parent_structure=structures[0], + parent_structure=parent or structures[0], vampire_output=vampire_output, ) jobs.append(doc) @@ -149,6 +150,7 @@ def make_from_ordering_doc(self, doc: MagneticOrderingsDocument) -> 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 @@ -157,4 +159,4 @@ def make_from_ordering_doc(self, doc: MagneticOrderingsDocument) -> Flow: structures.append(structure) energies.append(output.energy_per_atom) - return self.make(structures, energies) + return self.make(structures, energies, doc.parent_structure) diff --git a/src/atomate2/common/jobs/exchange.py b/src/atomate2/common/jobs/exchange.py index 1a2f0d74a9..6ea9dcb476 100644 --- a/src/atomate2/common/jobs/exchange.py +++ b/src/atomate2/common/jobs/exchange.py @@ -27,6 +27,7 @@ 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. @@ -52,7 +53,7 @@ def heisenberg_mapping( """ heisenberg_settings = heisenberg_settings or {} total_energies = [e * len(s) for s, e in zip(structures, energies, strict=True)] - hmapper = HeisenbergMapper(structures, total_energies, **heisenberg_settings) + hmapper = HeisenbergMapper(structures, total_energies, parent, **heisenberg_settings) return hmapper.get_heisenberg_model() @@ -69,9 +70,8 @@ def build_exchange_doc( heisenberg_model : HeisenbergModel The fitted Heisenberg model from :func:`heisenberg_mapping`. parent_structure : Structure or None - The full ground-state structure, used to populate the document's - parent_structure/formula fields. If None, the model's (magnetic-only) - structure is used instead. + 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. From b7ba08d10f33794d91a3be29e4d187f66bd821b4 Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Mon, 13 Jul 2026 10:54:25 +0200 Subject: [PATCH 05/11] Adapt vendored VampireCaller to the reworked HeisenbergModel API The pymatgen fork's HeisenbergModel replaced the ground-state unique_site_ids dict with per-ordering site_labels and dropped _get_j_exc, which crashed run_vampire (TypeError in _create_mat). _create_mat now groups sites into materials by (sublattice, spin sign) straight from site_labels[0], and _create_ucf reads the per-bond J_ij from the igraph edge weights. Also fixes the material-element lookup (it indexed the structure with a sublattice id instead of a site index) and drops the unused HeisenbergMapper construction path. Co-Authored-By: Claude Fable 5 --- src/atomate2/vampire/vampire_caller.py | 156 +++++++++---------------- 1 file changed, 52 insertions(+), 104 deletions(-) diff --git a/src/atomate2/vampire/vampire_caller.py b/src/atomate2/vampire/vampire_caller.py index efa9d307bc..8213c2fe5a 100644 --- a/src/atomate2/vampire/vampire_caller.py +++ b/src/atomate2/vampire/vampire_caller.py @@ -9,12 +9,14 @@ (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 is a verbatim copy of the file's final pre-removal state (pymatgen + 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). The original author is Nathan C. Frey (``ncfrey``); - ``HeisenbergMapper``, on which this depends, is still shipped by pymatgen. + pymatgen 2025.10.7) and has since been adapted to the reworked + ``HeisenbergModel`` API: per-ordering ``site_labels`` and the ``igraph`` + interaction graph replace the old ``unique_site_ids`` dict and + ``_get_j_exc`` lookup. The original author is Nathan C. Frey (``ncfrey``). It is excluded from ruff (see ``[tool.ruff] extend-exclude`` in - ``pyproject.toml``) so it stays a faithful copy of the upstream source. + ``pyproject.toml``). This module depends on a compiled vampire executable available in the PATH. Please download at https://vampire.york.ac.uk/download/ and @@ -37,7 +39,6 @@ from monty.dev import requires from monty.json import MSONable -from pymatgen.analysis.magnetism.heisenberg import HeisenbergMapper from atomate2.vampire.schemas.vampire_output import VampireOutput __author__ = "ncfrey" @@ -57,12 +58,12 @@ class VampireCaller: information to compute the critical temperature with classical Monte Carlo. Attributes: - sgraph (StructureGraph): Ground state graph. - unique_site_ids (dict): Maps each site to its unique identifier - nn_interactions (dict): {i: j} pairs of NN interactions - between unique sites. - ex_params (dict): Exchange parameter values (meV/atom) - mft_t (float): Mean field theory estimate of critical T + structure (Structure): Ground state structure (magnetic ions only). + site_labels (list[int]): Parent sublattice id of each site in the + ground state structure. + igraph (StructureGraph): Ground state graph with the fitted J_ij + exchange values (meV) as edge weights. + javg (float): average exchange parameter estimate (meV). mat_name (str): Formula unit label for input files mat_id_dict (dict): Maps sites to material id # for vampire indexing. @@ -75,8 +76,6 @@ class VampireCaller: ) def __init__( self, - ordered_structures=None, - energies=None, mc_box_size=4.0, equil_timesteps=2000, mc_timesteps=4000, @@ -91,8 +90,6 @@ def __init__( * temp_increment (int): Temp step size, defaults to 25 K. Args: - ordered_structures (list): Structure objects with magmoms. - energies (list): Energies of each relaxed magnetic structure. 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 @@ -117,21 +114,12 @@ def __init__( else: self.user_input_settings = user_input_settings - # Get exchange parameters and set instance variables - if not hm: - hmapper = HeisenbergMapper(ordered_structures, energies, cutoff=3.0, tol=0.02) - - hm = hmapper.get_heisenberg_model() - # Attributes from HeisenbergModel - self.hm = hm - self.structure = hm.structures[0] # ground state - self.sgraph = hm.sgraphs[0] # ground state graph - self.unique_site_ids = hm.unique_site_ids - self.nn_interactions = hm.nn_interactions - self.dists = hm.dists - self.tol = hm.tol - self.ex_params = hm.ex_params + if hm is None: + raise ValueError("A fitted HeisenbergModel (hm=...) is required.") + self.structure = hm.structures[0] # ground state (magnetic ions only) + self.site_labels = hm.site_labels[0] # site -> parent sublattice id + self.igraph = hm.igraph # ground state graph, J_ij edge weights in meV self.javg = hm.javg # Full structure name before reducing to only magnetic ions @@ -173,70 +161,35 @@ def _create_mat(self): mat_name = self.mat_name magmoms = structure.site_properties["magmom"] - # Maps sites to material id for vampire inputs - mat_id_dict = {} - - n_mats = 0 - for key in self.unique_site_ids: - spin_up, spin_down = False, False - n_mats += 1 # at least 1 mat for each unique site - - # Check which spin sublattices exist for this site id - for site in key: - if magmoms[site] > 0: - spin_up = True - if magmoms[site] < 0: - spin_down = True - - # Assign material id for each site - for site in key: - if spin_up and not spin_down: - mat_id_dict[site] = n_mats - if spin_down and not spin_up: - mat_id_dict[site] = n_mats - if spin_up and spin_down: - # Check if spin up or down shows up first - m0 = magmoms[key[0]] - if magmoms[site] > 0 and m0 > 0: - mat_id_dict[site] = n_mats - if magmoms[site] < 0 and m0 < 0: - mat_id_dict[site] = n_mats - if magmoms[site] > 0 > m0: - mat_id_dict[site] = n_mats + 1 - if magmoms[site] < 0 < m0: - mat_id_dict[site] = n_mats + 1 - - # Increment index if two sublattices - if spin_up and spin_down: - n_mats += 1 + # 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.site_labels, 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}"] - for key in self.unique_site_ids: - i = self.unique_site_ids[key] # unique site id + # 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 site in key: - mat_id = mat_id_dict[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 += [ + f"material[{mat_id}]:damping-constant=1.0", + f"material[{mat_id}]:uniaxial-anisotropy-constant=1.0e-24", # Only positive magmoms allowed - m_magnitude = abs(magmoms[site]) - - if magmoms[site] > 0: - spin = 1 - elif magmoms[site] < 0: - spin = -1 - else: - spin = 0 - - atom = structure[i].species.reduced_formula - - mat_file += [f"material[{mat_id}]:material-element={atom}"] - mat_file += [ - f"material[{mat_id}]:damping-constant=1.0", - f"material[{mat_id}]:uniaxial-anisotropy-constant=1.0e-24", # xx - do we need this? - f"material[{mat_id}]:atomic-spin-moment={m_magnitude:.2f} !muB", - f"material[{mat_id}]:initial-spin-direction=0,0,{spin}", - ] + 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" @@ -346,29 +299,24 @@ def _create_ucf(self): 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 - sgraph = self.sgraph + # J_ij exchange interaction matrix; the interaction graph carries the + # fitted J_ij (meV) of every bond as an edge weight. + igraph = self.igraph n_inter = 0 - for idx in range(len(sgraph.graph.nodes)): - n_inter += sgraph.get_coordination_of_site(idx) + 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(sgraph.graph.nodes)): - connections = sgraph.get_connected_sites(idx) - for c in connections: - jimage = c[1] # relative integer coordinates of atom j - dx = jimage[0] - dy = jimage[1] - dz = jimage[2] - j = c[2] # index of neighbor - dist = round(c[-1], 2) - - # Look up J_ij between the sites - # if case: Just use estimate - j_exc = self.hm.javg if self.avg is True else self.hm._get_j_exc(idx, j, dist) + 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 + + # Just use the estimate, or the fitted per-bond value + j_exc = self.javg if self.avg is True else conn.weight # Convert J_ij from meV to Joules j_exc *= 1.6021766e-22 From dad3e2cc5cf6ff91c730e79a7bd24e726781f6d9 Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Thu, 23 Jul 2026 12:24:36 +0200 Subject: [PATCH 06/11] Add each materials unit cell category to .mat file --- src/atomate2/vampire/vampire_caller.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/atomate2/vampire/vampire_caller.py b/src/atomate2/vampire/vampire_caller.py index 8213c2fe5a..ea366cc2c6 100644 --- a/src/atomate2/vampire/vampire_caller.py +++ b/src/atomate2/vampire/vampire_caller.py @@ -81,7 +81,7 @@ def __init__( mc_timesteps=4000, save_inputs=False, hm=None, - avg=True, + avg=False, user_input_settings=None, ): """user_input_settings is a dictionary that can contain: @@ -184,6 +184,16 @@ def _create_mat(self): 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 From 0a9792805b4770474c8cce20eca064e7ca2fbc54 Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Mon, 3 Aug 2026 19:39:04 +0200 Subject: [PATCH 07/11] Update the Exchange Flow and VampireCaller due to the refactor in pymatgen and update ExchangeDocument to include vampire_settings; add tests for VampireCaller input-file writers. --- src/atomate2/common/flows/exchange.py | 5 +- src/atomate2/common/jobs/exchange.py | 4 + src/atomate2/common/schemas/exchange.py | 20 +-- src/atomate2/vampire/vampire_caller.py | 26 +++- tests/vampire/test_vampire_caller.py | 169 ++++++++++++++++++++++++ 5 files changed, 205 insertions(+), 19 deletions(-) create mode 100644 tests/vampire/test_vampire_caller.py diff --git a/src/atomate2/common/flows/exchange.py b/src/atomate2/common/flows/exchange.py index d92cf3899a..42f4bb76fd 100644 --- a/src/atomate2/common/flows/exchange.py +++ b/src/atomate2/common/flows/exchange.py @@ -112,12 +112,11 @@ def make( jobs.append(vmc) vampire_output = vmc.output - # structures[0] is the full ground-state structure (Heisenberg strips - # non-magnetic atoms internally, so pass the original through for provenance) doc = build_exchange_doc( hmap.output, - parent_structure=parent or structures[0], + parent_structure=parent, vampire_output=vampire_output, + vampire_settings=self.mc_settings if self.run_vampire else None, ) jobs.append(doc) diff --git a/src/atomate2/common/jobs/exchange.py b/src/atomate2/common/jobs/exchange.py index 6ea9dcb476..d498f7a53a 100644 --- a/src/atomate2/common/jobs/exchange.py +++ b/src/atomate2/common/jobs/exchange.py @@ -62,6 +62,7 @@ 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. @@ -74,6 +75,8 @@ def build_exchange_doc( 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 ------- @@ -84,4 +87,5 @@ def build_exchange_doc( 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 index 86a1bf6996..7b7962a46b 100644 --- a/src/atomate2/common/schemas/exchange.py +++ b/src/atomate2/common/schemas/exchange.py @@ -37,11 +37,11 @@ class ExchangeDocument(BaseModel): None, description="The {cutoff, tol} settings used by the HeisenbergMapper.", ) - nn_cutoff: float | None = Field( - None, description="Nearest-neighbour cutoff radius (Angstrom) used in the fit." - ) - nn_tol: float | None = Field( - None, description="Tolerance for grouping near-equal bond distances." + 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, avg), if run. Unset keys fall " + "back to the VampireCaller defaults.", ) javg: float | None = Field( None, description="Estimated average exchange parameter in meV/atom (atom = magnetic ion) from the energy difference between the lowest energy FM and AFM orderings." @@ -69,6 +69,7 @@ def from_model( 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. @@ -78,10 +79,12 @@ def from_model( 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]``, which - is the magnetic sublattice only (non-magnetic atoms stripped). + 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] @@ -94,8 +97,7 @@ def from_model( "cutoff": heisenberg_model.cutoff, "tol": heisenberg_model.tol, }, - nn_cutoff=heisenberg_model.cutoff, - nn_tol=heisenberg_model.tol, + vampire_settings=vampire_settings if vampire_output else None, javg=heisenberg_model.javg, ex_params=heisenberg_model.ex_params, ex_mat=heisenberg_model.ex_mat.to_dict(), diff --git a/src/atomate2/vampire/vampire_caller.py b/src/atomate2/vampire/vampire_caller.py index ea366cc2c6..139f06d03a 100644 --- a/src/atomate2/vampire/vampire_caller.py +++ b/src/atomate2/vampire/vampire_caller.py @@ -12,9 +12,11 @@ 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 ``site_labels`` and the ``igraph`` + ``HeisenbergModel`` API: per-ordering ``sublattice_ids`` and the ``igraph`` interaction graph replace the old ``unique_site_ids`` dict and - ``_get_j_exc`` lookup. The original author is Nathan C. Frey (``ncfrey``). + ``_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 author is Nathan C. Frey (``ncfrey``). It is excluded from ruff (see ``[tool.ruff] extend-exclude`` in ``pyproject.toml``). @@ -58,8 +60,11 @@ class VampireCaller: information to compute the critical temperature with classical Monte Carlo. Attributes: - structure (Structure): Ground state structure (magnetic ions only). - site_labels (list[int]): Parent sublattice id of each site in the + 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 with the fitted J_ij exchange values (meV) as edge weights. @@ -117,11 +122,18 @@ def __init__( # Attributes from HeisenbergModel if hm is None: raise ValueError("A fitted HeisenbergModel (hm=...) is required.") - self.structure = hm.structures[0] # ground state (magnetic ions only) - self.site_labels = hm.site_labels[0] # site -> parent sublattice id + 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 self.javg = hm.javg + # 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 @@ -165,7 +177,7 @@ def _create_mat(self): # 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.site_labels, magmoms, strict=True)): + 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] diff --git a/tests/vampire/test_vampire_caller.py b/tests/vampire/test_vampire_caller.py new file mode 100644 index 0000000000..271734db6c --- /dev/null +++ b/tests/vampire/test_vampire_caller.py @@ -0,0 +1,169 @@ +"""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, + javg=5.0, + ) + + +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) From bffa13803d30cadc5ebd07fca9969af7a7674a9c Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Tue, 4 Aug 2026 14:40:31 +0200 Subject: [PATCH 08/11] Refactor exchange parameters handling in ExchangeMaker and VampireCaller; update ExchangeDocument schema to include residuals and remove avg option. --- src/atomate2/common/flows/exchange.py | 5 +++-- src/atomate2/common/jobs/exchange.py | 4 ++-- src/atomate2/common/schemas/exchange.py | 10 ++++----- src/atomate2/vampire/jobs/run_vampire.py | 2 +- src/atomate2/vampire/vampire_caller.py | 28 +++++++++++++----------- tests/vampire/test_vampire_caller.py | 1 - 6 files changed, 26 insertions(+), 24 deletions(-) diff --git a/src/atomate2/common/flows/exchange.py b/src/atomate2/common/flows/exchange.py index 42f4bb76fd..632e51d567 100644 --- a/src/atomate2/common/flows/exchange.py +++ b/src/atomate2/common/flows/exchange.py @@ -33,7 +33,8 @@ class ExchangeMaker(Maker): Given a set of magnetic structures and their energies (per atom), this fits a classical Heisenberg Hamiltonian via pymatgen's ``HeisenbergMapper`` to extract - exchange parameters ``J_ij`` and the average exchange ```` (``javg``). + 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 @@ -54,7 +55,7 @@ class ExchangeMaker(Maker): 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``, ``avg``). Only used if ``run_vampire``. + ``equil_timesteps``, ``mc_timesteps``). Only used if ``run_vampire``. """ name: str = "exchange" diff --git a/src/atomate2/common/jobs/exchange.py b/src/atomate2/common/jobs/exchange.py index d498f7a53a..e5ff72ebc8 100644 --- a/src/atomate2/common/jobs/exchange.py +++ b/src/atomate2/common/jobs/exchange.py @@ -32,8 +32,8 @@ def heisenberg_mapping( ) -> HeisenbergModel: """Fit a classical Heisenberg Hamiltonian to magnetic structures and energies. - This wraps pymatgen's ``HeisenbergMapper`` to extract exchange parameters - ``J_ij`` and the average exchange ```` (``javg``). + This wraps pymatgen's ``HeisenbergMapper`` to extract shell-resolved exchange + parameters ``J_ij``, fit by least squares over all supplied orderings. Parameters ---------- diff --git a/src/atomate2/common/schemas/exchange.py b/src/atomate2/common/schemas/exchange.py index 7b7962a46b..843a756e58 100644 --- a/src/atomate2/common/schemas/exchange.py +++ b/src/atomate2/common/schemas/exchange.py @@ -40,14 +40,14 @@ class ExchangeDocument(BaseModel): 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, avg), if run. Unset keys fall " + "mc_box_size, equil_timesteps, mc_timesteps), if run. Unset keys fall " "back to the VampireCaller defaults.", ) - javg: float | None = Field( - None, description="Estimated average exchange parameter in meV/atom (atom = magnetic ion) from the energy difference between the lowest energy FM and AFM orderings." + residual: float | None = Field( + None, description="Sum of squared residuals of the least-squares fit that produced ex_params, in (meV per magnetic ion)^2." ) ex_params: dict | None = Field( - None, description="Fitted exchange parameters J_ij keyed by interaction label in meV/atom (atom = magnetic ion)." + 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." @@ -98,7 +98,7 @@ def from_model( "tol": heisenberg_model.tol, }, vampire_settings=vampire_settings if vampire_output else None, - javg=heisenberg_model.javg, + residual=heisenberg_model.residual, ex_params=heisenberg_model.ex_params, ex_mat=heisenberg_model.ex_mat.to_dict(), heisenberg_model=heisenberg_model.as_dict(), diff --git a/src/atomate2/vampire/jobs/run_vampire.py b/src/atomate2/vampire/jobs/run_vampire.py index bff0028a8f..a036acd131 100644 --- a/src/atomate2/vampire/jobs/run_vampire.py +++ b/src/atomate2/vampire/jobs/run_vampire.py @@ -20,7 +20,7 @@ def run_vampire( 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``, ``avg``. + ``equil_timesteps``, ``mc_timesteps``. Returns ------- diff --git a/src/atomate2/vampire/vampire_caller.py b/src/atomate2/vampire/vampire_caller.py index 139f06d03a..7851b109ea 100644 --- a/src/atomate2/vampire/vampire_caller.py +++ b/src/atomate2/vampire/vampire_caller.py @@ -16,7 +16,11 @@ 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 author is Nathan C. Frey (``ncfrey``). + 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``). @@ -66,9 +70,10 @@ class VampireCaller: non-magnetic ions. sublattice_ids (list[int]): Parent sublattice id of each site in the ground state structure. - igraph (StructureGraph): Ground state graph with the fitted J_ij - exchange values (meV) as edge weights. - javg (float): average exchange parameter estimate (meV). + 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. @@ -86,7 +91,6 @@ def __init__( mc_timesteps=4000, save_inputs=False, hm=None, - avg=False, user_input_settings=None, ): """user_input_settings is a dictionary that can contain: @@ -101,8 +105,6 @@ def __init__( save_inputs (bool): if True, save scratch dir of vampire input files hm (HeisenbergModel): object already fit to low energy magnetic orderings. - avg (bool): If True, simply use exchange parameter estimate. - If False, attempt to use NN, NNN, etc. interactions. user_input_settings (dict): optional commands for VAMPIRE Monte Carlo Todo: @@ -112,7 +114,6 @@ def __init__( self.equil_timesteps = equil_timesteps self.mc_timesteps = mc_timesteps self.save_inputs = save_inputs - self.avg = avg if not user_input_settings: # set to empty dict self.user_input_settings = {} @@ -125,7 +126,6 @@ def __init__( 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 - self.javg = hm.javg # 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). @@ -321,8 +321,11 @@ def _create_ucf(self): 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 the - # fitted J_ij (meV) of every bond as an edge weight. + # 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)): @@ -337,8 +340,7 @@ def _create_ucf(self): dx, dy, dz = conn.jimage # relative integer coordinates of atom j j = conn.index # index of neighbor - # Just use the estimate, or the fitted per-bond value - j_exc = self.javg if self.avg is True else conn.weight + j_exc = conn.weight # Convert J_ij from meV to Joules j_exc *= 1.6021766e-22 diff --git a/tests/vampire/test_vampire_caller.py b/tests/vampire/test_vampire_caller.py index 271734db6c..0141fe8e77 100644 --- a/tests/vampire/test_vampire_caller.py +++ b/tests/vampire/test_vampire_caller.py @@ -63,7 +63,6 @@ def heisenberg_model(magnetic_structure, full_structure) -> HeisenbergModel: magnetic_structures=[magnetic_structure], sublattice_ids=[[0, 0]], igraph=igraph, - javg=5.0, ) From cd9a24a2be0deeccdf884c7dd6d382c20d059220 Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Tue, 18 Aug 2026 10:00:44 +0200 Subject: [PATCH 09/11] Update residual meaning --- src/atomate2/common/schemas/exchange.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/atomate2/common/schemas/exchange.py b/src/atomate2/common/schemas/exchange.py index 843a756e58..909676b0b2 100644 --- a/src/atomate2/common/schemas/exchange.py +++ b/src/atomate2/common/schemas/exchange.py @@ -44,7 +44,10 @@ class ExchangeDocument(BaseModel): "back to the VampireCaller defaults.", ) residual: float | None = Field( - None, description="Sum of squared residuals of the least-squares fit that produced ex_params, in (meV per magnetic ion)^2." + 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." From 8458981e466db89ac8133e7881e2dd0bffc391bf Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Tue, 18 Aug 2026 10:37:20 +0200 Subject: [PATCH 10/11] Apply ruff format and check --- src/atomate2/common/flows/exchange.py | 16 +++++------ src/atomate2/common/jobs/exchange.py | 10 ++++--- src/atomate2/common/schemas/exchange.py | 8 ++++-- src/atomate2/vampire/jobs/__init__.py | 1 + src/atomate2/vampire/jobs/run_vampire.py | 19 ++++++++++--- src/atomate2/vampire/schemas/__init__.py | 1 + .../vampire/schemas/vampire_output.py | 27 +++++++++++++++---- 7 files changed, 60 insertions(+), 22 deletions(-) create mode 100644 src/atomate2/vampire/jobs/__init__.py create mode 100644 src/atomate2/vampire/schemas/__init__.py diff --git a/src/atomate2/common/flows/exchange.py b/src/atomate2/common/flows/exchange.py index 632e51d567..4ee1205c68 100644 --- a/src/atomate2/common/flows/exchange.py +++ b/src/atomate2/common/flows/exchange.py @@ -1,8 +1,9 @@ """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. +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 @@ -12,10 +13,7 @@ from jobflow import Flow, Maker -from atomate2.common.jobs.exchange import ( - build_exchange_doc, - heisenberg_mapping, -) +from atomate2.common.jobs.exchange import build_exchange_doc, heisenberg_mapping from atomate2.vampire.jobs.run_vampire import run_vampire if TYPE_CHECKING: @@ -104,7 +102,9 @@ def make( "required to fit a Heisenberg model." ) - hmap = heisenberg_mapping(structures, energies, parent, self.heisenberg_settings) + hmap = heisenberg_mapping( + structures, energies, parent, self.heisenberg_settings + ) jobs = [hmap] vampire_output = None diff --git a/src/atomate2/common/jobs/exchange.py b/src/atomate2/common/jobs/exchange.py index e5ff72ebc8..c7070cd48a 100644 --- a/src/atomate2/common/jobs/exchange.py +++ b/src/atomate2/common/jobs/exchange.py @@ -18,6 +18,7 @@ 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__) @@ -53,7 +54,9 @@ def heisenberg_mapping( """ 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) + hmapper = HeisenbergMapper( + structures, total_energies, parent, **heisenberg_settings + ) return hmapper.get_heisenberg_model() @@ -71,8 +74,9 @@ def build_exchange_doc( 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. + 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 diff --git a/src/atomate2/common/schemas/exchange.py b/src/atomate2/common/schemas/exchange.py index 909676b0b2..407ae5b181 100644 --- a/src/atomate2/common/schemas/exchange.py +++ b/src/atomate2/common/schemas/exchange.py @@ -50,10 +50,14 @@ class ExchangeDocument(BaseModel): "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." + 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." + 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())." 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 index a036acd131..5ef6e5c48f 100644 --- a/src/atomate2/vampire/jobs/run_vampire.py +++ b/src/atomate2/vampire/jobs/run_vampire.py @@ -1,8 +1,19 @@ -from atomate2.vampire.vampire_caller import VampireCaller -from atomate2.vampire.schemas.vampire_output import VampireOutput -from pymatgen.analysis.magnetism.heisenberg import HeisenbergModel +"""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, @@ -29,4 +40,4 @@ def run_vampire( """ mc_settings = mc_settings or {} vampire_caller = VampireCaller(hm=heisenberg_model, **mc_settings) - return vampire_caller.output \ No newline at end of file + 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 index fd0172eda2..afa9fdeede 100644 --- a/src/atomate2/vampire/schemas/vampire_output.py +++ b/src/atomate2/vampire/schemas/vampire_output.py @@ -1,15 +1,32 @@ +"""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): - """This class processes results from a Vampire Monte Carlo simulation - and parses the critical temperature. + """Process results from a Vampire Monte Carlo simulation. + + Parses the critical temperature from the simulation output. """ - def __init__(self, parsed_out=None, nmats=None, critical_temp=None): - """ + 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). + 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 From 4f255d0bdd722994c0c29ae0fd36d6a7a3fc2bcb Mon Sep 17 00:00:00 2001 From: Luca Frey Date: Mon, 24 Aug 2026 16:38:35 +0200 Subject: [PATCH 11/11] Update VAMPIRE command handling in ExchangeMaker and VampireCaller; allow customization via ATOMATE2_VAMPIRE_CMD and improve error messaging. --- src/atomate2/common/flows/exchange.py | 6 ++++-- src/atomate2/vampire/jobs/run_vampire.py | 5 +++-- src/atomate2/vampire/vampire_caller.py | 23 +++++++++++++++-------- 3 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/atomate2/common/flows/exchange.py b/src/atomate2/common/flows/exchange.py index 4ee1205c68..3d08b6a26e 100644 --- a/src/atomate2/common/flows/exchange.py +++ b/src/atomate2/common/flows/exchange.py @@ -49,8 +49,10 @@ class ExchangeMaker(Maker): ``cutoff`` (Angstrom) and the distance-grouping ``tol``. run_vampire : bool Whether to run the Vampire Monte-Carlo step to estimate the critical - temperature. Requires the external ``vampire-serial`` binary on PATH; the step - raises a clear error if it is missing. Defaults to True (atomate1 parity). + 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``. diff --git a/src/atomate2/vampire/jobs/run_vampire.py b/src/atomate2/vampire/jobs/run_vampire.py index 5ef6e5c48f..08bc8dec73 100644 --- a/src/atomate2/vampire/jobs/run_vampire.py +++ b/src/atomate2/vampire/jobs/run_vampire.py @@ -22,8 +22,9 @@ def run_vampire( """Run Vampire Monte-Carlo to estimate the critical temperature. This wraps the (vendored) ``VampireCaller``, which shells out to the external - ``vampire-serial`` binary. A clear error is raised if the binary is not found - on PATH. + 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 ---------- diff --git a/src/atomate2/vampire/vampire_caller.py b/src/atomate2/vampire/vampire_caller.py index 7851b109ea..9b2f5c0da2 100644 --- a/src/atomate2/vampire/vampire_caller.py +++ b/src/atomate2/vampire/vampire_caller.py @@ -38,6 +38,8 @@ from __future__ import annotations import logging +import os +import shlex import subprocess from shutil import which @@ -47,16 +49,20 @@ from atomate2.vampire.schemas.vampire_output import VampireOutput -__author__ = "ncfrey" -__version__ = "0.1" -__maintainer__ = "Nathan C. Frey" -__email__ = "ncfrey@lbl.gov" +__author__ = "Luguza, ncfrey" +__version__ = "0.2" +__maintainer__ = "Luca Frey, Nathan C. Frey" +__email__ = "luca.frey@student.kit.edu, ncfrey@lbl.gov" __status__ = "Development" -__date__ = "June 2019" +__date__ = "August 2026" logger = logging.getLogger(__name__) -VAMP_EXE = which("vampire-serial") +# 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: @@ -81,7 +87,8 @@ class VampireCaller: @requires( VAMP_EXE is not None, - "VampireCaller requires vampire-serial to be in the path." + 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__( @@ -148,7 +155,7 @@ def __init__( self._create_ucf() # Call Vampire - with subprocess.Popen([VAMP_EXE], stdout=subprocess.PIPE, stderr=subprocess.PIPE) as process: + with subprocess.Popen([VAMP_EXE, *VAMP_CMD[1:]], stdout=subprocess.PIPE, stderr=subprocess.PIPE) as process: _stdout, stderr = process.communicate() stdout: str = _stdout.decode()