diff --git a/tardis/energy_input/util.py b/tardis/energy_input/util.py
index 061017158b7..1a1b4102a0a 100644
--- a/tardis/energy_input/util.py
+++ b/tardis/energy_input/util.py
@@ -70,7 +70,9 @@ def doppler_factor_3d(direction_vector, position_vector, time):
         Doppler factor
     """
     velocity_vector = position_vector / time
-    return 1 - (np.dot(direction_vector, velocity_vector) / C_CGS)
+    direction_vector_contiguous = np.ascontiguousarray(direction_vector)
+    velocity_vector_contiguous = np.ascontiguousarray(velocity_vector)
+    return 1 - (np.dot(direction_vector_contiguous, velocity_vector_contiguous) / C_CGS)
 
 
 @njit(**njit_dict_no_parallel)
diff --git a/tardis/io/model/parse_radiation_field_configuration.py b/tardis/io/model/parse_radiation_field_configuration.py
index d311ad94923..725f69c4f8c 100644
--- a/tardis/io/model/parse_radiation_field_configuration.py
+++ b/tardis/io/model/parse_radiation_field_configuration.py
@@ -167,15 +167,8 @@ def calculate_geometric_dilution_factor(geometry):
     np.array
         The dilution factors for the inout geometry
     """
-    return 0.5 * (
-        1
-        - np.sqrt(
-            1
-            - (
-                geometry.r_inner[geometry.v_inner_boundary_index] ** 2
-                / geometry.r_middle**2
-            )
-            .to(1)
-            .value
-        )
+    value = (
+        1 - (geometry.r_inner[geometry.v_inner_boundary_index] ** 2 / geometry.r_middle**2).to(1).value
     )
+    value = np.clip(value, 0, None)
+    return 0.5 * (1 - np.sqrt(value))
diff --git a/tardis/io/tests/test_atomic.py b/tardis/io/tests/test_atomic.py
index 4c810b8a7f2..7e8def3cc3e 100644
--- a/tardis/io/tests/test_atomic.py
+++ b/tardis/io/tests/test_atomic.py
@@ -48,6 +48,7 @@ def test_atom_data_levels(levels):
 
 
 def test_atom_data_lines(lines):
+    lines = lines.sort_index()
     assert_quantity_allclose(
         lines.loc[(2, 0, 0, 6), "wavelength_cm"].values[0] * u.Unit("cm"),
         584.335 * u.Unit("Angstrom"),
diff --git a/tardis/model/matter/composition.py b/tardis/model/matter/composition.py
index 67b0a0ec8ce..23f54882f42 100644
--- a/tardis/model/matter/composition.py
+++ b/tardis/model/matter/composition.py
@@ -93,7 +93,7 @@ def __init__(
 
     def assemble_isotope_masses(self):
         isotope_mass_df = pd.Series(
-            index=self.isotopic_mass_fraction.index, data=-1
+            index=self.isotopic_mass_fraction.index, data=-1,dtype="float64"
         )
         for isotope_tuple in self.isotopic_mass_fraction.index:
             isotope_symbol = int("{:03d}{:03d}0000".format(*isotope_tuple))
diff --git a/tardis/plasma/properties/atomic.py b/tardis/plasma/properties/atomic.py
index 47ac12eaafe..076738f7aae 100644
--- a/tardis/plasma/properties/atomic.py
+++ b/tardis/plasma/properties/atomic.py
@@ -787,7 +787,11 @@ def exp1_times_exp(x):
         array_like
             Output array.
         """
-        f = exp1(x) * np.exp(x)
+        x = np.asarray(x)
+        f = np.empty_like(x)
+        small_x = x <= 500
+        f[small_x] = exp1(x[small_x]) * np.exp(x[small_x])
+
         # Use Laurent series for large values to avoid infinite exponential
         mask = x > 500
         f[mask] = (x**-1 - x**-2 + 2 * x**-3 - 6 * x**-4)[mask]
diff --git a/tardis/plasma/properties/nlte_rate_equation_solver.py b/tardis/plasma/properties/nlte_rate_equation_solver.py
index c15a46bf196..ac2e7fe7a1c 100644
--- a/tardis/plasma/properties/nlte_rate_equation_solver.py
+++ b/tardis/plasma/properties/nlte_rate_equation_solver.py
@@ -554,7 +554,7 @@ def calculate_first_guess(
     """
     first_guess = pd.Series(0.0, index=rate_matrix_index)
     for atomic_number in atomic_numbers:
-        first_guess.loc[(atomic_number, 1)].iloc[0] = number_density.loc[
+        first_guess.at[(atomic_number, 1)] = number_density.loc[
             atomic_number
         ]
     # TODO: After the first iteration, the new guess can be the old solution.
@@ -937,7 +937,7 @@ def ion_matrix(ion_coefficients, atomic_number, ion_number):
     offdiag = np.zeros(atomic_number)
     index = ion_coefficients.index
     for i in index:
-        offdiag[i] = ion_coefficients.loc[i]
+        offdiag[i] = float(ion_coefficients.loc[i].iloc[0])
     diag = np.hstack([-offdiag, np.zeros(1)])
     return (np.diag(diag) + np.diag(offdiag, k=-1))[ion_number, :]
 
@@ -961,7 +961,7 @@ def recomb_matrix(recomb_coefficients, atomic_number, ion_number):
     offdiag = np.zeros(atomic_number)
     index = recomb_coefficients.index
     for i in index:
-        offdiag[i] = recomb_coefficients.loc[i]
+        offdiag[i] = float(recomb_coefficients.loc[i].values[0])
     diag = np.hstack([np.zeros(1), -offdiag])
     return (np.diag(diag) + np.diag(offdiag, k=1))[ion_number, :]