diff --git a/abipy/dfpt/ddb.py b/abipy/dfpt/ddb.py
index 6d01aa086..13b46ac79 100644
--- a/abipy/dfpt/ddb.py
+++ b/abipy/dfpt/ddb.py
@@ -1139,7 +1139,7 @@ def anaget_phmodes_at_qpoints(self, qpoints=None, asr=2, chneut=1, dipdip=1, dip
                 anaddbnc_path = task.outpath_from_ext("anaddb.nc")
                 ncfile.phbands.read_non_anal_from_file(anaddbnc_path)
 
-            print("Calculation completed.\nAnaddb results available in dir:", task.workdir)
+            print("Calculation completed.\nAnaddb results available in:", task.workdir)
             return ncfile.phbands if not return_input else (ncfile.phbands, inp)
 
     def anaget_phbst_and_phdos_files(self, nqsmall=10, qppa=None, ndivsm=20, line_density=None, asr=2, chneut=1,
diff --git a/abipy/eph/varpeq.py b/abipy/eph/varpeq.py
index dcb529d00..1049ebf41 100644
--- a/abipy/eph/varpeq.py
+++ b/abipy/eph/varpeq.py
@@ -800,19 +800,16 @@ def plot_bqnu_with_phbands(self, phbands_qpath,
             raise RuntimeError(f"{len(phbands_qmesh.qpoints)=} != {np.product(phdos_ngqpt)=}")
 
         #with_ibz_b2dos = False
-
         for pstate in range(self.nstates):
             # Compute B2(E) by looping over the full BZ.
             bqnu_dos = np.zeros(len(phdos_mesh))
             for iq_bz, qpoint in enumerate(phbands_qmesh.qpoints):
-                q_weight = qpoint.weight
-                if abs(q_weight - 1.0/phdos_nqbz) > 1e-6:
-                    raise RuntimeError(f"abs(q_weight - 1.0/phdos_nqbz) > 1e-6, {q_weight=}, {1.0/phdos_nqbz=}")
-
                 freqs_nu = phbands_qmesh.phfreqs[iq_bz]
                 for w, b2 in zip(freqs_nu, b2_interp_state[pstate].eval_kpoint(qpoint), strict=True):
-                    bqnu_dos += q_weight * b2 * gaussian(phdos_mesh, width, center=w)
+                    bqnu_dos += b2 * gaussian(phdos_mesh, width, center=w)
 
+            # NB: all the q-weights in PHBST.nc are set to 1.
+            bqnu_dos /= phdos_nqbz
             bqnu_dos = Function1D(phdos_mesh, bqnu_dos)
 
             ax = ax_mat[pstate, 1]
diff --git a/abipy/flowtk/flows.py b/abipy/flowtk/flows.py
index 4ffb8d590..932948c07 100644
--- a/abipy/flowtk/flows.py
+++ b/abipy/flowtk/flows.py
@@ -3202,8 +3202,8 @@ class PhononFlow(Flow):
     """
     @classmethod
     def from_scf_input(cls, workdir: str, scf_input: AbinitInput, ph_ngqpt,
-                       with_becs=True, with_quad=False, with_flexoe=False,
-                       manager=None, allocate=True) -> PhononFlow:
+                       qptopt=1, with_becs=True, with_quad=False, with_flexoe=False,
+                       manager=None, allocate=True, **kwargs) -> PhononFlow:
         """
         Create a `PhononFlow` for phonon calculations from an `AbinitInput` defining a ground-state run.
 
@@ -3213,6 +3213,7 @@ def from_scf_input(cls, workdir: str, scf_input: AbinitInput, ph_ngqpt,
             ph_ngqpt: q-mesh for phonons. Must be a sub-mesh of the k-mesh used for
                 electrons. e.g if ngkpt = (8, 8, 8). ph_ngqpt = (4, 4, 4) is a valid choice
                 whereas ph_ngqpt = (3, 3, 3) is not!
+            qptopt: Option for the generation of q-points. Default: 1 i.e. use nsym spatial symmetries + TR.
             with_becs: True if Born effective charges are wanted.
             with_quad: Activate calculation of dynamical quadrupoles.
                 Note that only selected features are compatible with dynamical quadrupoles.
@@ -3220,6 +3221,8 @@ def from_scf_input(cls, workdir: str, scf_input: AbinitInput, ph_ngqpt,
             with_flexoe: True to activate computation of flexoelectric tensor.
             manager: |TaskManager| object. Read from `manager.yml` if None.
             allocate: True if the flow should be allocated before returning.
+            kwargs: Extra optional arguments passed to PhononWork.from_scf_task.
+                See doc string for possible kwargs.
 
         Return:
             :class:`PhononFlow` object.
@@ -3236,8 +3239,8 @@ def from_scf_input(cls, workdir: str, scf_input: AbinitInput, ph_ngqpt,
         if any(scf_ngkpt % ph_ngqpt != 0):
             raise ValueError("ph_ngqpt %s should be a sub-mesh of scf_ngkpt %s" % (ph_ngqpt, scf_ngkpt))
 
-        ph_work = PhononWork.from_scf_task(scf_task, ph_ngqpt, is_ngqpt=True, with_becs=with_becs,
-                                           with_quad=with_quad, with_flexoe=with_flexoe)
+        ph_work = PhononWork.from_scf_task(scf_task, ph_ngqpt, is_ngqpt=True, qptopt=qptopt, with_becs=with_becs,
+                                           with_quad=with_quad, with_flexoe=with_flexoe, **kwargs)
         flow.register_work(ph_work)
 
         if allocate: flow.allocate()
diff --git a/abipy/flowtk/works.py b/abipy/flowtk/works.py
index e9cc91fa5..7f3590a4a 100644
--- a/abipy/flowtk/works.py
+++ b/abipy/flowtk/works.py
@@ -1733,7 +1733,7 @@ def from_scf_task(cls, scf_task: ScfTask,
                 This option is the recommended one if the k-path contains two high symmetry k-points that are very close
                 as ndivsm > 0 may produce a very large number of wavevectors.
                 if 0, deactivate band structure calculation.
-            qptopt: Option for the generation of q-points. Default: 1
+            qptopt: Option for the generation of q-points. Default: 1 i.e. use msym spatial symmetries + TR.
             prtwf: Controls the output of the first-order WFK.
                 By default we set it to -1 when q != 0 so that AbiPy is still able
                 to restart the DFPT task if the calculation is not converged (worst case scenario)