diff --git a/abipy/abio/factories.py b/abipy/abio/factories.py
index 81ad21f97..6ab55e0cf 100644
--- a/abipy/abio/factories.py
+++ b/abipy/abio/factories.py
@@ -154,6 +154,34 @@ def _find_scf_nband(structure, pseudos, electrons, spinat=None):
     return int(nband)
 
 
+def _find_nscf_nband_from_gsinput(gs_input: AbinitInput) -> int:
+    """
+    Find the value of nband for a NSCF calculation based on the input of a previous gs calculation.
+
+    Args:
+        gs_input (AbinitInput): the ground state input.
+
+    Returns:
+        the value of nband.
+    """
+    scf_nband = gs_input.get("nband")
+    if scf_nband is None:
+        occopt = gs_input.get("occopt", 1)
+        tsmear = gs_input.get("tsmear", 0.01)
+        smearing = aobj.Smearing(occopt, tsmear)
+
+        # only nsppol is used in _find_scf_nband, so just set that with a meaningful value
+        nsppol = gs_input.get("nsppol", 1)
+        spin_mode = aobj.SpinMode(mode="test", nsppol=nsppol, nspinor=1, nspden=nsppol)
+
+        charge = gs_input.get("charge", 0)
+        electrons = aobj.Electrons(spin_mode=spin_mode, smearing=smearing, charge=charge)
+        spinat = gs_input.get('spinat', None)
+        scf_nband = _find_scf_nband(gs_input.structure, gs_input.pseudos, electrons, spinat)
+
+    return scf_nband + 10
+
+
 def _get_shifts(shift_mode: str, structure: Structure):
     """
     Gives the shifts based on the selected shift mode and on the symmetry of the structure.
@@ -379,6 +407,42 @@ def ion_ioncell_relax_and_ebands_input(structure, pseudos,
     return relax_multi + ebands_multi
 
 
+def scr_from_nscfinput(nscf_input, nband=None, ecuteps=3.0, ecutwfn=None, inclvkb=2, w_type="RPA", sc_mode="one_shot", hilbert=None, accuracy="normal"):
+    """Return a screening input."""
+    scr_input = nscf_input.deepcopy()
+    scr_input.pop_irdvars()
+    if nband is None:
+        nband = nscf_input.get("nband")
+    screening = aobj.Screening(ecuteps, nband, w_type=w_type, sc_mode=sc_mode,
+                               hilbert=hilbert, ecutwfn=ecutwfn, inclvkb=inclvkb)
+
+    scr_input.set_vars(screening.to_abivars())
+    scr_input.set_vars(_stopping_criterion("screening", accuracy))  # Dummy
+
+    return scr_input
+
+
+def sigma_from_inputs(nscf_input, scr_input, nband=None, ecutwfn=None, ecuteps=None, ecutsigx=None, ppmodel="godby", gw_qprange=1, accuracy="normal"):
+    """Return a sigma input."""
+    self_input = nscf_input.deepcopy()
+    self_input.pop_irdvars()
+    if nband is None:
+        nband = nscf_input.get("nband")
+    screening = aobj.Screening(ecuteps=scr_input["ecuteps"], nband=scr_input["nband"],
+                               w_type="RPA",
+                               sc_mode="one_shot",
+                               hilbert=None,
+                               ecutwfn=scr_input["ecutwfn"],)
+    ecuteps = ecuteps if ecuteps is not None else screening.ecuteps
+    self_energy = aobj.SelfEnergy(se_type="gw", sc_mode="one_shot", nband=nband, ecutsigx=ecutsigx, screening=screening,
+                                  gw_qprange=gw_qprange, ppmodel=ppmodel, ecuteps=ecuteps, ecutwfn=ecutwfn, gwpara=2)
+
+    self_input.set_vars(self_energy.to_abivars())
+    self_input.set_vars(_stopping_criterion("sigma", accuracy))  # Dummy
+
+    return self_input
+
+
 def g0w0_with_ppmodel_inputs(structure, pseudos,
                              kppa, nscf_nband, ecuteps, ecutsigx,
                              ecut=None, pawecutdg=None, shifts=(0.0, 0.0, 0.0),
@@ -418,55 +482,34 @@ def g0w0_with_ppmodel_inputs(structure, pseudos,
     """
 
     structure = Structure.as_structure(structure)
-    multi = MultiDataset(structure, pseudos, ndtset=4)
-
-    # Set the cutoff energies.
-    multi.set_vars(_find_ecut_pawecutdg(ecut, pawecutdg, multi.pseudos, accuracy))
-
-    scf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=shifts)
+    # Scf input
+    # Note that kppa and shift_mode here are dummy as, they will be overwritten just after.
+    scf_inp = scf_input(structure=structure, pseudos=pseudos, kppa=kppa,
+                        ecut=ecut, pawecutdg=pawecutdg, nband=None, accuracy=accuracy,
+                        spin_mode=spin_mode, smearing=smearing, charge=charge, scf_algorithm=scf_algorithm,
+                        shift_mode="Monkhorst-Pack")
+    ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=shifts)
+    scf_inp.set_vars(ksampling.to_abivars())
+    scf_inp.set_vars(istwfk="*1")
+    # The following is to keep the previous behavior
+    # - The spinat is not set
+    # - The number of bands is not adapted for spinat
+    # TODO: Should we consider changing that and update the reference files accordingly ?
+    scf_inp.pop_vars('spinat')
     scf_electrons = aobj.Electrons(spin_mode=spin_mode, smearing=smearing, algorithm=scf_algorithm,
                                    charge=charge, nband=None, fband=None)
-
-    if scf_electrons.nband is None:
-        scf_electrons.nband = _find_scf_nband(structure, multi.pseudos, scf_electrons)
-
-    multi[0].set_vars(scf_ksampling.to_abivars())
-    multi[0].set_vars(scf_electrons.to_abivars())
-    multi[0].set_vars(_stopping_criterion("scf", accuracy))
-
-    nscf_ksampling = aobj.KSampling.automatic_density(structure, kppa, chksymbreak=0, shifts=shifts)
-    nscf_electrons = aobj.Electrons(spin_mode=spin_mode, smearing=smearing, algorithm={"iscf": -2},
-                                    charge=charge, nband=nscf_nband, fband=None)
-
-    multi[1].set_vars(nscf_ksampling.to_abivars())
-    multi[1].set_vars(nscf_electrons.to_abivars())
-    multi[1].set_vars(_stopping_criterion("nscf", accuracy))
-    # nbdbuf
-
-    # Screening.
-    if scr_nband is None: scr_nband = nscf_nband
-    screening = aobj.Screening(ecuteps, scr_nband, w_type="RPA", sc_mode="one_shot",
-                               hilbert=None, ecutwfn=None, inclvkb=inclvkb)
-
-    multi[2].set_vars(nscf_ksampling.to_abivars())
-    multi[2].set_vars(nscf_electrons.to_abivars())
-    multi[2].set_vars(screening.to_abivars())
-    multi[2].set_vars(_stopping_criterion("screening", accuracy)) # Dummy
-
-    # Sigma.
-    if sigma_nband is None: sigma_nband = nscf_nband
-    self_energy = aobj.SelfEnergy("gw", "one_shot", sigma_nband, ecutsigx, screening,
-                                  gw_qprange=gw_qprange, ppmodel=ppmodel)
-
-    multi[3].set_vars(nscf_ksampling.to_abivars())
-    multi[3].set_vars(nscf_electrons.to_abivars())
-    multi[3].set_vars(self_energy.to_abivars())
-    multi[3].set_vars(_stopping_criterion("sigma", accuracy)) # Dummy
-
-    # TODO: Cannot use istwfk != 1.
-    multi.set_vars(istwfk="*1")
-
-    return multi
+    nband = _find_scf_nband(structure, scf_inp.pseudos, scf_electrons)
+    scf_inp.set_vars(nband=nband)
+    # Non-Scf input
+    nscf_inp = nscf_from_gsinput(gs_input=scf_inp, nband=nscf_nband, accuracy=accuracy)
+    # Scr input
+    scr_inp = scr_from_nscfinput(nscf_input=nscf_inp, nband=scr_nband, ecuteps=ecuteps, ecutwfn=None,
+                                 inclvkb=inclvkb, w_type="RPA", sc_mode="one_shot", hilbert=None, accuracy="normal")
+    # Sigma input
+    sigma_inp = sigma_from_inputs(nscf_input=nscf_inp, scr_input=scr_inp, nband=sigma_nband,
+                                  ecutwfn=None, ecuteps=None, ecutsigx=ecutsigx,
+                                  ppmodel=ppmodel, gw_qprange=gw_qprange)
+    return MultiDataset.from_inputs([scf_inp, nscf_inp, scr_inp, sigma_inp])
 
 
 def g0w0_convergence_inputs(structure, pseudos, kppa, nscf_nband, ecuteps, ecutsigx, scf_nband, ecut,
@@ -1118,15 +1161,18 @@ def scf_input(structure, pseudos, kppa=None, ecut=None, pawecutdg=None, nband=No
     return abinit_input
 
 
-def ebands_from_gsinput(gs_input, nband=None, ndivsm=15, accuracy="normal") -> AbinitInput:
+def ebands_from_gsinput(gs_input, nband=None, ndivsm=15, accuracy="normal",
+                        projection=None) -> AbinitInput:
     """
     Return an |AbinitInput| object to compute a band structure from a GS SCF input.
 
     Args:
-        gs_input:
-        nband:
-        ndivsm:
-        accuracy:
+        gs_input: the |AbinitInput| that was used to calculated the charge density.
+        nband: the number of bands to be used for the calculation. If None it will be
+            automatically generated.
+        ndivsm: Number of divisions used to sample the smallest segment of the k-path.
+        accuracy: Accuracy of the calculation.
+        projection: which projection should be performed. If None no projection, otherwise "l" or "lm"
 
     Return: |AbinitInput|
     """
@@ -1137,29 +1183,111 @@ def ebands_from_gsinput(gs_input, nband=None, ndivsm=15, accuracy="normal") -> A
 
     nscf_ksampling = aobj.KSampling.path_from_structure(ndivsm, gs_input.structure)
     if nband is None:
-        nband = gs_input.get("nband", gs_input.structure.num_valence_electrons(gs_input.pseudos)) + 10
+        nband = _find_nscf_nband_from_gsinput(gs_input)
 
     bands_input.set_vars(nscf_ksampling.to_abivars())
     bands_input.set_vars(nband=nband, iscf=-2)
     bands_input.set_vars(_stopping_criterion("nscf", accuracy))
 
+    if projection is None:
+        pass
+    elif projection == "l":
+        bands_input.set_vars(prtdos=3)
+    elif projection == "lm":
+        bands_input.set_vars(prtdos=3, prtdosm=1)
+    else:
+        raise ValueError(f"Unrecognized value for projection: {projection}")
+
     return bands_input
 
 
-def dos_from_gsinput(gs_input, dos_kppa, nband=None, accuracy="normal", pdos=False):
+def nscf_from_gsinput(gs_input, kppa=None, nband=None, accuracy="normal",
+                      shift_mode="Monkhorst-Pack") -> AbinitInput:
+    """
+    Return an |AbinitInput| object to perform a NSCF calculation from a GS SCF input.
+
+    Args:
+        gs_input: the |AbinitInput| that was used to calculated the charge density.
+        kppa: defines the kpt sampling used for the NSCF run. If None the kpoint sampling and
+            shifts will be the same as in the SCF input.
+        nband: the number of bands to be used for the calculation. If None it will be
+            automatically generated.
+        accuracy: accuracy of the calculation.
+        shift_mode: the mode to be used for the shifts. Options are "Gamma", "Monkhorst-Pack",
+            "Symmetric", "OneSymmetric". See ShiftMode object for more details. Only used if kppa
+            is not None.
 
+    Return: |AbinitInput|
+    """
     # create a copy to avoid messing with the previous input
-    dos_input = gs_input.deepcopy()
-    dos_input.pop_irdvars()
+    nscf_input = gs_input.deepcopy()
+    nscf_input.pop_irdvars()
+
+    if kppa is not None:
+        shift_mode = ShiftMode.from_object(shift_mode)
+        shifts = _get_shifts(shift_mode, gs_input.structure)
+        dos_ksampling = aobj.KSampling.automatic_density(nscf_input.structure, kppa, chksymbreak=0, shifts=shifts)
+        nscf_input.set_vars(dos_ksampling.to_abivars())
+
+    if nband is None:
+        nband = _find_nscf_nband_from_gsinput(gs_input)
+
+    nscf_input.set_vars(nband=nband, iscf=-2)
+    nscf_input.set_vars(_stopping_criterion("nscf", accuracy))
 
-    dos_ksampling = aobj.KSampling.automatic_density(dos_input.structure, dos_kppa, chksymbreak=0)
-    dos_input.set_vars(dos_ksampling.to_abivars())
-    dos_input.set_vars(iscf=-2, ionmov=0)
-    dos_input.set_vars(_stopping_criterion("nscf", accuracy))
+    return nscf_input
 
-    if pdos:
-        # FIXME
-        raise NotImplementedError()
+
+def dos_from_gsinput(gs_input, kppa=None, nband=None, accuracy="normal", dos_method="tetra",
+                     projection="l", shift_mode="Monkhorst-Pack") -> AbinitInput:
+    """
+    Return an |AbinitInput| object to perform a DOS calculation from a GS SCF input.
+
+    Args:
+        gs_input: the |AbinitInput| that was used to calculated the charge density.
+        kppa: defines the kpt sampling used for the NSCF run. If None the kpoint sampling and
+            shifts will be the same as in the SCF input.
+        nband: the number of bands to be used for the calculation. If None it will be
+            automatically generated.
+        accuracy: accuracy of the calculation.
+        dos_method: method to calculate the DOS in abinit (NB: not the one used from postprocessing
+            in abipy). Set to "tetra" for the tetrahedron method (prtdos 2 or 3). If "smearing",
+            occopt and tsmear will be taken from gs_input else a "smearing-type: smearing value"
+            (prtdos 1 or 4).
+        projection: which projection should be performed. If None no projection, otherwise "l" or "lm"
+        shift_mode: the mode to be used for the shifts. Options are "Gamma", "Monkhorst-Pack",
+            "Symmetric", "OneSymmetric". See ShiftMode object for more details. Only used if kppa
+            is not None.
+
+    Return: |AbinitInput|
+    """
+    dos_input = nscf_from_gsinput(gs_input, kppa=kppa, nband=nband, accuracy=accuracy, shift_mode=shift_mode)
+
+    if dos_method == "tetra":
+        if projection is None:
+            dos_input.set_vars(prtdos=2)
+        elif projection == "l":
+            dos_input.set_vars(prtdos=3)
+        elif projection == "lm":
+            dos_input.set_vars(prtdos=3, prtdosm=1)
+        else:
+            ValueError(f"Unrecognized value for projection: {projection}")
+    else:
+        if dos_method != "smearing":
+            smear_obj = aobj.Smearing.as_smearing(dos_method)
+            dos_input.set_vars(smear_obj.to_abivars())
+
+        if projection is None:
+            dos_input.set_vars(prtdos=1)
+        elif projection == "l":
+            dos_input.set_vars(prtdos=4)
+        elif projection == "lm":
+            raise ValueError("lm projection is only allowed for dos_method 'tetra'")
+        else:
+            ValueError(f"Unrecognized value for projection: {projection}")
+
+    if projection is not None and "m" in projection.lower():
+        dos_input.set_vars(prtdosm=1)
 
     return dos_input
 
@@ -1219,13 +1347,17 @@ def scf_for_phonons(structure, pseudos, kppa=None, ecut=None, pawecutdg=None, nb
                     spin_mode="polarized", smearing="fermi_dirac:0.1 eV", charge=0.0, scf_algorithm=None,
                     shift_mode="Symmetric"):
 
+    # add the band for nbdbuf, if needed
+    nbdbuf = 4
+    if nband is not None:
+        nband += nbdbuf
+
     abiinput = scf_input(structure=structure, pseudos=pseudos, kppa=kppa, ecut=ecut, pawecutdg=pawecutdg, nband=nband,
                          accuracy=accuracy, spin_mode=spin_mode, smearing=smearing, charge=charge,
                          scf_algorithm=scf_algorithm, shift_mode=shift_mode)
 
-    nbdbuf = 4
-    # with no smearing set the minimum number of bands plus some nbdbuf
-    if smearing is None:
+    # with no bands set and no smearing the minimum number of bands plus some nbdbuf
+    if nband is None and smearing is None:
         nval = structure.num_valence_electrons(pseudos)
         nval -= abiinput['charge']
         nband = int(round(nval / 2) + nbdbuf)
diff --git a/abipy/abio/tests/test_factories.py b/abipy/abio/tests/test_factories.py
index 1144a4ba9..5685859d6 100644
--- a/abipy/abio/tests/test_factories.py
+++ b/abipy/abio/tests/test_factories.py
@@ -7,6 +7,7 @@
 from abipy.abio.factories import *
 from abipy.abio.factories import (BandsFromGsFactory, IoncellRelaxFromGsFactory, HybridOneShotFromGsFactory,
     ScfForPhononsFactory, PhononsFromGsFactory, PiezoElasticFactory, PiezoElasticFromGsFactory, ShiftMode)
+from abipy.abio.factories import _find_nscf_nband_from_gsinput
 import json
 
 write_inputs_to_json = False
@@ -23,6 +24,38 @@ def test_shiftmode(self):
             ShiftMode.from_object({})
 
 
+class HelperTest(AbipyTest):
+
+    @classmethod
+    def setUpClass(cls):
+        cls.si_structure = abidata.structure_from_cif("si.cif")
+        cls.si_pseudo = abidata.pseudos("14si.pspnc")
+
+    def test_find_nscf_nband_from_gsinput(self):
+        gsi = AbinitInput(self.si_structure, self.si_pseudo)
+        gsi.set_vars(nband=100)
+        assert _find_nscf_nband_from_gsinput(gsi) == 110
+
+        gsi = AbinitInput(self.si_structure, self.si_pseudo)
+        gsi.set_vars(occopt=1, nsppol=1)
+        assert _find_nscf_nband_from_gsinput(gsi) == 18
+
+        sc4 = self.si_structure.copy()
+        sc4.make_supercell([4, 4, 4])
+        gsi = AbinitInput(sc4, self.si_pseudo)
+        gsi.set_vars(occopt=1, nsppol=1)
+        assert _find_nscf_nband_from_gsinput(gsi) == 292
+
+        gsi.set_vars(occopt=3)
+        assert _find_nscf_nband_from_gsinput(gsi) == 318
+
+        gsi.set_vars(nsppol=2)
+        assert _find_nscf_nband_from_gsinput(gsi) == 318
+
+        gsi.set_vars(spinat=[[0, 0, 1]] * len(sc4))
+        assert _find_nscf_nband_from_gsinput(gsi) == 382
+
+
 class FactoryTest(AbipyTest):
 
     def setUp(self):
@@ -348,17 +381,37 @@ def test_scf_input(self):
         inp["ecut"] = 2
         self.abivalidate_input(inp)
 
-    def test_ebands_dos_from_gsinput(self):
+    def test_nscf_ebands_dos_from_gsinput(self):
         """Testing ebands_from_gsinput and dos_from_gsinput"""
-        from abipy.abio.factories import ebands_from_gsinput, dos_from_gsinput
+        from abipy.abio.factories import ebands_from_gsinput, dos_from_gsinput, nscf_from_gsinput
         gs_inp = gs_input(self.si_structure, self.si_pseudo, kppa=None, ecut=2, spin_mode="unpolarized")
+
+        nscf_inp = nscf_from_gsinput(gs_inp, kppa=None, nband=120)
+        self.assertArrayEqual(gs_inp["ngkpt"], nscf_inp["ngkpt"])
+        self.assertEqual(nscf_inp["nband"], 120)
+
         ebands_inp = ebands_from_gsinput(gs_inp, nband=None, ndivsm=15, accuracy="normal")
         self.abivalidate_input(ebands_inp)
 
+        ebands_inp = ebands_from_gsinput(gs_inp, nband=None, ndivsm=15, accuracy="normal", projection="lm")
+        self.assertEqual(ebands_inp["prtdos"], 3)
+        self.assertEqual(ebands_inp["prtdosm"], 1)
+
         dos_kppa = 3000
-        edos_inp = dos_from_gsinput(gs_inp, dos_kppa, nband=None, accuracy="normal", pdos=False)
+        edos_inp = dos_from_gsinput(gs_inp, dos_kppa, nband=None, accuracy="normal")
         self.abivalidate_input(edos_inp)
 
+        edos_inp = dos_from_gsinput(gs_inp, dos_method="smearing", projection="l")
+        self.assertEqual(gs_inp["occopt"], edos_inp["occopt"])
+        self.assertEqual(edos_inp["prtdos"], 4)
+
+        with self.assertRaises(ValueError):
+            edos_inp = dos_from_gsinput(gs_inp, dos_method="smearing", projection="lm")
+
+        edos_inp = dos_from_gsinput(gs_inp, dos_method="marzari5: 0.01 eV", projection="l")
+        self.assertEqual(edos_inp["occopt"], 5)
+        self.assertNotIn("prtdosm", edos_inp)
+
         factory_obj = BandsFromGsFactory(nband=None, ndivsm=15, accuracy="normal")
         self.assertMSONable(factory_obj)
         ebands_input_obj = factory_obj.build_input(gs_inp)
diff --git a/abipy/electrons/ebands.py b/abipy/electrons/ebands.py
index 26a1d6ccd..60c0fbf68 100644
--- a/abipy/electrons/ebands.py
+++ b/abipy/electrons/ebands.py
@@ -4741,6 +4741,17 @@ def plotly_up_minus_down(self, e0="fermie", fig=None, xlims=None, **kwargs):
 
         return fig
 
+    def to_pymatgen(self):
+        """
+        Return a pymatgen DOS object from an Abipy |ElectronDos| object.
+        """
+
+        from pymatgen.electronic_structure.dos import Dos
+        den = {s: d.values for d, s in zip(self.spin_dos, [PmgSpin.up, PmgSpin.down])}
+        pmg_dos = Dos(energies=self.spin_dos[0].mesh, densities=den, efermi=self.fermie)
+
+        return pmg_dos
+
 
 class ElectronDosPlotter(NotebookWriter):
     """
diff --git a/abipy/electrons/tests/test_ebands.py b/abipy/electrons/tests/test_ebands.py
index 9ed07a9ed..ddc9212ed 100644
--- a/abipy/electrons/tests/test_ebands.py
+++ b/abipy/electrons/tests/test_ebands.py
@@ -293,6 +293,8 @@ def test_silicon_ebands(self):
         with self.assertRaises(TypeError):
             ElectronDos.as_edos({}, {})
 
+        assert si_ebands_kmesh.to_pymatgen()
+
         mu = si_edos.find_mu(8)
         imu = si_edos.tot_idos.find_mesh_index(mu)
         self.assert_almost_equal(si_edos.tot_idos[imu][1], 8, decimal=2)