diff --git a/.gitignore b/.gitignore index 1038bc4..8df655b 100644 --- a/.gitignore +++ b/.gitignore @@ -14,11 +14,14 @@ DOSCAR EIGENVAL IBZKPT LOCPROJ +LOCPROJ_OPT OSZICAR OUTCAR PCDAT PROCAR +PROCAR_OPT PROJCAR +PROJCAR_OPT REPORT STOPCAR WAVECAR @@ -33,6 +36,10 @@ wannier90.win *.nnkp # Generated VASP converter test artifacts (regenerated from vaspout.h5 at test time) -test/python/vasp/converter/*.ctrl -test/python/vasp/converter/*.pg[0-9] -test/python/vasp/converter/*.test.h5 +test/python/vasp/converter/**/*.ctrl +test/python/vasp/converter/**/*.pg[0-9] +test/python/vasp/converter/**/*.test.h5 +# DMFT driver run artifacts +checkpoint.h5 +log +err diff --git a/doc/ChangeLog.md b/doc/ChangeLog.md index 3cffa7f..f70657d 100644 --- a/doc/ChangeLog.md +++ b/doc/ChangeLog.md @@ -25,6 +25,7 @@ Find below an itemized list of changes in this release. ### VASP * Add a VASP driver for charge self-consistent DFT+DMFT calculations +* 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 * Warn on misplaced or unknown tags in the PLOVASP configuration ### Wannier90 diff --git a/doc/examples/vasp_csc_svo/INCAR b/doc/examples/vasp_csc_svo/INCAR new file mode 100644 index 0000000..8e28350 --- /dev/null +++ b/doc/examples/vasp_csc_svo/INCAR @@ -0,0 +1,34 @@ +System = AFM NiO + +# better convergence for small kpt grids +ISMEAR = 0 +SIGMA = 0.1 + +KSPACING = 0.2 + +# converge wave functions +EDIFF = 1.E-7 + +# the energy window to optimize projector channels (absolute) +EMIN = -1.0 +EMAX = 10.1 +NBANDS = 32 + +LMAXMIX = 4 + +# project to Ni d +LORBIT = 14 +LOCPROJ = "2 : d : Pr" + +# CSC flags +ICHARG = 5 +NELM = 1000 +NELMIN = 1000 +NELMDL = -2 +IMIX=4 +BMIX=0.2 +AMIX=0.02 +LSYNCH5=True + +LWAVE = False +LCHG = False diff --git a/doc/examples/vasp_csc_svo/POSCAR b/doc/examples/vasp_csc_svo/POSCAR new file mode 100644 index 0000000..7424924 --- /dev/null +++ b/doc/examples/vasp_csc_svo/POSCAR @@ -0,0 +1,13 @@ +SrVO3 + 3.859374 + +1.00000000 +0.0000000000 +0.0000000000 + +0.0000000000 +1.00000000 +0.0000000000 + +0.0000000000 +0.0000000000 +1.0000000000 +Sr V O + 1 1 3 +Direct + +0.0000000000 +0.0000000000 +0.0000000000 + +0.5000000000 +0.5000000000 +0.5000000000 + +0.5000000000 +0.5000000000 +0.0000000000 + +0.5000000000 +0.0000000000 +0.5000000000 + +0.0000000000 +0.5000000000 +0.5000000000 diff --git a/doc/examples/vasp_csc_svo/POTCAR b/doc/examples/vasp_csc_svo/POTCAR new file mode 100644 index 0000000..4d8e7fb --- /dev/null +++ b/doc/examples/vasp_csc_svo/POTCAR @@ -0,0 +1,6 @@ + PAW_PBE Sr_sv 07Sep2000 + TITEL = PAW_PBE Sr_sv 07Sep2000 + PAW_PBE V_sv 02Aug2007 + TITEL = PAW_PBE V_sv 02Aug2007 + PAW_PBE O 08Apr2002 + TITEL = PAW_PBE O 08Apr2002 diff --git a/doc/examples/vasp_csc_svo/plo.cfg b/doc/examples/vasp_csc_svo/plo.cfg new file mode 100644 index 0000000..4270bdf --- /dev/null +++ b/doc/examples/vasp_csc_svo/plo.cfg @@ -0,0 +1,15 @@ +[General] + +[Group 1] +SHELLS = 1 +NORMALIZE = True +EWINDOW = -1.4 2.0 + +[Shell 1] +LSHELL = 2 +IONS = 2 + +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 0.0 1.0 0.0 + diff --git a/doc/examples/vasp_csc_svo/scf/INCAR b/doc/examples/vasp_csc_svo/scf/INCAR new file mode 100644 index 0000000..8976406 --- /dev/null +++ b/doc/examples/vasp_csc_svo/scf/INCAR @@ -0,0 +1,27 @@ +System = NiO + +NCORE = 1 +KPAR = 4 + +ISMEAR = 0 +SIGMA = 0.1 + +KSPACING = 0.2 + +# converge wave functions +EDIFF = 1.E-10 +NELMIN = 20 +NBANDS = 32 + + +# the energy window to optimize projector channels (absolute) +EMIN = -1.0 +EMAX = 10.1 +NEDOS = 5001 + +LMAXMIX = 4 + +# project to Ni d and O p states +LORBIT = 14 +LOCPROJ = "2 : d : Pr" + diff --git a/doc/examples/vasp_csc_svo/scf/POSCAR b/doc/examples/vasp_csc_svo/scf/POSCAR new file mode 120000 index 0000000..fbebe8c --- /dev/null +++ b/doc/examples/vasp_csc_svo/scf/POSCAR @@ -0,0 +1 @@ +../POSCAR \ No newline at end of file diff --git a/doc/examples/vasp_csc_svo/scf/converter.py b/doc/examples/vasp_csc_svo/scf/converter.py new file mode 100644 index 0000000..d667d18 --- /dev/null +++ b/doc/examples/vasp_csc_svo/scf/converter.py @@ -0,0 +1,9 @@ +from triqs_dftkit.vasp import Converter +import triqs_dftkit.vasp.plovasp.converter as plo_converter + +# Generate and store PLOs +plo_converter.generate_and_output_as_text('plo.cfg', vasp_dir='./') + +# run the converter +Converter = Converter(filename = 'vasp',proj_or_hk='proj') +Converter.convert_dft_input() diff --git a/doc/examples/vasp_csc_svo/scf/plo.cfg b/doc/examples/vasp_csc_svo/scf/plo.cfg new file mode 120000 index 0000000..fc003cc --- /dev/null +++ b/doc/examples/vasp_csc_svo/scf/plo.cfg @@ -0,0 +1 @@ +../plo.cfg \ No newline at end of file diff --git a/doc/examples/vasp_csc_svo/vasp_modest_csc.py b/doc/examples/vasp_csc_svo/vasp_modest_csc.py new file mode 100644 index 0000000..4fdf48d --- /dev/null +++ b/doc/examples/vasp_csc_svo/vasp_modest_csc.py @@ -0,0 +1,192 @@ +# ============================================================================ +# Charge self-consistent (CSC) DFT+DMFT for SrVO3 with VASP + TRIQS/modest +# ============================================================================ +# +# This script drives a full CSC DFT+DMFT calculation. The structure is two +# nested loops: +# +# * outer loop (n_total_loops): one DFT charge update per iteration. VASP is +# kept alive as a persistent background process by the dftkit VASP driver; +# each outer iteration feeds the DMFT charge-density correction back into +# VASP and reads the updated Kohn-Sham Hamiltonian / projectors. +# * inner loop (n_dmft_loops): the DMFT self-consistency cycle (build the +# hybridization, solve the impurity with CT-SEG, update the self-energy) +# for the current DFT Hamiltonian. +# +# The run always ends on a DMFT step: after the last outer iteration the DFT +# code is not invoked again (see the guard around the charge update below). +# +# Required inputs in this directory: INCAR, POSCAR, POTCAR, KPOINTS and the +# PLOVasp projector definition plo.cfg. + +import numpy as np + +import triqs.utility.mpi as mpi + +from triqs.gfs import BlockGf, MeshImFreq + +from triqs_ctseg.solve_generic import solve_generic + +import triqs_modest as M + +from triqs_dftkit.vasp import Driver as VaspDriver, MPIHandler +from triqs_modest.dft_driver import DftDriver + +from h5 import HDFArchive + + +# --- Physical and run parameters --------------------------------------------- +seedname = "vasp" +beta = 20.0 # inverse temperature (1/eV) +U = 4.50 # Kanamori intra-orbital interaction (eV) +J = 0.65 # Hund's coupling (eV) +Up = U -2*J # inter-orbital interaction (rotationally invariant choice) +n_iw = 1000 # number of Matsubara frequencies +n_total_loops = 10 # outer CSC (DFT charge) iterations +n_dmft_loops = 1 # inner DMFT iterations per outer loop +n_iter_dft = 2 # VASP SCF steps per outer loop (>1 for ALGO that needs + # several charge updates before DMFT; Sigma is held fixed + # and the charge correction is recomputed between steps) + +# --- DFT driver -------------------------------------------------------------- +# Wrap the dftkit VASP driver in the modest DftDriver. The VASP driver launches +# vasp_command under MPI and converts the output via PLOVasp (plo.cfg). +driver = DftDriver(VaspDriver(seedname=seedname, plo_cfg="plo.cfg", + mpi_handler=MPIHandler(mpi_exec="mpirun -np 16"), + vasp_command="vasp_std")) + +# Run the initial DFT, convert the output, and return the target electron count +# together with the one-body elements (Kohn-Sham Hamiltonian + projectors). +target_density, obe = driver.one_body_elements_from_dft() +mpi.report(obe) + +# Build the embedding (correlated subspace) from the projector space. +E = M.make_embedding(obe.C_space); mpi.report(E.description(True)) + +# Local Kanamori interaction Hamiltonian on the impurity. +h_int = M.make_kanamori(E.sigma_names, E.imp_decomposition(0), U, Up, J, False, False) + +# Double-counting correction (Held's fully-localized-limit flavour). +DcTerm = M.DcSolver("NonPolarized", "cHeld", U, J) + +mesh = MeshImFreq(beta=beta, S = "Fermion", n_iw=n_iw) + +# DFT-only chemical potential and local Green's function, used to detect the +# degenerate block structure that the impurity quantities are symmetrized over. +mu_dft = M.find_chemical_potential(target_density, obe, beta, verbosity=False) +Gdft = E.extract(M.local_gf.gloc(mesh, obe, mu_dft))[0] +mpi.report(f"Gdft density= {Gdft.total_density().real}") + +deg_blocks = M.analyze_degenerate_blocks(Gdft) +mpi.report(f"deg_blocks= {deg_blocks}") + +# Initialise the self-energy: zero dynamic part plus the static double counting. +Sigma_imp_dc = DcTerm.dc_self_energy(Gdft) + +Sigma_imp_dynamic, Sigma_imp_static = E.make_zero_imp_self_energies(mesh)[0] +for i in range(len(Sigma_imp_static)): Sigma_imp_static[i] += Sigma_imp_dc[i] + +# DFT + DMFT loop +try: + for n_iter in range(n_total_loops): + mpi.report(f"\n{'='*60}\n=== Global DFT+DMFT iteration {n_iter+1}/{n_total_loops} ===\n{'='*60}") + + Hloc0 = E.extract(M.atomic_levels_and_delta.impurity_levels(obe))[0] + mpi.report(f"Hloc0= {[h[0,0].real for h in Hloc0]}") + + # for first iteration converge Sigma imp first + if n_iter == 0: + n_dmft_loops_loc = 4 + else: + n_dmft_loops_loc = n_dmft_loops + # Begin DMFT loop + for n_dmft_iter in range(n_dmft_loops_loc): + mpi.report(f"\n--- DMFT iteration {n_dmft_iter+1}/{n_dmft_loops_loc} " + f"(global iter {n_iter+1}/{n_total_loops}) ---") + + Sigma_imp_static_minus_dc = [block-dc for (block, dc) in zip(Sigma_imp_static, Sigma_imp_dc)] + + Sigma_C_dynamic, Sigma_C_static = E.embed([Sigma_imp_dynamic], [Sigma_imp_static_minus_dc]) + + mu = M.find_chemical_potential(target_density, obe, Sigma_C_dynamic, Sigma_C_static) + + Gloc = E.extract(M.local_gf.gloc(obe, mu, Sigma_C_dynamic, Sigma_C_static))[0] + + Eimp = [h-mu*np.eye(h.shape[0])-dc for (h,dc) in zip(Hloc0, Sigma_imp_dc)] + + Delta = M.symmetrize(M.hybridization(Eimp, Gloc, Sigma_imp_dynamic, Sigma_imp_static), deg_blocks) + + solver_params = dict(n_iw=n_iw, n_tau=10*n_iw, length_cycle=50, + n_cycles = int(4e+6/mpi.size), + n_warmup_cycles = int(1e+4), + perform_tail_fit=True, + fit_min_w=10, fit_max_w=14, + ) + solver_results = solve_generic(Delta, Eimp, h_int, **solver_params) + + Sigma_imp_dynamic, Sigma_imp_static = solver_results.Sigma_dynamic, solver_results.Sigma_HartreeFock + + solver_results.G_iw << M.symmetrize(solver_results.G_iw, deg_blocks) + solver_results.Sigma_iw << M.symmetrize(solver_results.Sigma_iw, deg_blocks) + solver_results.Sigma_dynamic << M.symmetrize(solver_results.Sigma_dynamic, deg_blocks) + # End of DMFT loop + # Update Double-counting Term + Sigma_imp_dc = DcTerm.dc_self_energy(solver_results.G_iw) + Eimp_dc = DcTerm.dc_energy(solver_results.G_iw) + + Sigma_imp_static_minus_dc = [block-dc for (block, dc) in zip(Sigma_imp_static, Sigma_imp_dc)] + + # Re-embed the self-energy + Sigma_C_dynamic, Sigma_C_static = E.embed([Sigma_imp_dynamic], [Sigma_imp_static_minus_dc]) + + # Re-compute the chemical potential + mu = M.find_chemical_potential(target_density, obe, Sigma_C_dynamic, Sigma_C_static) + + # Compute the charge density correction + N_k = M.charge_density_correction(obe, mu, Sigma_C_dynamic, Sigma_C_static) + + # Impurity interaction energy + Eint = 0.5 * np.real((solver_results.Sigma_iw*solver_results.G_iw).total_density()) + mpi.report(f"Eint= {Eint}") + Eint_m_dc = (Eint - Eimp_dc) + mpi.report(f"Eint-Edc= {Eint_m_dc}") + + mpi.report("Saving DFT + DMFT iteration...") + if mpi.is_master_node(): + with HDFArchive("checkpoint.h5", "a") as ar: + path = f"it={len(ar) + 1}" + ar.create_group(path) + ar[path]["mu"] = mu + ar[path]["nloc"] = Gloc.total_density().real + ar[path]["nimp"] = solver_results.G_iw.total_density().real + ar[path]["Delta_iw"] = Delta + ar[path]["Eimp"] = Eimp + ar[path]["Gloc_iw"] = Gloc + ar[path]["Gimp_iw"] = solver_results.G_iw + ar[path]["Sigma_iw"] = solver_results.Sigma_iw + ar[path]["Sigma_iw_static"] = solver_results.Sigma_HartreeFock + ar[path]["Sigma_dc"] = Sigma_imp_dc + + # Update the one-body Hamiltonian with the charge density correction. + # Skip on the last outer iteration: the calculation must end on a DMFT + # step, so there is no point triggering another VASP charge update whose + # result would never be used. + if n_iter < n_total_loops - 1: + mpi.report(f"Calling VASP charge update / DFT driver " + f"(global iter {n_iter+1}/{n_total_loops})...") + # Run n_iter_dft VASP charge-update steps with the self-energy held + # fixed. Between steps the projectors change, so recompute mu and the + # charge density correction from the freshly converted one-body + # elements. The last step is followed directly by DMFT, which + # recomputes everything, so no recompute is needed there. + for i_dft in range(n_iter_dft): + obe = driver.update_one_body_elements_with_charge_correction(N_k, Eint_m_dc)[1] + if i_dft == n_iter_dft - 1: + break + Sigma_C_dynamic, Sigma_C_static = E.embed([Sigma_imp_dynamic], [Sigma_imp_static_minus_dc]) + mu = M.find_chemical_potential(target_density, obe, Sigma_C_dynamic, Sigma_C_static) + N_k = M.charge_density_correction(obe, mu, Sigma_C_dynamic, Sigma_C_static) + +finally: + # Ensure VASP is killed even if the script crashes + driver.driver.kill() diff --git a/python/triqs_dftkit/vasp/converter.py b/python/triqs_dftkit/vasp/converter.py index f4a6be2..a2f0b85 100644 --- a/python/triqs_dftkit/vasp/converter.py +++ b/python/triqs_dftkit/vasp/converter.py @@ -31,6 +31,7 @@ import numpy from h5 import * from ..converter_tools import * +import os import os.path try: import simplejson as json @@ -139,7 +140,10 @@ def read_header_and_data(self, filename): def convert_dft_input(self): """ - Reads the input files, and stores the data in the HDFfile + 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(). """ energy_unit = 1.0 # VASP interface always uses eV k_dep_projection = 1 @@ -409,11 +413,349 @@ def convert_dft_input(self): # 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) + # 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: + if not os.path.exists(candidate): + continue + try: + with HDFArchive(candidate, 'r') as ar: + _ = ar['results/electron_eigenvalues_kpoints_opt/eigenvalues'] + _ = ar['results/locproj_opt/data'] + kpoints_opt_found = True + break + except KeyError: + continue + + if kpoints_opt_found: + mpi.report("Detected KPOINTS_OPT band data in %s. Converting %s..." % (candidate, self.bands_subgrp)) + self.convert_bands_input() + # TODO: Implement misc_input # self.convert_misc_input(bandwin_file=self.bandwin_file,struct_file=self.struct_file,outputs_file=self.outputs_file, # misc_subgrp=self.misc_subgrp,SO=self.SO,SP=self.SP,n_k=self.n_k) + def convert_bands_input(self, cfg_filename=None): + """ + Reads KPOINTS_OPT band and projector data from vaspout.h5 and stores + it in the bands_subgrp in the hdf5 archive. + + The PLO config is required so shell transforms, energy windows, and + orthonormalization are applied via the same PLOVasp path as for + convert_dft_input(). + + This routine requires that convert_dft_input() has been called first, + because correlated shell definitions are taken from dft_subgrp. + """ + + def _decode_string(value): + if isinstance(value, bytes): + return value.decode('ascii').strip() + return str(value).strip() + + def _read_kpath_labels(vaspout_path, n_k): + """ + Read the high-symmetry k-path labels for a KPOINTS_OPT line-mode run + directly from vaspout.h5 (/input/kpoints_opt) and map them onto the + flattened band k-point index. + + Returns (labels, idx) where labels is a list of label strings and idx + is a 0-based numpy int array giving, for each label, the position of + that high-symmetry point in the n_k band path. Consecutive duplicate + labels at segment boundaries (e.g. the shared endpoint of two adjacent + segments) are collapsed into a single tick. Returns (None, None) if no + line-mode label data is available. + """ + try: + with HDFArchive(vaspout_path, 'r') as ar: + kopt = ar['input/kpoints_opt'] + mode = _decode_string(kopt['mode']) + raw_labels = kopt['labels_kpoints'] + nkps = int(kopt['number_kpoints']) + except KeyError: + return None, None + + if mode.lower() != 'l' or nkps <= 0: + return None, None + + labels = [_decode_string(lbl) for lbl in raw_labels] + n_seg = n_k // nkps + # KPOINTS_OPT line mode stores two labels (start, end) per segment. + if 2 * n_seg != len(labels): + mpi.report("convert_bands_input: KPOINTS_OPT label count (%i) inconsistent with %i segments; skipping k-path labels." % (len(labels), n_seg)) + return None, None + + merged_labels = [] + merged_idx = [] + for i, lab in enumerate(labels): + if not lab: + continue + seg = i // 2 + idx = seg * nkps if i % 2 == 0 else seg * nkps + nkps - 1 + # Collapse the shared endpoint of two adjacent segments. + if merged_labels and merged_labels[-1] == lab and idx - merged_idx[-1] == 1: + continue + merged_labels.append(lab) + merged_idx.append(idx) + + if not merged_labels: + return None, None + + return merged_labels, numpy.array(merged_idx, dtype=int) + + mpi.report("Processing VASP KPOINTS_OPT band/projector data...") + + def _to_complex(array): + arr = numpy.array(array) + if arr.shape and arr.shape[-1] == 2: + return arr[..., 0] + 1j * arr[..., 1] + return arr.astype(complex) + + def _locproj_to_canonical(locproj_data, locproj_format): + format_raw = locproj_format.strip() + if not (format_raw.startswith('[') and format_raw.endswith(']')): + raise IOError("convert_bands_input: Unexpected locproj_opt format string '%s'." % locproj_format) + + axis_labels = [item.strip() for item in format_raw[1:-1].split(',') if item.strip()] + axis_to_index = {label: i for i, label in enumerate(axis_labels)} + required_axes = ['proj_index', 'spin_index', 'kpts_index', 'band_index'] + if sorted(axis_to_index.keys()) != sorted(required_axes): + raise IOError("convert_bands_input: Unsupported locproj_opt axis labels '%s'." % axis_labels) + + locproj_complex = _to_complex(locproj_data) + if locproj_complex.ndim != len(axis_labels): + raise IOError("convert_bands_input: Unsupported locproj_opt data rank %i." % locproj_complex.ndim) + + perm = [axis_to_index[name] for name in required_axes] + return numpy.transpose(locproj_complex, axes=perm) + + def _find_cfg_path(user_cfg_filename): + if user_cfg_filename is not None: + return user_cfg_filename + cfg_candidates = [ + self.basename + '.cfg', + os.path.join(os.path.dirname(self.basename), 'plo.cfg') + ] + for candidate in cfg_candidates: + if os.path.exists(candidate): + return candidate + return None + + def _build_from_plovasp(cfg_path, eigvals_raw, fermiweights_raw, kpoint_coords, kpoint_weights, + locproj_raw, proj_sites_raw, proj_labels_raw, nc_flag, efermi): + from .plovasp.inpconf import ConfigParameters + from .plovasp.plotools import generate_plo + from .plovasp.vaspio import label_to_l_m + + class _BandsElStruct: + pass + + eigvals = numpy.array(eigvals_raw) + if eigvals.ndim == 2: + eigvals = eigvals[numpy.newaxis, :, :] + ferw = numpy.array(fermiweights_raw) + if ferw.ndim == 2: + ferw = ferw[numpy.newaxis, :, :] + + n_spin, n_k_loc, _ = eigvals.shape + + proj_params = [] + for ip in range(len(proj_labels_raw)): + l, m = label_to_l_m(proj_labels_raw[ip], ip, nc_flag) + proj_params.append({'isite': int(proj_sites_raw[ip]), 'l': l, 'm': m}) + + natom = max(int(max(proj_sites_raw)), 1) + qcoords = numpy.zeros((natom, 3), dtype=float) + for ip in range(len(proj_sites_raw)): + iat = int(proj_sites_raw[ip]) - 1 + if 0 <= iat < natom: + qcoords[iat, :] = 0.0 + + kweights = numpy.array(kpoint_weights, dtype=float) + if kweights.shape[0] != n_k_loc: + raise IOError("convert_bands_input: k-point weights have incompatible shape %s." % (kweights.shape,)) + wsum = kweights.sum() + if wsum > 0: + kweights = kweights / wsum + else: + kweights = numpy.full(n_k_loc, 1.0 / float(n_k_loc), dtype=float) + + pars = ConfigParameters(cfg_path, verbosity=0) + pars.parse_input() + if 'dosmesh' in pars.general: + del pars.general['dosmesh'] + efermi = pars.general.get('efermi', efermi) + + el_struct = _BandsElStruct() + el_struct.natom = natom + el_struct.type_of_ion = [0 for _ in range(natom)] + el_struct.kmesh = {'nktot': n_k_loc, 'nkibz': n_k_loc, 'kpoints': kpoint_coords, 'kweights': kweights} + el_struct.nc_flag = nc_flag + el_struct.efermi = efermi + # generate_plo expects eigvals with shape [nk, nb, ns] + el_struct.eigvals = numpy.transpose(eigvals, (1, 2, 0)) + # ferw is used as [spin, k, band] + el_struct.ferw = ferw + el_struct.proj_raw = locproj_raw + el_struct.proj_params = proj_params + el_struct.structure = {'qcoords': qcoords} + + pshells, pgroups = generate_plo(pars, el_struct, print_projector_diagnostics=False) + if len(pgroups) != 1: + raise IOError("convert_bands_input: Exactly one PLO group is supported, found %i in %s." % (len(pgroups), cfg_path)) + + pgroup = pgroups[0] + ib_win = pgroup.ib_win + nspin_ib = ib_win.shape[1] + + n_orbitals = numpy.zeros((n_k_loc, n_spin_blocs), dtype=int) + for isp in range(n_spin_blocs): + is_b = min(isp, nspin_ib - 1) + for ik in range(n_k_loc): + ib1, ib2 = int(ib_win[ik, is_b, 0]), int(ib_win[ik, is_b, 1]) + n_orbitals[ik, isp] = ib2 - ib1 + 1 + + nb_max = int(numpy.max(n_orbitals)) + hopping = numpy.zeros([n_k_loc, n_spin_blocs, nb_max, nb_max], complex) + eigvals_shifted = eigvals - efermi + for isp in range(n_spin_blocs): + is_b = min(isp, nspin_ib - 1) + is_e = min(isp, eigvals_shifted.shape[0] - 1) + for ik in range(n_k_loc): + ib1, ib2 = int(ib_win[ik, is_b, 0]), int(ib_win[ik, is_b, 1]) + nb = ib2 - ib1 + 1 + for ib in range(nb): + hopping[ik, isp, ib, ib] = eigvals_shifted[is_e, ik, ib1 + ib] + + max_corr_dim = max([crsh['dim'] for crsh in self.corr_shells]) + proj_mat = numpy.zeros([n_k_loc, n_spin_blocs, self.n_corr_shells, max_corr_dim, nb_max], complex) + + for icrsh, crsh in enumerate(self.corr_shells): + shell_atom = int(crsh['atom']) - 1 + shell_l = int(crsh['l']) + shell_dim = int(crsh['dim']) + + matches = [] + for ish, pshell in enumerate(pshells): + if not pshell.corr: + continue + if int(pshell.lorb) != shell_l: + continue + if shell_atom not in pshell.ion_list: + continue + io = pshell.ion_list.index(shell_atom) + if int(pshell.ndim) != shell_dim: + continue + matches.append((ish, io)) + + if len(matches) != 1: + raise IOError("convert_bands_input: Could not uniquely match correlated shell %i (atom=%i, l=%i, dim=%i) to transformed PLO shells from %s." % + (icrsh, shell_atom + 1, shell_l, shell_dim, cfg_path)) + + ish, io = matches[0] + pshell = pshells[ish] + ns_proj = pshell.proj_win.shape[1] + for isp in range(n_spin_blocs): + is_p = min(isp, ns_proj - 1) + is_b = min(isp, nspin_ib - 1) + for ik in range(n_k_loc): + ib1, ib2 = int(ib_win[ik, is_b, 0]), int(ib_win[ik, is_b, 1]) + nb = ib2 - ib1 + 1 + proj_mat[ik, isp, icrsh, :shell_dim, :nb] = pshell.proj_win[io, is_p, ik, :shell_dim, :nb] + + return n_orbitals, proj_mat, hopping + + if not (mpi.is_master_node()): + return + + # Read shell information from converter output + try: + with HDFArchive(self.hdf_file, 'r') as ar: + if not (self.dft_subgrp in ar): + raise IOError("convert_bands_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp) + things_to_read = ['SP', 'SO', 'n_corr_shells', 'corr_shells'] + for it in things_to_read: + if not hasattr(self, it): + setattr(self, it, ar[self.dft_subgrp][it]) + except KeyError: + raise IOError("convert_bands_input: Needed data not found in hdf file. Call convert_dft_input first.") + + 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_h5 = None + for candidate in vaspout_candidates: + if os.path.exists(candidate): + vaspout_h5 = candidate + break + if vaspout_h5 is None: + raise IOError("convert_bands_input: Could not find vaspout.h5. Tried: %s" % vaspout_candidates) + + try: + with HDFArchive(vaspout_h5, 'r') as ar: + eigvals = numpy.array(ar['results/electron_eigenvalues_kpoints_opt/eigenvalues']) + fermiweights = numpy.array(ar['results/electron_eigenvalues_kpoints_opt/fermiweights']) + kpoint_coords = numpy.array(ar['results/electron_eigenvalues_kpoints_opt/kpoint_coords']) + kpoint_weights = numpy.array(ar['results/electron_eigenvalues_kpoints_opt/kpoints_symmetry_weight']) + efermi = float(ar['results/electron_dos/efermi']) + + locproj_data = numpy.array(ar['results/locproj_opt/data']) + locproj_format = _decode_string(ar['results/locproj_opt/format']) + lnoncollinear = bool(int(ar['results/locproj_opt/parameters/lnoncollinear'])) + proj_sites = numpy.array(ar['results/locproj_opt/parameters/site'], dtype=int) + proj_labels_raw = numpy.array(ar['results/locproj_opt/parameters/ang_type']) + proj_labels = [_decode_string(lbl) for lbl in proj_labels_raw] + except KeyError as err: + raise IOError("convert_bands_input: Missing KPOINTS_OPT dataset in %s: %s" % (vaspout_h5, err)) + + locproj_canonical = _locproj_to_canonical(locproj_data, locproj_format) + + if fermiweights.shape != eigvals.shape: + raise IOError("convert_bands_input: Fermi weights shape %s does not match eigenvalues shape %s." % (fermiweights.shape, eigvals.shape)) + + cfg_path = _find_cfg_path(cfg_filename) + if cfg_path is None: + raise IOError("convert_bands_input: Could not find the PLO config needed to convert KPOINTS_OPT consistently. Provide cfg_filename or place plo.cfg next to %s." % vaspout_h5) + + n_orbitals, proj_mat, hopping = _build_from_plovasp( + cfg_path, + eigvals, + fermiweights, + kpoint_coords, + kpoint_weights, + locproj_canonical, + proj_sites, + proj_labels, + lnoncollinear, + efermi + ) + n_k = kpoint_coords.shape[0] + + n_parproj = numpy.array([0]) + proj_mat_all = numpy.array([0]) + + kpts_labels, kpts_labels_idx = _read_kpath_labels(vaspout_h5, n_k) + + with HDFArchive(self.hdf_file, 'a') as ar: + if not (self.bands_subgrp in ar): + ar.create_group(self.bands_subgrp) + things_to_save = ['n_k', 'n_orbitals', 'proj_mat', 'hopping', 'n_parproj', 'proj_mat_all'] + if kpts_labels is not None: + things_to_save += ['kpts_labels', 'kpts_labels_idx'] + mpi.report(" Stored %i high-symmetry k-path labels: %s" % (len(kpts_labels), ', '.join(kpts_labels))) + for it in things_to_save: + ar[self.bands_subgrp][it] = locals()[it] + + def convert_misc_input(self, bandwin_file, struct_file, outputs_file, misc_subgrp, SO, SP, n_k): """ Reads input for the band window from bandwin_file, which is case.oubwin, diff --git a/python/triqs_dftkit/vasp/plovasp/plotools.py b/python/triqs_dftkit/vasp/plovasp/plotools.py index abf13b6..a7b077f 100644 --- a/python/triqs_dftkit/vasp/plovasp/plotools.py +++ b/python/triqs_dftkit/vasp/plovasp/plotools.py @@ -107,7 +107,7 @@ def check_data_consistency(pars, el_struct): # generate_plo() # ################################################################################ -def generate_plo(conf_pars, el_struct): +def generate_plo(conf_pars, el_struct, print_projector_diagnostics=True): """ Parameters ---------- @@ -156,43 +156,44 @@ def generate_plo(conf_pars, el_struct): #from h5 import HDFArchive #with HDFArchive(testout, 'w') as h5test: # h5test['hk'] = pgroup.hk -# DEBUG output - print("Density matrix:") - nimp = 0.0 - ov_all = [] - for ish in pgroup.ishells: - if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): - print(" Shell %i"%(ish + 1)) - dm_all, ov_all_ = pshells[ish].density_matrix(el_struct) - ov_all.append(ov_all_[0]) - spin_fac = 2 if dm_all.shape[0] == 1 else 1 - for io in range(dm_all.shape[1]): - print(" Site %i"%(io + 1)) - dm = spin_fac * dm_all[:, io, : ,:].sum(0) - for row in dm: - print(''.join(map("{0:14.7f}".format, row))) - ndm = dm.trace() - if pshells[ish].corr: - nimp += ndm - print(" trace: ", ndm) - print() - print(" Impurity density:", nimp) - print() - print("Overlap:") - for io, ov in enumerate(ov_all): - print(" Site %i"%(io + 1)) - print(ov[0,...]) - print() - print("Local Hamiltonian:") - for ish in pgroup.ishells: - if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): - print(" Shell %i"%(ish + 1)) - loc_ham = pshells[pgroup.ishells[ish]].local_hamiltonian(el_struct) - for io in range(loc_ham.shape[1]): - print(" Site %i (real | complex part)"%(io + 1)) - for row in loc_ham[:, io, :, :].sum(0): - print(''.join(map("{0:14.7f}".format, row.real))+' |'+''.join(map("{0:14.7f}".format, row.imag))) -# END DEBUG output + if print_projector_diagnostics: + print("Density matrix:") + nimp = 0.0 + ov_all = [] + for ish in pgroup.ishells: + if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): + print(" Shell %i"%(ish + 1)) + dm_all, ov_all_ = pshells[ish].density_matrix(el_struct) + ov_all.append(ov_all_[0]) + spin_fac = 2 if dm_all.shape[0] == 1 else 1 + for io in range(dm_all.shape[1]): + print(" Site %i"%(io + 1)) + dm = spin_fac * dm_all[:, io, : ,:].sum(0) + for row in dm: + print(''.join(map("{0:14.7f}".format, row))) + ndm = dm.trace() + if pshells[ish].corr: + nimp += ndm + print(" trace: ", ndm) + print() + print(" Impurity density:", nimp) + print() + print("Overlap:") + for io, ov in enumerate(ov_all): + print(" Site %i"%(io + 1)) + print(ov[0,...]) + print() + print("Local Hamiltonian:") + for ish in pgroup.ishells: + if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): + print(" Shell %i"%(ish + 1)) + loc_ham = pshells[pgroup.ishells[ish]].local_hamiltonian(el_struct) + for io in range(loc_ham.shape[1]): + print(" Site %i (real | complex part)"%(io + 1)) + for row in loc_ham[:, io, :, :].sum(0): + print(''.join(map("{0:14.7f}".format, row.real))+' |'+''.join(map("{0:14.7f}".format, row.imag))) + else: + print(" Skipping density matrix/local Hamiltonian diagnostics for this projector path.") if 'dosmesh' in conf_pars.general: print() print("Evaluating DOS...") diff --git a/python/triqs_dftkit/vasp/plovasp/proj_group.py b/python/triqs_dftkit/vasp/plovasp/proj_group.py index 479dd72..7768b76 100644 --- a/python/triqs_dftkit/vasp/plovasp/proj_group.py +++ b/python/triqs_dftkit/vasp/plovasp/proj_group.py @@ -417,8 +417,11 @@ def orthogonalize_projector_matrix(self, p_matrix): overlap = np.dot(p_matrix, p_matrix.conj().T) # Calculate [O^{-1/2}]_{m m'} eig, eigv = np.linalg.eigh(overlap) - assert np.all(eig > 0.0), ("Negative eigenvalues of the overlap matrix:" - "projectors are ill-defined") + eps = 1e-12 + if np.any(eig < -eps): + raise AssertionError("Negative eigenvalues of the overlap matrix:projectors are ill-defined") + if np.any(eig <= eps): + raise AssertionError("Singular overlap matrix: projectors are linearly dependent") sqrt_eig = 1.0 / np.sqrt(eig) shalf = np.dot(eigv * sqrt_eig, eigv.conj().T) # Apply \tilde{P}_{m v} = \sum_{m'} [O^{-1/2}]_{m m'} P_{m' v} diff --git a/python/triqs_dftkit/vasp/plovasp/vaspio.py b/python/triqs_dftkit/vasp/plovasp/vaspio.py index e0d85d6..41c2465 100644 --- a/python/triqs_dftkit/vasp/plovasp/vaspio.py +++ b/python/triqs_dftkit/vasp/plovasp/vaspio.py @@ -45,6 +45,25 @@ log = logging.getLogger('plovasp.vaspio') +ORB_LABELS = ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2", + "fy(3x2-y2)", "fxyz", "fyz2", "fz3", "fxz2", "fz(x2-y2)", "fx(x2-3y2)"] + + +def label_to_l_m(label, iproj=None, nc_flag=False): + if isinstance(label, bytes): + label = label.decode('ascii') + label = str(label).strip() + if label not in ORB_LABELS: + raise ValueError("Unknown LOCPROJ orbital label '%s'" % label) + lm = ORB_LABELS.index(label) + l = int(np.sqrt(lm)) + m = lm - l * l + if nc_flag: + if iproj is None: + raise ValueError("Projector index is required for non-collinear LOCPROJ labels") + m = 2 * m if (iproj % 2) == 0 else 2 * m + 1 + return l, m + def read_lines(filename): r""" @@ -173,14 +192,6 @@ def locproj_parser(self, locproj_filename='LOCPROJ'): Returns projector parameters (site/orbital indices etc.) and an array with projectors. """ - orb_labels = ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2", - "fy(3x2-y2)", "fxyz", "fyz2", "fz3", "fxz2", "fz(x2-y2)", "fx(x2-3y2)"] - - def lm_to_l_m(lm): - l = int(np.sqrt(lm)) - m = lm - l * l - return l, m - # Read the first line of LOCPROJ to get the dimensions with open(locproj_filename, 'rt') as f: line = f.readline() @@ -221,20 +232,13 @@ def lm_to_l_m(lm): sline = line.split(':') isite = int(sline[1].split()[0]) label = sline[-1].strip() - lm = orb_labels.index(label) - l, m = lm_to_l_m(lm) + l, m = label_to_l_m(label, ip, self.nc_flag) # ip_new = iproj_site * norb + il # ip_prev = (iproj_site - 1) * norb + il proj_params[ip]['label'] = label proj_params[ip]['isite'] = isite proj_params[ip]['l'] = l - if self.nc_flag == True: - if (ip % 2) == 0: - proj_params[ip]['m'] = 2 * m - else: - proj_params[ip]['m'] = 2 * m + 1 - else: - proj_params[ip]['m'] = m + proj_params[ip]['m'] = m ip += 1 @@ -787,9 +791,6 @@ def __init__(self, h5path): # self.proj_params, self.plo = self.locproj_parser(locproj_filename=vasp_dir + "LOCPROJ") - orb_labels = ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2", - "fy(3x2-y2)", "fxyz", "fyz2", "fz3", "fxz2", "fz(x2-y2)", "fx(x2-3y2)"] - self.plo = np.zeros((self.nproj, self.nspin, nk, self.nband), dtype=complex) for proj in range(self.nproj): @@ -800,11 +801,6 @@ def __init__(self, h5path): imag_plo = plo[proj, spin, kpt, band, 1] # Imaginary part self.plo[proj, spin, kpt, band] = complex(real_plo, imag_plo) - def lm_to_l_m(lm): - l = int(np.sqrt(lm)) - m = lm - l * l - return l, m - self.proj_params = [{} for i in range(self.nproj)] with HDFArchive(h5path, 'a') as archive: for it in range(self.nproj): @@ -814,16 +810,9 @@ def lm_to_l_m(lm): self.proj_params[it]['coord'] = projectors['coordinates'][it] for it in range(self.nproj): - lm = orb_labels.index(self.proj_params[it]['label'].strip()) - l, m = lm_to_l_m(lm) + l, m = label_to_l_m(self.proj_params[it]['label'], it, self.nc_flag) self.proj_params[it]['l'] = l - if self.nc_flag == True: - if (it % 2) == 0: - self.proj_params[it]['m'] = 2 * m - else: - self.proj_params[it]['m'] = 2 * m + 1 - else: - self.proj_params[it]['m'] = m + self.proj_params[it]['m'] = m # assert ip == nproj, "Number of projectors in the header is wrong in LOCPROJ" print("Read parameters: LOCPROJ") diff --git a/test/python/vasp/converter/svo/kpoints_opt/INCAR b/test/python/vasp/converter/svo/kpoints_opt/INCAR new file mode 100644 index 0000000..ec61890 --- /dev/null +++ b/test/python/vasp/converter/svo/kpoints_opt/INCAR @@ -0,0 +1,35 @@ +SYSTEM = SrVO3 +NCORE = 1 +KPAR = 4 +LMAXMIX=6 +EDIFF = 1.E-10 +# make sure that the wavefunction is properly converged +NELMIN = 30 + +# Pin NBANDS so the fixture is reproducible independent of the MPI rank count. +# VASP rounds the default NBANDS up to a multiple of the band-group size, so at +# 4 ranks it picks 23, the topmost eg pair drops out of the [Group 1] EWINDOW +# and n_orbitals changes, which breaks test_converter_svo_bands.py. +NBANDS = 24 + +# DOS energy window +NEDOS = 3001 + +# Smearing procedure +ISMEAR = -5 + +# the energy window to optimize projector channels +EMIN = -1.0 +EMAX = 10.1 + +# use the PAW channel optimization +LORBIT=14 + +# project to V d +LOCPROJ = 2 : d : Pr + +# use exact diag for kpoints_opt +KPOINTS_OPT_MODE = 2 + +# POTCAR used Sr_sv, V_pv, O + diff --git a/test/python/vasp/converter/svo/kpoints_opt/KPOINTS b/test/python/vasp/converter/svo/kpoints_opt/KPOINTS new file mode 100644 index 0000000..584cac2 --- /dev/null +++ b/test/python/vasp/converter/svo/kpoints_opt/KPOINTS @@ -0,0 +1,6 @@ +Automatic Mesh +0 +Gamma + 15 15 15 + 0 0 0 + diff --git a/test/python/vasp/converter/svo/kpoints_opt/KPOINTS_OPT b/test/python/vasp/converter/svo/kpoints_opt/KPOINTS_OPT new file mode 100644 index 0000000..1c80b9c --- /dev/null +++ b/test/python/vasp/converter/svo/kpoints_opt/KPOINTS_OPT @@ -0,0 +1,15 @@ +k-points for bandstructure using seekpath GAMMA-X-X-M-M-GAMMA-GAMMA-R-R-X-R-M +50 +line +reciprocal + 0.00000000 0.00000000 0.00000000 GAMMA + 0.00000000 0.50000000 0.00000000 X + + 0.00000000 0.50000000 0.00000000 X + 0.50000000 0.50000000 0.00000000 M + + 0.50000000 0.50000000 0.00000000 M + 0.00000000 0.00000000 0.00000000 GAMMA + + 0.00000000 0.00000000 0.00000000 GAMMA + 0.50000000 0.50000000 0.50000000 R diff --git a/test/python/vasp/converter/svo/kpoints_opt/POSCAR b/test/python/vasp/converter/svo/kpoints_opt/POSCAR new file mode 100644 index 0000000..1bc468d --- /dev/null +++ b/test/python/vasp/converter/svo/kpoints_opt/POSCAR @@ -0,0 +1,13 @@ +SrVO3 +1.0 + 3.8420900000 0.0000000000 0.0000000000 + 0.0000000000 3.8420900000 0.0000000000 + 0.0000000000 0.0000000000 3.8420900000 + Sr V O + 1 1 3 +Direct + 0.000000000 0.000000000 0.000000000 + 0.500000000 0.500000000 0.500000000 + 0.500000000 0.500000000 0.000000000 + 0.000000000 0.500000000 0.500000000 + 0.500000000 0.000000000 0.500000000 diff --git a/test/python/vasp/converter/svo/kpoints_opt/plo.cfg b/test/python/vasp/converter/svo/kpoints_opt/plo.cfg new file mode 100644 index 0000000..00c1621 --- /dev/null +++ b/test/python/vasp/converter/svo/kpoints_opt/plo.cfg @@ -0,0 +1,16 @@ +[General] +BASENAME = converter/svo/kpoints_opt/plo_full +DOSMESH = -3.0 3.0 2001 + +[Group 1] +SHELLS = 1 +NORMALIZE = True +EWINDOW = -1.4 2.0 + +[Shell 1] +LSHELL = 2 +IONS = 2 + +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 0.0 1.0 0.0 diff --git a/test/python/vasp/converter/svo/kpoints_opt/vaspout.h5 b/test/python/vasp/converter/svo/kpoints_opt/vaspout.h5 new file mode 100644 index 0000000..d61788e Binary files /dev/null and b/test/python/vasp/converter/svo/kpoints_opt/vaspout.h5 differ diff --git a/test/python/vasp/converter/test_converter_svo_bands.py b/test/python/vasp/converter/test_converter_svo_bands.py new file mode 100644 index 0000000..3a52a87 --- /dev/null +++ b/test/python/vasp/converter/test_converter_svo_bands.py @@ -0,0 +1,177 @@ +import os +import rpath +_rpath = os.path.dirname(rpath.__file__) + '/' + +from h5 import HDFArchive +import numpy as np + +from triqs_dftkit.vasp.plovasp.converter import generate_and_output_as_text +from triqs_dftkit.vasp import Converter +import mytest + + +class TestConverterSVOBands(mytest.MyTestCase): + """ + Test conversion of KPOINTS_OPT + LOCPROJ_OPT into dft_bands_input. + """ + + def _check_bands_payload(self, test_file, vasp_dir): + with HDFArchive(test_file, 'r') as ar: + assert 'dft_bands_input' in ar, "Missing dft_bands_input group" + bands = ar['dft_bands_input'] + + things = ['n_k', 'n_orbitals', 'proj_mat', 'hopping', 'n_parproj', 'proj_mat_all'] + for it in things: + assert it in bands, "Missing key in dft_bands_input: %s" % it + + n_k = int(bands['n_k']) + n_orbitals = bands['n_orbitals'] + proj_mat = bands['proj_mat'] + hopping = bands['hopping'] + n_orb_min = int(np.min(n_orbitals)) + n_orb_max = int(np.max(n_orbitals)) + + assert n_k == 200, "Unexpected number of k-points in bands data" + self.assertEqual(n_orbitals.shape, (200, 1)) + self.assertEqual(n_orb_min, 3) + self.assertEqual(n_orb_max, 5) + self.assertEqual(proj_mat.shape[0:4], (200, 1, 1, 3)) + self.assertEqual(proj_mat.shape[4], n_orb_max) + self.assertEqual(hopping.shape[0:2], (200, 1)) + self.assertEqual(hopping.shape[2], n_orb_max) + self.assertEqual(hopping.shape[3], n_orb_max) + + # High-symmetry k-path labels read from vaspout.h5 (/input/kpoints_opt). + # The path is GAMMA-X-M-GAMMA-R with 50 points per segment; the shared + # endpoints of adjacent segments are collapsed into a single tick. + assert 'kpts_labels' in bands, "Missing kpts_labels in dft_bands_input" + assert 'kpts_labels_idx' in bands, "Missing kpts_labels_idx in dft_bands_input" + self.assertEqual(list(bands['kpts_labels']), ['GAMMA', 'X', 'M', 'GAMMA', 'R']) + np.testing.assert_array_equal(bands['kpts_labels_idx'], [0, 49, 99, 149, 199]) + + with HDFArchive(vasp_dir + 'vaspout.h5', 'r') as src: + eig = src['results/electron_eigenvalues_kpoints_opt/eigenvalues'] + efermi = float(src['results/electron_dos/efermi']) + with HDFArchive(test_file, 'r') as ar: + ib1 = int(ar['dft_misc_input']['band_window'][0][0, 0]) - 1 + expected_h = eig[0, 0, ib1] - efermi + + with HDFArchive(test_file, 'r') as ar: + h00 = ar['dft_bands_input']['hopping'][0, 0, 0, 0] + self.assertAlmostEqual(h00.real, expected_h) + self.assertAlmostEqual(h00.imag, 0.0) + + self._check_locproj_opt_scale(vasp_dir) + self._check_downfolded_t2g(test_file) + + def _check_locproj_opt_scale(self, vasp_dir): + """ + Compare the KPOINTS_OPT projectors with the regular-mesh ones at GAMMA, + the one k-point the line-mode path and the 15x15x15 mesh have in common. + + Both groups must express the same localized orbitals, so at a shared + k-point sum_orb sum_band ||^2 has to agree. That sum is + invariant under the arbitrary band phase and under any unitary mixing + inside a degenerate multiplet, hence independent of the MPI + decomposition (verified to 1e-8 relative between 4 and 8 ranks). + + PLOVasp orthonormalizes the projectors, so an overall factor on the raw + KPOINTS_OPT amplitudes divides out again and is invisible in + dft_bands_input. Two kinds of error are therefore only detectable here: + an amplitude that scales with the MPI decomposition, and a LORBIT=14 + "optimal" PAW channel rebuilt from the interpolated KPOINTS_OPT + k-points instead of kept from the ground-state mesh. Both have occurred + during development of the VASP side of this interface. + """ + def to_complex(raw): + arr = np.asarray(raw) + # (proj, spin, k, band, 2) -> (proj, k, band), single spin channel + return arr[:, 0, ..., 0] + 1j * arr[:, 0, ..., 1] + + with HDFArchive(vasp_dir + 'vaspout.h5', 'r') as src: + proj_opt = to_complex(src['results/locproj_opt/data']) + proj_mesh = to_complex(src['results/locproj/data']) + kpts_opt = np.asarray(src['results/electron_eigenvalues_kpoints_opt/kpoint_coords']) + kpts_mesh = np.asarray(src['results/projectors/kpoints']) + + self.assertEqual(proj_opt.shape, (5, 200, 24)) + + ik_opt = int(np.argmin(np.abs(kpts_opt).sum(axis=1))) + ik_mesh = int(np.argmin(np.abs(kpts_mesh).sum(axis=1))) + for label, kpt in (('path', kpts_opt[ik_opt]), ('mesh', kpts_mesh[ik_mesh])): + np.testing.assert_allclose(kpt, 0.0, atol=1e-8, + err_msg='GAMMA not found in the %s k-points' % label) + + weight_opt = (np.abs(proj_opt[:, ik_opt, :]) ** 2).sum() + weight_mesh = (np.abs(proj_mesh[:, ik_mesh, :]) ** 2).sum() + np.testing.assert_allclose(weight_opt, weight_mesh, rtol=1e-6) + + def _check_downfolded_t2g(self, test_file): + """ + Downfold hopping onto the t2g shell with proj_mat and check the SrVO3 + t2g dispersion along GAMMA-X-M-GAMMA-R. This validates the projectors + themselves (not just shapes): a wrong orbital character, a k/band + misalignment or a broken label mapping all show up here, while the + overall projector normalization does not. + """ + with HDFArchive(test_file, 'r') as ar: + bands = ar['dft_bands_input'] + proj_mat = bands['proj_mat'] + hopping = bands['hopping'] + n_orbitals = bands['n_orbitals'] + label_idx = np.asarray(bands['kpts_labels_idx']) + + eigs = [] + for ik in range(proj_mat.shape[0]): + nb = int(n_orbitals[ik, 0]) + p = proj_mat[ik, 0, 0, :, :nb] + h = hopping[ik, 0, :nb, :nb] + eigs.append(np.linalg.eigvalsh(p @ h @ p.conj().T)) + eigs = np.array(eigs) + + # t2g bandwidth of SrVO3. + np.testing.assert_allclose(eigs.max() - eigs.min(), 2.4866, atol=1e-3) + + # Cubic symmetry: threefold degenerate at GAMMA and at R. + for tick in (0, 3, 4): # GAMMA, GAMMA (second pass), R + ik = int(label_idx[tick]) + self.assertAlmostEqual(float(np.ptp(eigs[ik])), 0.0, places=5) + + # Both GAMMA ticks are the same k-point and must give the same bands. + np.testing.assert_allclose(eigs[int(label_idx[0])], eigs[int(label_idx[3])], atol=1e-6) + + # Band bottom at GAMMA, band top at R. + np.testing.assert_allclose(eigs[int(label_idx[0])].mean(), -0.9745, atol=1e-3) + np.testing.assert_allclose(eigs[int(label_idx[4])].mean(), 1.5121, atol=1e-3) + self.assertAlmostEqual(float(eigs.min()), float(eigs[int(label_idx[0])].min()), places=5) + self.assertAlmostEqual(float(eigs.max()), float(eigs[int(label_idx[4])].max()), places=5) + + def test_convert_svo_bands_auto(self): + vasp_dir = _rpath + 'svo/kpoints_opt/' + + generate_and_output_as_text(vasp_dir + 'plo.cfg', vasp_dir) + + test_file = _rpath + 'svo_bands_auto.test.h5' + converter = Converter(filename=vasp_dir + 'plo_full', hdf_filename=test_file) + + converter.convert_dft_input() + + self._check_bands_payload(test_file, vasp_dir) + + def test_convert_svo_bands_explicit(self): + vasp_dir = _rpath + 'svo/kpoints_opt/' + + generate_and_output_as_text(vasp_dir + 'plo.cfg', vasp_dir) + + test_file = _rpath + 'svo_bands_explicit.test.h5' + converter = Converter(filename=vasp_dir + 'plo_full', hdf_filename=test_file) + + converter.convert_dft_input() + converter.convert_bands_input(cfg_filename=vasp_dir + 'plo.cfg') + + self._check_bands_payload(test_file, vasp_dir) + + +if __name__ == '__main__': + import unittest + unittest.main(verbosity=2, buffer=False)