diff --git a/doc/guide/analysis.rst b/doc/guide/analysis.rst index 66620f58..50f6a19e 100644 --- a/doc/guide/analysis.rst +++ b/doc/guide/analysis.rst @@ -125,6 +125,8 @@ This spectral function is calculated by typing:: The figure above shows the DFT SrVO\ :sub:`3`\ spaghetti plot (generated using V t\ :sub:`2g`\ Wannier projectors generated within a correlated energy window of [-13.6, 13.6] eV). As before, the broadening input has been set to the temperature (i.e., 1/Beta). The left panel shows the total A(k, :math:`\omega`) whereas the right gives the Wannier A(k, :math:`\omega`), both generated from this SK.spaghettis(). +For VASP inputs containing band projectors, e.g. data converted from a VASP ``KPOINTS_OPT``/``LOCPROJ_OPT`` calculation, orbital-projected spaghetti plots are available with ``proj_type='vasp'``. + Energy contours of the k-resolved Spectral function --------------------------------------------------- @@ -161,4 +163,3 @@ which calculates the partial charges using the self energy, double counting, and defined in the list `SK.shells`. This list is constructed by the Wien2k converter routines and stored automatically in the hdf5 archive. For the structure of `dm`, see also :meth:`partial charges `. - diff --git a/doc/guide/conv_vasp.rst b/doc/guide/conv_vasp.rst index 4dd6d6d4..02a627e5 100644 --- a/doc/guide/conv_vasp.rst +++ b/doc/guide/conv_vasp.rst @@ -27,7 +27,7 @@ Here, we will present a guide how the interface `can` be used to create input fo For VASP version older than 6.5.0 there are a few limitation of the interface as it was not officially supported and only text file based. See `Remarks for VASP older than 6.5.0`_ for details. -Generation of projectors for k-point lines (option `Lines` in KPOINTS) needed for Bloch spectral function calculations is not possible at the moment. +For VASP versions that write ``LOCPROJ_OPT`` data for ``KPOINTS_OPT``, projectors for a separate k-point path or grid can be converted for Bloch spectral function calculations. This workflow uses the HDF5 output in ``vaspout.h5`` and requires VASP to be compiled with HDF5 support enabled. VASP: generating raw projectors =============================== @@ -188,6 +188,36 @@ in :class:`SumkDFT `, e.g.:: However, this should only be done after a careful study of the density matrix and the projected DOS in the localized basis. For the complete process for SrVO3 see the tutorial for the VASP interface `here <../tutorials/svo_vasp/svo_notebook.html>`_. +When VASP band projectors are available in the HDF5 archive, for example from a ``KPOINTS_OPT``/``LOCPROJ_OPT`` bands conversion, :meth:`SumkDFTTools.spaghettis ` supports VASP orbital-projected spectral functions with ``proj_type='vasp'``. + +Data for post-processing - Spectral functions +============================================= + +For momentum-resolved spectral functions along a separate VASP k-point path or grid, use ``KPOINTS_OPT`` together with ``LOCPROJ`` in the VASP calculation. Recent VASP versions write the corresponding optional-k-point projector output to: + +* ``LOCPROJ_OPT`` and ``PROJCAR_OPT`` for text output. +* ``vaspout.h5:/results/locproj_opt`` for HDF5 output. + +The HDF5 route is used by the converter. The same calculation also provides the optional-k-point eigenvalues and weights under ``vaspout.h5:/results/electron_eigenvalues_kpoints_opt``. The relevant datasets are: + +* ``results/electron_eigenvalues_kpoints_opt/eigenvalues`` +* ``results/electron_eigenvalues_kpoints_opt/fermiweights`` +* ``results/electron_eigenvalues_kpoints_opt/kpoint_coords`` +* ``results/electron_eigenvalues_kpoints_opt/kpoints_symmetry_weight`` +* ``results/locproj_opt/data`` +* ``results/locproj_opt/format`` +* ``results/locproj_opt/parameters/*`` + +The normal converter workflow will automatically detect the optional-k-point bands data into ``dft_bands_input``:: + + Converter.convert_dft_input() + +and will internally call the bands converter with the same config:: + + Converter.convert_bands_input() + +After this conversion, :meth:`SumkDFTTools.spaghettis ` can calculate total and VASP orbital-projected spectral functions. Use ``proj_type='vasp'`` for orbital-projected VASP spaghetti plots. + PLOVASP detailed guide ====================== diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py index 67d05b9a..3d1f07de 100644 --- a/python/triqs_dft_tools/sumk_dft_tools.py +++ b/python/triqs_dft_tools/sumk_dft_tools.py @@ -676,14 +676,15 @@ def spaghettis(self, mu=None, broadening=None, mesh=None, plot_shift=0.0, plot_r with_dc : boolean, optional If True the double counting correction is used. proj_type : string, optional - The type of projection used for the orbital-projected DOS. + The type of projection used for the orbital-projected spectral function. These projected spectral functions will be determined alongside the total spectral function. - By default, no projected DOS type will be calculated (the corresponding projected arrays will be empty). + By default, no projected spectral function will be calculated (the corresponding projected arrays will be empty). The following options are: - 'None' - Only total DOS calculated - 'wann' - Wannier DOS calculated from the Wannier projectors - 'wien2k' - Wien2k orbital-projected DOS from the wien2k theta projectors + 'None' - Only total A(k,w) calculated + 'wann' - Wannier A(k,w) calculated from the Wannier projectors + 'vasp' - Vasp orbital-projected A(k,w) from Vasp inputs + 'wien2k' - Wien2k orbital-projected A(k,w) from the wien2k theta projectors save_to_file : boolean, optional If True, text files with the calculated data will be created. @@ -705,7 +706,7 @@ def spaghettis(self, mu=None, broadening=None, mesh=None, plot_shift=0.0, plot_r # initialisation if (proj_type != None): - assert proj_type in ('wann', 'wien2k'), "'proj_type' must be either 'wann', 'wien2k'" + assert proj_type in ('wann', 'vasp', 'wien2k'), "'proj_type' must be either 'wann', 'vasp', 'wien2k'" if (proj_type != 'wann'): assert proj_type == self.dft_code, "proj_type must be from the corresponding dft inputs." things_to_read = ['n_k', 'n_orbitals', 'proj_mat', 'hopping'] @@ -825,6 +826,7 @@ def gen_Akw(self, mu, broadening, mesh, plot_shift, plot_range, shell_list, with proj_type : string Output the orbital-projected A(k,w) type from the following: 'wann' - Wannier A(k,w) calculated from the Wannier projectors + 'vasp' - Vasp orbital-projected A(k,w) from Vasp inputs 'wien2k' - Wien2k orbital-projected A(k,w) from the wien2k theta projectors Returns @@ -863,6 +865,11 @@ def gen_Akw(self, mu, broadening, mesh, plot_shift, plot_range, shell_list, with gf_struct = self.gf_struct_sumk.copy() dims = [self.corr_shells[ish]['dim'] for ish in range(n_shells)] shells_type = 'corr' + elif (proj_type == 'vasp'): + n_shells = self.n_corr_shells + gf_struct = self.gf_struct_sumk.copy() + dims = [self.corr_shells[ish]['dim'] for ish in range(n_shells)] + shells_type = 'corr' elif (proj_type == 'wien2k'): n_shells = self.n_shells gf_struct = [[(sp, self.shells[ish]['dim']) for sp in spn] @@ -914,10 +921,11 @@ def gen_Akw(self, mu, broadening, mesh, plot_shift, plot_range, shell_list, with G_loc[ish].zero() tmp = G_loc[ish].copy() tmp.zero() - tmp << self.proj_type_G_loc(G_latt_w, tmp, ik, ish, proj_type) + downfold_proj_type = 'wann' if proj_type == 'vasp' else proj_type + tmp << self.proj_type_G_loc(G_latt_w, tmp, ik, ish, downfold_proj_type) G_loc[ish] += tmp # Rotate to local frame - if (self.use_rotations): + if ((proj_type != 'vasp') and self.use_rotations): for ish in range(n_shells): jsh=shell_list[ish] for bname, gf in G_loc[ish]: @@ -928,7 +936,7 @@ def gen_Akw(self, mu, broadening, mesh, plot_shift, plot_range, shell_list, with pAkw_orb[ish][bname][ik,:,:,:] = -gf.data[numpy.where((mesh_val > om_minplot) & (mesh_val < om_maxplot)),:,:].imag/numpy.pi # shift pAkw_orb for plotting stacked k-resolved eps(k) curves - pAkw_orb[ish][sp][ik] += ik * plot_shift + pAkw_orb[ish][bname][ik] += ik * plot_shift # Collect data from mpi mpi.barrier() @@ -1057,5 +1065,3 @@ def print_hamiltonian(self): (ik, self.hopping[ik, 0, i, i].real)) f.write('\n') f.close() - - diff --git a/test/python/svo_vasp_bands.py b/test/python/svo_vasp_bands.py new file mode 100644 index 00000000..a2c5ef7b --- /dev/null +++ b/test/python/svo_vasp_bands.py @@ -0,0 +1,90 @@ +import numpy as np + +from h5 import HDFArchive +from triqs.gf import BlockGf, Gf +from triqs.utility.comparison_tests import assert_arrays_are_close + +from triqs_dft_tools.sumk_dft_tools import SumkDFTTools + + +def assert_spectral_data_is_valid(Akw, pAkw, pAkw_orb, n_k, n_om, dim): + for sp in ('up', 'down'): + assert Akw[sp].shape == (n_k, n_om) + assert pAkw[0][sp].shape == (n_k, n_om) + assert pAkw_orb[0][sp].shape == (n_k, n_om, dim, dim) + + assert np.all(np.isfinite(Akw[sp])) + assert np.all(np.isfinite(pAkw[0][sp])) + assert np.all(np.isfinite(pAkw_orb[0][sp])) + assert np.min(Akw[sp]) > -1.0e-12 + assert np.min(pAkw[0][sp]) > -1.0e-12 + + assert_arrays_are_close( + pAkw[0][sp], + pAkw_orb[0][sp].trace(axis1=2, axis2=3), + precision=1.0e-12, + ) + + +def make_diagonal_sigma_from_srvo3(): + with HDFArchive('SrVO3_Sigma.h5', 'r') as ar: + sigma_srvo3 = ar['dmft_transp_input']['Sigma_w'] + + block_list = [ + Gf(mesh=sigma_srvo3.mesh, target_shape=(3, 3)), + Gf(mesh=sigma_srvo3.mesh, target_shape=(3, 3)), + ] + sigma = BlockGf(name_list=['up', 'down'], block_list=block_list, make_copies=False) + sigma.zero() + + for sp in ('up', 'down'): + for orb in range(3): + sigma[sp].data[:, orb, orb] = sigma_srvo3[f'{sp}_{orb}'].data[:, 0, 0] + + return sigma + + +Sigma = make_diagonal_sigma_from_srvo3() +n_om = len(Sigma.mesh) + +SK = SumkDFTTools(hdf_file='svo_vasp_bands.test.h5', mesh=Sigma.mesh) + +Akw_vasp, pAkw_vasp, pAkw_orb_vasp = SK.spaghettis( + broadening=0.05, + mesh=Sigma.mesh, + with_Sigma=False, + with_dc=False, + proj_type='vasp', + save_to_file=False, +) + +Akw_wann, pAkw_wann, pAkw_orb_wann = SK.spaghettis( + broadening=0.05, + mesh=Sigma.mesh, + with_Sigma=False, + with_dc=False, + proj_type='wann', + save_to_file=False, +) + +assert_spectral_data_is_valid(Akw_vasp, pAkw_vasp, pAkw_orb_vasp, SK.n_k, n_om, 3) + +for sp in ('up', 'down'): + assert_arrays_are_close(Akw_vasp[sp], Akw_wann[sp], precision=1.0e-12) + assert_arrays_are_close(pAkw_vasp[0][sp], pAkw_wann[0][sp], precision=1.0e-12) + assert_arrays_are_close(pAkw_orb_vasp[0][sp], pAkw_orb_wann[0][sp], precision=1.0e-12) + +SK.set_Sigma([Sigma], transform_to_sumk_blocks=False) + +Akw_sigma, pAkw_sigma, pAkw_orb_sigma = SK.spaghettis( + with_Sigma=True, + with_dc=False, + proj_type='vasp', + save_to_file=False, +) + +assert_spectral_data_is_valid(Akw_sigma, pAkw_sigma, pAkw_orb_sigma, SK.n_k, n_om, 3) + +for sp in ('up', 'down'): + assert np.max(np.abs(Akw_sigma[sp] - Akw_vasp[sp])) > 1.0e-6 + assert np.max(np.abs(pAkw_sigma[0][sp] - pAkw_vasp[0][sp])) > 1.0e-6 diff --git a/test/python/svo_vasp_bands.test.h5 b/test/python/svo_vasp_bands.test.h5 new file mode 100644 index 00000000..a1fb4f2d Binary files /dev/null and b/test/python/svo_vasp_bands.test.h5 differ