From 2a4a91b87bdef67855d0ad8d957b1c13393724d5 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 09:09:04 +0200 Subject: [PATCH 01/52] move to larixite.fdmnes module --- CHANGELOG.md | 6 ++++++ larixite/fdmnes/__init__.py | 1 + larixite/{fdmnes.py => fdmnes/input.py} | 0 3 files changed, 7 insertions(+) create mode 100644 larixite/fdmnes/__init__.py rename larixite/{fdmnes.py => fdmnes/input.py} (100%) diff --git a/CHANGELOG.md b/CHANGELOG.md index da45ded..3705d7b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,12 @@ The versioning scheme followed is `YYYY.MAJOR.MINOR`: - `MAJOR`: major feature (= change in API, no backward compatibility) - `MINOR`: minor feature (= bug fixes, minor changes not affecting the API, backward compatible) +## [2026.2.0 - not released yet] + +### `fmdnes` + + - Create a `fdmnes` module + ## [2026.1.0 - 18-May-2026] ### web application diff --git a/larixite/fdmnes/__init__.py b/larixite/fdmnes/__init__.py new file mode 100644 index 0000000..196db76 --- /dev/null +++ b/larixite/fdmnes/__init__.py @@ -0,0 +1 @@ +from .input import FdmnesXasInput \ No newline at end of file diff --git a/larixite/fdmnes.py b/larixite/fdmnes/input.py similarity index 100% rename from larixite/fdmnes.py rename to larixite/fdmnes/input.py From 91b5ac46567da9a391b2908d18065827c14bec13 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 09:12:49 +0200 Subject: [PATCH 02/52] fix import --- larixite/fdmnes/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/larixite/fdmnes/__init__.py b/larixite/fdmnes/__init__.py index 196db76..40e694f 100644 --- a/larixite/fdmnes/__init__.py +++ b/larixite/fdmnes/__init__.py @@ -1 +1 @@ -from .input import FdmnesXasInput \ No newline at end of file +from .input import FdmnesXasInput, struct2fdmnes \ No newline at end of file From d6b73a113e70b6ee0dc42e8fb3d9c321078ac7e1 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 09:55:25 +0200 Subject: [PATCH 03/52] fix imports --- larixite/fdmnes/__init__.py | 2 +- larixite/fdmnes/input.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/larixite/fdmnes/__init__.py b/larixite/fdmnes/__init__.py index 40e694f..3321fb8 100644 --- a/larixite/fdmnes/__init__.py +++ b/larixite/fdmnes/__init__.py @@ -1 +1 @@ -from .input import FdmnesXasInput, struct2fdmnes \ No newline at end of file +from .input import FdmnesXasInput, struct2fdmnes diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index f463382..9fb60fa 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -18,9 +18,10 @@ from larixite.utils import get_logger, strict_ascii, isotime, read_textfile from larixite.version import __version__ as larixite_version + logger = get_logger("larixite.fdmnes") -TEMPLATE_FOLDER = Path(Path(__file__).parent, "templates") +TEMPLATE_FOLDER = Path(Path(__file__).parent.parent, "templates") FDMNES_DEFAULT_PARAMS = { #: "FDMNES key name": True/False "Energpho": False, From dafcc4e05080202ae01e5bcc6a948a6323a60d4a Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 10:18:06 +0200 Subject: [PATCH 04/52] add yaml store/restore --- larixite/fdmnes/input.py | 72 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 9fb60fa..24089ff 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -10,6 +10,7 @@ """ import time +import yaml from dataclasses import dataclass from pathlib import Path from pymatgen.core import __version__ as pymatgen_version, Element, Molecule @@ -442,6 +443,77 @@ def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | fp.write(tp.read().format(**kwargs)) logger.info(f"written {fp.name}") + def dump_params(self, yamlpath: str | Path | None = None) -> Path: + """Dump input parameters to a YAML file + + Parameters + ---------- + yamlpath : str | Path | None + path to the output YAML file + if None: writes `{jobdir}/{fileout_prefix}_params.yaml` + if jobdir doesn't exist yet, writes next to the structural file + + Returns + ------- + Path + path to the written YAML file + """ + params_dict = { + "structpath": str(self.structpath), + "absorber": self.absorber.symbol, + "edge": self.edge, + "struct_type": self.struct_type, + "frame": self.frame, + "radius": self.radius, + "erange": self.erange, + "rself": self.rself, + "nself": self.nself, + "pself": self.pself, + "vmax": self.vmax, + "fileout_prefix": self.fileout_prefix, + "params": self.params, + "optimize": self.optimize, + "spacer": self.spacer, + "outdir": str(self.outdir) if self.outdir else None, + } + + if yamlpath is None: + if self._jobs: + yamlpath = self._jobs[-1] / f"{self.fileout_prefix}_params.yaml" + else: + yamlpath = Path(self.structpath).parent / f"{self.fileout_prefix}_params.yaml" + + if isinstance(yamlpath, str): + yamlpath = Path(yamlpath) + + with open(yamlpath, "w") as fp: + yaml.dump(params_dict, fp, default_flow_style=False) + logger.info(f"written {yamlpath.name}") + return yamlpath + + @classmethod + def from_yaml(cls, yamlpath: str | Path) -> "FdmnesXasInput": + """Restore a FdmnesXasInput object from a YAML parameters file + + Parameters + ---------- + yamlpath : str | Path + path to the YAML file previously written by `dump_params()` + + Returns + ------- + FdmnesXasInput + reconstructed input object + """ + if isinstance(yamlpath, str): + yamlpath = Path(yamlpath) + + with open(yamlpath) as fp: + params_dict = yaml.safe_load(fp) + + kwargs = {k: v for k, v in params_dict.items() if k != "structpath"} + return cls(params_dict["structpath"], **kwargs) + def struct2fdmnes( inp: str | Path, From d29d1c5faf166c23ab0c1f47788a33d6e34e1130 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 10:58:50 +0200 Subject: [PATCH 05/52] add dump at input write --- larixite/fdmnes/input.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 24089ff..6b79df5 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -218,7 +218,7 @@ def optimize_params(self) -> dict: self.vmax = -6 logger.info(f"Vmax set to {self.vmax} for molecule") - #enable SCF + #: enable SCF params["SCF"] = True logger.info("SCF enabled") @@ -405,6 +405,7 @@ def write_input( fp.write(f"1\n{fileout_name}\n") logger.info(f"written {fp.name}") self._jobs.append(jobdir) + _ = self.dump_params(jobdir / f"{self.fileout_prefix}_params.yaml") return jobdir def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | None = None, ncpus: int = 8): From 5c14e666274673939cf9b722706b05ec309058d3 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 11:15:52 +0200 Subject: [PATCH 06/52] add visulize --- larixite/struct/__init__.py | 1 + larixite/struct/visualize.py | 142 +++++++++++++++++++++++++++++++++++ pyproject.toml | 1 + 3 files changed, 144 insertions(+) create mode 100644 larixite/struct/visualize.py diff --git a/larixite/struct/__init__.py b/larixite/struct/__init__.py index 2b7e24b..b6cf5cc 100644 --- a/larixite/struct/__init__.py +++ b/larixite/struct/__init__.py @@ -22,6 +22,7 @@ from larixite.struct.xas_xyz import XasStructureXyz from larixite.utils import get_logger, read_textfile from larixite.amcsd_utils import PMG_CIF_OPTS +from larixite.struct.visualize import visualize logger = get_logger("larixite.struct") diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py new file mode 100644 index 0000000..f302eec --- /dev/null +++ b/larixite/struct/visualize.py @@ -0,0 +1,142 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +3D visualization of XasStructure local environment using py3Dmol +=============================================================== +""" + +import math + +from larixite.struct.xas import XasStructure +from larixite.utils import get_logger + +logger = get_logger("larixite.struct") + +try: + import py3Dmol + + HAS_PY3DMOL = True +except ImportError: + HAS_PY3DMOL = False + + +def _round_up(x: float) -> float: + return math.ceil(x * 100) / 100 + + +def _cluster_to_xyz(mol) -> tuple[str, list[str]]: + """Convert a pymatgen Molecule to an XYZ format string and the list of unique elements.""" + coords = [] + elements = [] + for site in mol: + elem = site.species_string + cart = site.coords + coords.append((elem, cart)) + if elem not in elements: + elements.append(elem) + + elements = sorted(elements) + + output_str = f"{len(coords)}\n\n" + for elem, coord in coords: + coord_str = " ".join(f"{c:.6f}" for c in coord) + output_str += f"{elem} {coord_str}\n" + + return output_str, elements + + +def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False): + """ + Display a 3D visualization of the local structure around the absorber atom. + + Args: + struct: XasStructure object with the absorber defined. + radius: radius in Angstroms for the visualization sphere around the central atom. + unitcell: if True, displays the structure unit cell (CIF structures only). + + Returns: + py3Dmol.view object, or None if py3Dmol is not installed. + """ + if not HAS_PY3DMOL: + logger.error("py3Dmol not installed! -> run `pip install py3Dmol`") + return + + radius = _round_up(radius) + + mol = struct.build_cluster(radius=radius) + + xyz, elems = _cluster_to_xyz(mol) + + xyzview = py3Dmol.view(width=600, height=600) + xyzview.addModel(xyz, "xyz") + + if unitcell: + try: + lat = struct.struct.lattice + a, b, c = lat.a, lat.b, lat.c + alpha, beta, gamma = lat.alpha, lat.beta, lat.gamma + m = xyzview.getModel() + m.setCrystData(a, b, c, alpha, beta, gamma) + xyzview.addUnitCell() + except AttributeError: + logger.warning("Unit cell visualization not available for this structure type") + + colors = [ + "red", + "green", + "blue", + "orange", + "yellow", + "white", + "purple", + "pink", + "brown", + "black", + "gray", + "cyan", + "magenta", + "olive", + "navy", + "teal", + "maroon", + "turquoise", + "indigo", + "salmon", + ] + color_elems = {} + for idx, elem in enumerate(elems): + color_elems[elem] = colors[idx] + + for elem in elems: + color = color_elems[elem] + xyzview.setStyle( + {"elem": elem}, + { + "stick": { + "radius": 0.1, + "opacity": 1, + "hidden": False, + "color": color, + }, + "sphere": {"color": color, "radius": 0.4, "opacity": 1}, + }, + ) + + xyzview.addLabel( + "Abs", + { + "fontColor": "black", + "fontSize": 14, + "backgroundColor": "white", + "backgroundOpacity": 0.8, + "showBackground": True, + }, + {"index": 0}, + ) + + xyzview.zoomTo() + xyzview.show() + + logger.info(color_elems) + return xyzview diff --git a/pyproject.toml b/pyproject.toml index a90dc92..f455754 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,7 @@ dependencies = [ "mp-api>=0.45.8", "emmet-core>=0.84.9", "requests", + "py3Dmol", ] license = "MIT" From e116cbbf843313b91f41497f7c02b4320a8924a2 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 11:57:15 +0200 Subject: [PATCH 07/52] visualize with covalent radius --- larixite/struct/visualize.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py index f302eec..8acb672 100644 --- a/larixite/struct/visualize.py +++ b/larixite/struct/visualize.py @@ -8,6 +8,8 @@ import math +from pymatgen.analysis.local_env import CovalentRadius + from larixite.struct.xas import XasStructure from larixite.utils import get_logger @@ -108,8 +110,11 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) for idx, elem in enumerate(elems): color_elems[elem] = colors[idx] + covalent_r = CovalentRadius().radius + for elem in elems: color = color_elems[elem] + sphere_r = max(0.2, min(covalent_r.get(elem, 0.5) / 2, 0.75)) xyzview.setStyle( {"elem": elem}, { @@ -119,7 +124,7 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) "hidden": False, "color": color, }, - "sphere": {"color": color, "radius": 0.4, "opacity": 1}, + "sphere": {"color": color, "radius": sphere_r, "opacity": 1}, }, ) From 105a6d02a58d7331997bb86403a06b6bd10eec7d Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 12:55:34 +0200 Subject: [PATCH 08/52] better visualize central absorber --- larixite/struct/visualize.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py index 8acb672..313fc7c 100644 --- a/larixite/struct/visualize.py +++ b/larixite/struct/visualize.py @@ -94,8 +94,6 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) "purple", "pink", "brown", - "black", - "gray", "cyan", "magenta", "olive", @@ -128,16 +126,19 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) }, ) - xyzview.addLabel( - "Abs", + absorber_sym = struct.absorber.symbol + absorber_r = max(0.2, min(covalent_r.get(absorber_sym, 0.5) / 2, 0.75)) + xyzview.setStyle( + {"index": 0}, { - "fontColor": "black", - "fontSize": 14, - "backgroundColor": "white", - "backgroundOpacity": 0.8, - "showBackground": True, + "stick": { + "radius": 0.1, + "opacity": 1, + "hidden": False, + "color": "gray", + }, + "sphere": {"color": "gray", "radius": absorber_r, "opacity": 1}, }, - {"index": 0}, ) xyzview.zoomTo() From 5a4e97851c7e170a5e4f9fa85c1733a083b40124 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 13:15:21 +0200 Subject: [PATCH 09/52] use vesta color from pymatgen (hard coded) --- larixite/struct/visualize.py | 142 +++++++++++++++++++++++++++++------ 1 file changed, 118 insertions(+), 24 deletions(-) diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py index 313fc7c..791e6bc 100644 --- a/larixite/struct/visualize.py +++ b/larixite/struct/visualize.py @@ -23,6 +23,119 @@ HAS_PY3DMOL = False +VESTA_COLORS = { + "Ac": "#70abfa", + "Ag": "#c0c0c0", + "Al": "#81b2d6", + "Am": "#545cf2", + "Ar": "#cffec4", + "As": "#74d057", + "At": "#754f45", + "Au": "#ffd123", + "B": "#1fa20f", + "Ba": "#00c900", + "Be": "#5ed77b", + "Bh": "#e00038", + "Bi": "#9e4fb5", + "Bk": "#8a4fe3", + "Br": "#7e3102", + "C": "#4c4c4c", + "Ca": "#5a96bd", + "Cd": "#ffd98f", + "Ce": "#ffffc7", + "Cf": "#a136d4", + "Cl": "#31fc02", + "Cm": "#785ce3", + "Co": "#0000af", + "Cr": "#00009e", + "Cs": "#57178f", + "Cu": "#2247dc", + "Db": "#d1004f", + "Dy": "#1fffc7", + "Er": "#00e675", + "Es": "#b31fd4", + "Eu": "#61ffc7", + "F": "#b0b9e6", + "Fe": "#b57100", + "Fm": "#b31fba", + "Fr": "#420066", + "Ga": "#9ee373", + "Gd": "#45ffc7", + "Ge": "#7e6ea6", + "H": "#ffcccc", + "He": "#fce8ce", + "Hf": "#4dc2ff", + "Hg": "#b8b8d0", + "Ho": "#00ff9c", + "Hs": "#e6002e", + "I": "#940094", + "In": "#a67573", + "Ir": "#175487", + "K": "#a121f6", + "Kr": "#fac1f3", + "La": "#5ac449", + "Li": "#86df73", + "Lr": "#c70066", + "Lu": "#00ab24", + "Md": "#b30da6", + "Mg": "#fb7b15", + "Mn": "#a7089d", + "Mo": "#54b5b5", + "Mt": "#eb0026", + "N": "#b0b9e6", + "Na": "#f9dc3c", + "Nb": "#73c2c9", + "Nd": "#c7ffc7", + "Ne": "#fe37b5", + "Ni": "#b7bbbd", + "No": "#bd0d87", + "Np": "#0080ff", + "O": "#fe0300", + "Os": "#266696", + "P": "#c09cc2", + "Pa": "#00a1ff", + "Pb": "#575961", + "Pd": "#006985", + "Pm": "#a3ffc7", + "Po": "#ab5c00", + "Pr": "#d9ffc7", + "Pt": "#d0d0e0", + "Pu": "#006bff", + "Ra": "#007d00", + "Rb": "#702eb0", + "Re": "#267dab", + "Rf": "#cc0059", + "Rh": "#0a7d8c", + "Rn": "#428296", + "Ru": "#248f8f", + "S": "#fffa00", + "Sb": "#9e63b5", + "Sc": "#b563ab", + "Se": "#9aef0f", + "Sg": "#d90045", + "Si": "#1b3bfa", + "Sm": "#8fffc7", + "Sn": "#9a8eb9", + "Sr": "#00ff00", + "Ta": "#4da6ff", + "Tb": "#30ffc7", + "Tc": "#3b9e9e", + "Te": "#d47a00", + "Th": "#00baff", + "Ti": "#78caff", + "Tl": "#a6544d", + "Tm": "#00d452", + "U": "#008fff", + "V": "#e51900", + "W": "#2194d6", + "Xe": "#429eb0", + "Y": "#94ffff", + "Yb": "#00bf38", + "Zn": "#8f8f81", + "Zr": "#00ff00", +} + + def _round_up(x: float) -> float: return math.ceil(x * 100) / 100 @@ -84,29 +197,9 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) except AttributeError: logger.warning("Unit cell visualization not available for this structure type") - colors = [ - "red", - "green", - "blue", - "orange", - "yellow", - "white", - "purple", - "pink", - "brown", - "cyan", - "magenta", - "olive", - "navy", - "teal", - "maroon", - "turquoise", - "indigo", - "salmon", - ] color_elems = {} - for idx, elem in enumerate(elems): - color_elems[elem] = colors[idx] + for elem in elems: + color_elems[elem] = VESTA_COLORS.get(elem, "#888888") covalent_r = CovalentRadius().radius @@ -128,6 +221,7 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) absorber_sym = struct.absorber.symbol absorber_r = max(0.2, min(covalent_r.get(absorber_sym, 0.5) / 2, 0.75)) + absorber_color = color_elems.get(absorber_sym, "#888888") xyzview.setStyle( {"index": 0}, { @@ -135,9 +229,9 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) "radius": 0.1, "opacity": 1, "hidden": False, - "color": "gray", + "color": absorber_color, }, - "sphere": {"color": "gray", "radius": absorber_r, "opacity": 1}, + "sphere": {"color": absorber_color, "radius": absorber_r, "opacity": 0.3}, }, ) From 0a398dab991f7a76f0e9ad0a1babe6eac175890f Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 13:32:18 +0200 Subject: [PATCH 10/52] try lighter absorber --- larixite/struct/visualize.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py index 791e6bc..8d030bb 100644 --- a/larixite/struct/visualize.py +++ b/larixite/struct/visualize.py @@ -201,8 +201,19 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) for elem in elems: color_elems[elem] = VESTA_COLORS.get(elem, "#888888") + def _lighten_hex(hex_color: str, factor: float = 0.75) -> str: + r = int(int(hex_color[1:3], 16) * factor) + g = int(int(hex_color[3:5], 16) * factor) + b = int(int(hex_color[5:7], 16) * factor) + return "#{:02x}{:02x}{:02x}".format(r, g, b) + covalent_r = CovalentRadius().radius + absorber_sym = struct.absorber.symbol + absorber_color = color_elems.get(absorber_sym, "#888888") + absorber_r = max(0.2, min(covalent_r.get(absorber_sym, 0.5) / 2, 0.75)) + absorber_light = _lighten_hex(absorber_color, 0.8) + for elem in elems: color = color_elems[elem] sphere_r = max(0.2, min(covalent_r.get(elem, 0.5) / 2, 0.75)) @@ -219,9 +230,8 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) }, ) - absorber_sym = struct.absorber.symbol - absorber_r = max(0.2, min(covalent_r.get(absorber_sym, 0.5) / 2, 0.75)) - absorber_color = color_elems.get(absorber_sym, "#888888") + #: reset index 0 style to avoid duplication, then reapply with lighter color + xyzview.setStyle({"index": 0}, None) # clear xyzview.setStyle( {"index": 0}, { @@ -229,9 +239,9 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) "radius": 0.1, "opacity": 1, "hidden": False, - "color": absorber_color, + "color": absorber_light, }, - "sphere": {"color": absorber_color, "radius": absorber_r, "opacity": 0.3}, + "sphere": {"color": absorber_light, "radius": absorber_r, "opacity": 1}, }, ) From 9ebd5c7e5762e213c7054996ff67857765e19e51 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 13:58:40 +0200 Subject: [PATCH 11/52] start migrating coordination environment analysis from struct2xas --- larixite/struct/__init__.py | 1 + larixite/struct/analyze.py | 177 ++++++++++++++++++++++++++++++++++++ 2 files changed, 178 insertions(+) create mode 100644 larixite/struct/analyze.py diff --git a/larixite/struct/__init__.py b/larixite/struct/__init__.py index b6cf5cc..35750f0 100644 --- a/larixite/struct/__init__.py +++ b/larixite/struct/__init__.py @@ -23,6 +23,7 @@ from larixite.utils import get_logger, read_textfile from larixite.amcsd_utils import PMG_CIF_OPTS from larixite.struct.visualize import visualize +from larixite.struct.analyze import get_coord_env, show_coord_env logger = get_logger("larixite.struct") diff --git a/larixite/struct/analyze.py b/larixite/struct/analyze.py new file mode 100644 index 0000000..c3192ec --- /dev/null +++ b/larixite/struct/analyze.py @@ -0,0 +1,177 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +Coordination environment analysis for XAS structures +==================================================== +""" + +import numpy as np +from larixite.struct.xas import XasStructure +from larixite.utils import get_logger, pprint + +logger = get_logger("larixite.struct.analyze") + + +def get_coord_env(struct: XasStructure, radius: float = 3) -> list: + """Determine the coordination environment around the absorber atom. + + For crystal structures (CIF), uses pymatgen's chemenv toolkit + (LocalGeometryFinder, BVAnalyzer, MultiWeightsChemenvStrategy) to + identify coordination geometry symbols such as: + + - S:4 - Square Plane + - T:4 - Tetrahedral + - T:5 - Trigonal bipyramid + - S:5 - Square pyramidal + - O:6 - Octahedral + - T:6 - Trigonal prism + + For molecular structures (XYZ), uses pymatgen's CrystalNN to + identify neighboring atoms. + + Parameters + ---------- + struct : XasStructure + The XAS structure object. + radius : float, optional + Maximum distance [Angstrom] to search for neighboring atoms. + Default is 3. + + Returns + ------- + list + List of lists, each containing: + + [0]: Info string about which site is being analyzed. + [1]: Coordination environment dictionary with keys: + - ce_symbol: coordination geometry symbol + - ce_fraction: probability for given coordination env (0-1) + - csm: continuous symmetry measure (lower = more symmetric) + - permutation: atom arrangement order + [2]: List of neighbor site dictionaries with keys: + - site: pymatgen Site with element and coordinates + - site_index: structure index of the neighbor atom + - image: displacement vector from original site + - weight: significance of the neighbor to the coordination + """ + idx_abs_site = struct.absorber_idx + coord_env_list = [] + + if struct.struct_type == "crystal": + from pymatgen.analysis.bond_valence import BVAnalyzer + from pymatgen.analysis.chemenv.coordination_environments.coordination_geometry_finder import ( + LocalGeometryFinder, + ) + from pymatgen.analysis.chemenv.coordination_environments.structure_environments import ( + LightStructureEnvironments, + ) + from pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies import ( + MultiWeightsChemenvStrategy, + ) + + lgf = LocalGeometryFinder() + lgf.setup_structure(struct.struct) + + bva = BVAnalyzer() + try: + valences = bva.get_valences(structure=struct.struct) + except ValueError: + valences = "undefined" + + se = lgf.compute_structure_environments( + max_cn=6, + valences=valences, + only_indices=[idx_abs_site], + only_symbols=["S:4", "T:4", "T:5", "S:5", "O:6", "T:6"], + ) + + logger.debug(se.voronoi.neighbors_distances[idx_abs_site][0]["max"]) + + strategy = MultiWeightsChemenvStrategy.stats_article_weights_parameters() + lse = LightStructureEnvironments.from_structure_environments( + strategy=strategy, structure_environments=se + ) + coord_env_ce = lse.coordination_environments[idx_abs_site] + coord_env = lse._all_nbs_sites + + elif struct.struct_type == "molecule": + from pymatgen.analysis.local_env import CrystalNN + + obj = CrystalNN() + + coord_env = obj.get_nn_info(struct.struct, idx_abs_site) + for site in coord_env: + site["site"].cart_coords = struct.struct[site["site_index"]].coords + coord_dict = obj.get_cn_dict(struct.struct, idx_abs_site) + coord_env_ce = {"ce_symbol": f"Elements Dict = {coord_dict}"} + + else: + raise ValueError(f"Unknown structure type: {struct.struct_type}") + + coord_env_list.append( + [ + f"Coord. Env. for Site {struct.absorber_sites[0][0]}", + coord_env_ce, + coord_env, + ] + ) + return coord_env_list + + +def show_coord_env(struct: XasStructure, **kws): + """Pretty-print coordination environment information for the absorber atom. + + Calls get_coord_env() and formats the results as a table with element + species and distances from the absorber site. + + Parameters + ---------- + struct : XasStructure + The XAS structure object. + **kws + Additional keyword arguments passed to get_coord_env(). + """ + coord_env = get_coord_env(struct, **kws)[0] + abs_site_coord = struct.absorber_sites[0][4] + + elems_dist = [] + sites = coord_env[2] + + if len(sites) == 0: + logger.warning("Empty coordination environment") + return None + + if struct.struct_type == "crystal": + coord_sym = [ + coord_env[1][i]["ce_symbol"] for i in range(len(coord_env[1])) + ] + for site in sites: + elems_dist.append( + ( + site["site"].species, + round(np.linalg.norm(site["site"].coords - abs_site_coord), 5), + ) + ) + elif struct.struct_type == "molecule": + coord_sym = coord_env[1]["ce_symbol"] + for site in sites: + elems_dist.append( + ( + site["site"].species, + round( + np.linalg.norm(site["site"].cart_coords - abs_site_coord), 5 + ), + ) + ) + + elems_dist = sorted(elems_dist, key=lambda x: x[1]) + print( + f"Coord. Env. from absorber atom: {struct.absorber.symbol} " + f"at site {struct.absorber_sites[0][0]}" + ) + print(coord_sym) + header = ["Element", "Distance"] + matrix = [header] + matrix.extend(elems_dist) + pprint(matrix) From 15b4cd345b17f7dba503a5dcfaa6aab0648fd653 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 14:26:59 +0200 Subject: [PATCH 12/52] coordination environment analysis by shell --- larixite/struct/analyze.py | 284 +++++++++++++++++++++++++------------ 1 file changed, 197 insertions(+), 87 deletions(-) diff --git a/larixite/struct/analyze.py b/larixite/struct/analyze.py index c3192ec..f8ccf96 100644 --- a/larixite/struct/analyze.py +++ b/larixite/struct/analyze.py @@ -13,7 +13,40 @@ logger = get_logger("larixite.struct.analyze") -def get_coord_env(struct: XasStructure, radius: float = 3) -> list: +def _group_into_shells(dists, toler): + """Group sorted distances into shells. + + Distances within `toler` Angstrom of each other are considered + belonging to the same coordination shell. + + Parameters + ---------- + dists : list of float + Sorted distances from the absorber. + toler : float + Tolerance [Angstrom] to group distances into the same shell. + + Returns + ------- + list of list of float + Each inner list contains the distances belonging to one shell. + """ + if len(dists) == 0: + return [] + + shells = [[dists[0]]] + for d in dists[1:]: + if d - shells[-1][-1] <= toler: + shells[-1].append(d) + else: + shells.append([d]) + + return shells + + +def get_coord_env( + struct: XasStructure, shell: int = 1, toler: float = 0.1 +) -> list: """Determine the coordination environment around the absorber atom. For crystal structures (CIF), uses pymatgen's chemenv toolkit @@ -30,13 +63,20 @@ def get_coord_env(struct: XasStructure, radius: float = 3) -> list: For molecular structures (XYZ), uses pymatgen's CrystalNN to identify neighboring atoms. + Neighbors are grouped into coordination shells according to their + distances from the absorber: distances differing by less than + `toler` are considered the same shell. + Parameters ---------- struct : XasStructure The XAS structure object. - radius : float, optional - Maximum distance [Angstrom] to search for neighboring atoms. - Default is 3. + shell : int, optional + Number of coordination shells to return. Default is 1 (first + coordination shell only). + toler : float, optional + Distance tolerance [Angstrom] to group neighbors into the same + coordination shell. Default is 0.1. Returns ------- @@ -49,129 +89,199 @@ def get_coord_env(struct: XasStructure, radius: float = 3) -> list: - ce_fraction: probability for given coordination env (0-1) - csm: continuous symmetry measure (lower = more symmetric) - permutation: atom arrangement order - [2]: List of neighbor site dictionaries with keys: - - site: pymatgen Site with element and coordinates - - site_index: structure index of the neighbor atom - - image: displacement vector from original site - - weight: significance of the neighbor to the coordination + [2]: List of shell dictionaries, one per requested shell, with keys: + - shell_idx: shell index (1-based) + - avg_dist: average distance of atoms in the shell [Angstrom] + - cn: coordination number (number of atoms in the shell) + - elems: dict mapping element symbol to count in this shell + - neighbors: list of neighbor site dicts with keys: + - site: pymatgen Site with element and coordinates + - site_index: structure index of the neighbor atom + - dist: distance from absorber [Angstrom] + - image: displacement vector from original site (CIF) + - weight: significance of the neighbor (CIF) """ idx_abs_site = struct.absorber_idx - coord_env_list = [] + abs_site = struct.struct[idx_abs_site] + coord_env_ce = {} + # Determine geometry symbol for the 1st coordination shell if struct.struct_type == "crystal": - from pymatgen.analysis.bond_valence import BVAnalyzer - from pymatgen.analysis.chemenv.coordination_environments.coordination_geometry_finder import ( - LocalGeometryFinder, - ) - from pymatgen.analysis.chemenv.coordination_environments.structure_environments import ( - LightStructureEnvironments, - ) - from pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies import ( - MultiWeightsChemenvStrategy, - ) + try: + from pymatgen.analysis.bond_valence import BVAnalyzer + from pymatgen.analysis.chemenv.coordination_environments.coordination_geometry_finder import ( + LocalGeometryFinder, + ) + from pymatgen.analysis.chemenv.coordination_environments.structure_environments import ( + LightStructureEnvironments, + ) + from pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies import ( + MultiWeightsChemenvStrategy, + ) - lgf = LocalGeometryFinder() - lgf.setup_structure(struct.struct) + lgf = LocalGeometryFinder() + lgf.setup_structure(struct.struct) - bva = BVAnalyzer() - try: - valences = bva.get_valences(structure=struct.struct) - except ValueError: - valences = "undefined" - - se = lgf.compute_structure_environments( - max_cn=6, - valences=valences, - only_indices=[idx_abs_site], - only_symbols=["S:4", "T:4", "T:5", "S:5", "O:6", "T:6"], - ) + bva = BVAnalyzer() + try: + valences = bva.get_valences(structure=struct.struct) + except ValueError: + valences = "undefined" - logger.debug(se.voronoi.neighbors_distances[idx_abs_site][0]["max"]) + se = lgf.compute_structure_environments( + max_cn=6, + valences=valences, + only_indices=[idx_abs_site], + only_symbols=["S:4", "T:4", "T:5", "S:5", "O:6", "T:6"], + ) - strategy = MultiWeightsChemenvStrategy.stats_article_weights_parameters() - lse = LightStructureEnvironments.from_structure_environments( - strategy=strategy, structure_environments=se - ) - coord_env_ce = lse.coordination_environments[idx_abs_site] - coord_env = lse._all_nbs_sites + logger.debug(se.voronoi.neighbors_distances[idx_abs_site][0]["max"]) + + strategy = MultiWeightsChemenvStrategy.stats_article_weights_parameters() + lse = LightStructureEnvironments.from_structure_environments( + strategy=strategy, structure_environments=se + ) + coord_env_ce = lse.coordination_environments[idx_abs_site] + except Exception: + logger.warning( + "chemenv analysis failed, falling back to basic neighbor search" + ) + coord_env_ce = {"ce_symbol": "chemenv analysis failed"} elif struct.struct_type == "molecule": - from pymatgen.analysis.local_env import CrystalNN + try: + from pymatgen.analysis.local_env import CrystalNN - obj = CrystalNN() + obj = CrystalNN() + coord_dict = obj.get_cn_dict(struct.struct, idx_abs_site) + coord_env_ce = {"ce_symbol": f"Elements Dict = {coord_dict}"} + except Exception: + logger.warning( + "CrystalNN analysis failed, falling back to basic neighbor search" + ) + coord_env_ce = {"ce_symbol": "CrystalNN analysis failed"} - coord_env = obj.get_nn_info(struct.struct, idx_abs_site) - for site in coord_env: - site["site"].cart_coords = struct.struct[site["site_index"]].coords - coord_dict = obj.get_cn_dict(struct.struct, idx_abs_site) - coord_env_ce = {"ce_symbol": f"Elements Dict = {coord_dict}"} + # Build full neighbor list using struct.get_neighbors() up to cluster_size + # This works for arbitrary distances, unlike chemenv/CrystalNN which only + # cover the 1st coordination shell + search_radius = struct.cluster_size + sphere = struct.struct.get_neighbors(abs_site, search_radius) - else: - raise ValueError(f"Unknown structure type: {struct.struct_type}") + neighbors = [] + for item in sphere: + # Structure: (PeriodicNeighbor, distance, site_index, jimage) + # Molecule: (PeriodicNeighbor, distance, site_index) + site = item[0] + dist = round(float(item[1]), 5) + si = item[2] if len(item) > 2 else None + if dist < 1e-6: + continue + neighbors.append({"site": site, "dist": dist, "site_index": si}) + + neighbors.sort(key=lambda x: x["dist"]) + + # Group into shells by distance clustering + dists = [nb["dist"] for nb in neighbors] + dist_shells = _group_into_shells(dists, toler) - coord_env_list.append( + # Build shell dicts for requested number of shells + shell_dicts = [] + idx = 0 + for s in range(min(shell, len(dist_shells))): + shell_dists = dist_shells[s] + shell_len = len(shell_dists) + shell_nbs = neighbors[idx : idx + shell_len] + idx += shell_len + + avg_dist = round(float(np.mean(shell_dists)), 5) + elems = {} + for nb in shell_nbs: + sym = nb["site"].species_string + elems[sym] = elems.get(sym, 0) + 1 + + shell_dicts.append( + { + "shell_idx": s + 1, + "avg_dist": avg_dist, + "cn": shell_len, + "elems": elems, + "neighbors": shell_nbs, + } + ) + + coord_env_list = [ [ f"Coord. Env. for Site {struct.absorber_sites[0][0]}", coord_env_ce, - coord_env, + shell_dicts, ] - ) + ] return coord_env_list def show_coord_env(struct: XasStructure, **kws): """Pretty-print coordination environment information for the absorber atom. - Calls get_coord_env() and formats the results as a table with element - species and distances from the absorber site. + Calls get_coord_env() and formats the results as a table with shells, + distances, coordination numbers and element breakdown. Parameters ---------- struct : XasStructure The XAS structure object. **kws - Additional keyword arguments passed to get_coord_env(). + Additional keyword arguments passed to get_coord_env(): + shell : int, default 1 + toler : float, default 0.1 """ coord_env = get_coord_env(struct, **kws)[0] - abs_site_coord = struct.absorber_sites[0][4] - - elems_dist = [] - sites = coord_env[2] + shells = coord_env[2] - if len(sites) == 0: + if len(shells) == 0: logger.warning("Empty coordination environment") return None - if struct.struct_type == "crystal": - coord_sym = [ - coord_env[1][i]["ce_symbol"] for i in range(len(coord_env[1])) - ] - for site in sites: - elems_dist.append( - ( - site["site"].species, - round(np.linalg.norm(site["site"].coords - abs_site_coord), 5), - ) - ) - elif struct.struct_type == "molecule": - coord_sym = coord_env[1]["ce_symbol"] - for site in sites: - elems_dist.append( - ( - site["site"].species, - round( - np.linalg.norm(site["site"].cart_coords - abs_site_coord), 5 - ), - ) - ) + # Extract coordination geometry symbol(s) + ce_info = coord_env[1] + if isinstance(ce_info, list): + # chemenv returns a list, extract symbol from the first entry + coord_sym = ce_info[0]["ce_symbol"] if ce_info else "not determined" + else: + # fallback: direct dict with ce_symbol key + coord_sym = ce_info.get("ce_symbol", "not determined") - elems_dist = sorted(elems_dist, key=lambda x: x[1]) print( f"Coord. Env. from absorber atom: {struct.absorber.symbol} " f"at site {struct.absorber_sites[0][0]}" ) - print(coord_sym) - header = ["Element", "Distance"] + print(f"Geometry: {coord_sym}") + print() + + header = ["Shell", "Dist. (Ang)", "CN", "Elements"] matrix = [header] - matrix.extend(elems_dist) + + for sh in shells: + elems_str = ", ".join( + f"{sym}: {cnt}" for sym, cnt in sorted(sh["elems"].items()) + ) + matrix.append( + [ + sh["shell_idx"], + sh["avg_dist"], + sh["cn"], + elems_str, + ] + ) + pprint(matrix) + print() + + # Detailed neighbor table per shell + for sh in shells: + print(f"Shell {sh['shell_idx']} (avg dist: {sh['avg_dist']} A, CN: {sh['cn']})") + det_header = ["Element", "Distance (Ang)"] + det_matrix = [det_header] + for nb in sh["neighbors"]: + det_matrix.append([nb["site"].species_string, nb["dist"]]) + pprint(det_matrix) + print() From 1c6cdc0b539dfa1e83fe6171f442b00c1c8bbb5b Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 14:34:05 +0200 Subject: [PATCH 13/52] add ce_symbols --- larixite/struct/analyze.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/larixite/struct/analyze.py b/larixite/struct/analyze.py index f8ccf96..6fa9b34 100644 --- a/larixite/struct/analyze.py +++ b/larixite/struct/analyze.py @@ -12,6 +12,38 @@ logger = get_logger("larixite.struct.analyze") +# Mapping of chemenv coordination environment symbols to human-readable names +# Based on pymatgen chemenv ideal polyhedra definitions +CE_SYMBOL_DESCRIPTION = { + "L:2": "Linear", + "T:3": "Trigonal planar", + "B:3": "Bent", + "T:4": "Tetrahedral", + "S:4": "Square planar", + "F:4": "See-saw", + "T:5": "Trigonal bipyramidal", + "S:5": "Square pyramidal", + "B:5": "Bay", + "R:5": "Rocket", + "G:5": "Envelope", + "E:5": "Eclipsed pentagonal", + "S:6": "Hexagonal planar", + "A:6": "Hexagonal antiprismatic", + "T:6": "Trigonal prismatic", + "O:6": "Octahedral", + "Y:6": "Pentagonal pyramidal", + "M:6": "Pentagonal bisphenoidal", + "H:6": "Triangular cup", + "S:7": "Pentagonal bipyramidal", + "D:7": "Heptagonal planar", + "S:8": "Square antiprismatic", + "D:8": "Dodecahedral", + "A:8": "Hexagonal bipyramidal", + "A:9": "Monocapped square antiprismatic", + "A:10": "Bicapped square antiprismatic", + "S:12": "Icosahedral", +} + def _group_into_shells(dists, toler): """Group sorted distances into shells. From 25368f4cfa4710f27c3d6035e49843c71136a642 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 14:37:55 +0200 Subject: [PATCH 14/52] add ce description --- larixite/struct/analyze.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/larixite/struct/analyze.py b/larixite/struct/analyze.py index 6fa9b34..6a3e1a0 100644 --- a/larixite/struct/analyze.py +++ b/larixite/struct/analyze.py @@ -286,7 +286,11 @@ def show_coord_env(struct: XasStructure, **kws): f"Coord. Env. from absorber atom: {struct.absorber.symbol} " f"at site {struct.absorber_sites[0][0]}" ) - print(f"Geometry: {coord_sym}") + desc = CE_SYMBOL_DESCRIPTION.get(coord_sym, "") + if desc: + print(f"Geometry: {coord_sym} ({desc})") + else: + print(f"Geometry: {coord_sym}") print() header = ["Shell", "Dist. (Ang)", "CN", "Elements"] From a6fe670d7a9e32db5b43bb2b6257c7eed05fb2e1 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 14:42:16 +0200 Subject: [PATCH 15/52] add tests --- tests/test_analyze.py | 135 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 tests/test_analyze.py diff --git a/tests/test_analyze.py b/tests/test_analyze.py new file mode 100644 index 0000000..6891552 --- /dev/null +++ b/tests/test_analyze.py @@ -0,0 +1,135 @@ +"""Tests for larixite.struct.analyze — get_coord_env and show_coord_env.""" + +from pathlib import Path +from io import StringIO +import sys +import warnings + +import pytest +import numpy as np + +from larixite.struct import get_structure +from larixite.struct.analyze import get_coord_env, show_coord_env + +STRUCTS = Path(__file__).parent / "structs" + + +# ── helpers ─────────────────────────────────────────────────────── + + +def cap_show(output, struct, **kws): + """Capture the stdout produced by show_coord_env and return it.""" + old_out = sys.stdout + sys.stdout = buf = StringIO() + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + output = show_coord_env(struct, **kws) + sys.stdout = old_out + return buf.getvalue() + + +@pytest.fixture(scope="module") +def cif_zno(): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + return get_structure(STRUCTS / "ZnO_mp-2133.cif", "Zn") + + +@pytest.fixture(scope="module") +def xyz_cuo6(): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + return get_structure(STRUCTS / "CuO6_D4h.xyz", "Cu") + + +# ── CIF (crystal) tests ───────────────────────────────────────── + + +class TestCoordEnvCIF: + def test_shell1_returned(self, cif_zno): + ce = get_coord_env(cif_zno, shell=1) + shells = ce[0][2] + assert len(shells) == 1 + assert shells[0]["shell_idx"] == 1 + + def test_shell1_cn_positive(self, cif_zno): + ce = get_coord_env(cif_zno, shell=1) + shells = ce[0][2] + assert shells[0]["cn"] > 0 + + def test_shell1_has_elements(self, cif_zno): + ce = get_coord_env(cif_zno, shell=1) + shells = ce[0][2] + assert len(shells[0]["elems"]) > 0 + + def test_shell1_neighbors_sorted(self, cif_zno): + ce = get_coord_env(cif_zno, shell=1) + shells = ce[0][2] + dists = [nb["dist"] for nb in shells[0]["neighbors"]] + assert dists == sorted(dists) + + def test_multiple_shells_returned(self, cif_zno): + ce = get_coord_env(cif_zno, shell=3) + shells = ce[0][2] + assert len(shells) >= 2 + + def test_multiple_shells_increasing_distance(self, cif_zno): + ce = get_coord_env(cif_zno, shell=3) + shells = ce[0][2] + for i in range(1, len(shells)): + assert shells[i]["avg_dist"] > shells[i - 1]["avg_dist"] + + def test_cn_sum_matches_requested_shells(self, cif_zno): + shell = 3 + ce = get_coord_env(cif_zno, shell=shell) + shells = ce[0][2] + expected_cn = 0 + distances = sorted(nb["dist"] for s in shells for nb in s["neighbors"]) + assert len(distances) == sum(s["cn"] for s in shells) + + def test_show_coord_env_output(self, cif_zno): + output = cap_show(None, cif_zno, shell=3) + assert "Coord. Env. from absorber atom: Zn" in output + assert "Shell" in output + assert "CN" in output + + def test_toler_groups_closely_spaced_distances(self, cif_zno): + ce_strict = get_coord_env(cif_zno, shell=10, toler=0.01) + ce_loose = get_coord_env(cif_zno, shell=10, toler=1.0) + n_shells_strict = len(ce_strict[0][2]) + n_shells_loose = len(ce_loose[0][2]) + # Larger tolerance should group into fewer (or equal) shells + assert n_shells_loose <= n_shells_strict + + +# ── XYZ (molecule) tests ───────────────────────────────────────── + + +class TestCoordEnvXYZ: + def test_shell1_returned(self, xyz_cuo6): + ce = get_coord_env(xyz_cuo6, shell=1) + shells = ce[0][2] + assert len(shells) >= 1 + assert shells[0]["shell_idx"] == 1 + + def test_shell1_cn_positive(self, xyz_cuo6): + ce = get_coord_env(xyz_cuo6, shell=1) + shells = ce[0][2] + assert shells[0]["cn"] > 0 + + def test_shell1_has_elements(self, xyz_cuo6): + ce = get_coord_env(xyz_cuo6, shell=1) + shells = ce[0][2] + assert len(shells[0]["elems"]) > 0 + + def test_shell1_neighbors_sorted(self, xyz_cuo6): + ce = get_coord_env(xyz_cuo6, shell=1) + shells = ce[0][2] + dists = [nb["dist"] for nb in shells[0]["neighbors"]] + assert dists == sorted(dists) + + def test_show_coord_env_output(self, xyz_cuo6): + output = cap_show(None, xyz_cuo6, shell=2) + assert "Coord. Env. from absorber atom: Cu" in output + assert "Shell" in output + assert "CN" in output From ccc3ac49e8351efb667a8486d947d8d90b7690e2 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 14:47:57 +0200 Subject: [PATCH 16/52] add tests visualize --- larixite/struct/visualize.py | 14 +-- tests/test_visualize.py | 204 +++++++++++++++++++++++++++++++++++ 2 files changed, 212 insertions(+), 6 deletions(-) create mode 100644 tests/test_visualize.py diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py index 8d030bb..385e183 100644 --- a/larixite/struct/visualize.py +++ b/larixite/struct/visualize.py @@ -140,6 +140,14 @@ def _round_up(x: float) -> float: return math.ceil(x * 100) / 100 +def _lighten_hex(hex_color: str, factor: float = 0.75) -> str: + """Return a darker hex color obtained by scaling each RGB channel by *factor*.""" + r = int(int(hex_color[1:3], 16) * factor) + g = int(int(hex_color[3:5], 16) * factor) + b = int(int(hex_color[5:7], 16) * factor) + return "#{:02x}{:02x}{:02x}".format(r, g, b) + + def _cluster_to_xyz(mol) -> tuple[str, list[str]]: """Convert a pymatgen Molecule to an XYZ format string and the list of unique elements.""" coords = [] @@ -201,12 +209,6 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) for elem in elems: color_elems[elem] = VESTA_COLORS.get(elem, "#888888") - def _lighten_hex(hex_color: str, factor: float = 0.75) -> str: - r = int(int(hex_color[1:3], 16) * factor) - g = int(int(hex_color[3:5], 16) * factor) - b = int(int(hex_color[5:7], 16) * factor) - return "#{:02x}{:02x}{:02x}".format(r, g, b) - covalent_r = CovalentRadius().radius absorber_sym = struct.absorber.symbol diff --git a/tests/test_visualize.py b/tests/test_visualize.py new file mode 100644 index 0000000..721bcf0 --- /dev/null +++ b/tests/test_visualize.py @@ -0,0 +1,204 @@ +"""Tests for larixite.struct.visualize — _cluster_to_xyz, _lighten_hex, VESTA_COLORS, visualize().""" + +import math +from pathlib import Path + +import pytest +import numpy as np + +from larixite.struct import get_structure +from larixite.struct.visualize import ( + HAS_PY3DMOL, + VESTA_COLORS, + _round_up, + _cluster_to_xyz, + _lighten_hex, + visualize, +) + +STRUCTS = Path(__file__).parent / "structs" + +# ── fixtures ────────────────────────────────────────────────────── + + +@pytest.fixture(scope="module") +def cif_zno(): + return get_structure(STRUCTS / "ZnO_mp-2133.cif", "Zn") + + +@pytest.fixture(scope="module") +def xyz_cuo6(): + return get_structure(STRUCTS / "CuO6_D4h.xyz", "Cu") + + +# ── _round_up ─────────────────────────────────────────────────── + + +class TestRoundUp: + def test_integer(self): + assert _round_up(2.0) == 2.0 + + def test_rounds_up(self): + assert _round_up(2.31) == 2.31 + + def test_rounds_up_small(self): + assert _round_up(2.001) == 2.01 + + def test_rounds_up_boundary(self): + assert _round_up(2.501) == 2.51 + + +# ── _cluster_to_xyz ───────────────────────────────────────────── + + +class TestClusterToXyz: + def test_returns_strings_and_elems(self, cif_zno): + mol = cif_zno.build_cluster(radius=2.5) + xyz_str, elems = _cluster_to_xyz(mol) + assert isinstance(xyz_str, str) + assert isinstance(elems, list) + assert all(isinstance(e, str) for e in elems) + + def test_xyz_header_matches_atom_count(self, cif_zno): + mol = cif_zno.build_cluster(radius=2.5) + xyz_str, elems = _cluster_to_xyz(mol) + header = int(xyz_str.strip().split("\n")[0]) + assert header == len(mol) + + def test_absorber_first(self, cif_zno): + mol = cif_zno.build_cluster(radius=2.5) + xyz_str, _ = _cluster_to_xyz(mol) + lines = xyz_str.strip().split("\n") + absorber = lines[2].split()[0] + assert absorber == cif_zno.absorber.symbol + + def test_coords_zero_origin(self, cif_zno): + mol = cif_zno.build_cluster(radius=2.5) + xyz_str, _ = _cluster_to_xyz(mol) + lines = xyz_str.strip().split("\n") + coords = lines[2].split()[1:] + for c in coords: + assert float(c) == pytest.approx(0.0, abs=1e-5) + + def test_elements_sorted(self, cif_zno): + mol = cif_zno.build_cluster(radius=2.5) + _, elems = _cluster_to_xyz(mol) + assert elems == sorted(elems) + + def test_elements_match_structure(self, cif_zno): + mol = cif_zno.build_cluster(radius=7.0) + _, elems = _cluster_to_xyz(mol) + assert "Zn" in elems + assert "O" in elems + + def test_xyz_single_element(self, xyz_cuo6): + mol = xyz_cuo6.build_cluster(radius=0.5) + xyz_str, elems = _cluster_to_xyz(mol) + assert len(elems) == 1 + assert elems == ["Cu"] + + +# ── _lighten_hex ──────────────────────────────────────────────── + + +class TestLightenHex: + def test_returns_hex(self): + assert _lighten_hex("#ff0000", 0.5).startswith("#") + + def test_darkens_correctly(self): + result = _lighten_hex("#ff0000", 0.5) + assert result == "#800000" + + def test_identity(self): + assert _lighten_hex("#b57100", 1.0) == "#b57100" + + def test_factor_0_8(self): + result = _lighten_hex("#ff0000", 0.8) + r = int(int("ff", 16) * 0.8) + assert "#" in result and len(result) == 7 + + def test_fe_color(self): + fe_hex = VESTA_COLORS["Fe"] + light = _lighten_hex(fe_hex, 0.8) + assert light != fe_hex + assert light.startswith("#") + + def test_all_vesta_colors_valid(self): + for sym, hex_color in VESTA_COLORS.items(): + light = _lighten_hex(hex_color, 0.8) + assert light.startswith("#") + assert len(light) == 7 + + def test_invalid_hex_unchanged(self): + with pytest.raises(ValueError): + _lighten_hex("notahex", 0.5) + + +# ── VESTA_COLORS ──────────────────────────────────────────────── + + +class TestVestaColors: + COMMON = ["Fe", "O", "C", "N", "Si", "Al", "Ti", "Cu", "Zn", "S", "P", "Mg", "Ca"] + + def test_common_elements_present(self): + for sym in self.COMMON: + assert sym in VESTA_COLORS + + def test_colors_are_hex(self): + for v in VESTA_COLORS.values(): + assert v.startswith("#") and len(v) == 7 + + def test_all_vesta_values_valid_int(self): + for v in VESTA_COLORS.values(): + int(v[1:3], 16) + int(v[3:5], 16) + int(v[5:7], 16) + + def test_no_black_in_palette(self): + assert "#000000" not in VESTA_COLORS.values() + + def test_color_count(self): + assert len(VESTA_COLORS) > 100 + + +# ── visualize ─────────────────────────────────────────────────── + + +class TestVisualize: + @pytest.mark.skipif(not HAS_PY3DMOL, reason="py3Dmol not installed") + def test_returns_view_cif(self, cif_zno): + result = visualize(cif_zno, radius=2.5) + assert result is not None + + @pytest.mark.skipif(not HAS_PY3DMOL, reason="py3Dmol not installed") + def test_returns_view_xyz(self, xyz_cuo6): + result = visualize(xyz_cuo6, radius=2.0) + assert result is not None + + @pytest.mark.skipif(not HAS_PY3DMOL, reason="py3Dmol not installed") + def test_default_radius(self, cif_zno): + result = visualize(cif_zno) + assert result is not None + + @pytest.mark.skipif(not HAS_PY3DMOL, reason="py3Dmol not installed") + def test_unitcell_flag(self, cif_zno): + result = visualize(cif_zno, radius=2.5, unitcell=True) + assert result is not None + + @pytest.mark.skipif(not HAS_PY3DMOL, reason="py3Dmol not installed") + def test_unitcell_xyz_no_error(self, xyz_cuo6): + result = visualize(xyz_cuo6, radius=2.0, unitcell=True) + assert result is not None + + @pytest.mark.skipif(not HAS_PY3DMOL, reason="py3Dmol not installed") + def test_radius_affects_cluster(self, cif_zno): + result = visualize(cif_zno, radius=1.5) + assert result is not None + + @pytest.mark.skipif(HAS_PY3DMOL, reason="py3Dmol is installed") + def test_no_py3dmol_returns_none(self, cif_zno, monkeypatch): + import larixite.struct.visualize as viz + monkeypatch.setattr(viz, "HAS_PY3DMOL", False) + monkeypatch.setattr(viz, "py3Dmol", None) + result = viz.visualize(cif_zno) + assert result is None From b03f65aca50d4e7c3b02fd0579c7fbd4da60e594 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 15:19:48 +0200 Subject: [PATCH 17/52] improve colors output --- larixite/struct/visualize.py | 267 ++++++++++++++++++++--------------- tests/test_visualize.py | 31 ++-- 2 files changed, 173 insertions(+), 125 deletions(-) diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py index 385e183..cd3b080 100644 --- a/larixite/struct/visualize.py +++ b/larixite/struct/visualize.py @@ -7,6 +7,7 @@ """ import math +import sys from pymatgen.analysis.local_env import CovalentRadius @@ -24,115 +25,115 @@ VESTA_COLORS = { - "Ac": "#70abfa", - "Ag": "#c0c0c0", - "Al": "#81b2d6", - "Am": "#545cf2", - "Ar": "#cffec4", - "As": "#74d057", - "At": "#754f45", - "Au": "#ffd123", - "B": "#1fa20f", - "Ba": "#00c900", - "Be": "#5ed77b", - "Bh": "#e00038", - "Bi": "#9e4fb5", - "Bk": "#8a4fe3", - "Br": "#7e3102", - "C": "#4c4c4c", - "Ca": "#5a96bd", - "Cd": "#ffd98f", - "Ce": "#ffffc7", - "Cf": "#a136d4", - "Cl": "#31fc02", - "Cm": "#785ce3", - "Co": "#0000af", - "Cr": "#00009e", - "Cs": "#57178f", - "Cu": "#2247dc", - "Db": "#d1004f", - "Dy": "#1fffc7", - "Er": "#00e675", - "Es": "#b31fd4", - "Eu": "#61ffc7", - "F": "#b0b9e6", - "Fe": "#b57100", - "Fm": "#b31fba", - "Fr": "#420066", - "Ga": "#9ee373", - "Gd": "#45ffc7", - "Ge": "#7e6ea6", - "H": "#ffcccc", - "He": "#fce8ce", - "Hf": "#4dc2ff", - "Hg": "#b8b8d0", - "Ho": "#00ff9c", - "Hs": "#e6002e", - "I": "#940094", - "In": "#a67573", - "Ir": "#175487", - "K": "#a121f6", - "Kr": "#fac1f3", - "La": "#5ac449", - "Li": "#86df73", - "Lr": "#c70066", - "Lu": "#00ab24", - "Md": "#b30da6", - "Mg": "#fb7b15", - "Mn": "#a7089d", - "Mo": "#54b5b5", - "Mt": "#eb0026", - "N": "#b0b9e6", - "Na": "#f9dc3c", - "Nb": "#73c2c9", - "Nd": "#c7ffc7", - "Ne": "#fe37b5", - "Ni": "#b7bbbd", - "No": "#bd0d87", - "Np": "#0080ff", - "O": "#fe0300", - "Os": "#266696", - "P": "#c09cc2", - "Pa": "#00a1ff", - "Pb": "#575961", - "Pd": "#006985", - "Pm": "#a3ffc7", - "Po": "#ab5c00", - "Pr": "#d9ffc7", - "Pt": "#d0d0e0", - "Pu": "#006bff", - "Ra": "#007d00", - "Rb": "#702eb0", - "Re": "#267dab", - "Rf": "#cc0059", - "Rh": "#0a7d8c", - "Rn": "#428296", - "Ru": "#248f8f", - "S": "#fffa00", - "Sb": "#9e63b5", - "Sc": "#b563ab", - "Se": "#9aef0f", - "Sg": "#d90045", - "Si": "#1b3bfa", - "Sm": "#8fffc7", - "Sn": "#9a8eb9", - "Sr": "#00ff00", - "Ta": "#4da6ff", - "Tb": "#30ffc7", - "Tc": "#3b9e9e", - "Te": "#d47a00", - "Th": "#00baff", - "Ti": "#78caff", - "Tl": "#a6544d", - "Tm": "#00d452", - "U": "#008fff", - "V": "#e51900", - "W": "#2194d6", - "Xe": "#429eb0", - "Y": "#94ffff", - "Yb": "#00bf38", - "Zn": "#8f8f81", - "Zr": "#00ff00", + "Ac": ("#70abfa", "sky blue"), + "Ag": ("#c0c0c0", "silver"), + "Al": ("#81b2d6", "light grey blue"), + "Am": ("#545cf2", "warm blue"), + "Ar": ("#cffec4", "very light green"), + "As": ("#74d057", "fresh green"), + "At": ("#754f45", "brownish purple"), + "Au": ("#ffd123", "sun yellow"), + "B": ("#1fa20f", "green"), + "Ba": ("#00c900", "vibrant green"), + "Be": ("#5ed77b", "soft green"), + "Bh": ("#e00038", "cherry"), + "Bi": ("#9e4fb5", "purply"), + "Bk": ("#8a4fe3", "purple"), + "Br": ("#7e3102", "reddish brown"), + "C": ("#4c4c4c", "charcoal"), + "Ca": ("#5a96bd", "faded blue"), + "Cd": ("#ffd98f", "wheat"), + "Ce": ("#ffffc7", "ecru"), + "Cf": ("#a136d4", "dark orchid"), + "Cl": ("#31fc02", "vivid green"), + "Cm": ("#785ce3", "purple"), + "Co": ("#0000af", "royal blue"), + "Cr": ("#00009e", "royal blue"), + "Cs": ("#57178f", "indigo"), + "Cu": ("#2247dc", "true blue"), + "Db": ("#d1004f", "ruby"), + "Dy": ("#1fffc7", "greenish cyan"), + "Er": ("#00e675", "tealish green"), + "Es": ("#b31fd4", "magenta purple"), + "Eu": ("#61ffc7", "light greenish blue"), + "F": ("#b0b9e6", "light steel blue"), + "Fe": ("#b57100", "brown orange"), + "Fm": ("#b31fba", "magenta purple"), + "Fr": ("#420066", "royal purple"), + "Ga": ("#9ee373", "spring green"), + "Gd": ("#45ffc7", "greenish cyan"), + "Ge": ("#7e6ea6", "dark lavender"), + "H": ("#ffcccc", "light rose"), + "He": ("#fce8ce", "blanched almond"), + "Hf": ("#4dc2ff", "sky blue"), + "Hg": ("#b8b8d0", "cloudy blue"), + "Ho": ("#00ff9c", "medium spring green"), + "Hs": ("#e6002e", "cherry red"), + "I": ("#940094", "dark magenta"), + "In": ("#a67573", "reddish grey"), + "Ir": ("#175487", "light navy"), + "K": ("#a121f6", "electric purple"), + "Kr": ("#fac1f3", "light lavender"), + "La": ("#5ac449", "leafy green"), + "Li": ("#86df73", "light green"), + "Lr": ("#c70066", "deep pink"), + "Lu": ("#00ab24", "kelly green"), + "Md": ("#b30da6", "magenta"), + "Mg": ("#fb7b15", "pumpkin orange"), + "Mn": ("#a7089d", "barney purple"), + "Mo": ("#54b5b5", "medium aquamarine"), + "Mt": ("#eb0026", "cherry red"), + "N": ("#b0b9e6", "light steel blue"), + "Na": ("#f9dc3c", "off yellow"), + "Nb": ("#73c2c9", "seafoam blue"), + "Nd": ("#c7ffc7", "very pale green"), + "Ne": ("#fe37b5", "hot pink"), + "Ni": ("#b7bbbd", "silver"), + "No": ("#bd0d87", "medium violet red"), + "Np": ("#0080ff", "deep sky blue"), + "O": ("#fe0300", "red"), + "Os": ("#266696", "ugly blue"), + "P": ("#c09cc2", "pale purple"), + "Pa": ("#00a1ff", "azure"), + "Pb": ("#575961", "gunmetal"), + "Pd": ("#006985", "peacock blue"), + "Pm": ("#a3ffc7", "light seafoam"), + "Po": ("#ab5c00", "orange brown"), + "Pr": ("#d9ffc7", "very light green"), + "Pt": ("#d0d0e0", "light grey"), + "Pu": ("#006bff", "bright blue"), + "Ra": ("#007d00", "green"), + "Rb": ("#702eb0", "rebecca purple"), + "Re": ("#267dab", "bluish"), + "Rf": ("#cc0059", "deep pink"), + "Rh": ("#0a7d8c", "ocean"), + "Rn": ("#428296", "dirty blue"), + "Ru": ("#248f8f", "blue green"), + "S": ("#fffa00", "bright yellow"), + "Sb": ("#9e63b5", "amethyst"), + "Sc": ("#b563ab", "soft purple"), + "Se": ("#9aef0f", "acid green"), + "Sg": ("#d90045", "ruby"), + "Si": ("#1b3bfa", "vivid blue"), + "Sm": ("#8fffc7", "light seafoam"), + "Sn": ("#9a8eb9", "heather"), + "Sr": ("#00ff00", "lime"), + "Ta": ("#4da6ff", "cornflower blue"), + "Tb": ("#30ffc7", "greenish cyan"), + "Tc": ("#3b9e9e", "sea"), + "Te": ("#d47a00", "pumpkin"), + "Th": ("#00baff", "deep sky blue"), + "Ti": ("#78caff", "light blue"), + "Tl": ("#a6544d", "light maroon"), + "Tm": ("#00d452", "shamrock green"), + "U": ("#008fff", "azure"), + "V": ("#e51900", "tomato red"), + "W": ("#2194d6", "water blue"), + "Xe": ("#429eb0", "cool blue"), + "Y": ("#94ffff", "robin egg blue"), + "Yb": ("#00bf38", "shamrock green"), + "Zn": ("#8f8f81", "warm grey"), + "Zr": ("#00ff00", "lime"), } @@ -148,6 +149,48 @@ def _lighten_hex(hex_color: str, factor: float = 0.75) -> str: return "#{:02x}{:02x}{:02x}".format(r, g, b) +def _print_color_legend(color_elems: dict[str, str]) -> None: + """Print a human-readable color legend, adapting to the execution environment. + + In Jupyter: renders HTML color swatches. + In a terminal: uses ANSI 24-bit color escape codes. + Fallback: plain text with element names and color descriptions. + """ + try: + from IPython.display import display, HTML + + swatches = " ".join( + f'' + f'{e}' + for e, c in color_elems.items() + ) + display(HTML(f'
Colors: {swatches}
')) + return + except ImportError: + pass + + if hasattr(sys.stdout, "isatty") and sys.stdout.isatty() and sys.stdout.encoding == "UTF-8": + for elem, hex_color in color_elems.items(): + r = int(hex_color[1:3], 16) + g = int(hex_color[3:5], 16) + b = int(hex_color[5:7], 16) + sys.stdout.write("\x1b[38;2;{};{};{};m{}\x1b[0m ".format(r, g, b, elem)) + sys.stdout.write("\n") + return + + legend = ", ".join(f"{e}: {VESTA_COLORS.get(e, ('#888888', 'grey'))[1]}" for e in color_elems) + logger.info(f"Colors: {legend}") + + +def _adjust_font_color(hex_color: str, dark: str = "#222222", light: str = "#ffffff") -> str: + """Return dark or light text color for contrast against the given background hex color.""" + r = int(hex_color[1:3], 16) + g = int(hex_color[3:5], 16) + b = int(hex_color[5:7], 16) + luminance = (0.299 * r + 0.587 * g + 0.114 * b) / 255 + return dark if luminance > 0.55 else light + + def _cluster_to_xyz(mol) -> tuple[str, list[str]]: """Convert a pymatgen Molecule to an XYZ format string and the list of unique elements.""" coords = [] @@ -207,14 +250,14 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) color_elems = {} for elem in elems: - color_elems[elem] = VESTA_COLORS.get(elem, "#888888") + color_elems[elem] = VESTA_COLORS.get(elem, ("#888888", "grey"))[0] covalent_r = CovalentRadius().radius absorber_sym = struct.absorber.symbol absorber_color = color_elems.get(absorber_sym, "#888888") absorber_r = max(0.2, min(covalent_r.get(absorber_sym, 0.5) / 2, 0.75)) - absorber_light = _lighten_hex(absorber_color, 0.8) + absorber_light = _lighten_hex(absorber_color, 0.6) for elem in elems: color = color_elems[elem] @@ -250,5 +293,5 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) xyzview.zoomTo() xyzview.show() - logger.info(color_elems) + _print_color_legend(color_elems) return xyzview diff --git a/tests/test_visualize.py b/tests/test_visualize.py index 721bcf0..97fd01f 100644 --- a/tests/test_visualize.py +++ b/tests/test_visualize.py @@ -107,7 +107,8 @@ def test_returns_hex(self): def test_darkens_correctly(self): result = _lighten_hex("#ff0000", 0.5) - assert result == "#800000" + # int(255 * 0.5) == 127 (truncation), not 128 + assert result == "#7f0000" def test_identity(self): assert _lighten_hex("#b57100", 1.0) == "#b57100" @@ -118,16 +119,19 @@ def test_factor_0_8(self): assert "#" in result and len(result) == 7 def test_fe_color(self): - fe_hex = VESTA_COLORS["Fe"] + fe_hex, fe_name = VESTA_COLORS["Fe"] light = _lighten_hex(fe_hex, 0.8) assert light != fe_hex assert light.startswith("#") + assert fe_name == "brown orange" - def test_all_vesta_colors_valid(self): - for sym, hex_color in VESTA_COLORS.items(): + def test_all_vesta_colors_tuplified(self): + for sym, entry in VESTA_COLORS.items(): + assert isinstance(entry, tuple) and len(entry) == 2 + hex_color, name = entry light = _lighten_hex(hex_color, 0.8) - assert light.startswith("#") - assert len(light) == 7 + assert light.startswith("#") and len(light) == 7 + assert isinstance(name, str) and name def test_invalid_hex_unchanged(self): with pytest.raises(ValueError): @@ -145,17 +149,18 @@ def test_common_elements_present(self): assert sym in VESTA_COLORS def test_colors_are_hex(self): - for v in VESTA_COLORS.values(): - assert v.startswith("#") and len(v) == 7 + for hex_color, name in VESTA_COLORS.values(): + assert hex_color.startswith("#") and len(hex_color) == 7 + assert isinstance(name, str) and name def test_all_vesta_values_valid_int(self): - for v in VESTA_COLORS.values(): - int(v[1:3], 16) - int(v[3:5], 16) - int(v[5:7], 16) + for hex_color, name in VESTA_COLORS.values(): + int(hex_color[1:3], 16) + int(hex_color[3:5], 16) + int(hex_color[5:7], 16) def test_no_black_in_palette(self): - assert "#000000" not in VESTA_COLORS.values() + assert "#000000" not in [h for h, n in VESTA_COLORS.values()] def test_color_count(self): assert len(VESTA_COLORS) > 100 From 53f6201c2a086ba76b65baf41a69102c293de855 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 22 May 2026 15:47:47 +0200 Subject: [PATCH 18/52] return outsbatch --- larixite/fdmnes/input.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 6b79df5..5bef079 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -408,7 +408,7 @@ def write_input( _ = self.dump_params(jobdir / f"{self.fileout_prefix}_params.yaml") return jobdir - def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | None = None, ncpus: int = 8): + def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | None = None, ncpus: int = 8) -> Path: """Generates a SBATCH file (SLURM workload manager) using a template Arguments @@ -443,6 +443,7 @@ def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | with open(sbatchout, "w") as fp, open(template) as tp: fp.write(tp.read().format(**kwargs)) logger.info(f"written {fp.name}") + return sbatchout def dump_params(self, yamlpath: str | Path | None = None) -> Path: """Dump input parameters to a YAML file From a58765205e2c8523792af6839a78cfe95a81dd96 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Tue, 26 May 2026 13:22:18 +0200 Subject: [PATCH 19/52] start parsing output sims --- larixite/fdmnes/__init__.py | 1 + larixite/fdmnes/output.py | 39 +++++++++++++++++++++++++++++++++++++ 2 files changed, 40 insertions(+) create mode 100644 larixite/fdmnes/output.py diff --git a/larixite/fdmnes/__init__.py b/larixite/fdmnes/__init__.py index 3321fb8..3f750e7 100644 --- a/larixite/fdmnes/__init__.py +++ b/larixite/fdmnes/__init__.py @@ -1 +1,2 @@ from .input import FdmnesXasInput, struct2fdmnes +from .output import search_jobs \ No newline at end of file diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py new file mode 100644 index 0000000..24f4eab --- /dev/null +++ b/larixite/fdmnes/output.py @@ -0,0 +1,39 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +Search FDMNES job directories. +""" + +from pathlib import Path + +def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[Path]: + """Search for FDMNES job directories matching a glob pattern. + + Scans the given parent directory for subdirectories whose names match + the provided glob expression. + + Parameters + ---------- + globstr : str + Glob pattern applied to the *name* of each immediate subdirectory + (e.g. `"job_*"`, `"mi_*"`, `"*.done"`). + jobsdir : Path | str | None -> ~/.larixite/fdmnes + Parent directory containing job subdirectories. + + Returns + ------- + list[Path] + Sorted list of matching subdirectory paths. + """ + if jobsdir is None: + jobsdir = Path("~/.larixite/fdmnes").expanduser().resolve() + if isinstance(jobsdir, str): + jobsdir = Path(jobsdir).expanduser().resolve() + if not jobsdir.is_dir(): + return [] + + return sorted( + (entry for entry in jobsdir.glob(globstr) if entry.is_dir()), + key=lambda x: x.stat().st_ctime, + ) From 7fd4d02a5544519f33ddfb869892d078f082803b Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Tue, 26 May 2026 13:46:21 +0200 Subject: [PATCH 20/52] restore input --- larixite/fdmnes/output.py | 54 +++++++++++++++++++++++++++++++-------- 1 file changed, 44 insertions(+), 10 deletions(-) diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py index 24f4eab..e5d0175 100644 --- a/larixite/fdmnes/output.py +++ b/larixite/fdmnes/output.py @@ -5,26 +5,49 @@ Search FDMNES job directories. """ +from dataclasses import dataclass from pathlib import Path -def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[Path]: +from larixite.fdmnes.input import FdmnesXasInput +from larixite.utils import get_logger + +logger = get_logger("larixite.fdmnes.output") +logger.propagate = False + + +@dataclass +class FdmnesXasSim: + """A FDMNES XAS simulation associated with a job directory.""" + + jobdir: Path + input: FdmnesXasInput | None = None + + @property + def prefix(self) -> str | None: + if self.input is not None: + return self.input.fileout_prefix + return None + + +def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[FdmnesXasSim]: """Search for FDMNES job directories matching a glob pattern. - Scans the given parent directory for subdirectories whose names match - the provided glob expression. + For each matching subdirectory that contains a ``*_params.yaml`` file, + a :class:`FdmnesXasSim` instance is built via + :meth:`FdmnesXasInput.from_yaml`. Parameters ---------- globstr : str Glob pattern applied to the *name* of each immediate subdirectory - (e.g. `"job_*"`, `"mi_*"`, `"*.done"`). + (e.g. ``"job_*"``, ``"mi_*"``, ``"*.done"``). jobsdir : Path | str | None -> ~/.larixite/fdmnes Parent directory containing job subdirectories. Returns ------- - list[Path] - Sorted list of matching subdirectory paths. + list[FdmnesXasSim] + Sorted list of matching simulation objects. """ if jobsdir is None: jobsdir = Path("~/.larixite/fdmnes").expanduser().resolve() @@ -33,7 +56,18 @@ def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[Path]: if not jobsdir.is_dir(): return [] - return sorted( - (entry for entry in jobsdir.glob(globstr) if entry.is_dir()), - key=lambda x: x.stat().st_ctime, - ) + results: list[FdmnesXasSim] = [] + for entry in jobsdir.glob(globstr): + if not entry.is_dir(): + continue + + for param_file in entry.glob("*_params.yaml"): + try: + inp = FdmnesXasInput.from_yaml(param_file) + results.append(FdmnesXasSim(jobdir=entry, input=inp)) + except Exception: + continue + + results = sorted(results, key=lambda s: s.jobdir.stat().st_ctime) + logger.info(f"found {len(results)} job(s) in {jobsdir}") + return results From 0764fd6f15c1b4a8712ac87da92ccf8d2667d353 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Tue, 26 May 2026 14:45:18 +0200 Subject: [PATCH 21/52] avoid multiple messages --- larixite/utils.py | 1 + 1 file changed, 1 insertion(+) diff --git a/larixite/utils.py b/larixite/utils.py index 1789339..44b5afe 100644 --- a/larixite/utils.py +++ b/larixite/utils.py @@ -184,6 +184,7 @@ def get_logger(name="larixite", level="INFO"): """Get a custom logger for larixite""" logger = logging.getLogger(name) logger.setLevel(log_levels[level]) + logger.propagate = False formatter = ColorFormatter() handler = logging.StreamHandler() handler.setFormatter(formatter) From f17ac8e6eccd7335adcdd991acfa6a59836bd49d Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Tue, 26 May 2026 14:45:31 +0200 Subject: [PATCH 22/52] WIP: parse output --- larixite/fdmnes/output.py | 40 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py index e5d0175..0dda223 100644 --- a/larixite/fdmnes/output.py +++ b/larixite/fdmnes/output.py @@ -8,11 +8,12 @@ from dataclasses import dataclass from pathlib import Path +import yaml +import pandas as pd from larixite.fdmnes.input import FdmnesXasInput from larixite.utils import get_logger logger = get_logger("larixite.fdmnes.output") -logger.propagate = False @dataclass @@ -21,6 +22,8 @@ class FdmnesXasSim: jobdir: Path input: FdmnesXasInput | None = None + data: pd.DataFrame | None = None + metadata: dict | None = None @property def prefix(self) -> str | None: @@ -28,6 +31,41 @@ def prefix(self) -> str | None: return self.input.fileout_prefix return None + def load_data(self) -> pd.DataFrame | None: + datafile = self.jobdir / f"{self.prefix}.txt" + if not datafile.exists(): + error_msg = "cannot find the calculation output file" + logger.error(error_msg) + error_file = self.jobdir / "error.yaml" + error_file.write_text( + yaml.dump({"error": error_msg}) + ) + return + + lines = datafile.read_text().splitlines() + e_shift: float = 0.0 + skip = 0 + + for i, line in enumerate(lines): + stripped = line.strip() + if "=" in stripped: + eq_split = stripped.split(" = ", 1) + left_vals = [float(v) for v in eq_split[0].split() if v] + right_names = [n.strip() for n in eq_split[1].split(",")] + int_names = {"Z", "n_edge", "j_edge", "ninit", "ninit1", "natomsym_f"} + self.metadata = { + name: int(val) if name in int_names else val + for name, val in zip(right_names, left_vals) + } + e_shift = self.metadata.get("E_edge", 0.0) + skip = i + 2 + elif stripped.startswith("#"): + continue + else: + break + + self.data = pd.read_csv(datafile, sep=r"\s+", skiprows=skip) + self.data.iloc[:, 0] += e_shift def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[FdmnesXasSim]: """Search for FDMNES job directories matching a glob pattern. From 45a601f549883ea4ff5a9a9ccf3e1d8df3257a7a Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 06:13:13 +0200 Subject: [PATCH 23/52] fix parsing output --- larixite/fdmnes/output.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py index 0dda223..d7ac13e 100644 --- a/larixite/fdmnes/output.py +++ b/larixite/fdmnes/output.py @@ -31,7 +31,7 @@ def prefix(self) -> str | None: return self.input.fileout_prefix return None - def load_data(self) -> pd.DataFrame | None: + def load_data(self, shift_energy: bool = False) -> pd.DataFrame | None: datafile = self.jobdir / f"{self.prefix}.txt" if not datafile.exists(): error_msg = "cannot find the calculation output file" @@ -44,28 +44,31 @@ def load_data(self) -> pd.DataFrame | None: lines = datafile.read_text().splitlines() e_shift: float = 0.0 - skip = 0 + skiprows = 1 for i, line in enumerate(lines): stripped = line.strip() - if "=" in stripped: + if stripped.startswith("#"): + continue + elif "=" in stripped: eq_split = stripped.split(" = ", 1) left_vals = [float(v) for v in eq_split[0].split() if v] right_names = [n.strip() for n in eq_split[1].split(",")] + assert len(left_vals) == len(right_names), "Error parsing the metadata line" int_names = {"Z", "n_edge", "j_edge", "ninit", "ninit1", "natomsym_f"} self.metadata = { name: int(val) if name in int_names else val for name, val in zip(right_names, left_vals) } e_shift = self.metadata.get("E_edge", 0.0) - skip = i + 2 - elif stripped.startswith("#"): - continue + skiprows = i + 1 + print(skiprows) else: break - self.data = pd.read_csv(datafile, sep=r"\s+", skiprows=skip) - self.data.iloc[:, 0] += e_shift + self.data = pd.read_csv(datafile, sep=r"\s+", skiprows=skiprows) + if shift_energy: + self.data.iloc[:, 0] += e_shift def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[FdmnesXasSim]: """Search for FDMNES job directories matching a glob pattern. From 033977e1f2159e80762710464a539a1ee70d4b86 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 06:20:20 +0200 Subject: [PATCH 24/52] add e0/efermi proprierties --- larixite/fdmnes/output.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py index d7ac13e..0eff926 100644 --- a/larixite/fdmnes/output.py +++ b/larixite/fdmnes/output.py @@ -31,6 +31,20 @@ def prefix(self) -> str | None: return self.input.fileout_prefix return None + @property + def efermi(self) -> float: + if self.metadata is not None: + return self.metadata.get("VO_interstitial", 0) + logger.warning("metadata not available, run `load_data()` first") + return 0 + + @property + def e0(self) -> float: + if self.metadata is not None: + return self.metadata.get("E_edge", 0) + logger.warning("metadata not available, , run `load_data()` first") + return 0 + def load_data(self, shift_energy: bool = False) -> pd.DataFrame | None: datafile = self.jobdir / f"{self.prefix}.txt" if not datafile.exists(): From dcc1e1159978cf8180d42a16ca819ec418096f83 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 06:40:13 +0200 Subject: [PATCH 25/52] add energy/mu properties --- larixite/fdmnes/output.py | 48 ++++++++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 18 deletions(-) diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py index 0eff926..d1f6ea1 100644 --- a/larixite/fdmnes/output.py +++ b/larixite/fdmnes/output.py @@ -8,6 +8,7 @@ from dataclasses import dataclass from pathlib import Path +import numpy as np import yaml import pandas as pd from larixite.fdmnes.input import FdmnesXasInput @@ -33,31 +34,40 @@ def prefix(self) -> str | None: @property def efermi(self) -> float: - if self.metadata is not None: - return self.metadata.get("VO_interstitial", 0) - logger.warning("metadata not available, run `load_data()` first") - return 0 + if self.metadata is None: + self.load() + return self.metadata.get("VO_interstitial", 0) @property def e0(self) -> float: - if self.metadata is not None: - return self.metadata.get("E_edge", 0) - logger.warning("metadata not available, , run `load_data()` first") - return 0 + if self.metadata is None: + self.load() + return self.metadata.get("E_edge", 0) - def load_data(self, shift_energy: bool = False) -> pd.DataFrame | None: + @property + def energy(self, erel: bool = False) -> np.ndarray: + if self.data is None: + self.load() + eout = self.data["Energy"].values + return eout + self.e0 if erel else eout + + @property + def mu(self) -> np.ndarray: + if self.data is None: + self.load() + return self.data[""].values + + def load(self) -> pd.DataFrame | None: + """Load the FDMNES output data file and parse the metadata line.""" datafile = self.jobdir / f"{self.prefix}.txt" if not datafile.exists(): error_msg = "cannot find the calculation output file" logger.error(error_msg) error_file = self.jobdir / "error.yaml" - error_file.write_text( - yaml.dump({"error": error_msg}) - ) + error_file.write_text(yaml.dump({"error": error_msg})) return lines = datafile.read_text().splitlines() - e_shift: float = 0.0 skiprows = 1 for i, line in enumerate(lines): @@ -68,21 +78,23 @@ def load_data(self, shift_energy: bool = False) -> pd.DataFrame | None: eq_split = stripped.split(" = ", 1) left_vals = [float(v) for v in eq_split[0].split() if v] right_names = [n.strip() for n in eq_split[1].split(",")] - assert len(left_vals) == len(right_names), "Error parsing the metadata line" + assert len(left_vals) == len(right_names), ( + "Error parsing the metadata line" + ) int_names = {"Z", "n_edge", "j_edge", "ninit", "ninit1", "natomsym_f"} self.metadata = { name: int(val) if name in int_names else val for name, val in zip(right_names, left_vals) } - e_shift = self.metadata.get("E_edge", 0.0) skiprows = i + 1 - print(skiprows) else: break self.data = pd.read_csv(datafile, sep=r"\s+", skiprows=skiprows) - if shift_energy: - self.data.iloc[:, 0] += e_shift + assert len(self.data.columns) >= 2, "Output file must have at least two columns" + assert "Energy" in self.data.columns, "No 'Energy' column in the output file" + assert "" in self.data.columns, "No '' column in the output file" + def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[FdmnesXasSim]: """Search for FDMNES job directories matching a glob pattern. From 54aa5bfa37205dad77f494d17e5a956746e12a33 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 06:53:57 +0200 Subject: [PATCH 26/52] refactor search_jobs returning dict with only inputs from the _params.yaml --- larixite/fdmnes/output.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py index d1f6ea1..7685823 100644 --- a/larixite/fdmnes/output.py +++ b/larixite/fdmnes/output.py @@ -96,11 +96,11 @@ def load(self) -> pd.DataFrame | None: assert "" in self.data.columns, "No '' column in the output file" -def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[FdmnesXasSim]: +def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> dict[str, FdmnesXasSim]: """Search for FDMNES job directories matching a glob pattern. - For each matching subdirectory that contains a ``*_params.yaml`` file, - a :class:`FdmnesXasSim` instance is built via + For each matching subdirectory that contains a ``{fileout_prefix}_params.yaml`` + file, a :class:`FdmnesXasSim` instance is built via :meth:`FdmnesXasInput.from_yaml`. Parameters @@ -113,17 +113,17 @@ def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[FdmnesX Returns ------- - list[FdmnesXasSim] - Sorted list of matching simulation objects. + dict[str, FdmnesXasSim] + Dictionary mapping ``fileout_prefix`` to matching simulation objects. """ if jobsdir is None: jobsdir = Path("~/.larixite/fdmnes").expanduser().resolve() if isinstance(jobsdir, str): jobsdir = Path(jobsdir).expanduser().resolve() if not jobsdir.is_dir(): - return [] + return {} - results: list[FdmnesXasSim] = [] + results: dict[str, FdmnesXasSim] = {} for entry in jobsdir.glob(globstr): if not entry.is_dir(): continue @@ -131,10 +131,13 @@ def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> list[FdmnesX for param_file in entry.glob("*_params.yaml"): try: inp = FdmnesXasInput.from_yaml(param_file) - results.append(FdmnesXasSim(jobdir=entry, input=inp)) + expected = entry / f"{inp.fileout_prefix}_params.yaml" + if param_file.resolve() != expected.resolve(): + continue + sim = FdmnesXasSim(jobdir=entry, input=inp) + results[inp.fileout_prefix] = sim except Exception: continue - results = sorted(results, key=lambda s: s.jobdir.stat().st_ctime) logger.info(f"found {len(results)} job(s) in {jobsdir}") return results From 3c0915b17d41fd4766ed1b0dc3264185830131da Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 07:03:48 +0200 Subject: [PATCH 27/52] minor docstring --- larixite/fdmnes/output.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py index 7685823..531719c 100644 --- a/larixite/fdmnes/output.py +++ b/larixite/fdmnes/output.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ -Search FDMNES job directories. +Handle FDMNES output files. """ from dataclasses import dataclass From 918c690cc2d413db77f3376a8589d64e4fbfd832 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 07:11:25 +0200 Subject: [PATCH 28/52] start models --- larixite/models/__init__.py | 3 +++ larixite/models/xas.py | 38 +++++++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+) create mode 100644 larixite/models/__init__.py create mode 100644 larixite/models/xas.py diff --git a/larixite/models/__init__.py b/larixite/models/__init__.py new file mode 100644 index 0000000..14ae1c9 --- /dev/null +++ b/larixite/models/__init__.py @@ -0,0 +1,3 @@ +"""Models for data handled by Larixite""" + +from .xas import XasExpData \ No newline at end of file diff --git a/larixite/models/xas.py b/larixite/models/xas.py new file mode 100644 index 0000000..50621c0 --- /dev/null +++ b/larixite/models/xas.py @@ -0,0 +1,38 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +"""Experimental XAS data model.""" + +from dataclasses import dataclass +from pathlib import Path + +import pandas as pd +from larixite.utils import get_logger + +logger = get_logger("larixite.models.xas") + + +@dataclass +class XasExpData: + """Model for an Experimental XAS spectrum.""" + + data: pd.DataFrame | None = None + + @classmethod + def from_csv(cls, csvpath: Path | str) -> "XasExpData": + """Load experimental XAS data from a CSV file. + + Parameters + ---------- + csvpath : Path | str + Path to the CSV file containing the XAS data. + + Returns + ------- + XasExpData + Instance with the data loaded from the CSV. + """ + if isinstance(csvpath, str): + csvpath = Path(csvpath) + data = pd.read_csv(csvpath) + return cls(data=data) From d3a2c7ba0fe947e12bc05fa14bc2a98f8ee0c8fd Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 17:04:48 +0200 Subject: [PATCH 29/52] add note on polarization --- larixite/fdmnes/input.py | 2 +- larixite/templates/fdmnes_xas.tmpl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 5bef079..de9c314 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -27,7 +27,7 @@ FDMNES_DEFAULT_PARAMS = { #: "FDMNES key name": True/False "Energpho": False, "Quadrupole": False, - "Polarize": False, #: -> TODO + "Polarize": False, #: -> TODO: implement giving polarization vectors "Density": False, "Density_all": False, "Green": True, diff --git a/larixite/templates/fdmnes_xas.tmpl b/larixite/templates/fdmnes_xas.tmpl index 2e9d901..6a8cfec 100755 --- a/larixite/templates/fdmnes_xas.tmpl +++ b/larixite/templates/fdmnes_xas.tmpl @@ -22,8 +22,8 @@ Radius ! {Quadrupole} ! -{Polarize} -! =>TODO +{Polarize} ! without giving a vector it calculates non equivalent directions +! ! =>TODO: specify polarization and wave vectors {TDDFT} {Memory_save} {Relativism} From ce4062360bfa7ae95d189023e21f782fda6cdd2b Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 17:26:38 +0200 Subject: [PATCH 30/52] start math module, import convolution from Larch (TODO: move here to mantain only one version) --- larixite/math/__init__.py | 0 larixite/math/convolution.py | 216 +++++++++++++++++++++++++++++++++++ 2 files changed, 216 insertions(+) create mode 100644 larixite/math/__init__.py create mode 100644 larixite/math/convolution.py diff --git a/larixite/math/__init__.py b/larixite/math/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/larixite/math/convolution.py b/larixite/math/convolution.py new file mode 100644 index 0000000..a6d5bc0 --- /dev/null +++ b/larixite/math/convolution.py @@ -0,0 +1,216 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +Generic discrete (one-dimensional) convolution + +Description +----------- + + This is a manual (not optimized!) implementation of discrete one-dimensional + (1D) convolution intended for spectroscopy analysis. The difference with + commonly used methods is the possibility to adapt the width of the convolution + kernel for each convolution point, e.g. change the FWHM of the Gaussian kernel + as a function of the energy scale (x-axis). + +Resources +--------- + +.. [WPconv] +.. [Fisher] +.. [GP1202] + +""" +import numpy as np +from .lineshapes import gaussian, lorentzian + +DEBUG = 0 + +def polyfit(x, y, deg=1, reverse=False): + """ + simple emulation of deprecated numpy.polyfit, + including its ordering of coefficients + """ + pfit = np.polynomial.Polynomial.fit(x, y, deg=int(deg)) + coefs = pfit.convert().coef + if reverse: + coefs = list(reversed(coefs)) + return list(coefs) + +def get_ene_index(ene, cen, hwhm): + """returns the min/max indexes for array ene at (cen-hwhm) and (cen+hwhm) + very similar to index_of in larch + """ + try: + if (cen - hwhm) <= min(ene): + ene_imin = 0 + else: + ene_imin = max(np.where(ene < (cen - hwhm))[0]) + if (cen + hwhm) >= max(ene): + ene_imax = len(e) - 1 + else: + ene_imax = min(np.where(ene > (cen + hwhm))[0]) + return ene_imin, ene_imax + except Exception: + print("index not found for {0} +/- {1}".format(cen, hwhm)) + return None, None + + +def lin_gamma(ene, gamma_hole=0.5, linbroad=None): + """returns constant or linear energy-dependent broadening + + Parameters + ---------- + ene : energy array in eV + gamma_hole : initial full width at half maximum in eV + linbroad : list of 3-elements giving + 'final full width at half maximum' + 'start energy point of the linear increase' + 'end energy point of the linear increase' + """ + w = np.ones_like(ene) + if linbroad is None: + return w * gamma_hole + else: + try: + fwhm2 = linbroad[0] + e1 = linbroad[1] + e2 = linbroad[2] + except: + raise ValueError("wrong format for linbroad") + for en, ee in enumerate(ene): + if ee < e1: + w[en] *= gamma_hole + elif ee <= e2: + wlin = gamma_hole + (ee - e1) * (fwhm2 - gamma_hole) / (e2 - e1) + w[en] *= wlin + elif ee >= e2: + w[en] *= fwhm2 + return w + + +def atan_gamma(e, e_cut=0, e_cent=30, e_larg=30, gamma_hole=0.5, gamma_max=15): + """Arctangent energy dependent broadening as implemented in FDMNES (described in the "Convolution" section of the manual)""" + + f = (e - e_cut) / e_cent + a = np.pi * gamma_max * (f - 1 / f**2) / (3 * e_larg) + gammas = gamma_hole + gamma_max * (0.5 + 1 / np.pi * np.arctan(a)) + + # Set gamma to the gamma_hole below the cutoff energy. + mask = np.where(e < e_cut) + gammas[mask] = gamma_hole + + return gammas + + +def conv_fast(x, y, gammas, e_cut=None, kernel="gaussian", num=501, step=0.1): + """A significantly faster version of the `conv` function + + Parameters + ---------- + x : x-axis (energy) + y : f(x) to convolve with g(x) kernel, y(energy) + kernel : convolution kernel, g(x) + 'gaussian' + 'lorentzian' + gammas : the full width half maximum in eV for the kernel + broadening. It is an array of size 'e' with constants or + an energy-dependent values determined by a function as + 'lin_gamma()' or 'atan_gamma()' + """ + assert e_cut is not None, "starting energy for the convolution not given" + assert "lor" or "gauss" in kernel.lower(), ( + "the kernel should be either Lorentzian or Gaussian" + ) + + # Extend the X-axis array. + start = x[-1] + stop = start + num * step + x_ext = np.append(x, np.arange(start + step, stop, step)) + + ids = x_ext > e_cut + + x1 = x_ext[ids] + x1[0] = e_cut + + # Extend the intensity array by coping the last value `num - 1` times. + y = y[x > e_cut] + y = np.append(y, np.ones(num - 1) * y[-1]) + + y_conv = np.zeros_like(x) + for i, (xi, gamma) in enumerate(zip(x, gammas)): + gamma = gamma / 2.0 + if "gauss" in kernel.lower(): + ky = gaussian(x1, center=xi, sigma=gamma) + elif "lor" in kernel.lower(): + ky = lorentzian(x1, center=xi, sigma=gamma) + y_conv[i] = np.sum(ky * y) / np.pi + + return y_conv + + +def conv(x, y, gammas, e_cut=None, kernel="gaussian"): + """linear broadening + + Parameters + ---------- + x : x-axis (energy) + y : f(x) to convolve with g(x) kernel, y(energy) + kernel : convolution kernel, g(x) + 'gaussian' + 'lorentzian' + gammas : the full width half maximum in eV for the kernel + broadening. It is an array of size 'e' with constants or + an energy-dependent values determined by a function as + 'lin_gamma()' or 'atan_gamma()' + """ + assert e_cut is not None, "starting energy for the convolution not given" + + f = y[:] * 1.0 + z = np.zeros_like(f) + # ief = index_nearest(x, e_cut) + ief = np.argmin(np.abs(x - e_cut)) + f[0:ief] *= 0 + + if x.shape != gammas.shape: + print("Error: 'gammas' array does not have the same shape of 'x'") + return 0 + + # linear fit upper part of the spectrum to avoid border effects + # polyfit => pf + lpf = int(len(x) / 2.0) + cpf = polyfit(x[-lpf:], f[-lpf:], 1, reverse=False) + fpf = np.polynomial.Polynomial(cpf) + + # extend upper energy border to 3*fhwm_e[-1] + xstep = x[-1] - x[-2] + xup = np.append(x, np.arange(x[-1] + xstep, x[-1] + 3 * gammas[-1], xstep)) + + for n in range(len(f)): + # from now on I change e with eup + eimin, eimax = get_ene_index(xup, xup[n], 1.5 * gammas[n]) + if (eimin is None) or (eimax is None): + if DEBUG: + raise IndexError("e[{0}]".format(n)) + if len(range(eimin, eimax)) % 2 == 0: + kx = xup[eimin : eimax + 1] # odd range centered at the convolution point + else: + kx = xup[eimin:eimax] + ### kernel ### + hwhm = gammas[n] / 2.0 + if "gauss" in kernel.lower(): + ky = gaussian(kx, center=xup[n], sigma=hwhm) + elif "lor" in kernel.lower(): + ky = lorentzian(kx, center=xup[n], sigma=hwhm) + else: + raise ValueError("convolution kernel '{0}' not implemented".format(kernel)) + ky = ky / ky.sum() # normalize + zn = 0 + lk = len(kx) + for mf, mg in zip(range(-int(lk / 2), int(lk / 2) + 1), range(lk)): + if ((n + mf) >= 0) and ((n + mf) < len(f)): + zn += f[n + mf] * ky[mg] + elif (n + mf) >= 0: + zn += fpf(xup[n + mf]) * ky[mg] + z[n] = zn + return z From 0e79588437b6b22f88543d370c1ceab763165268 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 17:33:21 +0200 Subject: [PATCH 31/52] full local copy of convoluition --- larixite/math/convolution.py | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/larixite/math/convolution.py b/larixite/math/convolution.py index a6d5bc0..8cf1805 100644 --- a/larixite/math/convolution.py +++ b/larixite/math/convolution.py @@ -21,11 +21,29 @@ .. [GP1202] """ + import numpy as np -from .lineshapes import gaussian, lorentzian + +s2pi = np.sqrt(2 * np.pi) +tiny = 1.0e-15 DEBUG = 0 + +def gaussian(x, amplitude=1.0, center=0.0, sigma=1.0): + """Return a 1-dimensional Gaussian function.""" + return (amplitude / (max(tiny, s2pi * sigma))) * np.exp( + -((1.0 * x - center) ** 2) / max(tiny, (2 * sigma**2)) + ) + + +def lorentzian(x, amplitude=1.0, center=0.0, sigma=1.0): + """Return a 1-dimensional Lorentzian function.""" + return (amplitude / (1 + ((1.0 * x - center) / max(tiny, sigma)) ** 2)) / max( + tiny, (np.pi * sigma) + ) + + def polyfit(x, y, deg=1, reverse=False): """ simple emulation of deprecated numpy.polyfit, @@ -37,6 +55,7 @@ def polyfit(x, y, deg=1, reverse=False): coefs = list(reversed(coefs)) return list(coefs) + def get_ene_index(ene, cen, hwhm): """returns the min/max indexes for array ene at (cen-hwhm) and (cen+hwhm) very similar to index_of in larch From ba2bbbe9570851ca4cb76b9ed98e8305e97d0c45 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 27 May 2026 17:50:06 +0200 Subject: [PATCH 32/52] inheritance from base model for FdmnesXasSim --- larixite/fdmnes/output.py | 3 ++- larixite/models/__init__.py | 2 +- larixite/models/xas.py | 15 +++++++++++++-- 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py index 531719c..13c77d0 100644 --- a/larixite/fdmnes/output.py +++ b/larixite/fdmnes/output.py @@ -12,13 +12,14 @@ import yaml import pandas as pd from larixite.fdmnes.input import FdmnesXasInput +from larixite.models.xas import XasSimData from larixite.utils import get_logger logger = get_logger("larixite.fdmnes.output") @dataclass -class FdmnesXasSim: +class FdmnesXasSim(XasSimData): """A FDMNES XAS simulation associated with a job directory.""" jobdir: Path diff --git a/larixite/models/__init__.py b/larixite/models/__init__.py index 14ae1c9..a86de50 100644 --- a/larixite/models/__init__.py +++ b/larixite/models/__init__.py @@ -1,3 +1,3 @@ """Models for data handled by Larixite""" -from .xas import XasExpData \ No newline at end of file +from .xas import XasExpData, XasSimData \ No newline at end of file diff --git a/larixite/models/xas.py b/larixite/models/xas.py index 50621c0..fb70c20 100644 --- a/larixite/models/xas.py +++ b/larixite/models/xas.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -"""Experimental XAS data model.""" +"""Base model for XAS data.""" from dataclasses import dataclass from pathlib import Path @@ -12,12 +12,23 @@ logger = get_logger("larixite.models.xas") +@dataclass +class XasSimData: + """Base model for XAS simulation data.""" + + jobdir: Path + input: object | None = None + data: pd.DataFrame | None = None + metadata: dict | None = None + + @dataclass class XasExpData: """Model for an Experimental XAS spectrum.""" data: pd.DataFrame | None = None - + metadata: dict | None = None + @classmethod def from_csv(cls, csvpath: Path | str) -> "XasExpData": """Load experimental XAS data from a CSV file. From a6e12abe2c3a743f91e1d56e0d0e5a8ec32f8010 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Thu, 28 May 2026 10:47:51 +0200 Subject: [PATCH 33/52] revert back to a separate temporary project (a.k.a. melezia) --- larixite/fdmnes/__init__.py | 3 +- larixite/fdmnes/output.py | 144 --------------------- larixite/math/__init__.py | 0 larixite/math/convolution.py | 235 ----------------------------------- larixite/models/__init__.py | 3 - larixite/models/xas.py | 49 -------- 6 files changed, 1 insertion(+), 433 deletions(-) delete mode 100644 larixite/fdmnes/output.py delete mode 100644 larixite/math/__init__.py delete mode 100644 larixite/math/convolution.py delete mode 100644 larixite/models/__init__.py delete mode 100644 larixite/models/xas.py diff --git a/larixite/fdmnes/__init__.py b/larixite/fdmnes/__init__.py index 3f750e7..40e694f 100644 --- a/larixite/fdmnes/__init__.py +++ b/larixite/fdmnes/__init__.py @@ -1,2 +1 @@ -from .input import FdmnesXasInput, struct2fdmnes -from .output import search_jobs \ No newline at end of file +from .input import FdmnesXasInput, struct2fdmnes \ No newline at end of file diff --git a/larixite/fdmnes/output.py b/larixite/fdmnes/output.py deleted file mode 100644 index 13c77d0..0000000 --- a/larixite/fdmnes/output.py +++ /dev/null @@ -1,144 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -""" -Handle FDMNES output files. -""" - -from dataclasses import dataclass -from pathlib import Path - -import numpy as np -import yaml -import pandas as pd -from larixite.fdmnes.input import FdmnesXasInput -from larixite.models.xas import XasSimData -from larixite.utils import get_logger - -logger = get_logger("larixite.fdmnes.output") - - -@dataclass -class FdmnesXasSim(XasSimData): - """A FDMNES XAS simulation associated with a job directory.""" - - jobdir: Path - input: FdmnesXasInput | None = None - data: pd.DataFrame | None = None - metadata: dict | None = None - - @property - def prefix(self) -> str | None: - if self.input is not None: - return self.input.fileout_prefix - return None - - @property - def efermi(self) -> float: - if self.metadata is None: - self.load() - return self.metadata.get("VO_interstitial", 0) - - @property - def e0(self) -> float: - if self.metadata is None: - self.load() - return self.metadata.get("E_edge", 0) - - @property - def energy(self, erel: bool = False) -> np.ndarray: - if self.data is None: - self.load() - eout = self.data["Energy"].values - return eout + self.e0 if erel else eout - - @property - def mu(self) -> np.ndarray: - if self.data is None: - self.load() - return self.data[""].values - - def load(self) -> pd.DataFrame | None: - """Load the FDMNES output data file and parse the metadata line.""" - datafile = self.jobdir / f"{self.prefix}.txt" - if not datafile.exists(): - error_msg = "cannot find the calculation output file" - logger.error(error_msg) - error_file = self.jobdir / "error.yaml" - error_file.write_text(yaml.dump({"error": error_msg})) - return - - lines = datafile.read_text().splitlines() - skiprows = 1 - - for i, line in enumerate(lines): - stripped = line.strip() - if stripped.startswith("#"): - continue - elif "=" in stripped: - eq_split = stripped.split(" = ", 1) - left_vals = [float(v) for v in eq_split[0].split() if v] - right_names = [n.strip() for n in eq_split[1].split(",")] - assert len(left_vals) == len(right_names), ( - "Error parsing the metadata line" - ) - int_names = {"Z", "n_edge", "j_edge", "ninit", "ninit1", "natomsym_f"} - self.metadata = { - name: int(val) if name in int_names else val - for name, val in zip(right_names, left_vals) - } - skiprows = i + 1 - else: - break - - self.data = pd.read_csv(datafile, sep=r"\s+", skiprows=skiprows) - assert len(self.data.columns) >= 2, "Output file must have at least two columns" - assert "Energy" in self.data.columns, "No 'Energy' column in the output file" - assert "" in self.data.columns, "No '' column in the output file" - - -def search_jobs(globstr: str, jobsdir: Path | str | None = None) -> dict[str, FdmnesXasSim]: - """Search for FDMNES job directories matching a glob pattern. - - For each matching subdirectory that contains a ``{fileout_prefix}_params.yaml`` - file, a :class:`FdmnesXasSim` instance is built via - :meth:`FdmnesXasInput.from_yaml`. - - Parameters - ---------- - globstr : str - Glob pattern applied to the *name* of each immediate subdirectory - (e.g. ``"job_*"``, ``"mi_*"``, ``"*.done"``). - jobsdir : Path | str | None -> ~/.larixite/fdmnes - Parent directory containing job subdirectories. - - Returns - ------- - dict[str, FdmnesXasSim] - Dictionary mapping ``fileout_prefix`` to matching simulation objects. - """ - if jobsdir is None: - jobsdir = Path("~/.larixite/fdmnes").expanduser().resolve() - if isinstance(jobsdir, str): - jobsdir = Path(jobsdir).expanduser().resolve() - if not jobsdir.is_dir(): - return {} - - results: dict[str, FdmnesXasSim] = {} - for entry in jobsdir.glob(globstr): - if not entry.is_dir(): - continue - - for param_file in entry.glob("*_params.yaml"): - try: - inp = FdmnesXasInput.from_yaml(param_file) - expected = entry / f"{inp.fileout_prefix}_params.yaml" - if param_file.resolve() != expected.resolve(): - continue - sim = FdmnesXasSim(jobdir=entry, input=inp) - results[inp.fileout_prefix] = sim - except Exception: - continue - - logger.info(f"found {len(results)} job(s) in {jobsdir}") - return results diff --git a/larixite/math/__init__.py b/larixite/math/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/larixite/math/convolution.py b/larixite/math/convolution.py deleted file mode 100644 index 8cf1805..0000000 --- a/larixite/math/convolution.py +++ /dev/null @@ -1,235 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -""" -Generic discrete (one-dimensional) convolution - -Description ------------ - - This is a manual (not optimized!) implementation of discrete one-dimensional - (1D) convolution intended for spectroscopy analysis. The difference with - commonly used methods is the possibility to adapt the width of the convolution - kernel for each convolution point, e.g. change the FWHM of the Gaussian kernel - as a function of the energy scale (x-axis). - -Resources ---------- - -.. [WPconv] -.. [Fisher] -.. [GP1202] - -""" - -import numpy as np - -s2pi = np.sqrt(2 * np.pi) -tiny = 1.0e-15 - -DEBUG = 0 - - -def gaussian(x, amplitude=1.0, center=0.0, sigma=1.0): - """Return a 1-dimensional Gaussian function.""" - return (amplitude / (max(tiny, s2pi * sigma))) * np.exp( - -((1.0 * x - center) ** 2) / max(tiny, (2 * sigma**2)) - ) - - -def lorentzian(x, amplitude=1.0, center=0.0, sigma=1.0): - """Return a 1-dimensional Lorentzian function.""" - return (amplitude / (1 + ((1.0 * x - center) / max(tiny, sigma)) ** 2)) / max( - tiny, (np.pi * sigma) - ) - - -def polyfit(x, y, deg=1, reverse=False): - """ - simple emulation of deprecated numpy.polyfit, - including its ordering of coefficients - """ - pfit = np.polynomial.Polynomial.fit(x, y, deg=int(deg)) - coefs = pfit.convert().coef - if reverse: - coefs = list(reversed(coefs)) - return list(coefs) - - -def get_ene_index(ene, cen, hwhm): - """returns the min/max indexes for array ene at (cen-hwhm) and (cen+hwhm) - very similar to index_of in larch - """ - try: - if (cen - hwhm) <= min(ene): - ene_imin = 0 - else: - ene_imin = max(np.where(ene < (cen - hwhm))[0]) - if (cen + hwhm) >= max(ene): - ene_imax = len(e) - 1 - else: - ene_imax = min(np.where(ene > (cen + hwhm))[0]) - return ene_imin, ene_imax - except Exception: - print("index not found for {0} +/- {1}".format(cen, hwhm)) - return None, None - - -def lin_gamma(ene, gamma_hole=0.5, linbroad=None): - """returns constant or linear energy-dependent broadening - - Parameters - ---------- - ene : energy array in eV - gamma_hole : initial full width at half maximum in eV - linbroad : list of 3-elements giving - 'final full width at half maximum' - 'start energy point of the linear increase' - 'end energy point of the linear increase' - """ - w = np.ones_like(ene) - if linbroad is None: - return w * gamma_hole - else: - try: - fwhm2 = linbroad[0] - e1 = linbroad[1] - e2 = linbroad[2] - except: - raise ValueError("wrong format for linbroad") - for en, ee in enumerate(ene): - if ee < e1: - w[en] *= gamma_hole - elif ee <= e2: - wlin = gamma_hole + (ee - e1) * (fwhm2 - gamma_hole) / (e2 - e1) - w[en] *= wlin - elif ee >= e2: - w[en] *= fwhm2 - return w - - -def atan_gamma(e, e_cut=0, e_cent=30, e_larg=30, gamma_hole=0.5, gamma_max=15): - """Arctangent energy dependent broadening as implemented in FDMNES (described in the "Convolution" section of the manual)""" - - f = (e - e_cut) / e_cent - a = np.pi * gamma_max * (f - 1 / f**2) / (3 * e_larg) - gammas = gamma_hole + gamma_max * (0.5 + 1 / np.pi * np.arctan(a)) - - # Set gamma to the gamma_hole below the cutoff energy. - mask = np.where(e < e_cut) - gammas[mask] = gamma_hole - - return gammas - - -def conv_fast(x, y, gammas, e_cut=None, kernel="gaussian", num=501, step=0.1): - """A significantly faster version of the `conv` function - - Parameters - ---------- - x : x-axis (energy) - y : f(x) to convolve with g(x) kernel, y(energy) - kernel : convolution kernel, g(x) - 'gaussian' - 'lorentzian' - gammas : the full width half maximum in eV for the kernel - broadening. It is an array of size 'e' with constants or - an energy-dependent values determined by a function as - 'lin_gamma()' or 'atan_gamma()' - """ - assert e_cut is not None, "starting energy for the convolution not given" - assert "lor" or "gauss" in kernel.lower(), ( - "the kernel should be either Lorentzian or Gaussian" - ) - - # Extend the X-axis array. - start = x[-1] - stop = start + num * step - x_ext = np.append(x, np.arange(start + step, stop, step)) - - ids = x_ext > e_cut - - x1 = x_ext[ids] - x1[0] = e_cut - - # Extend the intensity array by coping the last value `num - 1` times. - y = y[x > e_cut] - y = np.append(y, np.ones(num - 1) * y[-1]) - - y_conv = np.zeros_like(x) - for i, (xi, gamma) in enumerate(zip(x, gammas)): - gamma = gamma / 2.0 - if "gauss" in kernel.lower(): - ky = gaussian(x1, center=xi, sigma=gamma) - elif "lor" in kernel.lower(): - ky = lorentzian(x1, center=xi, sigma=gamma) - y_conv[i] = np.sum(ky * y) / np.pi - - return y_conv - - -def conv(x, y, gammas, e_cut=None, kernel="gaussian"): - """linear broadening - - Parameters - ---------- - x : x-axis (energy) - y : f(x) to convolve with g(x) kernel, y(energy) - kernel : convolution kernel, g(x) - 'gaussian' - 'lorentzian' - gammas : the full width half maximum in eV for the kernel - broadening. It is an array of size 'e' with constants or - an energy-dependent values determined by a function as - 'lin_gamma()' or 'atan_gamma()' - """ - assert e_cut is not None, "starting energy for the convolution not given" - - f = y[:] * 1.0 - z = np.zeros_like(f) - # ief = index_nearest(x, e_cut) - ief = np.argmin(np.abs(x - e_cut)) - f[0:ief] *= 0 - - if x.shape != gammas.shape: - print("Error: 'gammas' array does not have the same shape of 'x'") - return 0 - - # linear fit upper part of the spectrum to avoid border effects - # polyfit => pf - lpf = int(len(x) / 2.0) - cpf = polyfit(x[-lpf:], f[-lpf:], 1, reverse=False) - fpf = np.polynomial.Polynomial(cpf) - - # extend upper energy border to 3*fhwm_e[-1] - xstep = x[-1] - x[-2] - xup = np.append(x, np.arange(x[-1] + xstep, x[-1] + 3 * gammas[-1], xstep)) - - for n in range(len(f)): - # from now on I change e with eup - eimin, eimax = get_ene_index(xup, xup[n], 1.5 * gammas[n]) - if (eimin is None) or (eimax is None): - if DEBUG: - raise IndexError("e[{0}]".format(n)) - if len(range(eimin, eimax)) % 2 == 0: - kx = xup[eimin : eimax + 1] # odd range centered at the convolution point - else: - kx = xup[eimin:eimax] - ### kernel ### - hwhm = gammas[n] / 2.0 - if "gauss" in kernel.lower(): - ky = gaussian(kx, center=xup[n], sigma=hwhm) - elif "lor" in kernel.lower(): - ky = lorentzian(kx, center=xup[n], sigma=hwhm) - else: - raise ValueError("convolution kernel '{0}' not implemented".format(kernel)) - ky = ky / ky.sum() # normalize - zn = 0 - lk = len(kx) - for mf, mg in zip(range(-int(lk / 2), int(lk / 2) + 1), range(lk)): - if ((n + mf) >= 0) and ((n + mf) < len(f)): - zn += f[n + mf] * ky[mg] - elif (n + mf) >= 0: - zn += fpf(xup[n + mf]) * ky[mg] - z[n] = zn - return z diff --git a/larixite/models/__init__.py b/larixite/models/__init__.py deleted file mode 100644 index a86de50..0000000 --- a/larixite/models/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -"""Models for data handled by Larixite""" - -from .xas import XasExpData, XasSimData \ No newline at end of file diff --git a/larixite/models/xas.py b/larixite/models/xas.py deleted file mode 100644 index fb70c20..0000000 --- a/larixite/models/xas.py +++ /dev/null @@ -1,49 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - -"""Base model for XAS data.""" - -from dataclasses import dataclass -from pathlib import Path - -import pandas as pd -from larixite.utils import get_logger - -logger = get_logger("larixite.models.xas") - - -@dataclass -class XasSimData: - """Base model for XAS simulation data.""" - - jobdir: Path - input: object | None = None - data: pd.DataFrame | None = None - metadata: dict | None = None - - -@dataclass -class XasExpData: - """Model for an Experimental XAS spectrum.""" - - data: pd.DataFrame | None = None - metadata: dict | None = None - - @classmethod - def from_csv(cls, csvpath: Path | str) -> "XasExpData": - """Load experimental XAS data from a CSV file. - - Parameters - ---------- - csvpath : Path | str - Path to the CSV file containing the XAS data. - - Returns - ------- - XasExpData - Instance with the data loaded from the CSV. - """ - if isinstance(csvpath, str): - csvpath = Path(csvpath) - data = pd.read_csv(csvpath) - return cls(data=data) From a0704cacc5ecb8fa06ce07b84b9de5c636a8a82e Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Thu, 28 May 2026 16:41:25 +0200 Subject: [PATCH 34/52] add total number of atoms, remove _round_up --- larixite/struct/visualize.py | 17 +++++------------ tests/test_visualize.py | 18 ------------------ 2 files changed, 5 insertions(+), 30 deletions(-) diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py index cd3b080..434a33d 100644 --- a/larixite/struct/visualize.py +++ b/larixite/struct/visualize.py @@ -6,7 +6,6 @@ =============================================================== """ -import math import sys from pymatgen.analysis.local_env import CovalentRadius @@ -137,10 +136,6 @@ } -def _round_up(x: float) -> float: - return math.ceil(x * 100) / 100 - - def _lighten_hex(hex_color: str, factor: float = 0.75) -> str: """Return a darker hex color obtained by scaling each RGB channel by *factor*.""" r = int(int(hex_color[1:3], 16) * factor) @@ -149,7 +144,7 @@ def _lighten_hex(hex_color: str, factor: float = 0.75) -> str: return "#{:02x}{:02x}{:02x}".format(r, g, b) -def _print_color_legend(color_elems: dict[str, str]) -> None: +def _print_color_legend(color_elems: dict[str, str], natoms: int) -> None: """Print a human-readable color legend, adapting to the execution environment. In Jupyter: renders HTML color swatches. @@ -164,7 +159,7 @@ def _print_color_legend(color_elems: dict[str, str]) -> None: f'{e}' for e, c in color_elems.items() ) - display(HTML(f'
Colors: {swatches}
')) + display(HTML(f'
{natoms} atoms — {swatches}
')) return except ImportError: pass @@ -175,11 +170,11 @@ def _print_color_legend(color_elems: dict[str, str]) -> None: g = int(hex_color[3:5], 16) b = int(hex_color[5:7], 16) sys.stdout.write("\x1b[38;2;{};{};{};m{}\x1b[0m ".format(r, g, b, elem)) - sys.stdout.write("\n") + sys.stdout.write(f" ({natoms} atoms)\n") return legend = ", ".join(f"{e}: {VESTA_COLORS.get(e, ('#888888', 'grey'))[1]}" for e in color_elems) - logger.info(f"Colors: {legend}") + logger.info(f"Colors ({natoms} atoms): {legend}") def _adjust_font_color(hex_color: str, dark: str = "#222222", light: str = "#ffffff") -> str: @@ -228,8 +223,6 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) logger.error("py3Dmol not installed! -> run `pip install py3Dmol`") return - radius = _round_up(radius) - mol = struct.build_cluster(radius=radius) xyz, elems = _cluster_to_xyz(mol) @@ -293,5 +286,5 @@ def visualize(struct: XasStructure, radius: float = 2.5, unitcell: bool = False) xyzview.zoomTo() xyzview.show() - _print_color_legend(color_elems) + _print_color_legend(color_elems, len(mol)) return xyzview diff --git a/tests/test_visualize.py b/tests/test_visualize.py index 97fd01f..4e5201e 100644 --- a/tests/test_visualize.py +++ b/tests/test_visualize.py @@ -1,16 +1,13 @@ """Tests for larixite.struct.visualize — _cluster_to_xyz, _lighten_hex, VESTA_COLORS, visualize().""" -import math from pathlib import Path import pytest -import numpy as np from larixite.struct import get_structure from larixite.struct.visualize import ( HAS_PY3DMOL, VESTA_COLORS, - _round_up, _cluster_to_xyz, _lighten_hex, visualize, @@ -31,21 +28,6 @@ def xyz_cuo6(): return get_structure(STRUCTS / "CuO6_D4h.xyz", "Cu") -# ── _round_up ─────────────────────────────────────────────────── - - -class TestRoundUp: - def test_integer(self): - assert _round_up(2.0) == 2.0 - - def test_rounds_up(self): - assert _round_up(2.31) == 2.31 - - def test_rounds_up_small(self): - assert _round_up(2.001) == 2.01 - - def test_rounds_up_boundary(self): - assert _round_up(2.501) == 2.51 # ── _cluster_to_xyz ───────────────────────────────────────────── From 95efa90ae41e136b6df24a0d66cce45b9bc4c4cc Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Thu, 28 May 2026 16:49:57 +0200 Subject: [PATCH 35/52] some notes for future devs --- larixite/struct/visualize.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/larixite/struct/visualize.py b/larixite/struct/visualize.py index 434a33d..8f41f83 100644 --- a/larixite/struct/visualize.py +++ b/larixite/struct/visualize.py @@ -4,6 +4,30 @@ """ 3D visualization of XasStructure local environment using py3Dmol =============================================================== + +TODO (future improvements): + +- Add coordination polyhedra around the absorber atom. py3Dmol does not support + polyhedron rendering directly. Possible approaches: + + 1. **nglview** — the `nglview` Jupyter widget supports `add_polyhedron()` and + works directly with pymatgen structures. Could be offered as an alternative + backend (e.g. `visualize(..., backend="nglview")`). + + 2. **py3Dmol workaround with triangles** — compute the convex hull of the first + coordination shell with `scipy.spatial.ConvexHull`, then draw the facets as + semi-transparent plane surfaces via `py3Dmol`'s `addSurface()` or by adding + dummy atoms at facet centroids with connecting sticks to form a wireframe. + + 3. **pymatgen IPyWidgets** — `pymatgen.vis.ipyvolumetric` or `pymatgen.vis` can + render structures with coordination environments in interactive HTML. + + 4. **Export for VESTA** — write the cluster as a `.vesta` file (VESTA's native + format supports polyhedra definitions) so the user can open it in VESTA. + +- Expose the pymatgen Molecule cluster for downstream rendering. +- Support partial occupancy visualization (transparency, multi-snapshot). +- Allow custom color schemes beyond VESTA. """ import sys From 241a9a0d93e6ea7ee5fb4c0f294a87f25d0d43fa Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 29 May 2026 09:27:01 +0200 Subject: [PATCH 36/52] more sbatch options --- larixite/fdmnes/input.py | 7 +++++-- larixite/templates/fdmnes_sbatch_esrf.tmpl | 4 ++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index de9c314..55d1014 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -408,7 +408,7 @@ def write_input( _ = self.dump_params(jobdir / f"{self.fileout_prefix}_params.yaml") return jobdir - def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | None = None, ncpus: int = 8) -> Path: + def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | None = None, ncpus: int = 8, nnodes: int = 1, mem_per_cpu: str = "16GB", **kwargs) -> Path: """Generates a SBATCH file (SLURM workload manager) using a template Arguments @@ -439,7 +439,10 @@ def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | jobdir = Path(jobdir) sbatchout = jobdir / f"{self.fileout_prefix}.sbatch" kwargs = {"jobname": self.fileout_prefix, - "ncpus": ncpus,} + "nnodes": nnodes, + "ncpus": ncpus, + "mem_per_cpu": mem_per_cpu, + **kwargs} with open(sbatchout, "w") as fp, open(template) as tp: fp.write(tp.read().format(**kwargs)) logger.info(f"written {fp.name}") diff --git a/larixite/templates/fdmnes_sbatch_esrf.tmpl b/larixite/templates/fdmnes_sbatch_esrf.tmpl index 2cce9e3..8f2d9a3 100644 --- a/larixite/templates/fdmnes_sbatch_esrf.tmpl +++ b/larixite/templates/fdmnes_sbatch_esrf.tmpl @@ -1,8 +1,8 @@ #!/bin/bash -l -#SBATCH --nodes=1 #: Number of nodes / if more nodes, add `#SBATCH --constraint="hpc7"` (type of machines) +#SBATCH --nodes={nnodes} #: Number of nodes / if more nodes, add `#SBATCH --constraint="hpc7"` (type of machines) #SBATCH --ntasks-per-node={ncpus} #: Total number of tasks (=CPUs cores) per node #SBATCH --cpus-per-task=1 #: CPUs for each task / DON'T CHANGE THIS - required by FDMNES -#SBATCH --mem-per-cpu=16G #: Total memory for each CPU +#SBATCH --mem-per-cpu={mem_per_cpu} #: Total memory for each CPU #SBATCH --time=120:00:00 #: Wall-time limit #SBATCH --partition=nice,nice-long #: ESRF cluster #SBATCH --job-name="{jobname}" #: Name used to set the output files From 88fb6fd7e9a08ee6b51400fc3022730b5cfc7a7a Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 29 May 2026 09:32:13 +0200 Subject: [PATCH 37/52] add run_sbatch --- larixite/fdmnes/input.py | 66 +++++++++++++++++++++++++++++++++++----- 1 file changed, 58 insertions(+), 8 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 55d1014..77267e0 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -9,10 +9,12 @@ spectroscopy (XAS, XES, RIXS) from the atomic structures """ + import time import yaml from dataclasses import dataclass from pathlib import Path +import subprocess from pymatgen.core import __version__ as pymatgen_version, Element, Molecule from larixite.struct import get_structure, get_structure_from_text from larixite.struct.xas import XasStructure @@ -408,12 +410,20 @@ def write_input( _ = self.dump_params(jobdir / f"{self.fileout_prefix}_params.yaml") return jobdir - def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | None = None, ncpus: int = 8, nnodes: int = 1, mem_per_cpu: str = "16GB", **kwargs) -> Path: + def write_sbatch( + self, + jobdir: str | Path | None = None, + template: str | Path | None = None, + ncpus: int = 8, + nnodes: int = 1, + mem_per_cpu: str = "16GB", + **kwargs, + ) -> Path: """Generates a SBATCH file (SLURM workload manager) using a template Arguments --------- - jobdir: str, Path or None + jobdir: str, Path or None path to the job directory to write the SBATCH file if None: uses last job directory template: str, Path or None @@ -438,16 +448,54 @@ def write_sbatch(self, jobdir: str | Path | None = None, template: str | Path | if isinstance(jobdir, str): jobdir = Path(jobdir) sbatchout = jobdir / f"{self.fileout_prefix}.sbatch" - kwargs = {"jobname": self.fileout_prefix, - "nnodes": nnodes, - "ncpus": ncpus, - "mem_per_cpu": mem_per_cpu, - **kwargs} + kwargs = { + "jobname": self.fileout_prefix, + "nnodes": nnodes, + "ncpus": ncpus, + "mem_per_cpu": mem_per_cpu, + **kwargs, + } with open(sbatchout, "w") as fp, open(template) as tp: fp.write(tp.read().format(**kwargs)) logger.info(f"written {fp.name}") return sbatchout + def run_sbatch(self, jobdir: Path) -> str | None: + """Submit the sbatch script via `sbatch --parsable`. + + Writes `status.yaml` with status=``submitted`` and the SLURM job ID. + Returns the job ID on success, None on failure. + """ + sbatch_script = jobdir / f"{self.fileout_prefix}.sbatch" + if not sbatch_script.exists(): + logger.warning(f"No sbatch script found for job {self.fileout_prefix}") + return None + + try: + result = subprocess.run( + ["sbatch", "--parsable", str(sbatch_script)], + capture_output=True, + text=True, + check=True, + ) + slurm_job_id = result.stdout.strip() + status_file = jobdir / "status.yaml" + status_file.write_text( + yaml.dump( + { + "status": "submitted", + "slurm_job_id": slurm_job_id, + } + ) + ) + logger.info( + f"SUBMITTED job {self.fileout_prefix} (slurm job id: {slurm_job_id})" + ) + return slurm_job_id + except subprocess.CalledProcessError as e: + logger.error(f"FAILED to submit job {self.fileout_prefix}: {e}") + return None + def dump_params(self, yamlpath: str | Path | None = None) -> Path: """Dump input parameters to a YAML file @@ -486,7 +534,9 @@ def dump_params(self, yamlpath: str | Path | None = None) -> Path: if self._jobs: yamlpath = self._jobs[-1] / f"{self.fileout_prefix}_params.yaml" else: - yamlpath = Path(self.structpath).parent / f"{self.fileout_prefix}_params.yaml" + yamlpath = ( + Path(self.structpath).parent / f"{self.fileout_prefix}_params.yaml" + ) if isinstance(yamlpath, str): yamlpath = Path(yamlpath) From 5acd3094dbce25397c881eec2cb770562a66b757 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 29 May 2026 09:50:00 +0200 Subject: [PATCH 38/52] update esrf sbatch --- larixite/fdmnes/input.py | 4 +++- larixite/templates/fdmnes_sbatch_esrf.tmpl | 5 +++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 77267e0..620f12d 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -414,9 +414,10 @@ def write_sbatch( self, jobdir: str | Path | None = None, template: str | Path | None = None, - ncpus: int = 8, + ncpus: int = 12, nnodes: int = 1, mem_per_cpu: str = "16GB", + constraint: str = "", **kwargs, ) -> Path: """Generates a SBATCH file (SLURM workload manager) using a template @@ -453,6 +454,7 @@ def write_sbatch( "nnodes": nnodes, "ncpus": ncpus, "mem_per_cpu": mem_per_cpu, + "constraint": constraint, **kwargs, } with open(sbatchout, "w") as fp, open(template) as tp: diff --git a/larixite/templates/fdmnes_sbatch_esrf.tmpl b/larixite/templates/fdmnes_sbatch_esrf.tmpl index 8f2d9a3..08ec360 100644 --- a/larixite/templates/fdmnes_sbatch_esrf.tmpl +++ b/larixite/templates/fdmnes_sbatch_esrf.tmpl @@ -5,9 +5,10 @@ #SBATCH --mem-per-cpu={mem_per_cpu} #: Total memory for each CPU #SBATCH --time=120:00:00 #: Wall-time limit #SBATCH --partition=nice,nice-long #: ESRF cluster +#SBATCH --constraint={constraint} #SBATCH --job-name="{jobname}" #: Name used to set the output files -#SBATCH --output=%x_output.txt -#SBATCH --error=%x_output.txt +#SBATCH --output=%x_slurm_stdout.txt +#SBATCH --error=%x_slurm_stderr.txt function compress {{ gzip -f {jobname}_bav.txt From 5d302d26e7be46d297a5422311400e43a06b16f1 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 29 May 2026 10:42:36 +0200 Subject: [PATCH 39/52] fix bug --- larixite/fdmnes/input.py | 1 + 1 file changed, 1 insertion(+) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 620f12d..cb40f7b 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -479,6 +479,7 @@ def run_sbatch(self, jobdir: Path) -> str | None: capture_output=True, text=True, check=True, + cwd=jobdir, ) slurm_job_id = result.stdout.strip() status_file = jobdir / "status.yaml" From f453bb77afacc73c85f33b29c8a6a50f9cb06dc8 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 29 May 2026 11:09:16 +0200 Subject: [PATCH 40/52] minor improvements --- larixite/templates/fdmnes_sbatch_esrf.tmpl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/larixite/templates/fdmnes_sbatch_esrf.tmpl b/larixite/templates/fdmnes_sbatch_esrf.tmpl index 08ec360..3d9f39b 100644 --- a/larixite/templates/fdmnes_sbatch_esrf.tmpl +++ b/larixite/templates/fdmnes_sbatch_esrf.tmpl @@ -5,10 +5,10 @@ #SBATCH --mem-per-cpu={mem_per_cpu} #: Total memory for each CPU #SBATCH --time=120:00:00 #: Wall-time limit #SBATCH --partition=nice,nice-long #: ESRF cluster -#SBATCH --constraint={constraint} +#SBATCH --constraint="{constraint}" #SBATCH --job-name="{jobname}" #: Name used to set the output files -#SBATCH --output=%x_slurm_stdout.txt -#SBATCH --error=%x_slurm_stderr.txt +#SBATCH --output="{jobname}_slurm.stdout" +#SBATCH --error="{jobname}_slurm.stderr" function compress {{ gzip -f {jobname}_bav.txt From f14ec6cf1bb43287bdee64d7990ff116adb56188 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Fri, 29 May 2026 11:25:36 +0200 Subject: [PATCH 41/52] fix running from jupyter-slurm --- larixite/fdmnes/input.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index cb40f7b..0add9a1 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -15,6 +15,7 @@ from dataclasses import dataclass from pathlib import Path import subprocess +import os from pymatgen.core import __version__ as pymatgen_version, Element, Molecule from larixite.struct import get_structure, get_structure_from_text from larixite.struct.xas import XasStructure @@ -49,6 +50,15 @@ } +def _clean_env(): + """Return environment without SLURM_ variables.""" + env = dict(os.environ) + for key in list(env): + if key.startswith("SLURM_"): + del env[key] + return env + + @dataclass class FdmnesXasInput: """Input generator for a XAS calculation with FDMNES""" @@ -480,6 +490,7 @@ def run_sbatch(self, jobdir: Path) -> str | None: text=True, check=True, cwd=jobdir, + env=_clean_env(), ) slurm_job_id = result.stdout.strip() status_file = jobdir / "status.yaml" From 8f89230dd5700c184d58e3f85e582424fae46162 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Sun, 31 May 2026 17:05:27 +0200 Subject: [PATCH 42/52] add convolution + refactor key logics --- larixite/fdmnes/input.py | 71 ++++++++++++++++++++---------- larixite/templates/fdmnes_xas.tmpl | 21 +++++++-- 2 files changed, 66 insertions(+), 26 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 0add9a1..c6239cd 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -39,6 +39,10 @@ "Spinorbit": False, "SCF": False, "SCF_exc": False, + "R_self": False, + "N_self": True, + "P_self": True, + "Vmax": False, "Screening": False, #: -> TODO "Dilatorb": False, #: -> TODO "TDDFT": False, @@ -47,6 +51,12 @@ "Atom_conf": False, #: preferred over `Atom` (permits to keep atomic number in the list of atoms) -> TODO "COOP": False, #: -> TODO "Convolution": True, + "E_cut": False, #: taken from E_fermi + "Dec": True, #: energy shift + "Ecent": True, #: arctan center + "Elarg": True, #: arctan width + "Gamma_hole": False, #: use core-hole broadening + "Gamma_max": True, } @@ -77,7 +87,7 @@ class FdmnesXasInput: rself: float | None = None #: radius for the SCF calculation nself: int = 100 #: number of maximum iterations for the SCF calculation pself: float = 0.025 #: initial weight for the SCF calculation - vmax: float | None = None #: maximum potential value for molecules + vmax: float | str | None = None #: maximum potential value for molecules fileout_prefix: str = ( "job" #: prefix of the output filename for the FDMNES job (extension: .inp) ) @@ -88,6 +98,12 @@ class FdmnesXasInput: optimize: bool = False #: optimize the input parameters spacer: str = " " #: spacer string for the FDMNES input text outdir: str | Path | None = None #: path to the output directory of the FDMNES jobs + ecut: float | str | None = None #: energy cutoff for the convolution + dec: float = 0 #: energy shift for the convolution + ecent: float = 30. #: energy center for the convolution + elarg: float = 30. #: energy width for the convolution + gamma_hole: float | str | None = None #: start width of the energy convolution (None -> core-hole broadening) + gamma_max: float = 8 #: maximum energy for the convolution def __post_init__(self): """Validate and optimize attributes""" @@ -109,6 +125,14 @@ def __post_init__(self): #: R_self if self.rself is None: self.set_rself() + else: + self.params["R_self"] = True + #: atictivate Vmax if given by the user + if self.vmax is None: + self.params["Vmax"] = False + self.vmax = "!" + else: + self.params["Vmax"] = True #: structure type if self.struct_type is None: self.struct_type = self._xs.struct_type @@ -132,7 +156,17 @@ def __post_init__(self): self.outdir = Path(self.outdir) if self.outdir is None: self.outdir = Path.home() / ".larixite" / "fdmnes" - + #: CONVOLUTION + if self.ecut is None: + self.params["E_cut"] = False + self.ecut = "!" + else: + self.params["E_cut"] = True + if self.gamma_hole is None: + self.params["Gamma_hole"] = False + self.gamma_hole = "!" + else: + self.params["Gamma_hole"] = True #: store a list of jobs self._jobs = [] @@ -227,11 +261,12 @@ def optimize_params(self) -> dict: ) if "mol" in self.struct_type.lower(): + params["Vmax"] = True self.vmax = -6 - logger.info(f"Vmax set to {self.vmax} for molecule") + logger.info(f"Molecule: Vmax enablend and set to {self.vmax}") #: enable SCF - params["SCF"] = True + self.scf = True logger.info("SCF enabled") return params @@ -337,22 +372,6 @@ def get_structure(self, struct_type: str | None = None) -> str: raise AttributeError(errmsg) return "\n".join(structout) - def get_vmax(self) -> str: - """Get the Vmax section of the input""" - if self.vmax is not None: - vmax = ["Vmax"] - vmax.append(f"{self.spacer}-6") - return "\n".join(vmax) - else: - return "! Vmax" - - def get_rself(self): - """Get the R_self section of the input""" - if self.rself is not None: - return f"{self.spacer}{self.rself}" - else: - return "!" - def get_input( self, comment: str = "", @@ -384,12 +403,18 @@ def get_input( "edge": f"{self.spacer}{self.edge}", "radius": f"{self.spacer}{self.radius:.2f}", "erange": self.spacer + self.erange, - "vmax": self.get_vmax(), - "rself": self.get_rself(), + "rself": f"{self.spacer}{self.rself}", "nself": f"{self.spacer}{self.nself}", - "pself": f"{self.spacer}{self.pself:.3f}", + "pself": f"{self.spacer}{self.pself}", + "vmax": f"{self.spacer}{self.vmax}", "struct_type": self.struct_type, "structure": self.get_structure(struct_type=struct_type), + "ecut": f"{self.spacer}{self.ecut}", + "dec": f"{self.spacer}{self.dec}", + "ecent": f"{self.spacer}{self.ecent}", + "elarg": f"{self.spacer}{self.elarg}", + "gamma_hole": f"{self.spacer}{self.gamma_hole}", + "gamma_max": f"{self.spacer}{self.gamma_max}", } for parkey, parval in params.items(): conf[parkey] = str(parkey) if parval is True else f"! {parkey}" diff --git a/larixite/templates/fdmnes_xas.tmpl b/larixite/templates/fdmnes_xas.tmpl index 6a8cfec..86a7a26 100755 --- a/larixite/templates/fdmnes_xas.tmpl +++ b/larixite/templates/fdmnes_xas.tmpl @@ -36,17 +36,32 @@ Radius {SpGr_atom} !`Full_atom` is default since March 2026 ! {SCF} -R_self +{R_self} {rself} -N_self +{N_self} {nself} -P_self +{P_self} {pself} {SCF_exc} +! +{Vmax} {vmax} +! {Atom_conf} ! {structure} ! {Convolution} +{E_cut} +{ecut} +{Dec} +{dec} +{Gamma_hole} +{gamma_hole} +{Ecent} +{ecent} +{Elarg} +{elarg} +{Gamma_max} +{gamma_max} End From a1a3d9679bed36611015e6001ecfe1932d2fe5b3 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Sun, 31 May 2026 17:09:29 +0200 Subject: [PATCH 43/52] init bug fix --- larixite/fdmnes/input.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index c6239cd..be3f8f3 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -120,6 +120,9 @@ def __post_init__(self): self.structpath = self._xs.filepath else: self._xs = get_structure(self.structpath, absorber=self.absorber) + #: parameters + if self.params is None: + self.params = FDMNES_DEFAULT_PARAMS #: radius self.set_radius(self.radius) #: R_self @@ -145,9 +148,6 @@ def __post_init__(self): self.tmplpath = Path(self.tmplpath) #: absorption edge self.validate_edge() - #: parameters - if self.params is None: - self.params = FDMNES_DEFAULT_PARAMS #: optimize params if self.optimize: self.params = self.optimize_params() From 69b3aa86b339606339a10fbe866dcec3c515deda Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Sun, 31 May 2026 17:32:11 +0200 Subject: [PATCH 44/52] implement cif file directly to fdmnes --- larixite/fdmnes/input.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index be3f8f3..cafac88 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -14,8 +14,9 @@ import yaml from dataclasses import dataclass from pathlib import Path -import subprocess import os +import shutil +import subprocess from pymatgen.core import __version__ as pymatgen_version, Element, Molecule from larixite.struct import get_structure, get_structure_from_text from larixite.struct.xas import XasStructure @@ -39,7 +40,7 @@ "Spinorbit": False, "SCF": False, "SCF_exc": False, - "R_self": False, + "R_self": True, "N_self": True, "P_self": True, "Vmax": False, @@ -286,12 +287,12 @@ def get_structure(self, struct_type: str | None = None) -> str: FDMNES supports various structure types: - Crystal -> Implemented (default) - Molecule -> TODO + - Cif_file -> Implemented - Film -> Not implemented yet - Surface -> Not implemented yet - Interface -> Not implemented yet - Pdb_file -> Not implemented yet - Film_Pdb_file -> Not implemented yet - - Cif_file -> Not implemented yet - Film_Cif_file -> Not implemented yet """ @@ -344,6 +345,11 @@ def get_structure(self, struct_type: str | None = None) -> str: ): sitestr = f"{elem.Z:>3d} {site.a:15.10f} {site.b:15.10f} {site.c:15.10f} {occupancy:>5.3f} !{site.label:>4s} {wyckoff:>4s} {elstr:>4s}" structout.append(sitestr) + elif "cif" in struct_type.lower(): + structout.append("Z_absorber") + structout.append(f"{self.spacer}{self.absorber.Z}") + structout.append("Cif_file") + structout.append(f"{self.spacer}{Path(self._xs.filepath).name}") elif "mol" in struct_type.lower(): #: build the cluster and map by distance from absorber at (0,0,0) mol = self._xs.cluster @@ -441,6 +447,10 @@ def write_input( with open(jobdir / "fdmfile.txt", "w") as fp: fp.write(f"1\n{fileout_name}\n") logger.info(f"written {fp.name}") + if self.struct_type is not None and "cif" in self.struct_type.lower(): + src = Path(self._xs.filepath) + shutil.copy2(src, jobdir) + logger.info(f"copied {src.name} to {jobdir}") self._jobs.append(jobdir) _ = self.dump_params(jobdir / f"{self.fileout_prefix}_params.yaml") return jobdir From 7a09abc398eb6d671553853f003dbc830da73f14 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Sun, 31 May 2026 17:47:56 +0200 Subject: [PATCH 45/52] green/scf at input --- larixite/fdmnes/input.py | 28 ++++++++++++---------------- tests/test_fdmnes.py | 19 +++++++++++++++++++ 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index cafac88..b54ba4d 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -84,6 +84,8 @@ class FdmnesXasInput: 0 #: index of the frame inside the structure (e.g. for multi-frame XYZ files) ) radius: float = 7 #: radius of the calculation + green: bool = True #: use muffin-tin approximation + scf: bool = False #: enable SCF erange: str = "-10.0 0.25 60.0 1.0 100.0" #: energy range rself: float | None = None #: radius for the SCF calculation nself: int = 100 #: number of maximum iterations for the SCF calculation @@ -124,6 +126,8 @@ def __post_init__(self): #: parameters if self.params is None: self.params = FDMNES_DEFAULT_PARAMS + self.params["Green"] = self.green + self.params["SCF"] = self.scf #: radius self.set_radius(self.radius) #: R_self @@ -171,26 +175,18 @@ def __post_init__(self): #: store a list of jobs self._jobs = [] + def __setattr__(self, name: str, value): + """Sync green/scf with params dict on assignment""" + super().__setattr__(name, value) + if name == "green" and hasattr(self, "params") and self.params is not None: + self.params["Green"] = value + elif name == "scf" and hasattr(self, "params") and self.params is not None: + self.params["SCF"] = value + @property def xs(self): return self._xs - @property - def green(self): - return self.params["Green"] - - @green.setter - def green(self, val): - self.params["Green"] = val - - @property - def scf(self): - return self.params["SCF"] - - @scf.setter - def scf(self, val): - self.params["SCF"] = val - def set_radius(self, value: float): self.radius = value self._xs.radius = value diff --git a/tests/test_fdmnes.py b/tests/test_fdmnes.py index a4d9f27..63c84ab 100644 --- a/tests/test_fdmnes.py +++ b/tests/test_fdmnes.py @@ -30,5 +30,24 @@ def test_fdmnes(): assert outdir.exists() +def test_green_scf_sync(): + db = get_amcsd() + cif = db.get_cif(4438) + outfile = cif.write_cif(verbose=True) + sg = get_structure(outfile, "Fe") + + f = FdmnesXasInput(sg, absorber="Fe", green=False, scf=True, optimize=False) + assert f.green is False + assert f.scf is True + assert f.params["Green"] is False + assert f.params["SCF"] is True + + f.green = True + f.scf = False + assert f.params["Green"] is True + assert f.params["SCF"] is False + + if __name__ == "__main__": test_fdmnes() + test_green_scf_sync() From 5f25605bf67c1a6e2480ffdaa042ba9572f7ee1f Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Sun, 31 May 2026 18:07:23 +0200 Subject: [PATCH 46/52] fix dec is not value --- larixite/fdmnes/input.py | 4 +--- larixite/templates/fdmnes_xas.tmpl | 1 - 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index b54ba4d..8852564 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -53,7 +53,7 @@ "COOP": False, #: -> TODO "Convolution": True, "E_cut": False, #: taken from E_fermi - "Dec": True, #: energy shift + "Dec": False, #: apply energy shift by E_cut "Ecent": True, #: arctan center "Elarg": True, #: arctan width "Gamma_hole": False, #: use core-hole broadening @@ -102,7 +102,6 @@ class FdmnesXasInput: spacer: str = " " #: spacer string for the FDMNES input text outdir: str | Path | None = None #: path to the output directory of the FDMNES jobs ecut: float | str | None = None #: energy cutoff for the convolution - dec: float = 0 #: energy shift for the convolution ecent: float = 30. #: energy center for the convolution elarg: float = 30. #: energy width for the convolution gamma_hole: float | str | None = None #: start width of the energy convolution (None -> core-hole broadening) @@ -412,7 +411,6 @@ def get_input( "struct_type": self.struct_type, "structure": self.get_structure(struct_type=struct_type), "ecut": f"{self.spacer}{self.ecut}", - "dec": f"{self.spacer}{self.dec}", "ecent": f"{self.spacer}{self.ecent}", "elarg": f"{self.spacer}{self.elarg}", "gamma_hole": f"{self.spacer}{self.gamma_hole}", diff --git a/larixite/templates/fdmnes_xas.tmpl b/larixite/templates/fdmnes_xas.tmpl index 86a7a26..8b344ed 100755 --- a/larixite/templates/fdmnes_xas.tmpl +++ b/larixite/templates/fdmnes_xas.tmpl @@ -55,7 +55,6 @@ Radius {E_cut} {ecut} {Dec} -{dec} {Gamma_hole} {gamma_hole} {Ecent} From f668c8390bdda3075a416827143b45fd659b7a60 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Tue, 2 Jun 2026 09:45:42 +0200 Subject: [PATCH 47/52] fix bug if fileout prefix none --- larixite/fdmnes/input.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index 8852564..a6dc2f0 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -91,12 +91,8 @@ class FdmnesXasInput: nself: int = 100 #: number of maximum iterations for the SCF calculation pself: float = 0.025 #: initial weight for the SCF calculation vmax: float | str | None = None #: maximum potential value for molecules - fileout_prefix: str = ( - "job" #: prefix of the output filename for the FDMNES job (extension: .inp) - ) - tmplpath: str | Path | None = ( - None #: path to the templates directory (FDMNES inputs and sbatch) - ) + fileout_prefix: str | None = None #: prefix for the FDMNES job (deafault: "job" / extension: .inp) + tmplpath: str | Path | None = None #: path to the templates directory params: dict[str, bool] | None = None #: enable/disable parameters for FDMNES optimize: bool = False #: optimize the input parameters spacer: str = " " #: spacer string for the FDMNES input text @@ -122,6 +118,8 @@ def __post_init__(self): self.structpath = self._xs.filepath else: self._xs = get_structure(self.structpath, absorber=self.absorber) + if self.fileout_prefix is None: + self.fileout_prefix = "job" #: parameters if self.params is None: self.params = FDMNES_DEFAULT_PARAMS From c82da5291644365b380a634ad4f47d1c820dc663 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Tue, 2 Jun 2026 10:29:01 +0200 Subject: [PATCH 48/52] add more convolution params --- larixite/fdmnes/input.py | 15 +++++++++++++++ larixite/templates/fdmnes_xas.tmpl | 4 ++++ 2 files changed, 19 insertions(+) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index a6dc2f0..b5539be 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -58,6 +58,7 @@ "Elarg": True, #: arctan width "Gamma_hole": False, #: use core-hole broadening "Gamma_max": True, + "Gaussian": False, } @@ -102,6 +103,8 @@ class FdmnesXasInput: elarg: float = 30. #: energy width for the convolution gamma_hole: float | str | None = None #: start width of the energy convolution (None -> core-hole broadening) gamma_max: float = 8 #: maximum energy for the convolution + gaussian: float | str | None = None #: gaussian broadening (= experimental braodening) + estart: float | str | None = None #: starting energy for the output spectrum def __post_init__(self): """Validate and optimize attributes""" @@ -169,6 +172,16 @@ def __post_init__(self): self.gamma_hole = "!" else: self.params["Gamma_hole"] = True + if self.gaussian is None: + self.params["Gaussian"] = False + self.gaussian = "!" + else: + self.params["Gaussian"] = True + if self.estart is None: + self.params["Estart"] = False + self.estart = "!" + else: + self.params["Estart"] = True #: store a list of jobs self._jobs = [] @@ -413,6 +426,8 @@ def get_input( "elarg": f"{self.spacer}{self.elarg}", "gamma_hole": f"{self.spacer}{self.gamma_hole}", "gamma_max": f"{self.spacer}{self.gamma_max}", + "gaussian": f"{self.spacer}{self.gaussian}", + "estart": f"{self.spacer}{self.estart}", } for parkey, parval in params.items(): conf[parkey] = str(parkey) if parval is True else f"! {parkey}" diff --git a/larixite/templates/fdmnes_xas.tmpl b/larixite/templates/fdmnes_xas.tmpl index 8b344ed..ef955e2 100755 --- a/larixite/templates/fdmnes_xas.tmpl +++ b/larixite/templates/fdmnes_xas.tmpl @@ -63,4 +63,8 @@ Radius {elarg} {Gamma_max} {gamma_max} +{Gaussian} +{gaussian} +{Estart} +{estart} End From b8904a0b207c4eeba04594b0f4fcb00b6c9399f5 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Tue, 2 Jun 2026 11:30:28 +0200 Subject: [PATCH 49/52] force conv_out --- larixite/templates/fdmnes_xas.tmpl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/larixite/templates/fdmnes_xas.tmpl b/larixite/templates/fdmnes_xas.tmpl index ef955e2..e0ed819 100755 --- a/larixite/templates/fdmnes_xas.tmpl +++ b/larixite/templates/fdmnes_xas.tmpl @@ -51,6 +51,8 @@ Radius ! {structure} ! +Conv_out +{filout}_conv {Convolution} {E_cut} {ecut} From cd8be36ceeb19f5cb61929e987e91fd99f483c57 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Wed, 3 Jun 2026 12:03:50 +0200 Subject: [PATCH 50/52] more control on sbatch inputs --- larixite/fdmnes/input.py | 4 ++++ larixite/templates/fdmnes_sbatch_esrf.tmpl | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index b5539be..c32ad5b 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -469,7 +469,9 @@ def write_sbatch( ncpus: int = 12, nnodes: int = 1, mem_per_cpu: str = "16GB", + walltime: str = "8:00:00", constraint: str = "", + partition: str = "nice,nice-long", **kwargs, ) -> Path: """Generates a SBATCH file (SLURM workload manager) using a template @@ -506,7 +508,9 @@ def write_sbatch( "nnodes": nnodes, "ncpus": ncpus, "mem_per_cpu": mem_per_cpu, + "walltime": walltime, "constraint": constraint, + "partition": partition, **kwargs, } with open(sbatchout, "w") as fp, open(template) as tp: diff --git a/larixite/templates/fdmnes_sbatch_esrf.tmpl b/larixite/templates/fdmnes_sbatch_esrf.tmpl index 3d9f39b..a621e3b 100644 --- a/larixite/templates/fdmnes_sbatch_esrf.tmpl +++ b/larixite/templates/fdmnes_sbatch_esrf.tmpl @@ -3,8 +3,8 @@ #SBATCH --ntasks-per-node={ncpus} #: Total number of tasks (=CPUs cores) per node #SBATCH --cpus-per-task=1 #: CPUs for each task / DON'T CHANGE THIS - required by FDMNES #SBATCH --mem-per-cpu={mem_per_cpu} #: Total memory for each CPU -#SBATCH --time=120:00:00 #: Wall-time limit -#SBATCH --partition=nice,nice-long #: ESRF cluster +#SBATCH --time={walltime} #: Wall-time limit +#SBATCH --partition="{partition}" #: ESRF cluster #SBATCH --constraint="{constraint}" #SBATCH --job-name="{jobname}" #: Name used to set the output files #SBATCH --output="{jobname}_slurm.stdout" From bef9f0ed858217ddb0452304eade05e517930926 Mon Sep 17 00:00:00 2001 From: Mauro Rovezzi Date: Tue, 23 Jun 2026 10:21:28 +0200 Subject: [PATCH 51/52] update changes --- CHANGELOG.md | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3705d7b..9509bb1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,9 +14,19 @@ The versioning scheme followed is `YYYY.MAJOR.MINOR`: ## [2026.2.0 - not released yet] -### `fmdnes` +### `fdmnes` + + - Create a `fdmnes` module with refactored input handling (`fdmnes/input.py`) + - Added output parsing and YAML store/restore for input configurations + - Extended sbatch options and improved SLURM job controls + +### `struct` + + - Migration of coordination environment analysis and visualization from `struct2xas` into separate modules `analyze` and `visualize` + +### Tests - - Create a `fdmnes` module + - Added tests for analysis, visualization, and fdmnes input modules ## [2026.1.0 - 18-May-2026] From 80bc4017993b5a25d7333e7b08955ad737bb8a4e Mon Sep 17 00:00:00 2001 From: "Mauro Rovezzi (on lfameangelo)" Date: Fri, 28 Aug 2026 15:54:34 +0200 Subject: [PATCH 52/52] reviewed and ready for merge --- CHANGELOG.md | 18 ++++++++++++++---- larixite/fdmnes/__init__.py | 2 +- larixite/fdmnes/input.py | 1 + 3 files changed, 16 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9509bb1..2627a1f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,13 +16,23 @@ The versioning scheme followed is `YYYY.MAJOR.MINOR`: ### `fdmnes` - - Create a `fdmnes` module with refactored input handling (`fdmnes/input.py`) - - Added output parsing and YAML store/restore for input configurations - - Extended sbatch options and improved SLURM job controls + - Split `fdmnes.py` into a `fdmnes` package (`fdmnes/input.py`) with refactored input handling; `FdmnesXasInput` and `struct2fdmnes` remain importable from `larixite.fdmnes` as before + - Direct `Cif_file` structure type support: a CIF can now be passed straight to FDMNES instead of only an expanded `Crystal` description + - Added convolution parameters (`E_cut`, `Dec`, `Ecent`, `Elarg`, `Gamma_hole`, `Gamma_max`, `Gaussian`, `Estart`) and automatic `Conv_out` section in the generated input + - `green`/`scf` are now plain attributes kept in sync with `params["Green"]`/`params["SCF"]` (previously computed properties) + - Added `run_sbatch()` to submit jobs via `sbatch --parsable` and track status in `status.yaml` + - Added `dump_params()` / `from_yaml()` to store and restore a `FdmnesXasInput` configuration as YAML + - Extended `write_sbatch()` with `nnodes`, `mem_per_cpu`, `walltime`, `constraint`, `partition` options (default `ncpus` raised from 8 to 12) + - Removed the `get_vmax()`/`get_rself()` helper methods; `vmax`/`rself` are now rendered directly from their attribute values (also fixes `Vmax` previously always being written as `-6` regardless of the value set) ### `struct` - - Migration of coordination environment analysis and visualization from `struct2xas` into separate modules `analyze` and `visualize` + - Migration of coordination environment analysis and visualization from `struct2xas` into separate modules `analyze` and `visualize` (new `py3Dmol`-based 3D visualization, new `pyproject.toml` dependency) + +### Bug fixes + + - Fix `FdmnesXasInput(optimize=True)` not actually enabling `SCF` in the generated input despite logging "SCF enabled" (regression introduced during the `fdmnes` refactor, caught before release; affects the web app's FDMNES "optimize" option) + - Restore `from larixite.fdmnes import logger`, broken by the module split (used in `examples/fdmnes_xas_workflow.ipynb`) ### Tests diff --git a/larixite/fdmnes/__init__.py b/larixite/fdmnes/__init__.py index 40e694f..e77f857 100644 --- a/larixite/fdmnes/__init__.py +++ b/larixite/fdmnes/__init__.py @@ -1 +1 @@ -from .input import FdmnesXasInput, struct2fdmnes \ No newline at end of file +from .input import FdmnesXasInput, struct2fdmnes, logger diff --git a/larixite/fdmnes/input.py b/larixite/fdmnes/input.py index c32ad5b..c1272f6 100644 --- a/larixite/fdmnes/input.py +++ b/larixite/fdmnes/input.py @@ -274,6 +274,7 @@ def optimize_params(self) -> dict: #: enable SCF self.scf = True + params["SCF"] = True logger.info("SCF enabled") return params