diff --git a/docs/source/postproc/tutorials/ephmat-spin.rst b/docs/source/postproc/tutorials/ephmat-spin.rst
new file mode 100644
index 00000000..2bb1dd12
--- /dev/null
+++ b/docs/source/postproc/tutorials/ephmat-spin.rst
@@ -0,0 +1,93 @@
+EphmatSpin tutorial
+===============
+
+In this section, we describe how to use Perturbopy to process a Perturbo ``'ephmat_spin'`` calculation. 
+
+The ``'ephmat_spin'`` calculation computes absolute values of the e-ph spin-flip matrix elements, summed over the number of electronic bands, given two lists of k-points and q-points. We first run the Perturbo calculation following the instructions on the Perturbo website and obtain the YAML file, *'diam_ephmat_spin.yml'*. For more information, please refer to the `Perturbo website <https://perturbo-code.github.io/mydoc_spin.html#e-ph-spin-flip-matrix-elementscalc_mode--ephmat_spin>`_. 
+
+Next, we create the :py:class:`.EphmatSpin` object using the YAML file as an input. This object contains all of the information from the YAML file.
+
+.. code-block :: python
+
+    import perturbopy.postproc as ppy
+
+    # Example using the ephmat_spin calculation mode.
+    diam_ephmat_spin = ppy.EphmatSpin.from_yaml('diam_ephmat_spin.yml')
+
+Accessing the data
+~~~~~~~~~~~~~~~~~~
+
+The attributes in an :py:class:`.EphmatSpin` object have the same name and format as an object from the :py:class:`.Ephmat`. The k-points are stored in :py:attr:`.EphmatSpin.kpt`, which is a :py:class:`RecipPtDB` object. The data for the phonons is stored analogously in :py:attr:`.EphmatSpin.qpt` (another :py:class:`RecipPtDB` object) and :py:attr:`EphmatSpin.phdisp`, a :py:class:`UnitsDict` object where the keys are the phonon mode and the values are the phonon energies computed across the q-points.
+
+.. code-block :: python
+    
+    # Access the k-point coordinates. There is only one in this calculation.
+    # The units are in crystal coordinates.
+    diam_ephmat_spin.kpt.points
+    >> array([[0.],
+              [0.],
+              [0.]])
+
+    diam_ephmat_spin.kpt.units
+    >> 'crystal'
+
+    # There are 196 q-point coordinates (we display the first two below).
+    # The units are in crystal coordinates.
+    diam_ephmat_spin.qpt.points.shape
+    >> (3, 196)
+
+    diam_ephmat_spin.qpt.points[:, :2]
+    >> array([[0.5   , 0.4902],
+              [0.5   , 0.4902],
+              [0.5   , 0.4902]])
+
+    diam_ephmat_spin.qpt.units
+    >> 'crystal'
+
+    # Access the phonon energies, which are a UnitsDict.
+    # There are 6 modes, which are the keys of the dictionary.
+    diam_ephmat_spin.phdisp.keys()
+    >> dict_keys([1, 2, 3, 4, 5, 6])
+
+    # Phonon energies of the first 2 q-points in phonon mode 3.
+    diam_ephmat_spin.phdisp[3][:2]
+    >> array([130.41105408, 130.31173133])
+
+    diam_ephmat_spin.phdisp.units
+    >> 'meV'
+
+Please see the section :ref:`handling_kpt_qpt` for more details on accessing information from :py:attr:`.EphmatSpin.kpt` and :py:attr:`.EphmatSpin.qpt`, such as labeling the k, q-points and converting to Cartesian coordinates.
+
+The ``'ephmat_spin'`` calculation interpolates the deformation potentials and e-ph elements from the spin-flip process which are stored in dictionaries :py:attr:`.EphmatSpin.defpot` and :py:attr:`.EphmatSpin.ephmat`, respectively. Both are :py:class:`UnitsDict` objects. The keys represent the phonon mode, and the values are (num_kpoints x num_qpoints) size arrays.
+
+.. code-block :: python
+
+    # There are 6 keys, one for each mode.
+    diam_ephmat_spin.ephmat.keys()
+    >> dict_keys([1, 2, 3, 4, 5, 6])
+
+    # There is 1 k-point and 196 q-points, so the e-ph matrix is 1 x 196.
+    diam_ephmat_spin.ephmat[1].shape
+    >> (1, 196)
+
+    # The e-ph spin-flip matrix elements corresponding to the first
+    # phonon mode, first (and only) k-point, and first two q-points.
+    diam_ephmat_spin.ephmat[1][0, :2]
+    >> array([[5.37973306e-06, 2.51372197e+00]])
+
+    # Units for the e-ph spin-flip matrix elements are in meV.
+    diam_ephmat_spin.ephmat.units
+    >> 'meV'
+
+    # We can extract analogous information from the deformation potential.
+    diam_ephmat_spin.defpot[1].shape
+    >> (1, 196)
+
+    # Units for the deformation potential are in eV/A.
+    diam_ephmat_spin.defpot.units
+    >> 'eV/A'
+
+Plotting the data
+-----------------
+
+Please refer to the :ref:`ephmat_tutorial` for details on plotting the data.
diff --git a/docs/source/postproc/tutorials/ephmat.rst b/docs/source/postproc/tutorials/ephmat.rst
index 64df1747..00e6ca7b 100644
--- a/docs/source/postproc/tutorials/ephmat.rst
+++ b/docs/source/postproc/tutorials/ephmat.rst
@@ -1,3 +1,5 @@
+.. _ephmat_tutorial:
+
 Ephmat tutorial
 ===============
 
diff --git a/docs/source/postproc/tutorials/figures/diam_spins.png b/docs/source/postproc/tutorials/figures/diam_spins.png
new file mode 100644
index 00000000..a1ffdf33
Binary files /dev/null and b/docs/source/postproc/tutorials/figures/diam_spins.png differ
diff --git a/docs/source/postproc/tutorials/figures/diam_spins_linear.png b/docs/source/postproc/tutorials/figures/diam_spins_linear.png
new file mode 100644
index 00000000..59dfb97b
Binary files /dev/null and b/docs/source/postproc/tutorials/figures/diam_spins_linear.png differ
diff --git a/docs/source/postproc/tutorials/imsigma-spin.rst b/docs/source/postproc/tutorials/imsigma-spin.rst
new file mode 100644
index 00000000..55991af6
--- /dev/null
+++ b/docs/source/postproc/tutorials/imsigma-spin.rst
@@ -0,0 +1,146 @@
+ImsigmaSpin tutorial
+====================
+
+In this section, we describe how to use Perturbopy to process a Perturbo ``'imsigma_spin'`` calculation.
+
+The ``'imsigma_spin'`` calculation mode computes the e-ph spin-flip self-energy for a list of electronic crystal momenta. We first run the Perturbo calculation following the instructions on the Perturbo website and obtain the YAML file, *'diam_imsigma_spin.yml'*. For more information, please refer to the `Perturbo website <https://perturbo-code.github.io/mydoc_spin.html#imaginary-part-of-e-ph-spin-flip-self-energycalc_mode--imsigma_spin>`_.
+
+Next, we create the :py:class:`.ImsigmaSpin` object using the YAML file as an input. This object contains all of the information from the YAML file.
+
+.. code-block :: python
+
+    import perturbopy.postproc as ppy
+
+    # Example using the imsigma_spin calculation mode.
+    diam_imsigma_spin = ppy.ImsigmaSpin.from_yaml('diam_imsigma_spin.yml')
+
+
+Accessing the data
+------------------
+
+The attributes in an :py:class:`.ImsigmaSpin` object have the same name and format as an object from the :py:class:`.Imsigma`. The main results of the results are categorized below: 
+
+* The k-points for which the e-ph self energies are computed and the corresponding band energies
+* Configuration information (temperature, chemical potential, carrier concentration)
+* The e-ph self energies for each electronic state, given both in total and by phonon mode
+
+See :ref:`exporting_data` to learn how to access data from ``diam_imsigma_spin`` that is general for all calculation modes, such as input parameters and the material's crystal structure.
+
+Bands data
+~~~~~~~~~~
+
+The k-points used for the ``'ephmat_spin'`` calculation are stored in the :py:attr:`.ImsigmaSpin.kpt` attribute, which is of type :py:class:`.RecipPtDB`. For example, to access the k-point coordinates and their units,
+
+.. code-block :: python
+  
+    # The first k-point (points are stored column-oriented).
+    diam_imsigma_spin.kpt.points[:, 0]
+    >> array([0.    , 0.25  , 0.2625])
+
+    # The points are in crystal coordinates.
+    diam_imsigma_spin.kpt.units
+    >> 'crystal'
+
+Please see the section :ref:`handling_kpt_qpt` for more details on accessing information from :py:attr:`.ImsigmaSpin.kpt`, such as labeling the k-points and converting to Cartesian coordinates.
+
+The  band energies are stored in the :py:attr:`.ImsigmaSpin.bands` attribute, which is a :py:class:`.UnitsDict` object. The keys represent the band index, and the values are arrays containing the band energies corresponding to each k-point. 
+
+.. code-block :: python
+
+    # There are two bands used in this calculation.
+    diam_imsigma_spin.bands.keys()
+    >> dict_keys([1, 2])
+
+    # Band energies of the first band corresponding to the first 10 k-points.
+    diam_imsigma_spin.bands[1][:10]
+    >> array([17.84019125, 17.75059409, 17.78470674, 17.69665692, 17.81793713,
+              17.76517492, 17.67960317, 17.83014117, 17.7809763 , 17.69861557])
+
+Please see the section :ref:`physical_quantities` for details on accessing the bands and their units.
+
+Configuration data
+~~~~~~~~~~~~~~~~~~
+
+The ``'imsigma_spin'`` calculations can be run for various system configurations, i.e. the temperature, chemical potential, and carrier concentration. Information about the configuration(s), are stored in the following attributes:
+
+* :py:attr:`.ImsigmaSpin.temper`
+* :py:attr:`.ImsigmaSpin.chem_pot`
+
+These attributes are :py:class:`.UnitsDict` objects, which are Python dictionaries with an additional attribute that stores the units. The keys of the dictionary represent the configuration number. The values are floats representing the temperature or chemical potential.
+
+For example, let's look at the temperatures.
+
+.. code-block :: python
+
+    # Keys are configuration number, values are temperatures.
+    diam_imsigma_spin.temper
+    >> {1: 25.85203}
+    
+    # Units are in meV.
+    diam_imsigma_spin.temper.units
+    >> 'meV'
+
+Please see the section :ref:`physical_quantities` for details on working with :py:class:`UnitsDict` objects.
+
+ImsigmaSpin results
+~~~~~~~~~~~~~~~~~~~
+
+The e-ph self energies are stored in the :py:attr:`.ImsigmaSpin.imsigma` object, which is of type :py:class:`.UnitsDict`. There are two levels in this dictionary. The first level gives the configuration number. The second level gives the band index. The values are arrays of the e-ph self energies computed along all the k-points, at that configuration and band index.
+
+.. code-block :: python
+
+    # The first key is the configuration number.
+    # Here we have one configuration.
+    diam_imsigma_spin.imsigma.keys()
+    >> dict_keys([1])
+
+    # The second key is the band index. Here we are looking at configuration 1,
+    # and we have 2 bands (matching the diam_imsigma_spin.bands attribute).
+    diam_imsigma_spin.imsigma[1].keys()
+    >> dict_keys([1, 2])
+
+    # The e-ph spin-flip self-energy array for configuration 1 and
+    # band index 2. The array size matches the number of k-points.
+    diam_imsigma_spin.imsigma[1][2].shape
+    >> (815,)
+
+    # The e-ph spin-flip self-energies for configuration 1, band index 2,
+    # and the first 4 k-points.
+    diam_imsigma_spin.imsigma[1][2][:4]
+    >> array([1.29425589e-05, 7.97155145e-06, 1.06093255e-05, 8.78246442e-06])
+
+    # The units are meV.
+    diam_imsigma_spin.imsigma.units
+    >> 'meV'
+
+We can also get the e-ph spin-flip self energies for each phonon mode through the :py:attr:`.ImsigmaSpin.imsigma_mode` object. This dictionary is similar, but there is an additional level that identifies the phonon mode.
+
+.. code-block :: python
+
+    # The first key is the configuration number.
+    diam_imsigma_spin.imsigma_mode.keys()
+    >> dict_keys([1])
+    
+    # The second key is the phonon mode. Here we have 6 modes.
+    diam_imsigma_spin.imsigma_mode[1].keys()
+    >> dict_keys([1, 2, 3, 4, 5, 6])
+
+    # The third key is the band index. Here we are looking at configuration 1,
+    # phonon mode 3, and we see we have 2 bands. Note that this matches
+    # the diam_imsigma_spin.bands attribute.
+    diam_imsigma_spin.imsigma_mode[1][3].keys()
+    >> dict_keys([1, 2])
+
+    # The e-ph spin-flip self-energy array for configuration 1, phonon mode 3,
+    # and band index 2. The array size matches the number of k-points.
+    diam_imsigma_spin.imsigma_mode[1][3][2].shape
+    >> (2445,)
+
+    # The e-ph spin-flip self-energies for configuration 1, phonon mode 3,
+    # band index 2, and the first 4 k-points.
+    diam_imsigma_spin.imsigma_mode[1][3][2][:4]
+    >> array([2.71039146e-06, 0.00000000e+00, 0.00000000e+00, 1.83641809e-06])
+
+    # The units are meV.
+    diam_imsigma_spin.imsigma_mode.units
+    >> 'meV'
diff --git a/docs/source/postproc/tutorials/imsigma.rst b/docs/source/postproc/tutorials/imsigma.rst
index 45484986..77e80b3a 100644
--- a/docs/source/postproc/tutorials/imsigma.rst
+++ b/docs/source/postproc/tutorials/imsigma.rst
@@ -1,9 +1,11 @@
+.. _imsigma_tutorial:
+
 Imsigma tutorial
 ================
 
-In this section, we describe how to use Perturbopy to process a Perturbo imsigma calculation.
+In this section, we describe how to use Perturbopy to process a Perturbo ``'imsigma'`` calculation.
 
-The imsigma calculation mode computes the e-ph self-energy for a list of electronic crystal momenta. We first run the Perturbo calculation following the instructions on the Perturbo website and obtain the YAML file, *'si_imsigma.yml'*. For more information, please refer to the `Perturbo website <https://perturbo-code.github.io/mydoc_scattering.html#imaginary-part-of-e-ph-self-energycalc_mode--imsigma>`_.
+The ``'imsigma'`` calculation mode computes the e-ph self-energy for a list of electronic crystal momenta. We first run the Perturbo calculation following the instructions on the Perturbo website and obtain the YAML file, *'si_imsigma.yml'*. For more information, please refer to the `Perturbo website <https://perturbo-code.github.io/mydoc_scattering.html#imaginary-part-of-e-ph-self-energycalc_mode--imsigma>`_.
 
 Next, we create the :py:class:`.Imsigma` object using the YAML file as an input. This object contains all of the information from the YAML file.
 
@@ -20,7 +22,7 @@ Accessing the data
 
 The main results of the results are categorized below: 
 
-* the k-points for which the e-ph self energies are computed and the corresponding band energies
+* The k-points for which the e-ph self energies are computed and the corresponding band energies
 * Configuration information (temperature, chemical potential, carrier concentration)
 * The e-ph self energies for each electronic state, given both in total and by phonon mode
 
@@ -140,4 +142,4 @@ We can also get the e-ph self energies for each phonon mode through the :py:attr
 
     # The units are meV
     si_imsigma.imsigma_mode.units
-    >> 'meV'
\ No newline at end of file
+    >> 'meV'
diff --git a/docs/source/postproc/tutorials/index.rst b/docs/source/postproc/tutorials/index.rst
index 3dbb5ee8..2dc07571 100644
--- a/docs/source/postproc/tutorials/index.rst
+++ b/docs/source/postproc/tutorials/index.rst
@@ -11,10 +11,13 @@ Below, we give general information on how to use Perturbopy. We also provide tut
     bands
     phdisp
     ephmat
+    ephmat-spin
     trans
     imsigma
+    imsigma-spin
     dynamics-run
     dynamics-pp
+    spins
 
 
 General information
diff --git a/docs/source/postproc/tutorials/spins.rst b/docs/source/postproc/tutorials/spins.rst
new file mode 100644
index 00000000..f6da20cd
--- /dev/null
+++ b/docs/source/postproc/tutorials/spins.rst
@@ -0,0 +1,117 @@
+Spins tutorial
+===============
+
+In this section, we describe how to use Perturbopy to process a Perturbo ``'spins'`` calculation. 
+
+The ``'spins'`` calculation mode interpolates the band structure using Wannier functions and computes the spin texture. We first run the Perturbo calculation following the instructions on the Perturbo website and obtain the YAML file, *'diam_spins.yml'*. For more information, please refer to the `Perturbo website <https://perturbo-code.github.io/mydoc_spin#spin-texturecalc_mode--spins>`_. 
+
+Next, we create the :py:class:`.Spins` object using the YAML file as an input. This object contains all of the information from the YAML file.
+
+.. code-block :: python
+
+    import perturbopy.postproc as ppy
+
+    # Example using the spins calculation mode
+    diam_spins = ppy.Spins.from_yaml('diam_spins.yml')
+
+Accessing the data
+~~~~~~~~~~~~~~~~~~
+
+The outputs of the calculation are stored in three objects:
+* :py:attr:`.Spins.kpt` stores the k-points used in the calculation
+* :py:attr:`.Spins.bands` stores the interpolated band energies
+* :py:attr:`.Spins.spins` stores the spin texture values
+
+Please see the :ref:`exporting_data` section for more information on accessing data general to all calculation modes, such as input parameters and material properties.
+
+K-points
+--------
+
+The k-points used for the spins calculation are stored in the :py:attr:`.Spins.kpt` attribute, which is of type :py:class:`.RecipPtDB`. For example, to access the k-point coordinates and their units (which are column-oriented):
+
+.. code-block :: python
+    
+    # Obtain the first k-point
+    diam_spins.kpt.points[:, 0]
+    >> array([0.5, 0.5, 0.5])
+
+    # There are 196 k-points
+    diam_spins.kpt.points.shape
+    >> (3, 196)
+
+    diam_spins.kpt.units
+    >> 'crystal'
+
+Please see the section :ref:`handling_kpt_qpt` for details on accessing the k-points through this attribute.
+
+Band energies
+-------------
+
+The interpolated band energies computed by the spins calculation are stored in the :py:attr:`.Spins.bands` attribute, which is a :py:class:`.UnitsDict` object. The keys represent the band index, and the values are arrays containing the band energies corresponding to each k-point. 
+
+.. code-block :: python
+
+    # The keys are the band indices, and here we have 16
+    diam_spins.bands.keys()
+    >> dict_keys([1, 2, 3, ..., 14, 15, 16])
+
+    # Band energies of the 8th band
+    diam_spins.bands[8]
+    >> array([10.67315828, 10.67472505, ..., 13.51506129, 13.52024087])
+
+Please see the section :ref:`physical_quantities` for details on accessing the bands and their units.
+
+Spin textures
+-------------
+
+The spin texture values computed by the spins calculation are stored in the :py:attr:`.Spins.spins` attribute, which is a :py:class:`.UnitsDict` object. The keys represent the band index, and the values are arrays containing the spin texture values corresponding to each k-point. 
+
+.. code-block :: python
+
+    # The keys are the band indices, and here we have 16
+    diam_spins.spins.keys()
+    >> dict_keys([1, 2, 3, ..., 14, 15, 16])
+
+    # Spin texture values of the 8th band
+    diam_spins.spins[8]
+    >> array([5.77350351e-01, 5.77348291e-01, ..., 8.64282329e-01, 1.00000000e+00])
+
+Please see the section :ref:`physical_quantities` for details on accessing the spin texture values and their units.
+
+Plotting the data
+-----------------
+
+We can quickly visualize the spin texture by plotting them as a colormap overlaid on the band structure. Below, we plot the spin texture along the k-point path.
+
+.. code-block :: python
+    
+    import matplotlib.pyplot as plt
+
+    plt.rcParams.update(ppy.plot_tools.plotparams)
+    diam_spins.kpt.add_labels(ppy.lattice.points_fcc)
+
+    fig, ax  = plt.subplots()
+    diam_spins.plot_spins(ax)
+    plt.show()
+
+.. image:: figures/diam_spins.png
+    :width: 600
+    :align: center
+
+
+We can choose whether or not we want to normalize the spin texture values on a log scale. For example, let's plot results on a linear scale. By default, the plot will normalize the values logarithmically.
+
+.. code-block :: python
+    
+    plt.rcParams.update(ppy.plot_tools.plotparams)
+    diam_spins.kpt.add_labels(ppy.lattice.points_fcc)
+
+    fig, ax  = plt.subplots()
+    diam_spins.plot_spins(ax, log=False)
+    plt.show()
+
+.. image:: figures/diam_spins_linear.png
+    :width: 600
+    :align: center
+
+Note that k-point labels can be removed from the plot by setting the ``show_kpoint_labels`` input to False.
diff --git a/setup.py b/setup.py
index afd9ab0c..47e0195a 100644
--- a/setup.py
+++ b/setup.py
@@ -5,7 +5,7 @@
     
 setup(
     name='perturbopy',
-    version='0.6.5',
+    version='0.6.9',
     description="Suite of Python scripts for Perturbo testing and postprocessing",
     long_description=long_description,
     long_description_content_type="text/markdown",
diff --git a/src/perturbopy/postproc/__init__.py b/src/perturbopy/postproc/__init__.py
index f610abb1..cf7cad88 100644
--- a/src/perturbopy/postproc/__init__.py
+++ b/src/perturbopy/postproc/__init__.py
@@ -3,6 +3,7 @@
 """
 from .calc_modes.calc_mode import CalcMode
 from .calc_modes.bands import Bands
+from .calc_modes.spectral_cumulant import SpectralCumulant
 from .calc_modes.spins import Spins
 from .calc_modes.phdisp import Phdisp
 from .calc_modes.ephmat import Ephmat
diff --git a/src/perturbopy/postproc/calc_modes/calc_mode.py b/src/perturbopy/postproc/calc_modes/calc_mode.py
index 63df9bdd..f4ad9170 100644
--- a/src/perturbopy/postproc/calc_modes/calc_mode.py
+++ b/src/perturbopy/postproc/calc_modes/calc_mode.py
@@ -45,8 +45,9 @@ def __init__(self, pert_dict):
 
         """
 
-        # Extract calculation mode name from pert_dict
+        # Extract calculation mode name and prefix from pert_dict
         self.calc_mode = pert_dict['input parameters']['after conversion'].pop('calc_mode')
+        self.prefix = pert_dict['input parameters']['after conversion'].pop('prefix')
 
         # Extract basic data from pert_dict
         self.alat = pert_dict['basic data']['alat']
diff --git a/src/perturbopy/postproc/calc_modes/imsigma_spin.py b/src/perturbopy/postproc/calc_modes/imsigma_spin.py
index a7c66a70..0dc332f6 100644
--- a/src/perturbopy/postproc/calc_modes/imsigma_spin.py
+++ b/src/perturbopy/postproc/calc_modes/imsigma_spin.py
@@ -28,12 +28,12 @@ class ImsigmaSpin(CalcMode):
         give the configuration number, and the values are floats giving the
         chemical potential (with units chem_pot.units)
     
-    imsigma_flip : UnitsDict
+    imsigma : UnitsDict
         Dictionary of spin flip imaginary self-energies computed for each configuration. The top level keys are the
         configuration number, and the second level keys are the band index. The values are arrays of length N giving the
         imaginary self-energies along all the k-points at that band index for the configuration. Units are in imsigma.units.
 
-    imsigma_flip_mode : UnitsDict
+    imsigma_mode : UnitsDict
         Dictionary of spin flip imaginary self-energies resolved by phonon mode computed. The top level keys are the
         configuration number, and the second level keys are the band index. The third level keys are
         the phonon mode. Finally,the values are arrays of length N giving the imaginary self-energies along all the k-points
@@ -72,22 +72,22 @@ def __init__(self, pert_dict):
 
         self.temper = UnitsDict(units=self._pert_dict['imsigma_spin'].pop('temperature units'))
         self.chem_pot = UnitsDict(units=self._pert_dict['imsigma_spin'].pop('chemical potential units'))
-        self.imsigma_flip = UnitsDict(units=self._pert_dict['imsigma_spin'].pop('Im(Sigma) units'))
-        self.imsigma_flip_mode = UnitsDict(units=self.imsigma_flip.units)
+        self.imsigma = UnitsDict(units=self._pert_dict['imsigma_spin'].pop('Im(Sigma) units'))
+        self.imsigma_mode = UnitsDict(units=self.imsigma.units)
         
         for config_idx in config_dat.keys():
-            self.imsigma_flip[config_idx] = {}
-            self.imsigma_flip_mode[config_idx] = {}
+            self.imsigma[config_idx] = {}
+            self.imsigma_mode[config_idx] = {}
             self.temper[config_idx] = config_dat[config_idx].pop('temperature')
             self.chem_pot[config_idx] = config_dat[config_idx].pop('chemical potential')
             
             imsigma_dat = config_dat[config_idx].pop('band index')
 
             for mode in np.arange(1, num_modes + 1):
-                self.imsigma_flip_mode[config_idx][mode] = {}
+                self.imsigma_mode[config_idx][mode] = {}
 
             for band_index in imsigma_dat.keys():
-                self.imsigma_flip[config_idx][band_index] = np.array(imsigma_dat[band_index]['Im(Sigma)']['total'])
+                self.imsigma[config_idx][band_index] = np.array(imsigma_dat[band_index]['Im(Sigma)']['total'])
 
                 for mode in np.arange(1, num_modes + 1):
-                    self.imsigma_flip_mode[config_idx][mode][band_index] = np.array(imsigma_dat[band_index]['Im(Sigma)']['phonon mode (total)'][mode])
+                    self.imsigma_mode[config_idx][mode][band_index] = np.array(imsigma_dat[band_index]['Im(Sigma)']['phonon mode'][mode])
diff --git a/src/perturbopy/postproc/calc_modes/spectral_cumulant.py b/src/perturbopy/postproc/calc_modes/spectral_cumulant.py
new file mode 100644
index 00000000..1d3e7ad4
--- /dev/null
+++ b/src/perturbopy/postproc/calc_modes/spectral_cumulant.py
@@ -0,0 +1,132 @@
+import numpy as np
+import h5py as h5
+import matplotlib.pyplot as plt
+from perturbopy.postproc.calc_modes.calc_mode import CalcMode
+from perturbopy.io_utils.io import open_yaml, open_hdf5, close_hdf5
+import os
+
+
+class SpectralCumulant(CalcMode):
+    """
+    Class representation of a Perturbo cumulant calculation
+
+    Parameters
+    ----------
+    prefix : str
+        Prefix for the HDF5 file containing spectral function data.
+
+    Attributes
+    ----------
+    temp_array : numpy.ndarray
+        Array of temperatures corresponding to the spectral data.
+    freq_array : numpy.ndarray
+        Array of energy value in electron volts (eV).
+    freq_step : float
+        Energy step size for the energy grid in eV.
+    Akw : numpy.ndarray
+        Spectral function data, indexed by k-point, band, temperature, and ω.
+    """
+    def __init__(self, spectral_file, pert_dict):
+        """
+        Constructor method
+
+        Parameters
+        ----------
+        spectral_file : dict
+            Dictionary for the prefix_spectral_cumulant.h5
+        pert_dict : dict
+            Dictionary containing the inputs from the spectral-cum calculation.
+
+        """
+        super().__init__(pert_dict)
+        if self.calc_mode != 'spectral-cum':
+            raise ValueError('Calculation mode for a SpectralCumulantCalcMode object should be "spectral-cum"')
+
+        Akw = spectral_file['spectral_functions']
+        w_lower = np.asarray(spectral_file['w_lower_index'])
+        w_upper = np.asarray(spectral_file['w_upper_index'])
+        freq_step = np.asarray(spectral_file['wfreq_step_eV'])
+        self.temp_array = np.asanyarray(spectral_file['temperatures'])
+        self.freq_array = np.arange(w_lower, w_upper + 1) * freq_step
+        self.freq_step = freq_step
+        Akw_np = []
+        for key in Akw.keys():
+            Akw_np.append(np.asarray(Akw[key]))
+        close_hdf5(spectral_file)
+
+        self.Akw = np.asarray(Akw_np)
+
+    @classmethod
+    def from_hdf5_yaml(cls, spectral_path, yaml_path='pert_output.yml'):
+        """
+        Class method to create a SpectralCumulantCalcMode object from the HDF5 file and YAML file
+        generated by a Perturbo calculation
+
+        Parameters
+        ----------
+        spectral_path : str
+           Path to the HDF5 file generated by a spectral-cum calculation
+        yaml_path : str, optional
+           Path to the YAML file generated by a spectral-cum calculation
+
+        Returns
+        -------
+
+        """
+
+        if not os.path.isfile(spectral_path):
+            raise FileNotFoundError(f'File {spectral_path} not found')
+        if not os.path.isfile(yaml_path):
+            raise FileNotFoundError(f'File {yaml_path} not found')
+
+        spectral_file = open_hdf5(spectral_path)
+        yaml_dict = open_yaml(yaml_path)
+
+        return cls(spectral_file, yaml_dict)
+
+    def plot_Aw(self, ax, ik=0, it=0, ib=0):
+        """
+        Plots the spectral function A(ω) for a given k-point, temperature, and band.
+
+        Parameters
+        ----------
+        ax : matplotlib.axes.Axes
+            Axis on which to plot the spectral function.
+        ik : int, optional
+            Index of the k-point in the grid. Default is 0.
+        it : int, optional
+            Index of the temperature in the temperature array. Default is 0.
+        ib : int, optional
+            Index of the band. Default is 0.
+
+        Returns
+        -------
+        matplotlib.axes.Axes
+            The axis object with the plotted spectral function.
+
+        Notes
+        -----
+        The spectral function is normalized before plotting. The plot includes
+        labels, legends, and adjusted aesthetics for better visualization.
+        """
+        if it > len(self.temp_array):
+            raise ValueError('Temperature index is out of range')
+        if ib > self.Akw[0].shape[0]:
+            raise ValueError('Band index is out of range')
+        if ik > len(self.Akw):
+            raise ValueError('k-point index is out of range')
+        A0w = self.Akw[ik][ib, it, :]
+        freq_step = self.freq_step
+        freq_array = self.freq_array
+        # normalize
+        A0w = A0w / (np.sum(A0w) * freq_step)
+        # plot
+        ax.plot(freq_array, A0w, lw=2, label=f'T={ int(self.temp_array[it])} K')
+        ax.legend(fontsize=18)
+        plt.ylim([0, np.max(A0w) * 1.1])
+        plt.xlabel(r'$\omega-\epsilon_{nk}$ (eV)', fontsize=20)
+        plt.ylabel(r'A($\omega$) (eV$^{-1}$)', fontsize=20)
+        plt.xticks(fontsize=20)
+        plt.yticks(fontsize=20)
+        plt.tight_layout()
+        return ax
diff --git a/src/perturbopy/postproc/calc_modes/spins.py b/src/perturbopy/postproc/calc_modes/spins.py
index b73faba9..5e538638 100644
--- a/src/perturbopy/postproc/calc_modes/spins.py
+++ b/src/perturbopy/postproc/calc_modes/spins.py
@@ -45,7 +45,7 @@ def __init__(self, pert_dict):
         energy_units = self._pert_dict['spins'].pop('band units')
 
         spins_dict = self._pert_dict['spins'].pop('band index (spins)')
-        spin_units = self._pert_dict['spins'].pop('<n|Sigma_z|n> units')
+        spin_units = self._pert_dict['spins'].pop('<n|sigma_z|n> units')
 
         self.kpt = RecipPtDB.from_lattice(kpoint, kpoint_units, self.lat, self.recip_lat, kpath, kpath_units)
         self.bands = UnitsDict.from_dict(energies_dict, energy_units)
@@ -53,7 +53,7 @@ def __init__(self, pert_dict):
 
     def plot_spins(self, ax, kpoint_idx=0, show_kpoint_labels=True, log=True, **kwargs):
         """
-        Method to plot the <n|Sigma_z|n> values over the band structure.
+        Method to plot the <n|sigma_z|n> values over the band structure.
 
         Parameters
         ----------
@@ -61,7 +61,7 @@ def plot_spins(self, ax, kpoint_idx=0, show_kpoint_labels=True, log=True, **kwar
            Axis on which to plot the bands.
         
         kpoint_idx : int, optional
-           Index of the k-point to plot the <n|Sigma_z|n> values for. <n|Sigma_z|n> elements will be plotted along k-points, at this k-point
+           Index of the k-point to plot the <n|sigma_z|n> values for. <n|sigma_z|n> elements will be plotted along k-points, at this k-point
            By default, it will be the first k-point.
 
         show_kpoint_labels : bool, optional
@@ -81,7 +81,7 @@ def plot_spins(self, ax, kpoint_idx=0, show_kpoint_labels=True, log=True, **kwar
 
         values = {key: 1 - val for key, val in self.spins.items()}
 
-        ax = plot_vals_on_bands(ax, self.kpt.path, self.bands, self.bands.units, values=values, log=log, label=r'$|<n|\Sigma_z|n>|$', **kwargs)
+        ax = plot_vals_on_bands(ax, self.kpt.path, self.bands, self.bands.units, values=values, log=log, label=r'$|\langle n|\sigma_z|n\rangle |$', **kwargs)
 
         if show_kpoint_labels:
             ax = plot_recip_pt_labels(ax, self.kpt.labels, self.kpt.points, self.kpt.path)
diff --git a/src/perturbopy/yml_files/input_template.yml b/src/perturbopy/yml_files/input_template.yml
index 6e3cb4aa..b2ed10e3 100644
--- a/src/perturbopy/yml_files/input_template.yml
+++ b/src/perturbopy/yml_files/input_template.yml
@@ -25,6 +25,8 @@ qe2pert:
       polar_alpha: convergence parameter used in the Ewald sum, only used for polar correction. 
       thickness_2d: for 2D polar e-ph correction. Only needed when system_2d=.true.
       eig_corr: usually used when one wants to use modified eigenvalues (e.g., from GW)
+      tdep: when =.true., read force constants from TDEP outputs, infile.forceconstant and infile.ucposcar
+      symmborn: Symmetrize Born effective charges read from TDEP files for polar materials. Only needed when tdep = .true. and the material is polar.
 
 
 bands:
diff --git a/tests/refs/gaas_spins.yml b/tests/refs/gaas_spins.yml
index e91964ce..1e5d2314 100644
--- a/tests/refs/gaas_spins.yml
+++ b/tests/refs/gaas_spins.yml
@@ -2163,7 +2163,7 @@ spins:
          -   13.0429498274
          -   13.0081597319
 
-   <n|Sigma_z|n> units: arbitrary
+   <n|sigma_z|n> units: arbitrary
    band index (spins):
 
       1:
diff --git a/tests/refs/sto_spectral-cum.yml b/tests/refs/sto_spectral-cum.yml
new file mode 100644
index 00000000..1a7191d7
--- /dev/null
+++ b/tests/refs/sto_spectral-cum.yml
@@ -0,0 +1,405 @@
+---
+program: perturbo
+
+start date and time:
+   date:  6 December 2024
+   time: "14:23:16"
+
+parallelization:
+   type: MPI+OpenMP
+   number of MPI tasks:   1
+   number of OpenMP threads:  12
+   number of nodes:   1
+
+input parameters:
+   before conversion (before reading from epr file):
+      band_max:                   3
+      band_min:                   1
+      boltz_acc_thr:           0.0000000000000000E+00
+      boltz_de:           1.0000000000000000E+00
+      boltz_efield:          [ 0.0000000000000000E+00,   0.0000000000000000E+00,   0.0000000000000000E+00,  ]
+      boltz_emax:           9.9990000000000000E+03
+      boltz_emin:          -9.9990000000000000E+03
+      boltz_init_dist:          
+      boltz_init_e0:          -9.9990000000000000E+03
+      boltz_init_smear:           2.0000000000000000E+01
+      boltz_kdim:          [         1,           1,           1,  ]
+      boltz_norm_dist:          False
+      boltz_nstep:                   0
+      boltz_nstep_min:                   1
+      boltz_qdim:          [         1,           1,           1,  ]
+      calc_mode:          spectral-cum
+      cauchy_scale:           5.0000000000000003E-02
+      scat_impl:          std
+      cum_de:           4.0000000000000002E-01
+      cum_inner_emax:           2.0000000000000001E-01
+      cum_inner_emin:          -2.0000000000000001E-01
+      cum_outer_emax:           5.9999999999999998E-01
+      cum_outer_emin:          -5.9999999999999998E-01
+      cum_outer_np:                   4
+      debug:          False
+      delta_smear:           1.0000000000000000E+01
+      find_efermi:          False
+      fklist:          band.kpt
+      fqlist:          
+      ftemper:          sto.temper
+      full_ite:          True 
+      hole:          False
+      lmagsym:          False
+      load_scatter_eph:          False
+      nsamples:              100000
+      output_nstep:                   1
+      phfreq_cutoff:           1.0000000000000000E+00
+      polar_split:          
+      prefix:          sto
+      sampling:          uniform
+      solver:          rk4
+      spectral_emax:           5.9999999999999998E-01
+      spectral_emin:          -5.9999999999999998E-01
+      spectral_np:                   1
+      time_step:           1.0000000000000000E+00
+      tmp_dir:          ./
+      trans_thr:           2.0000000000000000E-03
+      scat_impl:          std
+      use_mem:          True 
+      yaml_fname:          pert_output.yml
+   after conversion:
+      band_max:                   3
+      band_min:                   1
+      boltz_acc_thr:           0.0000000000000000E+00
+      boltz_de:           7.3498617648950546E-05
+      boltz_efield:          [ 0.0000000000000000E+00,   0.0000000000000000E+00,   0.0000000000000000E+00,  ]
+      boltz_emax:           7.3491267787185654E+02
+      boltz_emin:          -7.3491267787185654E+02
+      boltz_init_dist:          
+      boltz_init_e0:          -7.3491267787185654E+02
+      boltz_init_smear:           1.4699723529790110E-03
+      boltz_kdim:          [         1,           1,           1,  ]
+      boltz_norm_dist:          False
+      boltz_nstep:                   0
+      boltz_nstep_min:                   1
+      boltz_qdim:          [         1,           1,           1,  ]
+      calc_mode:          spectral-cum
+      cauchy_scale:           5.0000000000000003E-02
+      scat_impl:          std
+      cum_de:           2.9399447059580220E-05
+      cum_inner_emax:           1.4699723529790111E-02
+      cum_inner_emin:          -1.4699723529790111E-02
+      cum_outer_emax:           4.4099170589370330E-02
+      cum_outer_emin:          -4.4099170589370330E-02
+      cum_outer_np:                   4
+      debug:          False
+      delta_smear:           7.3498617648950549E-04
+      find_efermi:          False
+      fklist:          band.kpt
+      fqlist:          
+      ftemper:          sto.temper
+      full_ite:          True 
+      hole:          False
+      lmagsym:          False
+      load_scatter_eph:          False
+      nsamples:              100000
+      output_nstep:                   1
+      phfreq_cutoff:           7.3498617648950546E-05
+      polar_split:          
+      prefix:          sto
+      sampling:          uniform
+      solver:          rk4
+      spectral_emax:           4.4099170589370330E-02
+      spectral_emin:          -4.4099170589370330E-02
+      spectral_np:                   1
+      time_step:           2.0670686467503085E+01
+      tmp_dir:          ./
+      trans_thr:           2.0000000000000000E-03
+      scat_impl:          std
+      use_mem:          True 
+      yaml_fname:          sto_spectral-cum.yml
+
+
+basic data:
+   alat:           7.3699300000
+   alat units:       bohr
+   # a vector is a row
+   lattice vectors:
+         - [   1.00000,     0.00000,     0.00000,  ]
+         - [   0.00000,     1.00000,     0.00000,  ]
+         - [   0.00000,     0.00000,     1.00000,  ]
+   lattice vectors units:       alat
+   reciprocal lattice vectors:
+   # a vector is a row
+         - [   1.00000,     0.00000,     0.00000,  ]
+         - [   0.00000,     1.00000,     0.00000,  ]
+         - [   0.00000,     0.00000,     1.00000,  ]
+   reciprocal lattice vectors units:       2pi/alat
+   number of atoms in unit cell:                5
+   atomic positions:
+   # a vector is a row
+         - [   0.00000,     0.00000,     0.00000,  ]
+         - [   0.50000,     0.50000,     0.50000,  ]
+         - [   0.00000,     0.50000,     0.50000,  ]
+         - [   0.50000,     0.00000,     0.50000,  ]
+         - [   0.50000,     0.50000,     0.00000,  ]
+   atomic positions units:       alat
+   volume:         400.3041465593
+   volume units:       bohr^3
+   number of symmetry operations:               48
+   symop:
+      1:
+         - [    1,      0,      0,  ]
+         - [    0,      1,      0,  ]
+         - [    0,      0,      1,  ]
+      2:
+         - [   -1,      0,      0,  ]
+         - [    0,     -1,      0,  ]
+         - [    0,      0,      1,  ]
+      3:
+         - [   -1,      0,      0,  ]
+         - [    0,      1,      0,  ]
+         - [    0,      0,     -1,  ]
+      4:
+         - [    1,      0,      0,  ]
+         - [    0,     -1,      0,  ]
+         - [    0,      0,     -1,  ]
+      5:
+         - [    0,      1,      0,  ]
+         - [    1,      0,      0,  ]
+         - [    0,      0,     -1,  ]
+      6:
+         - [    0,     -1,      0,  ]
+         - [   -1,      0,      0,  ]
+         - [    0,      0,     -1,  ]
+      7:
+         - [    0,      1,      0,  ]
+         - [   -1,      0,      0,  ]
+         - [    0,      0,      1,  ]
+      8:
+         - [    0,     -1,      0,  ]
+         - [    1,      0,      0,  ]
+         - [    0,      0,      1,  ]
+      9:
+         - [    0,      0,      1,  ]
+         - [    0,     -1,      0,  ]
+         - [    1,      0,      0,  ]
+      10:
+         - [    0,      0,     -1,  ]
+         - [    0,     -1,      0,  ]
+         - [   -1,      0,      0,  ]
+      11:
+         - [    0,      0,      1,  ]
+         - [    0,      1,      0,  ]
+         - [   -1,      0,      0,  ]
+      12:
+         - [    0,      0,     -1,  ]
+         - [    0,      1,      0,  ]
+         - [    1,      0,      0,  ]
+      13:
+         - [   -1,      0,      0,  ]
+         - [    0,      0,      1,  ]
+         - [    0,      1,      0,  ]
+      14:
+         - [   -1,      0,      0,  ]
+         - [    0,      0,     -1,  ]
+         - [    0,     -1,      0,  ]
+      15:
+         - [    1,      0,      0,  ]
+         - [    0,      0,      1,  ]
+         - [    0,     -1,      0,  ]
+      16:
+         - [    1,      0,      0,  ]
+         - [    0,      0,     -1,  ]
+         - [    0,      1,      0,  ]
+      17:
+         - [    0,      1,      0,  ]
+         - [    0,      0,      1,  ]
+         - [    1,      0,      0,  ]
+      18:
+         - [    0,     -1,      0,  ]
+         - [    0,      0,      1,  ]
+         - [   -1,      0,      0,  ]
+      19:
+         - [    0,      1,      0,  ]
+         - [    0,      0,     -1,  ]
+         - [   -1,      0,      0,  ]
+      20:
+         - [    0,     -1,      0,  ]
+         - [    0,      0,     -1,  ]
+         - [    1,      0,      0,  ]
+      21:
+         - [    0,      0,      1,  ]
+         - [    1,      0,      0,  ]
+         - [    0,      1,      0,  ]
+      22:
+         - [    0,      0,      1,  ]
+         - [   -1,      0,      0,  ]
+         - [    0,     -1,      0,  ]
+      23:
+         - [    0,      0,     -1,  ]
+         - [   -1,      0,      0,  ]
+         - [    0,      1,      0,  ]
+      24:
+         - [    0,      0,     -1,  ]
+         - [    1,      0,      0,  ]
+         - [    0,     -1,      0,  ]
+      25:
+         - [   -1,      0,      0,  ]
+         - [    0,     -1,      0,  ]
+         - [    0,      0,     -1,  ]
+      26:
+         - [    1,      0,      0,  ]
+         - [    0,      1,      0,  ]
+         - [    0,      0,     -1,  ]
+      27:
+         - [    1,      0,      0,  ]
+         - [    0,     -1,      0,  ]
+         - [    0,      0,      1,  ]
+      28:
+         - [   -1,      0,      0,  ]
+         - [    0,      1,      0,  ]
+         - [    0,      0,      1,  ]
+      29:
+         - [    0,     -1,      0,  ]
+         - [   -1,      0,      0,  ]
+         - [    0,      0,      1,  ]
+      30:
+         - [    0,      1,      0,  ]
+         - [    1,      0,      0,  ]
+         - [    0,      0,      1,  ]
+      31:
+         - [    0,     -1,      0,  ]
+         - [    1,      0,      0,  ]
+         - [    0,      0,     -1,  ]
+      32:
+         - [    0,      1,      0,  ]
+         - [   -1,      0,      0,  ]
+         - [    0,      0,     -1,  ]
+      33:
+         - [    0,      0,     -1,  ]
+         - [    0,      1,      0,  ]
+         - [   -1,      0,      0,  ]
+      34:
+         - [    0,      0,      1,  ]
+         - [    0,      1,      0,  ]
+         - [    1,      0,      0,  ]
+      35:
+         - [    0,      0,     -1,  ]
+         - [    0,     -1,      0,  ]
+         - [    1,      0,      0,  ]
+      36:
+         - [    0,      0,      1,  ]
+         - [    0,     -1,      0,  ]
+         - [   -1,      0,      0,  ]
+      37:
+         - [    1,      0,      0,  ]
+         - [    0,      0,     -1,  ]
+         - [    0,     -1,      0,  ]
+      38:
+         - [    1,      0,      0,  ]
+         - [    0,      0,      1,  ]
+         - [    0,      1,      0,  ]
+      39:
+         - [   -1,      0,      0,  ]
+         - [    0,      0,     -1,  ]
+         - [    0,      1,      0,  ]
+      40:
+         - [   -1,      0,      0,  ]
+         - [    0,      0,      1,  ]
+         - [    0,     -1,      0,  ]
+      41:
+         - [    0,     -1,      0,  ]
+         - [    0,      0,     -1,  ]
+         - [   -1,      0,      0,  ]
+      42:
+         - [    0,      1,      0,  ]
+         - [    0,      0,     -1,  ]
+         - [    1,      0,      0,  ]
+      43:
+         - [    0,     -1,      0,  ]
+         - [    0,      0,      1,  ]
+         - [    1,      0,      0,  ]
+      44:
+         - [    0,      1,      0,  ]
+         - [    0,      0,      1,  ]
+         - [   -1,      0,      0,  ]
+      45:
+         - [    0,      0,     -1,  ]
+         - [   -1,      0,      0,  ]
+         - [    0,     -1,      0,  ]
+      46:
+         - [    0,      0,     -1,  ]
+         - [    1,      0,      0,  ]
+         - [    0,      1,      0,  ]
+      47:
+         - [    0,      0,      1,  ]
+         - [    1,      0,      0,  ]
+         - [    0,     -1,      0,  ]
+      48:
+         - [    0,      0,      1,  ]
+         - [   -1,      0,      0,  ]
+         - [    0,      1,      0,  ]
+   kc dimensions:
+         - [      4,        4,        4,  ]
+   spinor:       False
+   polar_alpha:           1.0000000000
+   epsil:
+   # a vector is a row
+         - [      6.25856984,        0.00000000,        0.00000000,  ]
+         - [      0.00000000,        6.25856984,        0.00000000,  ]
+         - [      0.00000000,        0.00000000,        6.25856984,  ]
+   qc dimensions:
+         - [      2,        2,        2,  ]
+   mass:
+      -       79860.74458118
+      -       43628.10158488
+      -       14582.19644550
+      -       14582.19644550
+      -       14582.19644550
+   mass units:       atomic
+   zstar:
+      1:
+         - [   2.55389,     0.00000,     0.00000,  ]
+         - [   0.00000,     2.55389,     0.00000,  ]
+         - [   0.00000,     0.00000,     2.55389,  ]
+      2:
+         - [   7.28719,     0.00000,     0.00000,  ]
+         - [   0.00000,     7.28719,     0.00000,  ]
+         - [   0.00000,     0.00000,     7.28719,  ]
+      3:
+         - [  -5.77445,     0.00000,     0.00000,  ]
+         - [   0.00000,    -2.03331,     0.00000,  ]
+         - [   0.00000,     0.00000,    -2.03331,  ]
+      4:
+         - [  -2.03331,     0.00000,     0.00000,  ]
+         - [   0.00000,    -5.77445,     0.00000,  ]
+         - [   0.00000,     0.00000,    -2.03331,  ]
+      5:
+         - [  -2.03331,     0.00000,     0.00000,  ]
+         - [   0.00000,    -2.03331,     0.00000,  ]
+         - [   0.00000,     0.00000,    -5.77445,  ]
+   system_2d:       False
+   number of Wannier functions:                3
+   wannier_center:
+   # a vector is a row
+         - [   0.50000,     0.50000,     0.50000,  ]
+         - [   0.50000,     0.50000,     0.50000,  ]
+         - [   0.50000,     0.50000,     0.50000,  ]
+   wannier_center units: cartesian
+   wannier_center_cryst:
+   # a vector is a row
+         - [   0.50000,     0.50000,     0.50000,  ]
+         - [   0.50000,     0.50000,     0.50000,  ]
+         - [   0.50000,     0.50000,     0.50000,  ]
+   wannier_center_cryst units: crystal
+
+memory:
+   physical: 65737300 kB                                                                     
+   max rss: 54812 kB                                                                        
+
+timings:
+
+   time units: s
+
+   PERTURBO:
+      cpu:         3.10515600
+      wall:        0.35491920
+      calls:         1
+...
diff --git a/tests/refs/sto_spectral_cumulant.h5 b/tests/refs/sto_spectral_cumulant.h5
new file mode 100644
index 00000000..850ef717
Binary files /dev/null and b/tests/refs/sto_spectral_cumulant.h5 differ
diff --git a/tests/test_spectral_cum.py b/tests/test_spectral_cum.py
new file mode 100644
index 00000000..8d1af09e
--- /dev/null
+++ b/tests/test_spectral_cum.py
@@ -0,0 +1,55 @@
+import numpy as np
+import pytest
+import os
+
+import perturbopy.postproc as ppy
+
+
+@pytest.fixture()
+def sto_spectral_cum():
+    """
+    Method to generate the spectral_cum object corresponding to spectralhdf5.
+
+    Returns
+    -------
+    phdisp : ppy.PhdispCalcMode
+
+    """
+    yml_path = os.path.join("refs", "sto_spectral-cum.yml")
+    spectral_path = os.path.join("refs", "sto_spectral_cumulant.h5")
+    return ppy.SpectralCumulant.from_hdf5_yaml(spectral_path, yml_path)
+
+def test_plot_spectral_cum(sto_spectral_cum, plt, with_plt):
+    """
+    Method to test SpectralCumulant.plot_Aw function
+
+    Method to plot the spectral function A(ω)
+
+    Parameters
+    ----------
+    show_qpoint_labels : bool, optional
+        If true, the q-point labels stored in the labels attribute will be shown on the plot. Default true.
+
+    """
+
+    if not with_plt:
+        pytest.skip("Test requires pytest-plt")
+
+    fig, ax = plt.subplots()
+    ppy.SpectralCumulant.plot_Aw(sto_spectral_cum, ax)
+
+def test_spectral_cum():
+    """
+    Method to test SpectralCumulant.plot_Aw function array shape
+
+    Parameters
+    ----------
+
+    """
+
+    yml_path = os.path.join("refs", "sto_spectral-cum.yml")
+    spectral_path = os.path.join("refs", "sto_spectral_cumulant.h5")
+    model = ppy.SpectralCumulant.from_hdf5_yaml(spectral_path, yml_path)
+    np.testing.assert_equal(model.Akw.shape, (1, 3, 2, 3001))
+    np.testing.assert_equal(model.freq_array.shape, (3001,))
+