diff --git a/doc/ChangeLog.md b/doc/ChangeLog.md index 64ebf96..be9ffd0 100644 --- a/doc/ChangeLog.md +++ b/doc/ChangeLog.md @@ -9,6 +9,7 @@ this section to the version number when the next release is cut. ### VASP * Add `KPOINTS_OPT` band conversion from `vaspout.h5`: when `LOCPROJ_OPT` data are available, the converter writes `dft_bands_input` for band/spectral workflows, applying the same PLO config settings (`EWINDOW`, `TRANSFORM`, `NORMALIZE`, and optional `EFERMI`) as the regular VASP conversion, and stores the high-symmetry k-path labels +* Run DMFT on the irreducible Brillouin zone: `convert_dft_input(use_ibz=...)` reduces the data to the irreducible k-points, sets `symm_op = 1` and writes the correlated-shell symmetry operations (real-harmonic rotations for `l = 0, 1, 2` and atom permutations, built from the space group in `vaspout.h5`) to `dft_symmcorr_input`, so that DFTTools and ModEST restore the full-BZ average by symmetrization. Enabled by default when VASP used symmetry and `vaspout.h5` holds the symmetry data, and only when it is known to work: correlated shells with `l <= 2`, a symmetry-closed orbital set with an orthonormal `TRANSFORM`, all equivalent atoms projected, no spin-orbit coupling, and the symmetrized IBZ occupation and local levels reproduce the full grid within `ibz_tol` (1e-3). Otherwise the converter falls back to the full grid with a note; `use_ibz=True` raises instead. The VASP `Driver` passes `use_ibz` through for CSC runs ### Wien2k * Read the high-symmetry k-path labels from the end of `case.outband` and store them as `kpts_labels` / `kpts_labels_idx` in `dft_bands_input`, matching the VASP band conversion diff --git a/python/triqs_dftkit/vasp/converter.py b/python/triqs_dftkit/vasp/converter.py index a2f0b85..4cd6848 100644 --- a/python/triqs_dftkit/vasp/converter.py +++ b/python/triqs_dftkit/vasp/converter.py @@ -138,18 +138,38 @@ def read_header_and_data(self, filename): return header, f_gen, fh - def convert_dft_input(self): + def convert_dft_input(self, use_ibz=None, vasp_h5=None, ibz_tol=1e-3): """ Reads the input files, and stores the data in the HDFfile. If KPOINTS_OPT projector data is detected in vaspout.h5, the bands input is converted automatically by calling convert_bands_input(). + + Parameters + ---------- + use_ibz : bool or None, optional + Run on the irreducible Brillouin zone (IBZ) with symmetrization + instead of unfolding to the full k-grid. Requires the VASP symmetry + data in ``vasp_h5``. If ``None`` (default) the IBZ path is enabled + automatically whenever symmetry was used (n_k_ibz < n_k) and the + symmetry data is available, all correlated shells are supported + (l <= 2, symmetry-closed orbital set, orthonormal TRANSFORM, all + equivalent atoms projected), there is no spin-orbit coupling, and + the symmetrized IBZ sums reproduce the full grid (see ``ibz_tol``); + otherwise the converter falls back to the full grid with a note. + ``True`` raises instead of falling back. In text-based mode (no + vaspout.h5) a note is printed suggesting the h5 interface. + vasp_h5 : string, optional + Path to vaspout.h5 holding the k-point symmetry data. By default it + is looked up next to the basename, as for the KPOINTS_OPT data. + ibz_tol : float, optional + Largest accepted deviation of the symmetrized IBZ occupation (and + local levels, in eV) of the correlated shells from the full-grid + ones. Default 1e-3. """ energy_unit = 1.0 # VASP interface always uses eV k_dep_projection = 1 -# Symmetries are switched off for the moment -# TODO: implement symmetries - symm_op = 0 # Use symmetry groups for the k-sum + symm_op = 0 # Use symmetry groups for the k-sum (set to 1 in IBZ mode) # Read and write only on the master node if not (mpi.is_master_node()): return @@ -222,6 +242,7 @@ def convert_dft_input(self): shells = [] corr_shells = [] + corr_transforms = [] # per corr-shell (dim x 2l+1) real-harmonic -> orbital transform shion_to_shell = [[] for ish in range(len(p_shells))] cr_shion_to_shell = [[] for ish in range(len(p_shells))] shorbs_to_globalorbs = [[] for ish in range(len(p_shells))] @@ -249,6 +270,16 @@ def convert_dft_input(self): shion_to_shell[ish].append(icsh) icsh += 1 corr_shells.append(pars) + # Transform built by PLOVASP for this ion, read straight + # from the .pg header (no plo.cfg re-parse). Falls back to + # the identity for legacy files written without it. + if 'tmatrix_re' in sh: + tm = (numpy.array(sh['tmatrix_re'][i], dtype=float) + + 1j * numpy.array(sh['tmatrix_im'][i], dtype=float)) + else: + nm = 2 * pars['l'] + 1 + tm = numpy.identity(nm, dtype=complex)[:pars['dim'], :] + corr_transforms.append(tm) # TODO: generalize this to the case of multiple shell groups @@ -389,6 +420,107 @@ def convert_dft_input(self): #used for certain routines within dft_tools if treating the inputs differently is required. dft_code = 'vasp' +# ----------------------------------------------------------------------------- +# IBZ mode: reduce the data to the irreducible BZ and build the correlated-shell +# symmetry operations, so that DMFT can run on the IBZ with symmetrization. +# ----------------------------------------------------------------------------- + from . import symmetry as _sym + symm_data = None + sym_was_on = (n_k_ibz is not None) and (n_k_ibz < n_k) + if vasp_h5 is None: + vasp_h5 = next((c for c in self._vaspout_candidates() if os.path.exists(c)), 'vaspout.h5') + h5_sym = _sym.read_vasp_symmetry(vasp_h5) if os.path.exists(vasp_h5) else None + h5_available = h5_sym is not None + + auto_ibz = use_ibz is None + if auto_ibz: + use_ibz = sym_was_on and h5_available and not SO + if sym_was_on and h5_available and SO: + mpi.report(" Note: IBZ symmetrization is not implemented with spin-orbit " + "coupling; unfolding to the full k-grid.") + + if use_ibz and not h5_available: + mpi.report(" WARNING: use_ibz requested but no VASP symmetry data found in " + f"'{vasp_h5}'. Falling back to the full k-grid (use_ibz=False).") + use_ibz = False + elif not h5_available: + mpi.report(f" Note: running in text-based mode (no '{vasp_h5}'); using the full " + "k-grid. For symmetry-reduced (IBZ) runs use the vaspout.h5 interface.") + elif (not use_ibz) and sym_was_on: + mpi.report(" Note: symmetry was used in VASP but use_ibz=False; unfolding to the " + "full k-grid.") + + if use_ibz: + # In auto mode a projector set that cannot be symmetrized falls back to + # the full grid, so that setups that worked before keep working. + try: + # Guard against under-determined correlated shells. If VASP emitted + # only a subset of the 2l+1 real harmonics (e.g. a partial LOCPROJ + # selection) without a TRANSFORM defining a proper subspace, the + # absent orbitals are zero-filled in proj_mat. Their symmetry + # partners are missing, so symmetrizing on the IBZ would leak the + # kept orbitals into the dead ones (the invariance guard cannot see + # this: an identity transform is trivially invariant). Detect and + # refuse before building the symmetry operations. + for icrsh, csh in enumerate(corr_shells): + dim = csh['dim'] + orb_weight = numpy.max(numpy.abs(proj_mat[:, :, icrsh, :dim, :]), + axis=(0, 1, 3)) # (dim,) + scale = orb_weight.max() + dead = numpy.where(orb_weight <= max(1e-10, 1e-8 * scale))[0] + if scale > 0 and len(dead) > 0: + raise RuntimeError( + f"IBZ symmetrization: correlated shell {icrsh} (atom " + f"{csh.get('atom')}, l={csh['l']}, dim={dim}) has zero-weight " + f"projector orbital(s) {list(map(int, dead))}. This happens when " + "only a subset of the 2l+1 real harmonics was projected (e.g. a " + "partial LOCPROJ selection) without a TRANSFORM defining a " + "symmetry-adapted subspace. The missing symmetry partners cannot " + "be restored, so the k-sum cannot be symmetrized on the " + "irreducible BZ. Project the full shell (or a complete irrep via " + "TRANSFORM), or run on the full grid (use_ibz=False).") + + symm_data = _sym.build_symmcorr(vasp_h5, corr_shells, corr_transforms, SP, SO, sym=h5_sym) + + # The symmetrized IBZ sums must reproduce the full grid, which + # is still at hand here, for the occupation and the local levels. + err_occ, err_hloc = _sym.symmetrization_error(symm_data, proj_mat, hopping, f_weights, + bz_weights, n_orbitals, corr_shells) + mpi.report(f" IBZ check: max deviation from the full grid {err_occ:.1e} " + f"(occupation), {err_hloc:.1e} eV (local levels).") + if max(err_occ, err_hloc) > ibz_tol: + raise RuntimeError( + f"IBZ symmetrization: the symmetrized IBZ sums deviate from the full " + f"grid by more than ibz_tol={ibz_tol:.0e}. Run on the full grid " + "(use_ibz=False).") + except (RuntimeError, NotImplementedError) as err: + if not auto_ibz: raise + mpi.report(f" Note: {err}\n Unfolding to the full k-grid instead.") + use_ibz, symm_data = False, None + + if use_ibz: + nkibz = int(symm_data['n_k_ibz']) + mpi.report(f" IBZ mode: reducing {n_k} k-points to {nkibz} irreducible " + f"k-points with {symm_data['n_symm']} symmetry operations.") + + # The first nkibz points of the (plotools-ordered) full grid are the IBZ + # representatives, so slicing recovers the irreducible data. + sl = slice(0, nkibz) + hopping = hopping[sl] + f_weights = f_weights[sl] + n_orbitals = n_orbitals[sl] + proj_mat = proj_mat[sl] + kpts = kpts[sl] + kpts_cart = kpts_cart[sl] + band_window = [bw[sl] for bw in band_window] + if self.proj_or_hk == 'hk' or self.proj_or_hk == True: + proj_mat_csc = proj_mat_csc[sl] + n_k = nkibz + self.n_k = nkibz + bz_weights = numpy.asarray(symm_data['ibz_weights']).copy() + kpt_weights = bz_weights.copy() + symm_op = 1 + # Save it to the HDF: with HDFArchive(self.hdf_file,'a') as ar: if not (self.dft_subgrp in ar): ar.create_group(self.dft_subgrp) @@ -410,16 +542,15 @@ def convert_dft_input(self): if n_k_ibz is not None: ar[self.misc_subgrp]['n_k_ibz'] = n_k_ibz - # Symmetries are used, so now convert symmetry information for *correlated* orbitals: - self.convert_symmetry_input(ctrl_head, orbits=self.corr_shells, symm_subgrp=self.symmcorr_subgrp) + # Symmetry information for *correlated* orbitals: + if use_ibz and symm_data is not None: + self.write_symmetry_input(symm_data, symm_subgrp=self.symmcorr_subgrp) + else: + self.convert_symmetry_input(ctrl_head, orbits=self.corr_shells, symm_subgrp=self.symmcorr_subgrp) # Auto-convert KPOINTS_OPT band/projector data when available. - vaspout_candidates = [ - os.path.join(self.basename, 'vaspout.h5'), - os.path.join(os.path.dirname(self.basename), 'vaspout.h5') - ] kpoints_opt_found = False - for candidate in vaspout_candidates: + for candidate in self._vaspout_candidates(): if not os.path.exists(candidate): continue try: @@ -688,10 +819,7 @@ class _BandsElStruct: n_spin_blocs = 1 if int(self.SO) == 1 else int(self.SP) + 1 # Read KPOINTS_OPT data directly from vaspout.h5 - vaspout_candidates = [ - os.path.join(self.basename, 'vaspout.h5'), - os.path.join(os.path.dirname(self.basename), 'vaspout.h5') - ] + vaspout_candidates = self._vaspout_candidates() vaspout_h5 = None for candidate in vaspout_candidates: if os.path.exists(candidate): @@ -870,6 +998,31 @@ def convert_misc_input(self, bandwin_file, struct_file, outputs_file, misc_subgr for it in things_to_save: ar[misc_subgrp][it] = locals()[it] + def _vaspout_candidates(self): + """Locations where vaspout.h5 is looked up, relative to the basename.""" + return [os.path.join(self.basename, 'vaspout.h5'), + os.path.join(os.path.dirname(self.basename), 'vaspout.h5')] + + def write_symmetry_input(self, symm_data, symm_subgrp): + """ + Write the correlated-shell symmetry operations (built from the VASP + symmetry data) to the symmetry subgroup of the HDF5 archive. + + Parameters + ---------- + symm_data : dict + Output of ``triqs_dftkit.vasp.symmetry.build_symmcorr`` holding + n_symm, n_atoms, perm, orbits, SO, SP, time_inv, mat, mat_tinv. + symm_subgrp : string + Name of the symmetry group in the HDF5 archive. + """ + with HDFArchive(self.hdf_file, 'a') as ar: + if not (symm_subgrp in ar): ar.create_group(symm_subgrp) + things_to_save = ['n_symm', 'n_atoms', 'perm', 'orbits', 'SO', 'SP', + 'time_inv', 'mat', 'mat_tinv'] + for it in things_to_save: + ar[symm_subgrp][it] = symm_data[it] + def convert_symmetry_input(self, ctrl_head, orbits, symm_subgrp): """ Reads input for the symmetrisations from symm_file, which is case.sympar or case.symqmc. diff --git a/python/triqs_dftkit/vasp/driver.py b/python/triqs_dftkit/vasp/driver.py index fa5517f..0a3a369 100644 --- a/python/triqs_dftkit/vasp/driver.py +++ b/python/triqs_dftkit/vasp/driver.py @@ -87,18 +87,24 @@ class Driver: VASP executable name (e.g., "vasp_std"). mpi_handler : MPIHandler MPI configuration handler. + use_ibz : bool or None + Passed to Converter.convert_dft_input(): run on the irreducible BZ + with symmetrization. None (default) enables it whenever VASP used + symmetry and vaspout.h5 holds the symmetry data. """ def __init__(self, seedname: str, plo_cfg: str, mpi_handler: MPIHandler, - vasp_command: str = "vasp_std") -> None: + vasp_command: str = "vasp_std", use_ibz: bool | None = None) -> None: self.seedname = seedname self.plo_cfg = plo_cfg self.vasp_command = vasp_command self.mpi_handler = mpi_handler + self.use_ibz = use_ibz self._vasp_process_id = None def __repr__(self): - return f"VaspDriver(seedname={self.seedname}, plo_cfg={self.plo_cfg}, vasp_command={self.vasp_command})" + return (f"VaspDriver(seedname={self.seedname}, plo_cfg={self.plo_cfg}, " + f"vasp_command={self.vasp_command}, use_ibz={self.use_ibz})") __str__ = __repr__ @@ -194,7 +200,7 @@ def _run_plo_converter(self): generate_and_output_as_text(self.plo_cfg, vasp_dir='./') self.mpi_handler.report(f"[{datetime.now()}] Running VASP HDF5 converter") - Converter(filename=self.seedname).convert_dft_input() + Converter(filename=self.seedname).convert_dft_input(use_ibz=self.use_ibz) mpi.barrier(poll_msec=100) diff --git a/python/triqs_dftkit/vasp/plovasp/plotools.py b/python/triqs_dftkit/vasp/plovasp/plotools.py index a7b077f..fa1a4df 100644 --- a/python/triqs_dftkit/vasp/plovasp/plotools.py +++ b/python/triqs_dftkit/vasp/plovasp/plotools.py @@ -418,7 +418,12 @@ def plo_output(conf_pars, el_struct, pshells, pgroups): sh_dict['ion_list'] = ion_output sh_dict['ion_sort'] = shell.ion_sort -# TODO: add the output of transformation matrices +# Per-ion transformation matrices (real harmonics -> correlated orbitals), +# shape (nion, ndim, 2l+1). Stored split into real/imag because JSON has no +# complex type. Consumed by the converter to build the correlated-shell +# symmetry operations for IBZ runs (no need to re-parse the PLO config). + sh_dict['tmatrix_re'] = shell.tmatrices.real.tolist() + sh_dict['tmatrix_im'] = shell.tmatrices.imag.tolist() head_shells.append(sh_dict) diff --git a/python/triqs_dftkit/vasp/symmetry.py b/python/triqs_dftkit/vasp/symmetry.py new file mode 100644 index 0000000..54a3616 --- /dev/null +++ b/python/triqs_dftkit/vasp/symmetry.py @@ -0,0 +1,480 @@ +################################################################################ +# +# TRIQS / dftkit +# +# Construction of correlated-shell symmetry operations for VASP, enabling DMFT +# to run on the irreducible Brillouin zone (IBZ) with symmetrization instead of +# unfolding to the full k-grid. +# +# VASP exposes, in vaspout.h5, the space-group operations that relate every +# full-grid k-point to its IBZ representative (as integer matrices in the basis +# of the reciprocal lattice). From these we build the rotation matrices acting +# on the correlated local orbitals (the `mat` entries of `dft_symmcorr_input`) +# that triqs_modest / dft_tools use to restore the full-BZ average of any +# k-summed local quantity. +# +################################################################################ +""" +VASP symmetry tools for IBZ-based DMFT. +""" +import numpy as np + + +# --------------------------------------------------------------------------- +# Real cubic-harmonic rotation matrices D^l(R) +# +# The local orbitals are real (cubic) harmonics in the VASP/PLO ordering. Under +# a Cartesian point-group operation R (proper or improper, det = +-1) they +# transform with a (2l+1)x(2l+1) orthogonal matrix D^l(R), built here without +# any external dependency: +# +# l = 1 : real harmonics are the Cartesian components; D^1 is R itself, +# re-ordered to the VASP basis [y, z, x]. +# l = 2 : each real d harmonic maps to a symmetric traceless 3x3 tensor T_i; +# under R the tensor transforms as T -> R T R^T, and +# D^2[i, j] = (Frobenius product on the orthonormal +# tensor basis). Correct for improper R (d harmonics are parity-even). +# +# Basis orderings match modest's dft_tools/vasp/spherical_rotation.cpp: +# l = 1 : "y", "z", "x" (m = -1, 0, +1) +# l = 2 : "xy","yz","z2","xz","x2-y2" +# --------------------------------------------------------------------------- + +# l = 1: real-harmonic basis vectors in VASP order [y, z, x] +_P_BASIS = [np.array([0.0, 1.0, 0.0]), # y + np.array([0.0, 0.0, 1.0]), # z + np.array([1.0, 0.0, 0.0])] # x + + +def _d_tensor_basis(): + """Symmetric traceless 3x3 tensors for the 5 real d harmonics (VASP order), + Frobenius-orthonormal: = delta_ij.""" + s2 = np.sqrt(2.0) + T = [] + m = np.zeros((3, 3)); m[0, 1] = m[1, 0] = 1 / s2; T.append(m) # xy + m = np.zeros((3, 3)); m[1, 2] = m[2, 1] = 1 / s2; T.append(m) # yz + T.append(np.diag([-1.0, -1.0, 2.0]) / np.sqrt(6.0)) # z^2 + m = np.zeros((3, 3)); m[0, 2] = m[2, 0] = 1 / s2; T.append(m) # xz + T.append(np.diag([1.0, -1.0, 0.0]) / s2) # x^2-y^2 + return T + + +_D_BASIS = _d_tensor_basis() + +# Angular momenta for which real_harmonic_rotation is implemented +SUPPORTED_L = (0, 1, 2) + + +def real_harmonic_rotation(R, l): + """ + Rotation matrix of the real (cubic) harmonics of angular momentum ``l`` + under the Cartesian orthogonal operation ``R`` (det = +-1). + + Parameters + ---------- + R : (3, 3) array + Cartesian point-group operation (proper or improper). + l : int + Angular momentum (0, 1 and 2 are implemented). + + Returns + ------- + (2l+1, 2l+1) ndarray + Orthogonal rotation matrix in the VASP real-harmonic basis. + """ + if l == 0: + return np.ones((1, 1)) + if l == 1: + D = np.zeros((3, 3)) + for j, bj in enumerate(_P_BASIS): + Rbj = R @ bj + for i, bi in enumerate(_P_BASIS): + D[i, j] = bi @ Rbj + return D + if l == 2: + D = np.zeros((5, 5)) + for j, Tj in enumerate(_D_BASIS): + RTjRt = R @ Tj @ R.T + for i, Ti in enumerate(_D_BASIS): + D[i, j] = np.tensordot(Ti, RTjRt) + return D + raise NotImplementedError( + f"real_harmonic_rotation is only implemented for l = 0, 1, 2 (got l = {l}). " + "Higher l (f electrons) needs the Wigner-D / rank-l tensor extension.") + + +def _subspace_invariance(T, R_cart, l): + """ + Maximum 'leakage' of the correlated subspace out of itself under the real- + harmonic representation D(R). + + The correlated orbitals span the row space of ``T`` (shape ``(dim, 2l+1)``) + inside the (2l+1)-dimensional real-harmonic space. The returned value is + ``max_R || (1 - P) D(R) P ||_2`` with ``P`` the orthogonal projector onto + that subspace. It is zero iff the subspace is invariant under every ``R`` + (the condition for ``Q = T D T^dag`` to be a group representation). The + measure is gauge-free: it depends only on the span of ``T``, not on its + particular rows, so a non-orthonormal or rotated transform gives the same + answer. + """ + nm = 2 * l + 1 + T = np.asarray(T, dtype=complex) + # Orthonormal basis of the row span via SVD (rows of Vh with non-zero sv). + _, s, Vh = np.linalg.svd(T, full_matrices=False) + Q = Vh[s > 1e-10] # (r, nm), orthonormal rows + P = Q.conj().T @ Q # (nm, nm) Hermitian projector + eye = np.eye(nm) + viol = 0.0 + for R in R_cart: + D = real_harmonic_rotation(R, l).astype(complex) + viol = max(viol, np.linalg.norm((eye - P) @ D @ P, ord=2)) + return viol + + +# --------------------------------------------------------------------------- +# Reading the VASP symmetry data from vaspout.h5 +# --------------------------------------------------------------------------- +def _close_group(ops, max_size=200): + """ + Close a set of integer matrices under multiplication to recover the full + group. ``ops`` is an (n, 3, 3) integer array; returns an (m, 3, 3) array + with the complete group. Exact (integer arithmetic). + """ + group = [np.asarray(o, dtype=int) for o in ops] + seen = {o.tobytes() for o in group} + changed = True + while changed: + changed = False + for a in list(group): + for b in list(group): + c = a @ b + key = c.tobytes() + if key not in seen: + seen.add(key) + group.append(c) + changed = True + if len(group) > max_size: + raise RuntimeError("symmetry group closure exceeded %d operations; " + "stored symmetry data may be inconsistent" % max_size) + return np.array(group) + + +def read_vasp_symmetry(vasp_h5): + """ + Read the symmetry-relevant data from a vaspout.h5 file. + + Returns a dict with: + n_k_ibz : number of irreducible k-points + ibz_weights : (n_k_ibz,) IBZ k-point weights (sum to 1) + ksymmap : (n_k,) index of the IBZ representative of each k-point + R_cart : list of unique Cartesian point-group operations (3x3) + positions : (natom, 3) fractional ion positions + lattice : (3, 3) Cartesian lattice vectors (rows) + type_of_ion : (natom,) integer ion-type index + or ``None`` if the file does not contain k-point symmetry data (no symmetry + was used / not a symmetry-reduced run). + """ + from h5 import HDFArchive + with HDFArchive(vasp_h5, 'r') as ar: + if 'results' not in ar or 'electron_eigenvalues' not in ar['results']: + return None + ee = ar['results']['electron_eigenvalues'] + if 'kpoints_symmetry_symop' not in ee: + return None + symop = np.asarray(ee['kpoints_symmetry_symop']) # (nktot, 3, 3) int, reciprocal-fractional + ibz_weights = np.asarray(ee['kpoints_symmetry_weight']) # (nkibz,) sum = 1 + n_k_ibz = int(ee['kpoints']) if 'kpoints' in ee else len(ibz_weights) + ksymmap = np.asarray(ee['kpoints_symmetry_mapping']) - 1 # full grid -> IBZ index + positions_grp = ar['results']['positions'] + lattice = np.asarray(positions_grp['lattice_vectors']) # rows = Cartesian lattice vectors + positions = np.asarray(positions_grp['position_ions']) # (natom, 3) fractional + number_ion_types = np.asarray(positions_grp['number_ion_types']) + + # type index per ion + type_of_ion = np.concatenate([[it] * n for it, n in enumerate(number_ion_types)]) + + # reciprocal basis (columns = reciprocal vectors): B = 2pi (A^-1)^T, A rows = lattice vectors + B = (2 * np.pi * np.linalg.inv(lattice).T).T + Binv = np.linalg.inv(B) + # VASP stores, per full-grid k-point, the single operation mapping the IBZ + # representative to that point (a coset representative). The set of stored + # operations is therefore complete only if some k-point has a full-size star. + # On coarse / special meshes it is a proper subset, so we close it under + # multiplication to recover the full point group (exact: integer matrices). + uniq = _close_group(np.unique(symop.reshape(-1, 9), axis=0).reshape(-1, 3, 3)) + R_cart = [B @ s.astype(float) @ Binv for s in uniq] + + return dict(n_k_ibz=n_k_ibz, ibz_weights=ibz_weights, ksymmap=ksymmap, R_cart=R_cart, + positions=positions, lattice=lattice, type_of_ion=type_of_ion) + + +# --------------------------------------------------------------------------- +# Atom permutations under the space-group operations +# --------------------------------------------------------------------------- +def _frac_mod(x): + """Fold fractional coordinates into [0, 1).""" + return x - np.floor(x + 1e-8) + + +def atom_permutation(R_cart, positions, lattice, type_of_ion, tol=1e-4): + """ + Permutation of the crystal atoms under the Cartesian operation ``R_cart``. + + A global fractional translation ``t`` (space-group translation) is detected + so the routine works for symmorphic operations regardless of cell origin. + Returns a 1-indexed permutation array ``perm`` such that atom ``j`` maps to + atom ``perm[j]`` (1-indexed), or ``None`` if no consistent permutation is + found (e.g. a genuinely non-symmorphic operation needing its own glide/screw + translation, not handled here). + """ + natom = len(positions) + pos_cart = positions @ lattice # (natom, 3) Cartesian + rot_cart = pos_cart @ R_cart.T # apply R to each position + Linv = np.linalg.inv(lattice) + rot_frac = rot_cart @ Linv # back to fractional + + # candidate global translations from mapping atom 0 -> any same-type atom + for i0 in range(natom): + if type_of_ion[i0] != type_of_ion[0]: + continue + t = positions[i0] - rot_frac[0] + perm = np.full(natom, -1, dtype=int) + ok = True + for j in range(natom): + img = _frac_mod(rot_frac[j] + t) + # find matching atom of the same type + match = -1 + for i in range(natom): + if type_of_ion[i] != type_of_ion[j]: + continue + d = _frac_mod(positions[i] - img) + d = np.minimum(d, 1.0 - d) + if np.linalg.norm(d) < tol: + match = i + break + if match < 0: + ok = False + break + perm[j] = match + 1 # 1-indexed + if ok and len(set(perm)) == natom: + return perm + return None + + +# --------------------------------------------------------------------------- +# Top-level: build the dft_symmcorr_input payload +# --------------------------------------------------------------------------- +def build_symmcorr(vasp_h5, corr_shells, transforms, SP, SO, sym=None): + """ + Build the ``dft_symmcorr_input`` data for the correlated shells. + + Parameters + ---------- + vasp_h5 : str + Path to vaspout.h5. + corr_shells : list of dict + Correlated shells (as built by the converter); each has 'atom', 'l', + 'dim'. + transforms : list of (dim, 2l+1) arrays + Per-correlated-shell transformation matrix T mapping the real harmonics + onto the correlated orbitals (the PLO TRANSFORM). One per corr shell. + SP, SO : int + Spin-polarization and spin-orbit flags. Spin-orbit coupling (SO = 1) is + not supported: the operations would have to act on the spinors too. + sym : dict, optional + Output of ``read_vasp_symmetry``, to avoid reading ``vasp_h5`` again. + + Returns + ------- + dict with keys n_symm, n_atoms, perm, orbits, SO, SP, time_inv, mat, mat_tinv + and additionally n_k_ibz, ibz_weights for the IBZ slicing. + Returns ``None`` if no symmetry data is available. + """ + if SO: + raise NotImplementedError( + "IBZ symmetrization is not implemented with spin-orbit coupling (SO = 1): the " + "symmetry operations do not act on the spin. Run on the full grid (use_ibz=False).") + + if sym is None: + sym = read_vasp_symmetry(vasp_h5) + if sym is None: + return None + + # The converter recovers the IBZ data by slicing the first n_k_ibz points + # of the full grid, which requires them to be the IBZ representatives. + n_k_ibz = sym['n_k_ibz'] + if not np.array_equal(sym['ksymmap'][:n_k_ibz], np.arange(n_k_ibz)): + raise RuntimeError("IBZ symmetrization: the first n_k_ibz points of the full k-grid in " + "vaspout.h5 are not the irreducible k-points.") + + R_cart = sym['R_cart'] + n_symm = len(R_cart) + n_corr = len(corr_shells) + + # Only shells with known real-harmonic rotation matrices, and transforms + # with orthonormal rows (otherwise T D T^dag is not unitary). + for ish, csh in enumerate(corr_shells): + if csh['l'] not in SUPPORTED_L: + raise NotImplementedError( + f"IBZ symmetrization: correlated shell {ish} (atom {csh.get('atom')}) has " + f"l={csh['l']}; the symmetry operations are only implemented for l = " + f"{', '.join(map(str, SUPPORTED_L))}. Run on the full grid (use_ibz=False).") + T = np.asarray(transforms[ish], dtype=complex) + if np.abs(T @ T.conj().T - np.eye(T.shape[0])).max() > 1e-8: + raise RuntimeError( + f"IBZ symmetrization: the TRANSFORM of correlated shell {ish} (atom " + f"{csh.get('atom')}) does not have orthonormal rows. Run on the full grid " + "(use_ibz=False).") + + # The projected operations Q = T D(R) T^dag form a representation of the + # point group on the correlated orbitals only if the orbital subspace (the + # row span of T) is invariant under every D(R). If it is not (e.g. an + # arbitrary LOCPROJ subset that splits a degenerate multiplet such as t2g), + # the symmetrization of any k-summed quantity on the IBZ is invalid: the + # missing partners would be needed to close the star. Refuse rather than + # silently produce a wrong full-BZ average. + for ish, csh in enumerate(corr_shells): + viol = _subspace_invariance(transforms[ish], R_cart, csh['l']) + if viol > 1e-6: + raise RuntimeError( + f"IBZ symmetrization: the correlated subspace of shell {ish} " + f"(atom {csh.get('atom')}, l={csh['l']}, dim={csh['dim']}) is not " + f"invariant under the crystal point group (leakage {viol:.2e}). The " + "chosen orbital set is not symmetry-closed, so k-summed quantities " + "cannot be symmetrized on the irreducible BZ. Use a symmetry-adapted " + "orbital set (the full 2l+1 shell, or a complete irrep such as t2g or " + "eg in cubic symmetry), or run on the full grid (use_ibz=False).") + + # symmetrization matrices: mat[g] = list over corr shells of (dim x dim) + mat = [] + for R in R_cart: + per_shell = [] + for ish, csh in enumerate(corr_shells): + l = csh['l'] + T = np.asarray(transforms[ish], dtype=complex) # (dim, 2l+1) + D = real_harmonic_rotation(R, l).astype(complex) # (2l+1, 2l+1) + Q = T @ D @ T.conj().T # (dim, dim) in the correlated basis + per_shell.append(Q) + mat.append(per_shell) + + # atom permutations (over all crystal atoms, 1-indexed) + natom = len(sym['positions']) + perm = [] + for R in R_cart: + p = atom_permutation(R, sym['positions'], sym['lattice'], sym['type_of_ion']) + if p is None: + raise RuntimeError( + "IBZ symmetrization: could not determine the atom permutation of a VASP " + "symmetry operation. Run on the full grid (use_ibz=False).") + perm.append(list(map(int, p))) + + # Every symmetry image of a correlated shell must itself be a correlated + # shell (same l and dim), otherwise the star of the local quantity cannot be + # closed, e.g. when only one of two equivalent atoms is treated as correlated. + shell_key = {(csh['atom'], csh['l'], csh['dim']): ish for ish, csh in enumerate(corr_shells)} + for p in perm: + for csh in corr_shells: + if (p[csh['atom'] - 1], csh['l'], csh['dim']) not in shell_key: + raise RuntimeError( + f"IBZ symmetrization: atom {csh['atom']} (l={csh['l']}) is mapped onto " + f"atom {p[csh['atom'] - 1]} by a symmetry operation, but the latter is not " + "a correlated shell with the same l and dim. Project all symmetry-equivalent " + "atoms, or run on the full grid (use_ibz=False).") + + # dft_tools' Symmetry finds the image of an orbit by comparing the whole + # orbit dict, 'sort' included, while the VASP converter gives every site its + # own sort. Shells mapped onto each other must therefore share one sort. + # Union-find over the shell images; each class takes its smallest sort. + parent = list(range(n_corr)) + def find(i): + while parent[i] != i: + i = parent[i] + return i + for p in perm: + for ish, csh in enumerate(corr_shells): + jsh = shell_key[(p[csh['atom'] - 1], csh['l'], csh['dim'])] + ri, rj = find(ish), find(jsh) + if ri != rj: + parent[max(ri, rj)] = min(ri, rj) + orbits = [dict(csh) for csh in corr_shells] + for ish, orb in enumerate(orbits): + orb['sort'] = corr_shells[find(ish)]['sort'] + + time_inv = [0] * n_symm + mat_tinv = [np.identity(csh['dim'], complex) for csh in corr_shells] + + return dict(n_symm=n_symm, n_atoms=natom, perm=perm, orbits=orbits, + SO=SO, SP=SP, time_inv=time_inv, mat=mat, mat_tinv=mat_tinv, + n_k_ibz=sym['n_k_ibz'], ibz_weights=sym['ibz_weights']) + + +# --------------------------------------------------------------------------- +# Check of the symmetry operations against the full-grid data +# --------------------------------------------------------------------------- +def symmetrize_local(symm_data, local): + """ + Symmetrize per-shell matrices ``local[icrsh][spin]`` (summed on the IBZ) + with the operations of ``symm_data``, following + triqs_dft_tools.symmetry.Symmetry, including how the image of a shell is + found (comparison of the whole orbit dict). + """ + orbits = symm_data['orbits'] + n_symm = symm_data['n_symm'] + out = [np.zeros_like(x) for x in local] + for isym in range(n_symm): + for iorb, orb in enumerate(orbits): + jorb = orbits.index(dict(orb, atom=symm_data['perm'][isym][orb['atom'] - 1])) + mat = symm_data['mat'][isym][iorb] + out[jorb] += mat @ local[iorb] @ mat.conj().T / n_symm + return out + + +def local_sums(proj_mat, hopping, f_weights, weights, n_orbitals, corr_shells, n_k): + """ + Occupation sum_k w_k P f P^dag and local Hamiltonian sum_k w_k P H P^dag of + every correlated shell over the first ``n_k`` k-points, shape + (n_spin, dim, dim) per shell. + """ + n_spin = proj_mat.shape[1] + occ, hloc = [], [] + for icrsh, csh in enumerate(corr_shells): + dim = csh['dim'] + o = np.zeros((n_spin, dim, dim), dtype=complex) + h = np.zeros((n_spin, dim, dim), dtype=complex) + for ik in range(n_k): + for isp in range(n_spin): + nb = n_orbitals[ik, isp] + P = proj_mat[ik, isp, icrsh, :dim, :nb] + o[isp] += weights[ik] * (P * f_weights[ik, isp, :nb]) @ P.conj().T + h[isp] += weights[ik] * P @ hopping[ik, isp, :nb, :nb] @ P.conj().T + occ.append(o) + hloc.append(h) + return occ, hloc + + +def symmetrization_error(symm_data, proj_mat, hopping, f_weights, bz_weights, n_orbitals, corr_shells): + """ + Largest deviation of the symmetrized IBZ sums from the full-grid sums of the + occupation and of the local Hamiltonian, for full-grid input arrays. + + This checks the symmetry operations on the data themselves, independent of + how they were constructed: wrong rotations, permutations or weights show up + here. What remains for correct operations is the symmetry error of the VASP + projectors (~1e-4 for the NiO test data). + + Returns + ------- + (float, float) + Maximum absolute deviation of the occupation matrices and of the local + Hamiltonian (eV). + """ + n_k_ibz = int(symm_data['n_k_ibz']) + occ_full, hloc_full = local_sums(proj_mat, hopping, f_weights, bz_weights, n_orbitals, + corr_shells, len(bz_weights)) + occ_ibz, hloc_ibz = local_sums(proj_mat, hopping, f_weights, symm_data['ibz_weights'], + n_orbitals, corr_shells, n_k_ibz) + occ_ibz = symmetrize_local(symm_data, occ_ibz) + hloc_ibz = symmetrize_local(symm_data, hloc_ibz) + err_occ = max(np.abs(a - b).max() for a, b in zip(occ_ibz, occ_full)) + err_hloc = max(np.abs(a - b).max() for a, b in zip(hloc_ibz, hloc_full)) + return err_occ, err_hloc diff --git a/test/python/vasp/converter/nio_afm.cfg b/test/python/vasp/converter/nio_afm.cfg new file mode 100644 index 0000000..e488552 --- /dev/null +++ b/test/python/vasp/converter/nio_afm.cfg @@ -0,0 +1,27 @@ +[General] +BASENAME = converter/nio_afm + +[Group 1] +SHELLS = 1 2 +EWINDOW = -10 10 +BANDS = 9 28 +NORMION = False +NORMALIZE = True + +[Shell 1] # Ni d shell (2 equivalent Ni, full d, l=2) +LSHELL = 2 +IONS = 1 2 +CORR = True +TRANSFORM = 1.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 + 0.0 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 0.0 1.0 + +[Shell 2] # O p shell (2 O, full p, l=1) +LSHELL = 1 +IONS = 3 4 +CORR = True +TRANSFORM = 1.0 0.0 0.0 + 0.0 1.0 0.0 + 0.0 0.0 1.0 diff --git a/test/python/vasp/converter/nio_afm/vaspout.h5 b/test/python/vasp/converter/nio_afm/vaspout.h5 new file mode 100644 index 0000000..ada2ea1 Binary files /dev/null and b/test/python/vasp/converter/nio_afm/vaspout.h5 differ diff --git a/test/python/vasp/converter/nio_afm_ibz.ref.h5 b/test/python/vasp/converter/nio_afm_ibz.ref.h5 new file mode 100644 index 0000000..6a4fe3b Binary files /dev/null and b/test/python/vasp/converter/nio_afm_ibz.ref.h5 differ diff --git a/test/python/vasp/converter/nio_afm_one_o.cfg b/test/python/vasp/converter/nio_afm_one_o.cfg new file mode 100644 index 0000000..a3439e2 --- /dev/null +++ b/test/python/vasp/converter/nio_afm_one_o.cfg @@ -0,0 +1,27 @@ +[General] +BASENAME = converter/nio_afm_one_o + +[Group 1] +SHELLS = 1 2 +EWINDOW = -10 10 +BANDS = 9 28 +NORMION = False +NORMALIZE = True + +[Shell 1] # Ni d shell (2 equivalent Ni, full d, l=2) +LSHELL = 2 +IONS = 1 2 +CORR = True +TRANSFORM = 1.0 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 0.0 + 0.0 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 0.0 1.0 + +[Shell 2] # O p shell (only one of the two equivalent O) +LSHELL = 1 +IONS = 3 +CORR = True +TRANSFORM = 1.0 0.0 0.0 + 0.0 1.0 0.0 + 0.0 0.0 1.0 diff --git a/test/python/vasp/converter/svo_ibz.ref.h5 b/test/python/vasp/converter/svo_ibz.ref.h5 new file mode 100644 index 0000000..196d421 Binary files /dev/null and b/test/python/vasp/converter/svo_ibz.ref.h5 differ diff --git a/test/python/vasp/converter/test_converter_ibz.py b/test/python/vasp/converter/test_converter_ibz.py new file mode 100644 index 0000000..f14cc56 --- /dev/null +++ b/test/python/vasp/converter/test_converter_ibz.py @@ -0,0 +1,173 @@ + +import os +import unittest.mock +import numpy as np +import rpath +_rpath = os.path.dirname(rpath.__file__) + '/' + +from h5 import HDFArchive +from triqs_dftkit.vasp.plovasp.converter import generate_and_output_as_text +from triqs_dftkit.vasp import Converter +from triqs_dftkit.vasp import symmetry +import mytest + + +def local_occupations(h5_file): + """ + DFT occupation matrix of every correlated shell, sum_k w_k P_k f_k P_k^dag. + """ + with HDFArchive(h5_file, 'r') as ar: + dft = ar['dft_input'] + proj_mat, bz_weights, n_orbitals = dft['proj_mat'], dft['bz_weights'], dft['n_orbitals'] + corr_shells = dft['corr_shells'] + f_weights = ar['dft_misc_input']['dft_fermi_weights'] + + n_k, n_spin = proj_mat.shape[:2] + occ = [] + for icrsh, csh in enumerate(corr_shells): + o = np.zeros((n_spin, csh['dim'], csh['dim']), dtype=complex) + for ik in range(n_k): + for isp in range(n_spin): + nb = n_orbitals[ik, isp] + P = proj_mat[ik, isp, icrsh, :csh['dim'], :nb] + o[isp] += bz_weights[ik] * (P * f_weights[ik, isp, :nb]) @ P.conj().T + occ.append(o) + return occ + + +def symmetrize(occ, h5_file): + """ + Symmetrize per-shell matrices with the dft_symmcorr_input of ``h5_file``, + following triqs_dft_tools.symmetry.Symmetry, including how the image of a + shell is found (comparison of the whole orbit dict, 'sort' included), so the + archive is checked in the form a DMFT code consumes it. + """ + with HDFArchive(h5_file, 'r') as ar: + symm = ar['dft_symmcorr_input'] + orbits = symm['orbits'] + symm_occ = [np.zeros_like(o) for o in occ] + for isym in range(symm['n_symm']): + for iorb, orb in enumerate(orbits): + jorb = orbits.index(dict(orb, atom=symm['perm'][isym][orb['atom'] - 1])) + mat = symm['mat'][isym][iorb] + symm_occ[jorb] += mat @ occ[iorb] @ mat.conj().T / symm['n_symm'] + return symm_occ + + +################################################################################ +# +# TestConverterIBZ +# +# Converts VASP output to the irreducible-BZ representation with symmetrization +# (use_ibz=True): the data are reduced to the irreducible k-points and the +# correlated-shell symmetry operations are written to dft_symmcorr_input so that +# DMFT can run on the IBZ instead of the unfolded full grid. +# +# Besides the comparison with the reference archive, the symmetrized IBZ +# occupation matrices must reproduce the symmetrized full-grid ones, which only +# holds if the weights, the orbital rotations and the atom permutations are all +# correct. The unsymmetrized full grid agrees only to the symmetry of the VASP +# projectors themselves, hence the looser tolerance for that comparison. +# +################################################################################ +class TestConverterIBZ(mytest.MyTestCase): + + def convert(self, name, use_ibz): + generate_and_output_as_text(_rpath + name + '.cfg', _rpath + name + '/') + test_file = _rpath + name + ('_ibz' if use_ibz else '_full') + '.test.h5' + if os.path.exists(test_file): os.remove(test_file) + converter = Converter(filename=_rpath + name, hdf_filename=test_file) + converter.convert_dft_input(use_ibz=use_ibz, vasp_h5=_rpath + name + '/vaspout.h5') + return test_file + + def check_ibz(self, name, tol): + full_file = self.convert(name, use_ibz=False) + ibz_file = self.convert(name, use_ibz=True) + + self.assertH5FileEqual(ibz_file, _rpath + name + '_ibz.ref.h5') + + occ_full = local_occupations(full_file) + occ_ibz = symmetrize(local_occupations(ibz_file), ibz_file) + occ_full_symm = symmetrize(occ_full, ibz_file) + self.assertEqual(len(occ_full), len(occ_ibz)) + for icrsh in range(len(occ_full)): + np.testing.assert_allclose(occ_ibz[icrsh], occ_full_symm[icrsh], rtol=0, atol=1e-8, + err_msg=f"correlated shell {icrsh}") + np.testing.assert_allclose(occ_ibz[icrsh], occ_full[icrsh], rtol=0, atol=tol, + err_msg=f"correlated shell {icrsh}") + + def test_convert_svo_ibz(self): + """SrVO3: single V, cubic O_h, t2g (l=2) projection.""" + self.check_ibz('svo', tol=1e-10) + + def test_convert_nio_afm_ibz(self): + """ + NiO in the AFM cell (non-spin-polarized): full Ni d (l=2) and O p (l=1), + non-cubic cell, operations that exchange the two O. The full-grid VASP + projectors are only symmetric to ~1e-5 here. + """ + self.check_ibz('nio_afm', tol=1e-4) + + + def test_rotation_representation(self): + """D^l(R) is orthogonal and a representation, D(R1 R2) = D(R1) D(R2).""" + rng = np.random.default_rng(1) + def random_op(): + q, r = np.linalg.qr(rng.normal(size=(3, 3))) + return q * np.sign(np.diag(r)) # random O(3) element, proper or improper + for l in symmetry.SUPPORTED_L: + for _ in range(5): + R1, R2 = random_op(), random_op() + D1, D2 = (symmetry.real_harmonic_rotation(R, l) for R in (R1, R2)) + np.testing.assert_allclose(D1 @ D1.T, np.eye(2 * l + 1), atol=1e-12) + np.testing.assert_allclose(symmetry.real_harmonic_rotation(R1 @ R2, l), D1 @ D2, atol=1e-12) + + def test_unsupported_shells(self): + """f shells and non-orthonormal transforms are refused.""" + vasp_h5 = _rpath + 'svo/vaspout.h5' + f_shell = [{'atom': 2, 'sort': 2, 'l': 3, 'dim': 7, 'SO': 0, 'irep': 0}] + with self.assertRaises(NotImplementedError): + symmetry.build_symmcorr(vasp_h5, f_shell, [np.eye(7)], 0, 0) + d_shell = [{'atom': 2, 'sort': 2, 'l': 2, 'dim': 5, 'SO': 0, 'irep': 0}] + with self.assertRaises(RuntimeError): + symmetry.build_symmcorr(vasp_h5, d_shell, [2 * np.eye(5)], 0, 0) + with self.assertRaises(NotImplementedError): + symmetry.build_symmcorr(vasp_h5, d_shell, [np.eye(5)], 0, 1) + + def check_fallback(self, name, cfg, **kwargs): + """use_ibz=None falls back to the full grid, use_ibz=True raises.""" + generate_and_output_as_text(_rpath + cfg + '.cfg', _rpath + name + '/') + test_file = _rpath + cfg + '_fallback.test.h5' + if os.path.exists(test_file): os.remove(test_file) + converter = Converter(filename=_rpath + cfg, hdf_filename=test_file) + converter.convert_dft_input(vasp_h5=_rpath + name + '/vaspout.h5', **kwargs) + with HDFArchive(test_file, 'r') as ar: + self.assertEqual(ar['dft_input']['symm_op'], 0) + self.assertGreater(ar['dft_input']['n_k'], ar['dft_misc_input']['n_k_ibz']) + with self.assertRaises(RuntimeError): + Converter(filename=_rpath + cfg, hdf_filename=test_file).convert_dft_input( + use_ibz=True, vasp_h5=_rpath + name + '/vaspout.h5', **kwargs) + + def test_fallback_missing_equivalent_atom(self): + """Only one of the two equivalent O is projected: no IBZ.""" + self.check_fallback('nio_afm', 'nio_afm_one_o') + + def test_fallback_wrong_operations(self): + """ + Wrong d rotation matrices (z^2 and x^2-y^2 exchanged) pass all structural + checks but not the comparison with the full grid. + """ + rotation = symmetry.real_harmonic_rotation + def wrong_rotation(R, l): + D = rotation(R, l) + if l == 2: + idx = [0, 1, 4, 3, 2] + D = D[np.ix_(idx, idx)] + return D + with unittest.mock.patch.object(symmetry, 'real_harmonic_rotation', wrong_rotation): + self.check_fallback('nio_afm', 'nio_afm') + + +if __name__ == '__main__': + import unittest + unittest.main(verbosity=2, buffer=False) diff --git a/test/python/vasp/converter/test_converter_svo.py b/test/python/vasp/converter/test_converter_svo.py index de35962..118871d 100644 --- a/test/python/vasp/converter/test_converter_svo.py +++ b/test/python/vasp/converter/test_converter_svo.py @@ -32,7 +32,8 @@ def test_convert_svo(self): converter = Converter(filename=_rpath + 'svo', hdf_filename=test_file) - converter.convert_dft_input() + # full-grid reference; the IBZ conversion is covered by test_converter_ibz + converter.convert_dft_input(use_ibz=False) expected_file = _rpath + 'svo.ref.h5' self.assertH5FileEqual(test_file, expected_file)