diff --git a/docs/examples/Magnetized_OCP/Magnetized_Plasma.ipynb b/docs/examples/Magnetized_OCP/Magnetized_Plasma.ipynb
index ee694ab2..b23d0e8a 100644
--- a/docs/examples/Magnetized_OCP/Magnetized_Plasma.ipynb
+++ b/docs/examples/Magnetized_OCP/Magnetized_Plasma.ipynb
@@ -441,46 +441,29 @@
     "\tsnapshot interval time = 5.5000e-17 [s] = 9.2153e-02 w_p = 1.4667e-02 plasma periods\n",
     "\tTotal number of snapshots = 14000\n",
     "\n",
-    "Bernu ran a simulation for 244 plasma periods. We decided to run for a time 5 times longer. This is because we want to show the time slicing option in the Postprocessing phase.\n",
-    "\n",
-    "The estimated times, on this machine, are \n",
-    "\n",
-    "    ----------------------- Total Estimated Times ------------------------ \n",
-    "\n",
-    "\n",
-    "    Equilibration Time: 0 hrs 2 min 20 sec\n",
-    "\n",
-    "    Magnetization Time: 0 sec 0 msec 0 usec 0 nsec\n",
-    "\n",
-    "    Production Time: 0 hrs 33 min 11 sec\n",
-    "\n",
-    "    Total Run Time: 0 hrs 35 min 31 sec\n",
-    "\n",
-    "and the entire run requires\n",
-    "\n",
-    "    Total minimum required space: 1 GB 43 MB 1009 KB 192 bytes\n",
-    "\n",
-    "We are ready to start the simulation."
+    "Bernu ran a simulation for 244 plasma periods. We decided to run for a time 5 times longer. This is because we want to show the time slicing option in the Postprocessing phase."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 1,
    "id": "d585f2ca",
-   "metadata": {},
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2024-03-22T20:31:15.612765Z",
+     "iopub.status.busy": "2024-03-22T20:31:15.611857Z",
+     "iopub.status.idle": "2024-03-22T20:31:15.632210Z",
+     "shell.execute_reply": "2024-03-22T20:31:15.630891Z",
+     "shell.execute_reply.started": "2024-03-22T20:31:15.612697Z"
+    }
+   },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "{'platform': 'Darwin',\n",
-       " 'platform-release': '22.6.0',\n",
-       " 'architecture': 'arm64',\n",
-       " 'processor': 'arm',\n",
-       " 'ram': '32 GB'}"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Running on an arm processor with 32 GB GB of RAM\n"
+     ]
     }
    ],
    "source": [
@@ -493,7 +476,32 @@
     "info['processor']=platform.processor()\n",
     "info['ram']=str(round(psutil.virtual_memory().total / (1024.0 **3)))+\" GB\"\n",
     "\n",
-    "display(info)\n"
+    "print(f\"Running on an {info['processor']} processor with {info['ram']} GB of RAM\")\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "efabcbd6-07c6-43fc-8bf0-c2bd7dc34549",
+   "metadata": {},
+   "source": [
+    "The estimated times, on this machine, are \n",
+    "\n",
+    "    ----------------------- Total Estimated Times ------------------------ \n",
+    "\n",
+    "\n",
+    "    Equilibration Time: 0 hrs 2 min 20 sec\n",
+    "\n",
+    "    Magnetization Time: 0 sec 0 msec 0 usec 0 nsec\n",
+    "\n",
+    "    Production Time: 0 hrs 33 min 11 sec\n",
+    "\n",
+    "    Total Run Time: 0 hrs 35 min 31 sec\n",
+    "\n",
+    "and the entire run requires\n",
+    "\n",
+    "    Total minimum required space: 1 GB 43 MB 1009 KB 192 bytes\n",
+    "\n",
+    "We are ready to start the simulation."
    ]
   },
   {
@@ -1336,7 +1344,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
diff --git a/sarkas/core.py b/sarkas/core.py
index 2532d853..5d5c2ce5 100644
--- a/sarkas/core.py
+++ b/sarkas/core.py
@@ -769,6 +769,11 @@ def set_species_attributes(self, species: list):
             List of :class:`sarkas.plasma.Species`.
 
         """
+
+        # DEV NOTE: This method could be defined in the Particles class, however, it is more convenient to have it here because it uses the correct unit system. 
+        # If I were to move it to the Particles class, I would need to pass the units as an argument. 
+        # TODO: It could be moved in the future. 
+
         # Loop over species and assign missing attributes
         # Collect species properties in single arrays
 
diff --git a/sarkas/particles.py b/sarkas/particles.py
index 4e3042bd..980d3a3f 100644
--- a/sarkas/particles.py
+++ b/sarkas/particles.py
@@ -6,7 +6,7 @@
 from copy import deepcopy
 from h5py import File as h5File
 from numba import float64, int64, jit, njit, void
-from numpy import arange, array, empty, floor, int64
+from numpy import arange, array, empty, floor, full, int64
 from numpy import load as np_load
 from numpy import (
     loadtxt,
@@ -536,6 +536,9 @@ def initialize_arrays(self):
 
         self.species_initial_velocity = zeros((self.num_species, 3))
         self.species_thermal_velocity = zeros((self.num_species, 3))
+        
+        self.species_initial_spatial_distribution = empty((self.num_species, 3), dtype=str)
+        self.species_initial_velocity_distribution = empty((self.num_species, 3), dtype=str)
 
         self.species_kinetic_energy = zeros(self.num_species)
         self.species_potential_energy = zeros(self.num_species)
@@ -569,25 +572,25 @@ def initialize_arrays(self):
             self.heat_flux_species_tensor = zeros((3, self.num_species, self.num_species))
             self.species_heat_flux = zeros((self.num_species, 3))
 
-    def initialize_positions(self):
+    def initialize_positions(self, species: list = None):
         """
         Initialize particles' positions based on the load method.
         """
 
-        # position distribution.
+        # TODO: Break this method into smaller methods. It is too long.
         if self.load_method == "lattice":
             self.lattice(self.load_perturb)
 
         elif self.load_method == "random_reject":
             # check
             if not hasattr(self, "load_rejection_radius"):
-                raise AttributeError("Rejection radius not defined. " "Please define Parameters.load_rejection_radius.")
+                raise AttributeError("Rejection radius not defined. Please define Parameters.load_rejection_radius.")
             self.random_reject(self.load_rejection_radius)
 
         elif self.load_method == "halton_reject":
             # check
             if not hasattr(self, "load_rejection_radius"):
-                raise AttributeError("Rejection radius not defined. " "Please define Parameters.load_rejection_radius.")
+                raise AttributeError("Rejection radius not defined. Please define Parameters.load_rejection_radius.")
             self.halton_reject(self.load_halton_bases, self.load_rejection_radius)
 
         elif self.load_method in ["uniform", "random_no_reject"]:
@@ -596,6 +599,10 @@ def initialize_positions(self):
             )
 
         elif self.load_method == "gaussian":
+
+            if self.load_gauss_sigma is None:
+                raise AttributeError("Gaussian sigma not defined. Please define Parameters.load_gauss_sigma.")
+            
             sp_start = 0
             sp_end = 0
             for sp, sp_num in enumerate(self.species_num):
@@ -604,6 +611,7 @@ def initialize_positions(self):
                     self.box_lengths[0] / 2.0, self.load_gauss_sigma[sp], (sp_num, 3)
                 )
                 sp_start += sp_num
+
         elif self.load_method in ["qmc", "quasi_monte_carlo"]:
             # Ensure that we are not making silly typos
             self.qmc_sequence = self.qmc_sequence.lower()
@@ -614,12 +622,81 @@ def initialize_positions(self):
                 )
 
             self.pos = self.quasi_monte_carlo(self.qmc_sequence, self.total_num_ptcls, self.dimensions, self.box_lengths)
+        
+        elif self.load_method == "species_specific":
+            sp_start = 0
+            sp_end = 0
+            for ic, sp in enumerate(species):
+                # TODO: Add the read_from_file option
+                if sp.name != "electron_background":
+                    sp_end += sp.num
+                    if sp.initial_spatial_distribution in ["uniform", "random_no_reject"]:
+                        self.pos[sp_start:sp_end, :] = self.uniform_no_reject(
+                            0.5 * self.box_lengths - 0.5 * self.pbox_lengths, 0.5 * self.box_lengths + 0.5 * self.pbox_lengths
+                        )
+                    elif sp.initial_spatial_distribution == "gaussian":
+                        
+                        if sp.gaussian_sigma is None:
+                            raise AttributeError("Gaussian sigma not defined. Please define Species.gaussian_sigma.")
+                        
+                        if sp.gaussian_mean is None:
+                            raise AttributeError("Gaussian mean not defined. Please define Species.gaussian_mean.")
+
+                        # Set the mean and sigma for the species
+                        if sp.gaussian_mean == "center":
+                            sp_mean = self.box_lengths / 2.0
+                        else:
+                            if isinstance(sp.gaussian_mean, (int, float)):
+                                sp_mean = full(self.dimensions, sp.gaussian_mean)
+                            elif isinstance(sp.gaussian_mean, (list, ndarray)):
+                                # Check that it has length 3
+                                if len(sp.gaussian_mean) != 3 and self.dimensions == 3:
+                                    raise AttributeError("Gaussian mean must be a list or array of length 3.")
+                                elif len(sp.gaussian_mean) != 2 and self.dimensions == 2:
+                                    raise AttributeError("Gaussian mean must be a list or array of length 2.")
+                                else:
+                                    sp_mean = sp.gaussian_mean
+
+                        if isinstance(sp.gaussian_sigma, (int, float)):
+                            sp_sigma = full(self.dimensions, sp.gaussian_sigma)
+                        else:  
+                            # Check that it has the correct dimensions
+                            if len(sp.gaussian_sigma) != 3 and self.dimensions == 3:
+                                raise AttributeError("Gaussian sigma must be a list or array of length 3.")
+                            elif len(sp.gaussian_sigma) != 2 and self.dimensions == 2:
+                                raise AttributeError("Gaussian sigma must be a list or array of length 2.")
+                            else:
+                                sp_sigma = sp.gaussian_sigma
+                        
+                        for dim in range(self.dimensions):
+                            self.pos[sp_start:sp_end, dim] = self.gaussian(sp_mean[dim], sp_sigma[dim], (sp.num, 1))
+                    elif sp.initial_spatial_distribution in ["quasi_monte_carlo", "qmc"]:
+                        if sp.qmc_sequence not in self.available_qmc_sequences:
+                            raise AttributeError(
+                                f"Quasi Monte Carlo sequence not recognized. Please choose from {self.available_qmc_sequences}"
+                            )
+                        self.pos[sp_start:sp_end, :] = self.quasi_monte_carlo(
+                            sp.qmc_sequence, sp.num, self.dimensions, self.box_lengths
+                        )
+                    elif sp.initial_spatial_distribution in ["quasi_monte_carlo_rejection", "qmc_reject"]:
+                        if self.qmc_sequence not in self.available_qmc_sequences:
+                            raise AttributeError(
+                                f"Quasi Monte Carlo sequence not recognized. Please choose from {self.available_qmc_sequences}"
+                            )
+                        self.pos[sp_start:sp_end, :] = self.quasi_monte_carlo_rejection(
+                                sp.qmc_sequence, sp_num, self.dimensions, self.load_rejection_radius, self.box_lengths
+                            )
+                    sp_start += sp.num
         else:
             raise AttributeError("Incorrect particle placement scheme specified.")
 
     @staticmethod
     def quasi_monte_carlo(
-        sequence: str = "sobol", num_ptcls: int = 2**10, dimensions: int = 3, box_lengths: ndarray = None
+        sequence: str = "sobol",
+        num_ptcls: int = 2**10,
+        dimensions: int = 3,
+        box_lengths: ndarray = None,
+        **kwargs
     ):
         """
         Place particles according to a quasi-Monte Carlo sequence.
@@ -642,6 +719,8 @@ def quasi_monte_carlo(
         box_lengths : array-like, optional
             Lengths of the simulation box in each dimension. Default is None.
 
+        kwargs : dict, optional
+            Additional keyword arguments to pass for the initialization of the `scipy.stats.qmc` class.
         Returns
         -------
         array-like
@@ -664,22 +743,75 @@ def quasi_monte_carlo(
 
         """
         if sequence == "sobol":
-            qmc_sampler = qmc.Sobol(dimensions)
+            qmc_sampler = qmc.Sobol(d= dimensions, **kwargs)
         elif sequence == "latinhypercube":
-            qmc_sampler = qmc.LatinHypercube(dimensions)
-        elif sequence == "poissondisk":
-            qmc_sampler = qmc.PoissonDisk(dimensions)
+            qmc_sampler = qmc.LatinHypercube(d =dimensions, **kwargs)
         elif sequence == "halton":
-            qmc_sampler = qmc.Halton(dimensions)
+            qmc_sampler = qmc.Halton(d=dimensions, **kwargs)
         else:
             raise ValueError("Invalid sequence name.")
 
         positions = qmc_sampler.random(num_ptcls)
-        if box_lengths is not None:
-            positions *= box_lengths
+        positions = qmc.scale(positions, l_bounds=[0.0, 0.0, 0.0], u_bounds=box_lengths)
 
         return positions
 
+    @staticmethod
+    def quasi_monte_carlo_rejection(sequence: str = "poissondisk", num_ptcls: int = 2**10, dimensions: int = 3, rejection_radius: float = 0.5,  box_lengths=None, **kwargs):
+        """
+        Place particles according to a quasi-Monte Carlo sequence with rejection.
+
+        Parameters
+        ----------
+        sequence : str, optional
+            Name of the sequence to use. Options: PoissonDisk. Default is "PoissonDisk".
+        
+        num_ptcls : int, optional
+            Number of particles to place. Default is 2**10.
+        
+        dimensions : int, optional
+            Number of dimensions. Default is 3.
+        
+        rejection_radius : float, optional
+            Value of rejection radius. Default is 0.5. 
+        
+        box_lengths : array-like, optional
+            Lengths of the simulation box in each dimension. Default is None.
+        
+        kwargs : dict, optional
+            Additional keyword arguments to pass for the initialization of the `scipy.stats.qmc.PoissonDisk` class.
+
+        Returns
+        -------
+        array-like
+            Randomly generated positions of the particles, scaled by the box lengths.
+        
+        Raises
+        ------
+        ValueError
+            If an invalid sequence name is provided.
+        
+        Notes
+        -----
+        The `quasi_monte_carlo_rejection` method initializes a quasi-Monte Carlo sequence generator based on the provided sequence
+        name. It then generates random positions for the specified number of particles in the specified number of
+        dimensions. The generated positions are scaled by the box lengths of the simulation.
+
+        Examples
+        --------
+        >>> positions = quasi_monte_carlo_rejection("poissondisk", 100, 3, 0.5, [10, 10, 10])
+
+        """
+        if sequence == "poissondisk":
+            qmc_sampler = qmc.poisson(d=dimensions, radius = rejection_radius, **kwargs)
+        else:
+            raise ValueError("Invalid sequence name.")
+
+        positions = qmc_sampler.random(num_ptcls)
+        positions = qmc.scale(positions, l_bounds=[0.0, 0.0, 0.0], u_bounds=box_lengths)
+
+        return positions
+    
     def initialize_velocities(self, species):
         """
         Initialize particles' velocities based on the species input values. The velocities can be initialized from a
@@ -706,13 +838,13 @@ def initialize_velocities(self, species):
 
                     self.species_thermal_velocity[ic] = sqrt(self.kB * sp_temperature / sp.mass)
                     # Note gaussian(0.0, 0.0, N) = array of zeros
-                    self.vel[species_start:species_end, :] = self.gaussian(
-                        sp.initial_velocity, self.species_thermal_velocity[ic], (sp.num, 3)
+                    self.vel[species_start:species_end, :self.dimensions] = self.gaussian(
+                        sp.initial_velocity, self.species_thermal_velocity[ic], (sp.num, self.dimensions)
                     )
 
                 elif sp.initial_velocity_distribution == "monochromatic":
                     vrms = sqrt(self.dimensions * self.kB * sp.temperature / sp.mass)
-                    self.vel[species_start:species_end, :] = vrms * self.random_unit_vectors(sp.num, self.dimensions)
+                    self.vel[species_start:species_end, :self.dimensions] = vrms * self.random_unit_vectors(sp.num, self.dimensions)
 
                 species_start += sp.num
 
@@ -1057,23 +1189,6 @@ def calculate_observables(self):
         for i in self.observables_list:
             self.species_observables_calculator_dict[i]()
 
-    # def calculate_pressure(self):
-    #     """Calculate the pressure of each particle.
-
-    #     Return
-    #     ------
-    #     pressure : numpy.ndarray
-    #         Pressure of each particle. Shape = (:attr:`total_num_ptcls`)
-
-    #     Notes
-    #     -----
-    #     It does not calculate the tensor since that could lead to very large arrays.
-
-    #     """
-
-    #     self.pressure = 2.0 * self.kinetic_energy + self.virial_species_tensor[0, 0] + self.virial_species_tensor[1, 1] + self.virial_species_tensor[2, 2]
-    #     self.pressure /= self.box_volume
-
     def calculate_species_electric_current(self):
         """Calculate the electric current of each species from :attr:`vel` and stores it into :attr:`species_electric_current`."""
         self.species_electric_current = self.species_charges * vector_species_loop(self.vel, self.species_num)
@@ -1425,7 +1540,7 @@ def setup(self, params, species):
         ]:
             # checks
             if self.restart_step is None:
-                raise AttributeError("Restart step not defined." "Please define Parameters.restart_step.")
+                raise AttributeError("Restart step not defined. Please define Parameters.restart_step.")
 
             if type(self.restart_step) is not int:
                 self.restart_step = int(self.restart_step)
@@ -1447,8 +1562,8 @@ def setup(self, params, species):
             else:
                 self.load_from_file(self.particles_input_file)
         else:
-            self.initialize_positions()
-            self.initialize_velocities(species)
+            self.initialize_positions(species=species)
+            self.initialize_velocities(species=species)
             self.initialize_accelerations()
 
         if len(self.observables_list) > 3:
diff --git a/sarkas/potentials/coulomb.py b/sarkas/potentials/coulomb.py
index f027b817..bb6d08dd 100644
--- a/sarkas/potentials/coulomb.py
+++ b/sarkas/potentials/coulomb.py
@@ -173,11 +173,13 @@ def pretty_print_info(potential):
 
     """
 
+    msg = f"Effective coupling constant: Gamma_eff = {potential.coupling_constant:.2f}\n"
+
     if potential.method != "fmm":
         fmm_msg = f"Short-range cutoff radius: a_rs = {potential.a_rs:.6e} {potential.units_dict['length']}"
-    msg = f"Effective coupling constant: Gamma_eff = {potential.coupling_constant:.2f}\n"
+        msg += fmm_msg
 
-    print(msg + fmm_msg)
+    print(msg)
 
 
 def update_params(potential):
diff --git a/sarkas/processes.py b/sarkas/processes.py
index fb8e8a10..83d1e742 100644
--- a/sarkas/processes.py
+++ b/sarkas/processes.py
@@ -84,42 +84,89 @@ class Process:
 
     """
 
-    def __init__(self, input_file: str = None):
-        self.potential = Potential()
-        self.integrator = Integrator()
-        self.parameters = Parameters()
-        self.particles = Particles()
+    from .potentials.core import Potential
+
+    def __init__(self, 
+                 input_file: str = None, 
+                 potential_class: Potential = None,
+                 integrator_class: Integrator = None,
+                 particles_class: Particles = None,
+                 parameters_class: Parameters = None,
+                 io_class: InputOutput = None,
+                 species: list = None,
+                 
+                 ):
+        
+        if potential_class is not None:
+            self.potential = potential_class
+        else:
+            self.potential = Potential()
 
-        self.species = []  # Deprecated
-        self.species = []
+        if integrator_class is not None:
+            self.integrator = integrator_class
+        else:
+            self.integrator = Integrator()
+
+        if particles_class is not None:
+            self.particles = particles_class
+        else:   
+            self.particles = Particles()
+        
+        if parameters_class is not None:
+            self.parameters = parameters_class
+        else:
+            self.parameters = Parameters()
+        
+        if io_class is not None:
+            self.io = io_class
+        else:
+            self.io = InputOutput(process=self.__name__)
+
+        if species is not None:
+            self.species = species
+        else:
+            self.species = []
+    
         self.threads_ls = []
         self.observables_dict = {}
         self.transport_dict = {}
 
         self.input_file = input_file
         self.timer = SarkasTimer()
-        self.io = InputOutput(process=self.__name__)
 
     def common_parser(self, filename: str = None):
         """
-        Parse simulation parameters from YAML file.
+        Parse simulation parameters from a YAML file.
 
         Parameters
         ----------
-        filename: str
-            Input YAML file
+        filename : str, optional
+            Path to the YAML input file. If not provided, the input file path specified during object initialization will be used.
 
-        Return
-        ------
-        dics : dict
-            Nested dictionary from reading the YAML input file.
+        Returns
+        -------
+        dict
+            A nested dictionary containing the parsed simulation parameters.
+
+        Notes
+        -----
+        This method reads the simulation parameters from a YAML file and returns them as a nested dictionary. 
+        It uses the :meth:`sarkas.utilities.io.InputOutput.from_yaml` to read the YAML file.
+
+        If the `filename` parameter is provided, it will override the input file path specified during object initialization.
+
+        Examples
+        --------
+        >>> process = Process(input_file='/path/to/input.yaml')
+        >>> params_dict = process.common_parser()
+        
         """
         if filename:
             self.input_file = filename
 
-        dics = self.io.from_yaml(self.input_file)
+        params_dict = self.io.from_yaml(self.input_file)
 
-        return dics
+        return params_dict
 
     def update_subclasses_from_dict(self, nested_dict: dict):
         """Update the subclasses parameters using a dictionary.
@@ -256,7 +303,7 @@ def directory_sizes(self):
             )
 
         # Grab one file from the dump directory and get the size of it.
-        prod_dump_size = os_stat(join(self.io.eq_dump_dir, listdir(self.io.eq_dump_dir)[0])).st_size
+        prod_dump_size = os_stat(join(self.io.prod_dump_dir, listdir(self.io.prod_dump_dir)[0])).st_size
         prod_dump_fldr_size = prod_dump_size * (self.parameters.production_steps / self.parameters.prod_dump_step)
         # Prepare arguments to pass for print out
         sizes = array([[eq_dump_size, eq_dump_fldr_size], [prod_dump_size, prod_dump_fldr_size]])
@@ -380,8 +427,8 @@ def initialization(self):
         self.io.setup()
 
         # Copy relevant subsclasses attributes into parameters class. This is needed for post-processing.
+        # it updates parameters' dictionary with filenames and directories
         self.parameters.copy_io_attrs(self.io)
-        # Update parameters' dictionary with filenames and directories
 
         self.parameters.potential_type = self.potential.type.lower()
         self.parameters.setup(self.species)
diff --git a/sarkas/time_evolution/integrators.py b/sarkas/time_evolution/integrators.py
index 5698008d..65c68d97 100644
--- a/sarkas/time_evolution/integrators.py
+++ b/sarkas/time_evolution/integrators.py
@@ -1091,13 +1091,13 @@ def pretty_print(self):
             low_wc = abs(self.species_cyclotron_frequencies).min() 
             t_wc_high = 2.0 * pi/high_wc
             t_wc_low = 2.0 * pi/low_wc
-
+            mag_msg = "\nMagnetic Timescales:\n"
             if high_wc > low_wc:
-                mag_msg = (
+                mag_msg += (
                     f"Largest cyclotron frequency (w_c) = {high_wc:.6e} {self.units_dict['frequency']}\n"
                     f"w_c dt = {high_wc * self.dt:.2e} [rad]\n"
                     f"Cyclotron period (T) = {t_wc_high:.6e} {self.units_dict['time']}\n"
-                    f"dt/T = = {self.dt/t_wc_high:.2e}"
+                    f"dt/T = = {self.dt/t_wc_high:.2e}\n"
                     f"Timesteps per plasma cycle ~ {int(t_wc_high/self.dt)}\n"
                     f"Smallest cyclotron frequency (w_c) = {low_wc:.6e} {self.units_dict['frequency']}\n"
                     f"w_c dt = {low_wc * self.dt:.2e} [rad]\n"
@@ -1106,7 +1106,7 @@ def pretty_print(self):
                     f"Timesteps per plasma cycle ~ {int(t_wc_low/self.dt)}\n"
                         )
             else:
-                mag_msg = (
+                mag_msg += (
                     f"Cyclotron frequency (w_c) = {high_wc:.6e} {self.units_dict['frequency']}\n"
                     f"w_c dt = {high_wc * self.dt:.2e} [rad]\n"
                     f"Cyclotron period (T) = {t_wc_high:.6e} {self.units_dict['time']}\n"