diff --git a/src/mechanalyzer/calculator/bf.py b/src/mechanalyzer/calculator/bf.py
index 5ce0af7..512894a 100644
--- a/src/mechanalyzer/calculator/bf.py
+++ b/src/mechanalyzer/calculator/bf.py
@@ -8,6 +8,12 @@
 import numpy as np
 import pandas as pd
 from mechanalyzer import calculator
+# AVC: Backward compatibility for numpy < 2.0, where the function is called `trapz`
+try:
+    from numpy import trapezoid
+except ImportError:
+    from numpy import trapz as trapezoid
+
 
 ######################## wrapper functions #########################################################
 
@@ -160,7 +166,7 @@ def bf_tp_df_full(ped_df, hotbf_df):
                         kind='cubic', fill_value=(hoten_spc.values[0], hoten_spc.values[-1]))
                     hoten_vect = f_hoten(ene_vect)
                     # recompute in an appropriate range
-                    bf_series[spc] = np.trapz(ped_vect*hoten_vect, x=ene_vect)
+                    bf_series[spc] = trapezoid(ped_vect*hoten_vect, x=ene_vect)
                 # renormalize for all species and put in dataframe
                 if any(bf_series.values < 0):
                     print('Warning: found negative BFs at {:1.0f} K and {:1.1e} atm'
diff --git a/src/mechanalyzer/calculator/ene_partition.py b/src/mechanalyzer/calculator/ene_partition.py
index 0e40f7d..0ca7c7b 100644
--- a/src/mechanalyzer/calculator/ene_partition.py
+++ b/src/mechanalyzer/calculator/ene_partition.py
@@ -10,6 +10,11 @@
 from phydat import phycon
 from mechanalyzer import calculator
 from mechanalyzer.calculator.spinfo_frommess import get_dof_info_fromspcdct
+# AVC: Backward compatibility for numpy < 2.0, where the function is called `trapz`
+try:
+    from numpy import trapezoid
+except ImportError:
+    from numpy import trapz as trapezoid
 
 
 #################### wrapper function that calls the class ################################
@@ -144,7 +149,7 @@ def ped_df_rescale(starthot_df, ped_df_fromhot, save=False, name=''):
                     kernelsize)/kernelsize, mode='same')
                 cycles -= 1
             # renormalize and put in dataframe
-            prob_vect /= np.trapz(prob_vect, x=ene_vect)
+            prob_vect /= trapezoid(prob_vect, x=ene_vect)
             ped_df.at[temp, pressure] = pd.Series(prob_vect, index=ene_vect)
             if ped_df[pressure][temp].empty:
                 T_del.append(temp)
@@ -280,7 +285,7 @@ def equip_simple(self):
         for pressure in self.ped_df.columns:
             for temp in self.ped_df.sort_index().index:
                 idx_new = self.ped_df[pressure][temp].index * beta_prod
-                norm_factor = np.trapz(
+                norm_factor = trapezoid(
                     self.ped_df[pressure][temp].values, x=idx_new)
                 vals = self.ped_df[pressure][temp].values/norm_factor
                 ped_df_prod.at[temp, pressure] = pd.Series(vals, index=idx_new)
@@ -348,7 +353,7 @@ def init_dos(pressure, temp):
                     rho_rovib_prod2[idx_ene_int] *
                     rho_trasl[idx_ene_minus_ene_int]
                 )
-                rho_non1.append(np.trapz(rho_non1_integrand,
+                rho_non1.append(trapezoid(rho_non1_integrand,
                                         x=self.ene1_vect[idx_ene_int]))
 
             rho_non1 = np.array(rho_non1)
@@ -378,7 +383,7 @@ def dos(idx_ene1, idx_ene_new):
                 # rhonon1(ene-ene1) with ene1<ene (fixed ene)
                 rho_non1_array = self.rho_non1[idx_ene_minus_ene1_array]
                 num = rho1_ene1 * rho_non1
-                den = np.trapz(rho1_ene1_array*rho_non1_array,
+                den = trapezoid(rho1_ene1_array*rho_non1_array,
                                x=self.ene1_vect[idx_ene1_array])
                 # if for some reason you get den=0: append 0 as value
                 if den == 0:
@@ -408,13 +413,13 @@ def dostherm_rhovibtrasl(idx_ene_vect, temp):
                 rho1_ene1_array = self.rho_rovib_prod1[idx_ene1_array]
                 rho_non1_array = self.rho_non1[idx_ene_minus_ene1_array]
 
-                rhotot_E1 = np.trapz(rho1_ene1_array*rho_non1_array,
+                rhotot_E1 = trapezoid(rho1_ene1_array*rho_non1_array,
                                      x=self.ene1_vect[idx_ene1_array])
 
                 ene = self.ene1_vect[idx_ene1]
                 f_Etot_num.append(rhotot_E1*np.exp(-ene/phycon.RC_KCAL/temp))
 
-            den = np.trapz(f_Etot_num, x=self.ene1_vect[idx_ene_vect])
+            den = trapezoid(f_Etot_num, x=self.ene1_vect[idx_ene_vect])
             f_Etot = pd.Series(
                 f_Etot_num/den, index=self.ene1_vect[idx_ene_vect])
 
@@ -425,7 +430,7 @@ def dostherm_rhovib(temp):
             # rovib only
             rho1_ene1_array = self.rho_rovib_prod1 * \
                 np.exp(-self.ene1_vect/phycon.RC_KCAL/temp)
-            den = np.trapz(rho1_ene1_array, x=self.ene1_vect)
+            den = trapezoid(rho1_ene1_array, x=self.ene1_vect)
             f_Etot = pd.Series(rho1_ene1_array/den, index=self.ene1_vect)
 
             return f_Etot
@@ -489,11 +494,11 @@ def dostherm_rhovib(temp):
                             ped_series.values[idx_ene_vect >= idx_ene1]
                         )
 
-                        prob_ene1 = np.trapz(
+                        prob_ene1 = trapezoid(
                             prob_ene1ene_tot_pressure_ped, ene_new)
                         prob_ene1_vect.append(prob_ene1)
 
-                    norm_factor_prob_ene1 = np.trapz(
+                    norm_factor_prob_ene1 = trapezoid(
                         prob_ene1_vect, x=self.ene1_vect)
                     prob_ene1_norm = prob_ene1_vect/norm_factor_prob_ene1
 
diff --git a/src/mechanalyzer/tests/test__prompt.py b/src/mechanalyzer/tests/test__prompt.py
index 4629b3d..c375313 100644
--- a/src/mechanalyzer/tests/test__prompt.py
+++ b/src/mechanalyzer/tests/test__prompt.py
@@ -9,6 +9,11 @@
 from mechanalyzer.calculator import spinfo_frommess
 from mechanalyzer.calculator import ene_partition
 from mechanalyzer.calculator import bf
+# AVC: Backward compatibility for numpy < 2.0, where the function is called `trapz`
+try:
+    from numpy import trapezoid
+except ImportError:
+    from numpy import trapz as trapezoid
 
 PATH = os.path.dirname(os.path.realpath(__file__))
 INP_PATH = os.path.join(PATH, 'data', 'prompt', 'C3H8_OH')
@@ -70,8 +75,8 @@ def test_equip_simple():
 
     assert np.isclose((ped_600.iloc[100]), 0.1043, atol=1e-4, rtol=1e-4)
     assert np.isclose((ped_1200.iloc[100]), 0.02105, atol=1e-4, rtol=1e-4)
-    assert np.isclose(np.trapz(ped_600.values, x=ped_600.index), 1)
-    assert np.isclose(np.trapz(ped_1200.values, x=ped_1200.index), 1)
+    assert np.isclose(trapezoid(ped_600.values, x=ped_600.index), 1)
+    assert np.isclose(trapezoid(ped_1200.values, x=ped_1200.index), 1)
 
 
 def test_equip_phi():
@@ -89,8 +94,8 @@ def test_equip_phi():
 
     assert np.isclose((ped_600.iloc[166]), 0.045, atol=1e-3, rtol=1e-2)
     assert np.isclose((ped_1200.iloc[166]), 0.0376, atol=1e-3, rtol=1e-2)
-    assert np.isclose(np.trapz(ped_600.values, x=ped_600.index), 1)
-    assert np.isclose(np.trapz(ped_1200.values, x=ped_1200.index), 1)
+    assert np.isclose(trapezoid(ped_600.values, x=ped_600.index), 1)
+    assert np.isclose(trapezoid(ped_1200.values, x=ped_1200.index), 1)
 
 
 def test_beta_phi1a():
@@ -108,8 +113,8 @@ def test_beta_phi1a():
 
     assert np.isclose((ped_400.iloc[166]), 0.001, atol=1e-3, rtol=1e-2)
     assert np.isclose((ped_800.iloc[166]), 0.0589, atol=1e-3, rtol=1e-2)
-    assert np.isclose(np.trapz(ped_400.values, x=ped_400.index), 1)
-    assert np.isclose(np.trapz(ped_800.values, x=ped_800.index), 1)
+    assert np.isclose(trapezoid(ped_400.values, x=ped_400.index), 1)
+    assert np.isclose(trapezoid(ped_800.values, x=ped_800.index), 1)
 
 
 def test_beta_phi2a():
@@ -127,8 +132,8 @@ def test_beta_phi2a():
 
     assert np.isclose((ped_400.iloc[166]), 0.00045, atol=1e-4, rtol=1e-2)
     assert np.isclose((ped_800.iloc[166]), 0.0522, atol=1e-3, rtol=1e-2)
-    assert np.isclose(np.trapz(ped_400.values, x=ped_400.index), 1)
-    assert np.isclose(np.trapz(ped_800.values, x=ped_800.index), 1)
+    assert np.isclose(trapezoid(ped_400.values, x=ped_400.index), 1)
+    assert np.isclose(trapezoid(ped_800.values, x=ped_800.index), 1)
 
 
 def test_beta_phi3a():
@@ -146,8 +151,8 @@ def test_beta_phi3a():
 
     assert np.isclose((ped_400.iloc[166]), 0.001, atol=1e-3, rtol=1e-2)
     assert np.isclose((ped_800.iloc[166]), 0.0589, atol=1e-3, rtol=1e-2)
-    assert np.isclose(np.trapz(ped_400.values, x=ped_400.index), 1)
-    assert np.isclose(np.trapz(ped_800.values, x=ped_800.index), 1)
+    assert np.isclose(trapezoid(ped_400.values, x=ped_400.index), 1)
+    assert np.isclose(trapezoid(ped_800.values, x=ped_800.index), 1)
 
 
 def test_rovib_dos():
@@ -167,8 +172,8 @@ def test_rovib_dos():
 
     assert np.isclose((ped_1800.iloc[166]), 0.01985, atol=1e-3, rtol=1e-2)
     assert np.isclose((ped_2000.iloc[166]), 0.01724, atol=1e-3, rtol=1e-2)
-    assert np.isclose(np.trapz(ped_1800.values, x=ped_1800.index), 1)
-    assert np.isclose(np.trapz(ped_2000.values, x=ped_2000.index), 1)
+    assert np.isclose(trapezoid(ped_1800.values, x=ped_1800.index), 1)
+    assert np.isclose(trapezoid(ped_2000.values, x=ped_2000.index), 1)
 
 
 def test_thermal():
@@ -187,8 +192,8 @@ def test_thermal():
 
     assert np.isclose((ped_1800.iloc[166]), 0.01924792, atol=1e-5, rtol=1e-5)
     assert np.isclose((ped_2000.iloc[166]), 0.01132143, atol=1e-5, rtol=1e-5)
-    assert np.isclose(np.trapz(ped_1800.values, x=ped_1800.index), 1)
-    assert np.isclose(np.trapz(ped_2000.values, x=ped_2000.index), 1)
+    assert np.isclose(trapezoid(ped_1800.values, x=ped_1800.index), 1)
+    assert np.isclose(trapezoid(ped_2000.values, x=ped_2000.index), 1)
 
 
 def test_bf_from_phi1a():