[vasp] run DMFT on the irreducible BZ with symmetrization

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

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.

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)
 

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)