[vasp] add feat to read dft_bands_input for VASP calculations

see TRIQS/dftkit#3 for more information.

- add test for spaghettis for vasp
- add documentation for new vasp bands converter feature

Author: the-hampel

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 <dft.sumk_dft_tools.SumkDFTTools.partial_charges>`.
 
-

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 <dft.sumk_dft.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 <dft.sumk_dft_tools.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 <dft.sumk_dft_tools.SumkDFTTools.spaghettis>` can calculate total and VASP orbital-projected spectral functions. Use ``proj_type='vasp'`` for orbital-projected VASP spaghetti plots.
+
 PLOVASP detailed guide
 ======================
 

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()
-
-

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