diff --git a/pypl/config_coord_1d_solver.py b/pypl/config_coord_1d_solver.py
index 5b575fd..7918527 100644
--- a/pypl/config_coord_1d_solver.py
+++ b/pypl/config_coord_1d_solver.py
@@ -404,7 +404,7 @@ def compute_spectrum(self, eneaxis=None, linshape=None, tdm=1.0, zpl=0.0, spectr
     ax.tick_params(direction="in")
     ax.xaxis.set_ticks_position("both")
     ax.yaxis.set_ticks_position("both")
-    ax.set_xlabel("$\hbar\omega$ (eV)")
+    ax.set_xlabel(r"$\hbar\omega$ (eV)")
     ax.set_ylabel("PL/Abs (arb. unit.)")
 
     plt.show()
diff --git a/pypl/hr_factors.py b/pypl/hr_factors.py
index f0e0406..f48a3a6 100644
--- a/pypl/hr_factors.py
+++ b/pypl/hr_factors.py
@@ -220,7 +220,7 @@ def compute_hrf_forces(self, forces):
         deltaq = deltaq / self.freqs[self.nom_translational + self.nom_imag_freq :] ** 2
 
         print(
-            "Total \Delta Q is % .12e amu^{0.5} \AA"
+            r"Total \Delta Q is % .12e amu^{0.5} \AA"
             % (np.linalg.norm(deltaq) / constants.physical_constants["atomic mass constant"][0] ** 0.5 * 1e10)
         )
 
@@ -276,7 +276,7 @@ def compute_hrf_dis(self, gs_coord, es_coord, cell_parameters):
         deltaq = np.dot(modes[self.nom_translational + self.nom_imag_freq :, :], mass_dis)
 
         print(
-            "Total \Delta Q is % .12e amu^{0.5} \AA"
+            r"Total \Delta Q is % .12e amu^{0.5} \AA"
             % (np.linalg.norm(deltaq) / constants.physical_constants["atomic mass constant"][0] ** 0.5 * 1e10)
         )