diff --git a/CHANGELOG.rst b/CHANGELOG.rst index dacf92095..d1d5f9fa3 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -24,6 +24,17 @@ aware this didn't cause problems *yet*, but could misbehave in some environment. +- New features + + - Added VASP HDF5 input support (`vaspout.h5`) for importing force + constants, Born effective charges, and dielectric tensors via + ``ForceConstants.from_vasp``, as well as precalculated phonon mode + data via ``QpointPhononModes.from_vasp`` and + ``QpointFrequencies.from_vasp``. When primitive cell data is + present, a supercell-to-primitive force constant transformation is + attempted: less data is available is available than the equivalent + phonopy scenario, which may impact reliability. + - Compatibility fixes - The [brille] optional dependency group will no longer attempt to diff --git a/euphonic/force_constants.py b/euphonic/force_constants.py index 065312283..fcd5785c3 100644 --- a/euphonic/force_constants.py +++ b/euphonic/force_constants.py @@ -26,7 +26,7 @@ ) from euphonic.qpoint_frequencies import QpointFrequencies from euphonic.qpoint_phonon_modes import QpointPhononModes -from euphonic.readers import castep, phonopy +from euphonic.readers import castep, phonopy, vasp from euphonic.ureg import ureg from euphonic.util import ( _get_supercell_relative_idx, @@ -1875,6 +1875,28 @@ def from_castep(cls, filename: Path | str) -> Self: data = castep.read_interpolation_data(filename) return cls.from_dict(data) + @classmethod + def from_vasp(cls, filename: Path | str) -> Self: + """ + Reads force constants data from a VASP HDF5 file (e.g. vaspout.h5). + + Parameters + ---------- + filename + The path and name of the VASP HDF5 file to read + + Returns + ------- + forceconstants + """ + data = vasp.read_interpolation_data(Path(filename)) + fc = cls.from_dict(data) + if fc.born is not None: + fc = cls.from_total_fc_with_dipole( + fc.crystal, fc.force_constants, fc.sc_matrix, fc.cell_origins, + born=fc.born, dielectric=fc.dielectric) + return fc + @classmethod def from_phonopy(cls, *, diff --git a/euphonic/qpoint_frequencies.py b/euphonic/qpoint_frequencies.py index a1afe8afc..abae6e61b 100644 --- a/euphonic/qpoint_frequencies.py +++ b/euphonic/qpoint_frequencies.py @@ -15,7 +15,7 @@ _obj_to_json_file, _process_dict, ) -from euphonic.readers import castep, phonopy +from euphonic.readers import castep, phonopy, vasp from euphonic.spectra import Spectrum1D, Spectrum1DCollection, Spectrum2D from euphonic.ureg import Quantity, ureg from euphonic.util import ( @@ -464,3 +464,19 @@ def from_phonopy(cls, path=path, phonon_name=phonon_name, phonon_format=phonon_format, summary_name=summary_name, read_eigenvectors=False) return cls.from_dict(data) + + @classmethod + def from_vasp( + cls, + filename: Path | str, + ) -> Self: + """ + Reads phonon frequency data from a VASP HDF5 file (e.g. vaspout.h5) + + Parameters + ---------- + filename + The path and name of the VASP HDF5 file to read + """ + data = vasp.read_phonon_data(Path(filename)) + return cls.from_dict(data) diff --git a/euphonic/qpoint_phonon_modes.py b/euphonic/qpoint_phonon_modes.py index 2298e34ef..c3367bcdb 100644 --- a/euphonic/qpoint_phonon_modes.py +++ b/euphonic/qpoint_phonon_modes.py @@ -21,7 +21,7 @@ AdaptiveMethod, QpointFrequencies, ) -from euphonic.readers import castep, phonopy +from euphonic.readers import castep, phonopy, vasp from euphonic.spectra import Spectrum1DCollection from euphonic.structure_factor import StructureFactor from euphonic.ureg import Quantity, ureg @@ -722,6 +722,22 @@ def from_phonopy( summary_name=summary_name) return cls.from_dict(data) + @classmethod + def from_vasp( + cls, + filename: Path | str, + ) -> Self: + """ + Reads phonon mode data from a VASP HDF5 file (e.g. vaspout.h5) + + Parameters + ---------- + filename + The path and name of the VASP HDF5 file to read + """ + data = vasp.read_phonon_data(Path(filename)) + return cls.from_dict(data) + def _get_isotope_data(scattering_lengths: IsotopeDataset) -> IsotopeData: """Get dataset with coherent_scattering_length for coherent S(q, ω)""" diff --git a/euphonic/readers/phonopy.py b/euphonic/readers/phonopy.py index 27bbeb8b5..75d6403d3 100644 --- a/euphonic/readers/phonopy.py +++ b/euphonic/readers/phonopy.py @@ -17,11 +17,13 @@ class ImportPhonopyReaderError(ModuleNotFoundError): def __init__(self): - self.message = ( - '\n\nCannot import yaml, h5py to read Phonopy files, maybe ' - 'they are not installed. To install the optional ' - "dependencies for Euphonic's Phonopy reader, try:\n\n" - 'pip install euphonic[phonopy-reader]\n') + self.message = format_error( + 'Cannot import yaml and h5py to read Phonopy files.', + fix=( + 'To install optional dependencies for Phonopy reader, try: ' + 'pip install euphonic[phonopy-reader]' + ), + ) def __str__(self): return self.message diff --git a/euphonic/readers/vasp.py b/euphonic/readers/vasp.py new file mode 100644 index 000000000..9b44a4929 --- /dev/null +++ b/euphonic/readers/vasp.py @@ -0,0 +1,577 @@ +from contextlib import contextmanager +from pathlib import Path +import re +from typing import TYPE_CHECKING, TypedDict + +import numpy as np + +from euphonic.util import convert_fc_phases, format_error + +if TYPE_CHECKING: + import h5py + +BOUNDARY_TOLERANCE: float = 1e-12 + + +class CrystalDict(TypedDict): + cell_vectors: np.ndarray + cell_vectors_unit: str + atom_r: np.ndarray + atom_type: np.ndarray + atom_mass: np.ndarray + atom_mass_unit: str + + +class PhononDataDict(TypedDict, total=False): + crystal: CrystalDict + qpts: np.ndarray + frequencies: np.ndarray + frequencies_unit: str + weights: np.ndarray + eigenvectors: np.ndarray + + +class InterpolationDataDict(TypedDict, total=False): + crystal: CrystalDict + force_constants: np.ndarray + force_constants_unit: str + sc_matrix: np.ndarray + cell_origins: np.ndarray + born: np.ndarray + born_unit: str + dielectric: np.ndarray + dielectric_unit: str + + +class BornDict(TypedDict, total=False): + born_raw: np.ndarray + dielectric: np.ndarray + + +class ImportVaspReaderError(ModuleNotFoundError): + """ + Error raised when h5py is required to read VASP HDF5 files but is missing. + """ + + def __init__(self) -> None: + self.message = format_error( + 'Cannot import h5py to read VASP HDF5 files.', + fix=( + 'To install optional HDF5 dependencies for Euphonic, try: ' + 'pip install euphonic[phonopy-reader]' + ), + ) + + def __str__(self) -> str: + return self.message + + +class MissingPhononModesError(KeyError): + """ + Error raised when precalculated phonon modes/frequencies are missing. + """ + + +class MissingPrimitiveCellError(KeyError): + """ + Error raised when primitive cell structure is missing in VASP HDF5 file. + """ + + +def _normalize_fractional_coords(pos: np.ndarray) -> np.ndarray: + """ + Normalizes fractional atomic positions to [0.0, 1.0), snapping boundary + values within BOUNDARY_TOLERANCE of 1.0 or 0.0 to 0.0. + """ + atom_r = pos % 1.0 + near_boundary = np.isclose( + atom_r, 1.0, atol=BOUNDARY_TOLERANCE + ) | np.isclose(atom_r, 0.0, atol=BOUNDARY_TOLERANCE) + atom_r[near_boundary] = 0.0 + return atom_r + + +@contextmanager +def _open_vasp_h5(filename: Path): + """ + Context manager to open a VASP HDF5 file with error handling for h5py. + """ + try: + import h5py + except (ModuleNotFoundError, ImportError) as err: + raise ImportVaspReaderError from err + + filepath = Path(filename) + if not filepath.exists(): + msg = format_error( + f'VASP file not found at {filepath}.', + fix='Provide a valid path to an existing VASP HDF5 file.', + ) + raise FileNotFoundError(msg) + + with h5py.File(filepath, 'r') as h5_file: + yield h5_file + + +def _extract_pomass(h5_file: 'h5py.File') -> list[float]: + """ + Extracts atomic masses (POMASS) per species from INCAR or POTCAR + content datasets stored inside the VASP HDF5 file. + + Parameters + ---------- + h5_file + Opened h5py.File object representing the VASP HDF5 container + + Returns + ------- + masses_per_type + List of float atomic masses in amu per species + + Raises + ------ + ValueError + If atomic masses (POMASS) cannot be found in INCAR or POTCAR datasets. + """ + filename = Path(h5_file.filename) + + # 1. Try active input/incar/POMASS + if 'input/incar/POMASS' in h5_file: + val = h5_file['input/incar/POMASS'].asstr()[()] + raw_vals = re.findall(r'[0-9.]+', str(val)) + if not raw_vals: + msg = format_error( + f'POMASS found in input/incar/POMASS but could not parse ' + f'numeric values from: {val!r}', + fix='Ensure POMASS contains valid numeric values.', + ) + raise ValueError(msg) + return [float(mass) for mass in raw_vals] + + # 2. Try original/incar/content + if 'original/incar/content' in h5_file: + incar_content = h5_file['original/incar/content'].asstr()[()] + match = re.search( + r'POMASS\s*=\s*(?P[0-9.\s,]+)', incar_content + ) + if match: + raw_vals = re.findall(r'[0-9.]+', match.group('masses')) + if not raw_vals: + msg = format_error( + f'POMASS found in original/incar/content but could not ' + f'parse numeric values from: {match.group("masses")!r}', + fix='Ensure POMASS contains valid numeric values.', + ) + raise ValueError(msg) + return [float(mass) for mass in raw_vals] + + # 3. Try POTCAR content stored inside the HDF5 file + if 'input/potcar/content' in h5_file: + potcar_content = h5_file['input/potcar/content'].asstr()[()] + matches = re.findall(r'POMASS\s*=\s*(?P[0-9.]+)', potcar_content) + if matches: + return [float(mass) for mass in matches] + + # 4. If missing from all, raise error + msg = format_error( + f'Could not find atomic masses (POMASS) in {filename}.', + fix=( + 'Ensure POMASS is set in INCAR, or that POTCAR contains ' + 'POMASS data.' + ), + ) + raise ValueError(msg) + + +def _read_cell_from_group( + h5_file: 'h5py.File', group_path: str +) -> CrystalDict: + """ + Helper function to parse crystal structure from a specific HDF5 group. + """ + pos_group = h5_file[group_path] + latt = pos_group['lattice_vectors'][()] + pos = pos_group['position_ions'][()] + species_counts = pos_group['number_ion_types'][()] + species_types = pos_group['ion_types'].asstr()[()] + + species_masses = _extract_pomass(h5_file) + + atom_type = np.repeat(species_types, species_counts) + atom_mass = np.repeat(species_masses, species_counts) + + atom_r = _normalize_fractional_coords(pos) + + return { + 'cell_vectors': latt, + 'cell_vectors_unit': 'angstrom', + 'atom_r': atom_r, + 'atom_type': atom_type, + 'atom_mass': atom_mass, + 'atom_mass_unit': 'amu', + } + + +def read_cell(filename: Path) -> CrystalDict: + """ + Reads calculation cell structure from input/poscar in VASP HDF5 file. + + Parameters + ---------- + filename + Path to the VASP HDF5 output file + + Returns + ------- + crystal_dict + A CrystalDict for the calculation cell + + Raises + ------ + FileNotFoundError + If the file does not exist at filename. + ImportVaspReaderError + If h5py is not installed. + KeyError + If input/poscar group is missing from the file. + ValueError + If atomic masses (POMASS) cannot be found in INCAR or POTCAR. + """ + with _open_vasp_h5(filename) as h5_file: + if 'input/poscar' in h5_file: + return _read_cell_from_group(h5_file, 'input/poscar') + + msg = format_error( + f'Crystal position data not found in {filename}.', + fix='Ensure the file contains input/poscar group.', + ) + raise KeyError(msg) + + +def read_primitive_cell(filename: Path) -> CrystalDict: + """ + Reads primitive cell structure from results/phonons/primitive. + + Parameters + ---------- + filename + Path to the VASP HDF5 output file + + Returns + ------- + crystal_dict + A CrystalDict for the primitive cell + + Raises + ------ + MissingPrimitiveCellError + If primitive cell structure is not found in the file. + FileNotFoundError + If the file does not exist at filename. + ImportVaspReaderError + If h5py is not installed. + ValueError + If atomic masses (POMASS) cannot be found in INCAR or POTCAR. + """ + with _open_vasp_h5(filename) as h5_file: + if 'results/phonons/primitive' in h5_file: + return _read_cell_from_group(h5_file, 'results/phonons/primitive') + + msg = format_error( + f'Primitive cell structure not found in {filename}.', + fix='Ensure the file contains results/phonons/primitive data.', + ) + raise MissingPrimitiveCellError(msg) + + +def read_phonon_data(filename: Path) -> PhononDataDict: + """ + Reads precalculated phonon mode/band data from a VASP HDF5 file in native + THz units. + + Parameters + ---------- + filename + Path to the VASP HDF5 file + + Returns + ------- + data_dict + A PhononDataDict with keys: 'crystal', 'qpts', 'frequencies', + 'frequencies_unit', 'eigenvectors', 'weights' + + Raises + ------ + MissingPhononModesError + If precalculated phonon mode/band data is not found in the file. + FileNotFoundError + If the file does not exist at filename. + ImportVaspReaderError + If h5py is not installed. + """ + with _open_vasp_h5(filename) as h5_file: + if 'results/phonons/frequencies' not in h5_file: + msg = format_error( + f'Pre-calculated phonon band data not found in {filename}.', + fix='Use ForceConstants.from_vasp to read force constants.', + ) + raise MissingPhononModesError(msg) + + phonon_group = h5_file['results/phonons'] + + try: + crystal_dict = read_primitive_cell(filename) + except MissingPrimitiveCellError: + crystal_dict = read_cell(filename) + + n_atoms = len(crystal_dict['atom_r']) + + qpts = phonon_group['qpoint_coords'][()] + freqs = phonon_group['frequencies'][()] + weights = phonon_group['qpoints_symmetry_weight'][()] + evecs_raw = phonon_group['eigenvectors'][()] + + expected_shape = (len(qpts), 3 * n_atoms, n_atoms, 3, 2) + if evecs_raw.shape != expected_shape: + msg = format_error( + f'Unexpected eigenvector array shape {evecs_raw.shape} ' + f'in {filename} (expected {expected_shape}).', + fix='Ensure the file contains valid VASP 6 eigenvectors.', + ) + raise ValueError(msg) + + evecs_complex = ( + np.ascontiguousarray(evecs_raw).view(dtype=complex).squeeze(axis=-1) + ) + + return { + 'crystal': crystal_dict, + 'qpts': qpts, + 'frequencies': freqs, + 'frequencies_unit': 'THz', + 'weights': weights, + 'eigenvectors': evecs_complex, + } + + +def _find_fc_key(h5_file: 'h5py.File') -> str: + """ + Finds dataset key for force constants or Hessian matrix in HDF5 container. + + Raises + ------ + KeyError + If results/linear_response force constants group is missing. + """ + if 'results/linear_response/force_constants' in h5_file: + return 'results/linear_response/force_constants' + if 'results/linear_response/hessian' in h5_file: + return 'results/linear_response/hessian' + + filename = Path(h5_file.filename) + msg = format_error( + f'Force constants not found in {filename}.', + fix=( + 'Ensure the file contains results/linear_response ' + 'force constants.' + ), + ) + raise KeyError(msg) + + +def _extract_born_and_dielectric(h5_file: 'h5py.File') -> BornDict: + """ + Extracts Born charges and electronic dielectric tensor from HDF5 container. + """ + born_dict = {} + if 'results/linear_response/born_charges' in h5_file: + born_dict['born_raw'] = h5_file[ + 'results/linear_response/born_charges' + ][()] + + if 'results/linear_response/electron_dielectric_tensor' in h5_file: + born_dict['dielectric'] = h5_file[ + 'results/linear_response/electron_dielectric_tensor' + ][()] + + return born_dict + + +def _build_supercell_data( + crystal_dict: CrystalDict, + sc_hessian: np.ndarray, + born_dict: BornDict, +) -> InterpolationDataDict: + """ + Builds interpolation data dictionary for supercell-as-unit-cell fallback. + + Parameters + ---------- + crystal_dict + Crystal dictionary for the supercell (treated as unit cell) + sc_hessian + Raw supercell Hessian matrix of shape + (3 * n_atoms, 3 * n_atoms) in eV/Angstrom^2 + born_dict + Dictionary containing optional born charges and dielectric tensor + """ + uc_n_atoms = len(crystal_dict['atom_r']) + assert sc_hessian.shape == (3 * uc_n_atoms, 3 * uc_n_atoms), ( + f'Expected sc_hessian shape ({3 * uc_n_atoms}, {3 * uc_n_atoms}), ' + f'got {sc_hessian.shape}' + ) + fc = -sc_hessian.reshape(1, 3 * uc_n_atoms, 3 * uc_n_atoms) + + result: InterpolationDataDict = { + 'crystal': crystal_dict, + 'force_constants': fc, + 'force_constants_unit': 'eV/angstrom**2', + 'sc_matrix': np.eye(3, dtype=int), + 'cell_origins': np.zeros((1, 3), dtype=int), + } + + if 'born_raw' in born_dict: + result['born'] = born_dict['born_raw'] + result['born_unit'] = 'e' + + if 'dielectric' in born_dict: + result['dielectric'] = born_dict['dielectric'] + result['dielectric_unit'] = '(e**2)/(bohr*hartree)' + + return result + + +def _build_primitive_data( + h5_file: 'h5py.File', + crystal_dict: CrystalDict, + sc_hessian: np.ndarray, + born_dict: BornDict, +) -> InterpolationDataDict: + """ + Builds interpolation data dictionary by transforming supercell + force constants to primitive cell coordinates and cell origins. + + Parameters + ---------- + h5_file + Opened h5py.File object representing the VASP HDF5 container + crystal_dict + Crystal dictionary for the primitive cell + sc_hessian + Raw supercell Hessian matrix of shape + (3 * sc_n_atoms, 3 * sc_n_atoms) in eV/Angstrom^2 + born_dict + Dictionary containing optional born charges and dielectric tensor + """ + uc_n_atoms = len(crystal_dict['atom_r']) + + prim_group = h5_file['results/phonons/primitive'] + prim_lat = prim_group['lattice_vectors'][()] + prim_pos = prim_group['position_ions'][()] + atom_r = _normalize_fractional_coords(prim_pos) + + # Use input/poscar for equilibrium structure (results/positions may contain + # displaced positions from finite-difference Hessian calculation) + pos_group = h5_file['input/poscar'] + sc_lat = pos_group['lattice_vectors'][()] + sc_pos = pos_group['position_ions'][()] + sc_n_atoms = len(sc_pos) + + assert sc_hessian.shape == (3 * sc_n_atoms, 3 * sc_n_atoms), ( + f'Expected sc_hessian shape ({3 * sc_n_atoms}, {3 * sc_n_atoms}), ' + f'got {sc_hessian.shape}' + ) + + exact_sc_matrix = sc_lat @ np.linalg.inv(prim_lat) + sc_matrix = np.rint(exact_sc_matrix).astype(int) + sc_atom_r = sc_pos @ exact_sc_matrix + + cell_origins_per_atom = np.floor(sc_atom_r + 1e-5).astype(int) + r_in_p = _normalize_fractional_coords( + sc_atom_r - cell_origins_per_atom + ) + + # Map each supercell atom index (0 to sc_n_atoms - 1) to its + # equivalent primitive unit cell atom index (0 to uc_n_atoms - 1) + sc_to_uc_atom_idx = np.zeros(sc_n_atoms, dtype=int) + for i, pos in enumerate(r_in_p): + diffs = np.linalg.norm(atom_r - pos, axis=1) + sc_to_uc_atom_idx[i] = np.argmin(diffs) + + # Map each primitive unit cell atom index to one corresponding + # supercell atom index (used for indexing Born charges) + uc_to_sc_atom_idx = np.zeros(uc_n_atoms, dtype=int) + for k in range(uc_n_atoms): + uc_to_sc_atom_idx[k] = np.where(sc_to_uc_atom_idx == k)[0][0] + + fc_4d = -sc_hessian.reshape(sc_n_atoms, 3, sc_n_atoms, 3).transpose( + 0, 2, 1, 3 + ) + + fc_converted, cell_origins = convert_fc_phases( + fc_4d, + atom_r, + sc_atom_r, + uc_to_sc_atom_idx, + sc_to_uc_atom_idx, + sc_matrix, + ) + + result: InterpolationDataDict = { + 'crystal': crystal_dict, + 'force_constants': fc_converted, + 'force_constants_unit': 'eV/angstrom**2', + 'sc_matrix': sc_matrix, + 'cell_origins': cell_origins, + } + + if 'born_raw' in born_dict: + result['born'] = born_dict['born_raw'][uc_to_sc_atom_idx] + result['born_unit'] = 'e' + + if 'dielectric' in born_dict: + result['dielectric'] = born_dict['dielectric'] + result['dielectric_unit'] = '(e**2)/(bohr*hartree)' + + return result + + +def read_interpolation_data(filename: Path) -> InterpolationDataDict: + """ + Reads force constants, Born charges, dielectric tensor, and crystal + structure data from a VASP HDF5 file in native VASP units. + + Parameters + ---------- + filename + Path to the VASP HDF5 file + + Returns + ------- + data_dict + A dict with keys: 'crystal', 'force_constants', 'force_constants_unit', + 'sc_matrix', 'cell_origins'. Also optionally contains 'born', + 'born_unit', 'dielectric', and 'dielectric_unit' if present. + + Raises + ------ + KeyError + If results/linear_response force constants group is missing. + FileNotFoundError + If the file does not exist at filename. + ImportVaspReaderError + If h5py is not installed. + """ + with _open_vasp_h5(filename) as h5_file: + fc_key = _find_fc_key(h5_file) + sc_hessian = h5_file[fc_key][()] + born_dict = _extract_born_and_dielectric(h5_file) + + try: + crystal_dict = read_primitive_cell(filename) + return _build_primitive_data( + h5_file, crystal_dict, sc_hessian, born_dict + ) + except MissingPrimitiveCellError: + crystal_dict = read_cell(filename) + return _build_supercell_data(crystal_dict, sc_hessian, born_dict) diff --git a/meson.build b/meson.build index ad45cc551..90bb3da9a 100644 --- a/meson.build +++ b/meson.build @@ -1,15 +1,16 @@ project('euphonic', - 'c', version: run_command('python', 'build_utils/version.py', check: true).stdout().strip(), meson_version: '>=1.6', ) +add_languages('c', required: false) + build = get_option('python_only') ? disabler() : [] fs = import('fs') py = import('python').find_installation(pure: false) -py_dep = py.dependency() +py_dep = get_option('python_only') ? disabler() : py.dependency() py_src = { 'euphonic': ['__init__.py', 'brille.py', 'broadening.py', @@ -26,7 +27,7 @@ py_src = { '_grids.py', '_kwargs.py', '_loaders.py', '_pdos.py', '_plotting.py'], 'euphonic/isotopes': ['__init__.py', '_core.py', '_csv.py', '_legacy.py'], 'euphonic/isotopes/data': ['__init__.py', 'bluebook.json', 'sears-1992.json', 'sears-1992.csv'], - 'euphonic/readers': ['__init__.py', 'castep.py', 'phonopy.py'], + 'euphonic/readers': ['__init__.py', 'castep.py', 'phonopy.py', 'vasp.py'], 'euphonic/spectra': ['__init__.py', 'base.py', 'collections.py'], 'euphonic/styles': ['__init__.py', 'base.mplstyle', 'intensity_widget.mplstyle'], 'euphonic/ureg': ['__init__.py'], @@ -59,12 +60,14 @@ if not np.found() # Try default np = declare_dependency(include_directories: np_inc) endif -openmp = dependency('openmp', required: true, language: 'c') +if not get_option('python_only') + openmp = dependency('openmp', required: true, language: 'c') -py.extension_module( - '_euphonic', - src, - dependencies: [build, py_dep, np, openmp], - install: true, - subdir: 'euphonic', -) + py.extension_module( + '_euphonic', + src, + dependencies: [build, py_dep, np, openmp], + install: true, + subdir: 'euphonic', + ) +endif diff --git a/pyproject.toml b/pyproject.toml index da2a53431..b5e726ccd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ requires = ["meson-python", "numpy>=1.24.0"] [tool.meson-python.args] -setup = ['--vsenv'] +setup = ['--vsenv', '-Dpython_only=false'] [project] name = "Euphonic" @@ -84,7 +84,7 @@ dev = [ [project.optional-dependencies] matplotlib = ["matplotlib>=3.8.0"] -phonopy_reader = ["h5py>=3.6.0", "PyYAML>=6.0"] # Deprecated, will be removed in future versions. +# phonopy_reader = ["h5py>=3.6.0", "PyYAML>=6.0"] # Deprecated, will be removed in future versions. phonopy-reader = ["h5py>=3.6.0", "PyYAML>=6.0"] # brille build on ARM-Linux is currently a bad time; skip it so tests run properly brille = ["brille>=0.7.0; sys_platform != 'linux' or platform_machine not in 'aarch64 arm64'"] @@ -98,6 +98,13 @@ euphonic-show-sampling = "euphonic.cli.show_sampling:main" euphonic-intensity-map = "euphonic.cli.intensity_map:main" euphonic-powder-map = "euphonic.cli.powder_map:main" +[tool.uv] +default-groups = ["test"] +conflicts = [ + [{group = 'min_reqs'}, {group = 'docs'}], + [{group = 'min_reqs'}, {group = 'dev'}], +] + [tool.pytest.ini_options] markers = [ "brille: test requires 'brille' extra", @@ -105,6 +112,7 @@ markers = [ "matplotlib: test requires 'matplotlib' extra", "multiple_extras: test requires multiple extras to be installed", "phonopy_reader: test requires 'phonopy-reader' extra", + "vasp_reader: test requires 'phonopy-reader' extra", ] [tool.ruff] @@ -217,12 +225,6 @@ extend-ignore-names = ["k_B", "H_ab", "Q"] # Prefer single quotes over double quotes. quote-style = "single" -[tool.uv] -conflicts = [ - [{group = 'min_reqs'}, {group = 'docs'}], - [{group = 'min_reqs'}, {group = 'dev'}], -] - [tool.coverage.run] branch = true source = ["euphonic"] diff --git a/tests_and_analysis/test/data/crystal/crystal_vasp_fc_and_qpts_cell.json b/tests_and_analysis/test/data/crystal/crystal_vasp_fc_and_qpts_cell.json new file mode 100644 index 000000000..d83d13cd8 --- /dev/null +++ b/tests_and_analysis/test/data/crystal/crystal_vasp_fc_and_qpts_cell.json @@ -0,0 +1,140 @@ +{ + "cell_vectors": [ + [ + 0.0, + 5.7626884882564005, + 5.7626884882564005 + ], + [ + 5.7626884882564005, + 0.0, + 5.7626884882564005 + ], + [ + 5.7626884882564005, + 5.7626884882564005, + 0.0 + ] + ], + "cell_vectors_unit": "angstrom", + "n_atoms": 16, + "atom_r": [ + [ + 0.0, + 0.0, + 0.0 + ], + [ + 0.5, + 0.0, + 0.0 + ], + [ + 0.0, + 0.5, + 0.0 + ], + [ + 0.5, + 0.5, + 0.0 + ], + [ + 0.0, + 0.0, + 0.5 + ], + [ + 0.5, + 0.0, + 0.5 + ], + [ + 0.0, + 0.5, + 0.5 + ], + [ + 0.5, + 0.5, + 0.5 + ], + [ + 0.375, + 0.375, + 0.375 + ], + [ + 0.875, + 0.375, + 0.375 + ], + [ + 0.375, + 0.875, + 0.375 + ], + [ + 0.875, + 0.875, + 0.375 + ], + [ + 0.375, + 0.375, + 0.875 + ], + [ + 0.875, + 0.375, + 0.875 + ], + [ + 0.375, + 0.875, + 0.875 + ], + [ + 0.875, + 0.875, + 0.875 + ] + ], + "atom_type": [ + "Ga", + "Ga", + "Ga", + "Ga", + "Ga", + "Ga", + "Ga", + "Ga", + "As", + "As", + "As", + "As", + "As", + "As", + "As", + "As" + ], + "atom_mass": [ + 69.723, + 69.723, + 69.723, + 69.723, + 69.723, + 69.723, + 69.723, + 69.723, + 74.922, + 74.922, + 74.922, + 74.922, + 74.922, + 74.922, + 74.922, + 74.922 + ], + "atom_mass_unit": "amu" +} \ No newline at end of file diff --git a/tests_and_analysis/test/data/crystal/crystal_vasp_fc_and_qpts_prim.json b/tests_and_analysis/test/data/crystal/crystal_vasp_fc_and_qpts_prim.json new file mode 100644 index 000000000..94a1c0cff --- /dev/null +++ b/tests_and_analysis/test/data/crystal/crystal_vasp_fc_and_qpts_prim.json @@ -0,0 +1,42 @@ +{ + "cell_vectors": [ + [ + 0.0, + 2.8813442441282002, + 2.8813442441282002 + ], + [ + 2.8813442441282002, + 0.0, + 2.8813442441282002 + ], + [ + 2.8813442441282002, + 2.8813442441282002, + 0.0 + ] + ], + "cell_vectors_unit": "angstrom", + "n_atoms": 2, + "atom_r": [ + [ + 0.0, + 0.0, + 0.0 + ], + [ + 0.7499999999999998, + 0.7499999999999998, + 0.7499999999999998 + ] + ], + "atom_type": [ + "Ga", + "As" + ], + "atom_mass": [ + 69.723, + 74.922 + ], + "atom_mass_unit": "amu" +} \ No newline at end of file diff --git a/tests_and_analysis/test/data/crystal/crystal_vasp_fc_no_qpts.json b/tests_and_analysis/test/data/crystal/crystal_vasp_fc_no_qpts.json new file mode 100644 index 000000000..d83d13cd8 --- /dev/null +++ b/tests_and_analysis/test/data/crystal/crystal_vasp_fc_no_qpts.json @@ -0,0 +1,140 @@ +{ + "cell_vectors": [ + [ + 0.0, + 5.7626884882564005, + 5.7626884882564005 + ], + [ + 5.7626884882564005, + 0.0, + 5.7626884882564005 + ], + [ + 5.7626884882564005, + 5.7626884882564005, + 0.0 + ] + ], + "cell_vectors_unit": "angstrom", + "n_atoms": 16, + "atom_r": [ + [ + 0.0, + 0.0, + 0.0 + ], + [ + 0.5, + 0.0, + 0.0 + ], + [ + 0.0, + 0.5, + 0.0 + ], + [ + 0.5, + 0.5, + 0.0 + ], + [ + 0.0, + 0.0, + 0.5 + ], + [ + 0.5, + 0.0, + 0.5 + ], + [ + 0.0, + 0.5, + 0.5 + ], + [ + 0.5, + 0.5, + 0.5 + ], + [ + 0.375, + 0.375, + 0.375 + ], + [ + 0.875, + 0.375, + 0.375 + ], + [ + 0.375, + 0.875, + 0.375 + ], + [ + 0.875, + 0.875, + 0.375 + ], + [ + 0.375, + 0.375, + 0.875 + ], + [ + 0.875, + 0.375, + 0.875 + ], + [ + 0.375, + 0.875, + 0.875 + ], + [ + 0.875, + 0.875, + 0.875 + ] + ], + "atom_type": [ + "Ga", + "Ga", + "Ga", + "Ga", + "Ga", + "Ga", + "Ga", + "Ga", + "As", + "As", + "As", + "As", + "As", + "As", + "As", + "As" + ], + "atom_mass": [ + 69.723, + 69.723, + 69.723, + 69.723, + 69.723, + 69.723, + 69.723, + 69.723, + 74.922, + 74.922, + 74.922, + 74.922, + 74.922, + 74.922, + 74.922, + 74.922 + ], + "atom_mass_unit": "amu" +} \ No newline at end of file diff --git a/tests_and_analysis/test/data/crystal/crystal_vasp_only_qpts_prim.json b/tests_and_analysis/test/data/crystal/crystal_vasp_only_qpts_prim.json new file mode 100644 index 000000000..94a1c0cff --- /dev/null +++ b/tests_and_analysis/test/data/crystal/crystal_vasp_only_qpts_prim.json @@ -0,0 +1,42 @@ +{ + "cell_vectors": [ + [ + 0.0, + 2.8813442441282002, + 2.8813442441282002 + ], + [ + 2.8813442441282002, + 0.0, + 2.8813442441282002 + ], + [ + 2.8813442441282002, + 2.8813442441282002, + 0.0 + ] + ], + "cell_vectors_unit": "angstrom", + "n_atoms": 2, + "atom_r": [ + [ + 0.0, + 0.0, + 0.0 + ], + [ + 0.7499999999999998, + 0.7499999999999998, + 0.7499999999999998 + ] + ], + "atom_type": [ + "Ga", + "As" + ], + "atom_mass": [ + 69.723, + 74.922 + ], + "atom_mass_unit": "amu" +} \ No newline at end of file diff --git a/tests_and_analysis/test/data/vasp_files/vaspout_al_no_born.h5 b/tests_and_analysis/test/data/vasp_files/vaspout_al_no_born.h5 new file mode 100644 index 000000000..c602226f6 Binary files /dev/null and b/tests_and_analysis/test/data/vasp_files/vaspout_al_no_born.h5 differ diff --git a/tests_and_analysis/test/data/vasp_files/vaspout_dos_rerun_sanitized.h5 b/tests_and_analysis/test/data/vasp_files/vaspout_dos_rerun_sanitized.h5 new file mode 100644 index 000000000..b05daf447 Binary files /dev/null and b/tests_and_analysis/test/data/vasp_files/vaspout_dos_rerun_sanitized.h5 differ diff --git a/tests_and_analysis/test/data/vasp_files/vaspout_dos_sanitized.h5 b/tests_and_analysis/test/data/vasp_files/vaspout_dos_sanitized.h5 new file mode 100644 index 000000000..3cfdf437e Binary files /dev/null and b/tests_and_analysis/test/data/vasp_files/vaspout_dos_sanitized.h5 differ diff --git a/tests_and_analysis/test/data/vasp_files/vaspout_sanitized.h5 b/tests_and_analysis/test/data/vasp_files/vaspout_sanitized.h5 new file mode 100644 index 000000000..7f1afe83f Binary files /dev/null and b/tests_and_analysis/test/data/vasp_files/vaspout_sanitized.h5 differ diff --git a/tests_and_analysis/test/euphonic_test/test_vasp_reader.py b/tests_and_analysis/test/euphonic_test/test_vasp_reader.py new file mode 100644 index 000000000..9b4ea6e09 --- /dev/null +++ b/tests_and_analysis/test/euphonic_test/test_vasp_reader.py @@ -0,0 +1,544 @@ +import builtins +import json + +import numpy as np +from numpy.testing import assert_allclose +import pytest + +from euphonic import ForceConstants, QpointFrequencies, QpointPhononModes +from euphonic.readers.vasp import ( + ImportVaspReaderError, + MissingPhononModesError, + MissingPrimitiveCellError, + read_cell, + read_interpolation_data, + read_phonon_data, + read_primitive_cell, +) +from tests_and_analysis.test.euphonic_test.test_crystal import ( + ExpectedCrystal, + check_crystal, +) +from tests_and_analysis.test.utils import get_data_path + +FC_NO_QPTS_H5 = get_data_path('vasp_files', 'vaspout_sanitized.h5') +ONLY_QPTS_H5 = get_data_path('vasp_files', 'vaspout_dos_sanitized.h5') +FC_AND_QPTS_H5 = get_data_path( + 'vasp_files', 'vaspout_dos_rerun_sanitized.h5' +) +# Non-diagonal supercell without Born charges (Al FCC, 32 primitive cells) +AL_NO_BORN_H5 = get_data_path('vasp_files', 'vaspout_al_no_born.h5') + + +def get_crystal_path(*subpaths): + return get_data_path('crystal', *subpaths) + + +@pytest.mark.vasp_reader +class TestVaspReaderCell: + + def test_read_cell_fc_no_qpts(self): + cell_dict = read_cell(FC_NO_QPTS_H5) + cell_data = ExpectedCrystal( + {**cell_dict, 'n_atoms': len(cell_dict['atom_r'])} + ) + with open(get_crystal_path('crystal_vasp_fc_no_qpts.json')) as fp: + expected = ExpectedCrystal(json.load(fp)) + check_crystal(cell_data, expected) + + def test_read_primitive_cell_only_qpts(self): + prim_dict = read_primitive_cell(ONLY_QPTS_H5) + prim_data = ExpectedCrystal( + {**prim_dict, 'n_atoms': len(prim_dict['atom_r'])} + ) + with open(get_crystal_path('crystal_vasp_only_qpts_prim.json')) as fp: + expected = ExpectedCrystal(json.load(fp)) + check_crystal(prim_data, expected) + + def test_read_combined_cells(self): + cell_dict = read_cell(FC_AND_QPTS_H5) + cell_data = ExpectedCrystal( + {**cell_dict, 'n_atoms': len(cell_dict['atom_r'])} + ) + with open( + get_crystal_path('crystal_vasp_fc_and_qpts_cell.json') + ) as fp: + exp_cell = ExpectedCrystal(json.load(fp)) + check_crystal(cell_data, exp_cell) + + prim_dict = read_primitive_cell(FC_AND_QPTS_H5) + prim_data = ExpectedCrystal( + {**prim_dict, 'n_atoms': len(prim_dict['atom_r'])} + ) + with open( + get_crystal_path('crystal_vasp_fc_and_qpts_prim.json') + ) as fp: + exp_prim = ExpectedCrystal(json.load(fp)) + check_crystal(prim_data, exp_prim) + + def test_read_cell_from_incar_override(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_vaspout_incar.h5' + with h5py.File(dummy_h5, 'w') as f: + pos_group = f.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset('position_ions', data=np.zeros((1, 3))) + pos_group.create_dataset('number_ion_types', data=np.array([1])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + incar_group = f.create_group('original/incar') + incar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855\nISMEAR = 0') + ) + + data = read_cell(dummy_h5) + assert_allclose(data['atom_mass'], 28.0855) + + def test_read_cell_incar_overrides_potcar(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_vaspout_priority.h5' + with h5py.File(dummy_h5, 'w') as f: + pos_group = f.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset('position_ions', data=np.zeros((1, 3))) + pos_group.create_dataset('number_ion_types', data=np.array([1])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + potcar_group = f.create_group('input/potcar') + potcar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 10.0') + ) + + incar_group = f.create_group('original/incar') + incar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855') + ) + + data = read_cell(dummy_h5) + assert_allclose(data['atom_mass'], 28.0855) + + def test_read_cell_input_incar_overrides_original_incar(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_vaspout_input_incar.h5' + with h5py.File(dummy_h5, 'w') as f: + pos_group = f.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset('position_ions', data=np.zeros((1, 3))) + pos_group.create_dataset('number_ion_types', data=np.array([1])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + potcar_group = f.create_group('input/potcar') + potcar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 10.0') + ) + + incar_group = f.create_group('original/incar') + incar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 20.0') + ) + + input_incar = f.create_group('input/incar') + input_incar.create_dataset('POMASS', data=np.bytes_(b'30.0')) + + data = read_cell(dummy_h5) + assert_allclose(data['atom_mass'], 30.0) + + def test_read_cell_negative_positions(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_vaspout_neg.h5' + with h5py.File(dummy_h5, 'w') as f: + pos_group = f.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset( + 'position_ions', + data=np.array([[-0.1, -0.5, -1.0], [1.1, -1e-15, 0.25]]), + ) + pos_group.create_dataset('number_ion_types', data=np.array([2])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + incar_group = f.create_group('original/incar') + incar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855') + ) + + data = read_cell(dummy_h5) + assert_allclose(data['atom_r'][0], [0.9, 0.5, 0.0]) + assert_allclose(data['atom_r'][1], [0.1, 0.0, 0.25]) + + def test_read_primitive_cell_missing_raises_error(self): + with pytest.raises(MissingPrimitiveCellError): + read_primitive_cell(FC_NO_QPTS_H5) + + def test_read_cell_missing_pomass_raises_error(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_vaspout_empty.h5' + with h5py.File(dummy_h5, 'w') as f: + pos_group = f.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset('position_ions', data=np.zeros((1, 3))) + pos_group.create_dataset('number_ion_types', data=np.array([1])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + with pytest.raises(ValueError, match='Could not find atomic masses'): + read_cell(dummy_h5) + + def test_read_cell_input_incar_unparseable_pomass(self, tmp_path): + """POMASS in input/incar/POMASS exists but can't parse numeric values.""" + import h5py + + dummy_h5 = tmp_path / 'dummy_bad_input_incar_pomass.h5' + with h5py.File(dummy_h5, 'w') as f: + self._create_minimal_poscar(f) + f.create_dataset('input/incar/POMASS', data=np.bytes_(b'invalid')) + + with pytest.raises(ValueError, match='could not parse.*numeric'): + read_cell(dummy_h5) + + def test_read_cell_original_incar_unparseable_pomass(self, tmp_path): + """POMASS in original/incar/content exists but can't parse numeric values.""" + import h5py + + dummy_h5 = tmp_path / 'dummy_bad_original_incar_pomass.h5' + with h5py.File(dummy_h5, 'w') as f: + self._create_minimal_poscar(f) + f.create_dataset( + 'original/incar/content', data=np.bytes_(b'POMASS = invalid') + ) + + with pytest.raises(ValueError, match='could not parse.*numeric'): + read_cell(dummy_h5) + + def test_read_cell_potcar_no_pomass(self, tmp_path): + """POTCAR content exists but contains no POMASS entries.""" + import h5py + + dummy_h5 = tmp_path / 'dummy_potcar_no_pomass.h5' + with h5py.File(dummy_h5, 'w') as f: + self._create_minimal_poscar(f) + f.create_dataset( + 'input/potcar/content', data=np.bytes_(b'PAW_PBE Si') + ) + + with pytest.raises(ValueError, match='Could not find atomic masses'): + read_cell(dummy_h5) + + @staticmethod + def _create_minimal_poscar(h5_file): + """Create minimal input/poscar group required for read_cell.""" + pos_group = h5_file.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset('position_ions', data=np.zeros((1, 3))) + pos_group.create_dataset('number_ion_types', data=np.array([1])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + +@pytest.mark.vasp_reader +class TestVaspReaderPhononData: + + def test_read_phonon_data_missing_precalculated_raises_error(self): + with pytest.raises(MissingPhononModesError): + read_phonon_data(FC_NO_QPTS_H5) + + def test_read_phonon_data_from_only_qpts(self): + phonon_data = read_phonon_data(ONLY_QPTS_H5) + assert 'crystal' in phonon_data + assert phonon_data['crystal']['atom_r'].shape == ( + 2, + 3, + ) # Primitive cell + assert phonon_data['qpts'].shape == (3, 3) + assert phonon_data['frequencies'].shape == (3, 6) + assert phonon_data['frequencies_unit'] == 'THz' + assert phonon_data['eigenvectors'].shape == (3, 6, 2, 3) + + # Check Gamma-point optical max frequency ~ 7.6315 THz + max_freq = np.max(phonon_data['frequencies']) + assert_allclose(max_freq, 7.6315, rtol=1e-3) + + +@pytest.mark.vasp_reader +class TestQpointPhononModesFromVasp: + + def test_from_vasp_modes_missing_data_raises_error(self): + with pytest.raises(MissingPhononModesError): + QpointPhononModes.from_vasp(FC_NO_QPTS_H5) + + def test_from_vasp_modes_from_only_qpts(self): + modes = QpointPhononModes.from_vasp(ONLY_QPTS_H5) + assert modes.crystal.n_atoms == 2 # Primitive GaAs cell + assert modes.frequencies.shape == (3, 6) + assert modes.eigenvectors.shape == (3, 6, 2, 3) + + max_freq_mev = np.max(modes.frequencies.to('meV').magnitude) + assert_allclose(max_freq_mev, 31.561, rtol=1e-3) + + def test_from_vasp_frequencies_from_only_qpts(self): + freqs = QpointFrequencies.from_vasp(ONLY_QPTS_H5) + assert freqs.crystal.n_atoms == 2 + assert freqs.frequencies.shape == (3, 6) + + max_freq_mev = np.max(freqs.frequencies.to('meV').magnitude) + assert_allclose(max_freq_mev, 31.561, rtol=1e-3) + + +@pytest.mark.vasp_reader +class TestForceConstantsFromVasp: + + def test_read_interpolation_data(self): + import h5py + + data = read_interpolation_data(FC_NO_QPTS_H5) + assert 'crystal' in data + assert data['force_constants'].shape == (1, 48, 48) + assert data['sc_matrix'].shape == (3, 3) + assert 'born' in data + assert data['born'].shape == (16, 3, 3) + assert 'dielectric' in data + assert data['dielectric'].shape == (3, 3) + + # Assert raw force constants dimensions (3 * n_atoms_sc, 3 * n_atoms_sc) + with h5py.File(FC_NO_QPTS_H5, 'r') as f: + fc_raw = f['results/linear_response/force_constants'][()] + n_atoms_sc = len(f['results/positions/position_ions'][()]) + assert fc_raw.shape == (3 * n_atoms_sc, 3 * n_atoms_sc) + + def test_fc_from_vasp_and_fallback_calculation(self): + # 1. ForceConstants.from_vasp loads Hessian/force_constants + fc = ForceConstants.from_vasp(FC_NO_QPTS_H5) + assert fc.crystal.n_atoms == 16 + assert fc.n_cells_in_sc == 1 + assert fc.force_constants.shape == (1, 48, 48) + + # 2. Outer caller calculates modes explicitly from ForceConstants + q_freqs = fc.calculate_qpoint_frequencies(np.array([[0.0, 0.0, 0.0]])) + + # 3. Compare with QpointFrequencies loaded from precalculated QPOINTS file + precalc_freqs = QpointFrequencies.from_vasp(ONLY_QPTS_H5) + + fc_freqs_mev = q_freqs.frequencies.to('meV').magnitude[0] + precalc_gamma_freqs_mev = ( + precalc_freqs.frequencies.to('meV').magnitude[0] + ) + + # Max optical frequency at Gamma (31.561 meV) must match + assert_allclose( + np.max(fc_freqs_mev), np.max(precalc_gamma_freqs_mev), rtol=1e-3 + ) + + @pytest.mark.vasp_reader + def test_fc_from_vasp_without_born(self): + # Test with no Born charges/dielectric tensor + # Also a non-diagonal supercell: FCC primitive in cubic supercell + fc = ForceConstants.from_vasp(AL_NO_BORN_H5) + + # Primitive cell has 1 atom (FCC) + assert fc.crystal.n_atoms == 1 + # Supercell contains 32 primitive cells + assert fc.n_cells_in_sc == 32 + # Force constants shape: (n_cells_in_sc, 3*n_atoms, 3*n_atoms) + assert fc.force_constants.shape == (32, 3, 3) + + # No Born charges or dielectric tensor + assert fc.born is None + assert fc.dielectric is None + + @pytest.mark.vasp_reader + def test_read_interpolation_data_without_born(self): + data = read_interpolation_data(AL_NO_BORN_H5) + + assert 'crystal' in data + assert data['force_constants'].shape == (32, 3, 3) + assert data['sc_matrix'].shape == (3, 3) + # Non-diagonal supercell matrix for FCC primitive in cubic supercell + # Determinant = 32 (32 primitive cells in supercell) + assert_allclose(np.linalg.det(data['sc_matrix']), 32.0) + + # No Born charges or dielectric tensor + assert 'born' not in data + assert 'dielectric' not in data + + +@pytest.mark.vasp_reader +class TestVaspReaderCombined: + + def test_combined_fc_and_modes(self): + # 1. ForceConstants reads primitive cell force constants (8, 6, 6) + fc = ForceConstants.from_vasp(FC_AND_QPTS_H5) + assert fc.crystal.n_atoms == 2 + assert fc.n_cells_in_sc == 8 + assert fc.force_constants.shape == (8, 6, 6) + assert fc.born is not None + assert fc.born.shape == (2, 3, 3) + + # 2. QpointPhononModes reads primitive precalculated modes (3 qpts, 6 branches) + modes = QpointPhononModes.from_vasp(FC_AND_QPTS_H5) + assert modes.crystal.n_atoms == 2 + assert modes.frequencies.shape == (3, 6) + assert modes.eigenvectors.shape == (3, 6, 2, 3) + + # 3. Compare calculated frequencies from primitive FC vs precalculated modes + q_freqs = fc.calculate_qpoint_frequencies(modes.qpts) + assert_allclose( + q_freqs.frequencies.to('meV').magnitude, + modes.frequencies.to('meV').magnitude, + rtol=1e-2, + atol=1e-2, + ) + + +class TestVaspReaderEdgeCases: + """Tests for edge cases and error handling. + + Note: test_missing_h5py_import_error does NOT have @pytest.mark.vasp_reader + because it tests behavior when h5py is not installed. + """ + + def test_missing_h5py_import_error(self, mocker, tmp_path): + """Test that a helpful error is raised when h5py is not installed.""" + dummy_h5 = tmp_path / 'dummy.h5' + dummy_h5.write_text('dummy') + + real_import = builtins.__import__ + + def mocked_import(name, *args, **kwargs): + if name == 'h5py': + raise ModuleNotFoundError + return real_import(name, *args, **kwargs) + + mocker.patch('builtins.__import__', side_effect=mocked_import) + with pytest.raises(ImportVaspReaderError) as exc_info: + read_cell(dummy_h5) + assert 'Cannot import h5py' in str(exc_info.value) + + @pytest.mark.vasp_reader + def test_unexpected_eigenvector_shape_raises_error(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_bad_evec.h5' + with h5py.File(dummy_h5, 'w') as f: + pos_group = f.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset('position_ions', data=np.zeros((1, 3))) + pos_group.create_dataset('number_ion_types', data=np.array([1])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + potcar_group = f.create_group('input/potcar') + potcar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855') + ) + + ph_group = f.create_group('results/phonons') + ph_group.create_dataset('qpoint_coords', data=np.zeros((1, 3))) + ph_group.create_dataset('frequencies', data=np.zeros((1, 3))) + ph_group.create_dataset('qpoints_symmetry_weight', data=np.ones(1)) + # Bad eigenvector shape + ph_group.create_dataset('eigenvectors', data=np.zeros((1, 1, 1))) + + with pytest.raises( + ValueError, match='Unexpected eigenvector array shape' + ): + read_phonon_data(dummy_h5) + + @pytest.mark.vasp_reader + def test_read_cell_file_not_found_raises_error(self, tmp_path): + non_existent = tmp_path / 'non_existent.h5' + with pytest.raises(FileNotFoundError, match='VASP file not found'): + read_cell(non_existent) + + @pytest.mark.vasp_reader + def test_read_cell_missing_group_raises_key_error(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_nogroup.h5' + with h5py.File(dummy_h5, 'w') as f: + potcar_group = f.create_group('input/potcar') + potcar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855') + ) + + with pytest.raises(KeyError, match='Crystal position data not found'): + read_cell(dummy_h5) + + @pytest.mark.vasp_reader + def test_read_cell_falls_back_to_poscar(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_poscar.h5' + with h5py.File(dummy_h5, 'w') as f: + poscar_group = f.create_group('input/poscar') + poscar_group.create_dataset('lattice_vectors', data=np.eye(3)) + poscar_group.create_dataset('position_ions', data=np.zeros((1, 3))) + poscar_group.create_dataset('number_ion_types', data=np.array([1])) + poscar_group.create_dataset( + 'ion_types', data=np.array([b'Si']) + ) + + potcar_group = f.create_group('input/potcar') + potcar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855') + ) + + data = read_cell(dummy_h5) + assert data['cell_vectors'].shape == (3, 3) + + @pytest.mark.vasp_reader + def test_read_cell_no_positions_or_poscar_raises_key_error(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_nopos.h5' + with h5py.File(dummy_h5, 'w') as f: + potcar_group = f.create_group('input/potcar') + potcar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855') + ) + + with pytest.raises(KeyError, match='Crystal position data not found'): + read_cell(dummy_h5) + + @pytest.mark.vasp_reader + def test_find_fc_key_missing_raises_key_error(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_nofc.h5' + with h5py.File(dummy_h5, 'w') as f: + pos_group = f.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset('position_ions', data=np.zeros((1, 3))) + pos_group.create_dataset('number_ion_types', data=np.array([1])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + potcar_group = f.create_group('input/potcar') + potcar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855') + ) + + with pytest.raises(KeyError, match='Force constants not found'): + read_interpolation_data(dummy_h5) + + @pytest.mark.vasp_reader + def test_find_fc_key_hessian(self, tmp_path): + import h5py + + dummy_h5 = tmp_path / 'dummy_hessian.h5' + with h5py.File(dummy_h5, 'w') as f: + pos_group = f.create_group('input/poscar') + pos_group.create_dataset('lattice_vectors', data=np.eye(3)) + pos_group.create_dataset('position_ions', data=np.zeros((1, 3))) + pos_group.create_dataset('number_ion_types', data=np.array([1])) + pos_group.create_dataset('ion_types', data=np.array([b'Si'])) + + potcar_group = f.create_group('input/potcar') + potcar_group.create_dataset( + 'content', data=np.bytes_(b'POMASS = 28.0855') + ) + + lin_group = f.create_group('results/linear_response') + lin_group.create_dataset('hessian', data=np.zeros((3, 3))) + + data = read_interpolation_data(dummy_h5) + assert 'force_constants' in data diff --git a/tox.ini b/tox.ini index 3877b37eb..725671adb 100644 --- a/tox.ini +++ b/tox.ini @@ -40,11 +40,11 @@ commands = # Py3.15 is pre-release, test without extras [testenv:py315] -commands = {[testenv]test_command} {posargs} -m "not (phonopy_reader or matplotlib or brille)" +commands = {[testenv]test_command} {posargs} -m "not (phonopy_reader or vasp_reader or matplotlib or brille)" # Test with no extras [testenv:py310-base] -commands = {[testenv]test_command} {posargs} -m "not (phonopy_reader or matplotlib or brille)" +commands = {[testenv]test_command} {posargs} -m "not (phonopy_reader or vasp_reader or matplotlib or brille)" # Test with matplotlib extra only [testenv:py310-matplotlib] @@ -52,11 +52,11 @@ extras = matplotlib commands = {[testenv]test_command} {posargs} -m "matplotlib and not multiple_extras" -# Test with phonopy-reader extra only +# Test with hdf5/yaml readers only [testenv:py310-phonopy-reader] extras = phonopy-reader -commands = {[testenv]test_command} {posargs} -m "phonopy_reader and not multiple_extras" +commands = {[testenv]test_command} {posargs} -m "(phonopy_reader or vasp_reader) and not multiple_extras" # Test with brille extra only [testenv:py310-brille]