diff --git a/statmechcrack/isometric.py b/statmechcrack/isometric.py
index e1a6422..ec5d2de 100644
--- a/statmechcrack/isometric.py
+++ b/statmechcrack/isometric.py
@@ -130,13 +130,14 @@ def beta_A_abs_isometric(self, v, approach='asymptotic'):
             numpy.ndarray: The absolute nondimensional Helmholtz free energy.
 
         """
-        s_hat = self.minimize_beta_U(v)[2]
-        lambda_hat = s_hat[-self.M:]
-        return 0.5*self.M*np.log(
-            self.alpha**2*self.varepsilon/np.pi
-        ) + self.beta_A_0_isometric(v, lambda_hat)
+        if approach == 'asymptotic':
+            s_hat = self.minimize_beta_U(v)[2]
+            lambda_hat = s_hat[-self.M:]
+            return 0.5*self.M*np.log(
+                self.alpha**2*self.varepsilon/np.pi
+            ) + self.beta_A_0_isometric(v, lambda_hat)
 
-    def beta_A_isometric(self, v, approach='asymptotic'):
+    def beta_A_isometric(self, v, approach='asymptotic', **kwargs):
         r"""The relative nondimensional Helmholtz free energy
         as a function of the nondimensional end separation,
 
@@ -147,6 +148,8 @@ def beta_A_isometric(self, v, approach='asymptotic'):
             v (array_like): The nondimensional end separation.
             approach (str, optional, default='asymptotic'):
                 The calculation approach.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`~.beta_A_isometric_monte_carlo`.
 
         Returns:
             numpy.ndarray: The relative nondimensional Helmholtz free energy.
@@ -180,8 +183,13 @@ def beta_A_isometric(self, v, approach='asymptotic'):
                 >>> plt.show()
 
         """
-        return self.beta_A_abs_isometric(v, approach=approach) \
-            - self.beta_A_abs_isometric(1, approach=approach)
+        if approach == 'asymptotic':
+            beta_A_isometric = \
+                self.beta_A_abs_isometric(v, approach=approach) - \
+                self.beta_A_abs_isometric(1, approach=approach)
+        elif approach == 'monte carlo':
+            beta_A_isometric = self.beta_A_isometric_monte_carlo(v, **kwargs)
+        return beta_A_isometric
 
     def p_isometric(self, v, approach='asymptotic', **kwargs):
         r"""The nondimensional end force
@@ -237,7 +245,7 @@ def p_isometric(self, v, approach='asymptotic', **kwargs):
             p_isometric = self.p_isometric_monte_carlo(v, **kwargs)
         return p_isometric
 
-    def k_isometric(self, v, approach='asymptotic'):
+    def k_isometric(self, v, approach='asymptotic', **kwargs):
         r"""The nondimensional forward reaction rate coefficient
         as a function of the nondimensional end separation
         in the isometric ensemble.
@@ -246,6 +254,8 @@ def k_isometric(self, v, approach='asymptotic'):
             v (array_like): The nondimensional end separation.
             approach (str, optional, default='asymptotic'):
                 The calculation approach.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`~.k_isometric_monte_carlo`.
 
         Returns:
             numpy.ndarray: The nondimensional forward reaction rate.
@@ -277,13 +287,15 @@ def k_isometric(self, v, approach='asymptotic'):
                 >>> plt.show()
 
         """
-        return (
-            self.Q_isometric(v, approach=approach, transition_state=True) /
-            self.Q_isometric(v, approach=approach)
-        ) / (
-            self.Q_isometric(1, approach=approach, transition_state=True) /
-            self.Q_isometric(1, approach=approach)
-        )
+        if approach == 'asymptotic':
+            k_isometric = (
+                self.Q_isometric(v, transition_state=True)/self.Q_isometric(v)
+            ) / (
+                self.Q_isometric(1, transition_state=True)/self.Q_isometric(1)
+            )
+        elif approach == 'monte carlo':
+            k_isometric = self.k_isometric_monte_carlo(v, **kwargs)
+        return k_isometric
 
     def Q_0_isometric(self, v, lambda_):
         r"""The nondimensional isometric partition function
diff --git a/statmechcrack/isotensional.py b/statmechcrack/isotensional.py
index 587cde0..b8da785 100644
--- a/statmechcrack/isotensional.py
+++ b/statmechcrack/isotensional.py
@@ -125,11 +125,12 @@ def beta_G_abs_isotensional(self, p, approach='asymptotic'):
             numpy.ndarray: The absolute nondimensional Gibbs free energy.
 
         """
-        s_hat = self.minimize_beta_Pi(p)[2]
-        lambda_hat = s_hat[-self.M:]
-        return self.beta_G_0_abs_isotensional(p, lambda_hat)
+        if approach == 'asymptotic':
+            s_hat = self.minimize_beta_Pi(p)[2]
+            lambda_hat = s_hat[-self.M:]
+            return self.beta_G_0_abs_isotensional(p, lambda_hat)
 
-    def beta_G_isotensional(self, p, approach='asymptotic'):
+    def beta_G_isotensional(self, p, approach='asymptotic', **kwargs):
         r"""The relative nondimensional Gibbs free energy
         as a function of the nondimensional end force,
 
@@ -140,6 +141,8 @@ def beta_G_isotensional(self, p, approach='asymptotic'):
             p (array_like): The nondimensional end force.
             approach (str, optional, default='asymptotic'):
                 The calculation approach.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`~.beta_G_isotensional_monte_carlo`.
 
         Returns:
             numpy.ndarray: The relative nondimensional Gibbs free energy.
@@ -174,8 +177,13 @@ def beta_G_isotensional(self, p, approach='asymptotic'):
                 >>> plt.show()
 
         """
-        return self.beta_G_abs_isotensional(p, approach=approach) \
-            - self.beta_G_abs_isotensional(0, approach=approach)
+        if approach == 'asymptotic':
+            G_isotensional = \
+                self.beta_G_abs_isotensional(p, approach=approach) - \
+                self.beta_G_abs_isotensional(0, approach=approach)
+        elif approach == 'monte carlo':
+            G_isotensional = self.beta_G_isotensional_monte_carlo(p, **kwargs)
+        return G_isotensional
 
     def v_isotensional(self, p, approach='asymptotic', **kwargs):
         r"""The nondimensional end separation
@@ -237,7 +245,7 @@ def v_isotensional(self, p, approach='asymptotic', **kwargs):
         elif approach == 'monte carlo':
             return self.v_isotensional_monte_carlo(p, **kwargs)
 
-    def k_isotensional(self, p, approach='asymptotic'):
+    def k_isotensional(self, p, approach='asymptotic', **kwargs):
         r"""The nondimensional forward reaction rate coefficient
         as a function of the nondimensional end force
         in the isotensional ensemble.
@@ -278,13 +286,17 @@ def k_isotensional(self, p, approach='asymptotic'):
                 >>> plt.show()
 
         """
-        return (
-            self.Z_isotensional(p, approach=approach, transition_state=True) /
-            self.Z_isotensional(p, approach=approach)
-        ) / (
-            self.Z_isotensional(0, approach=approach, transition_state=True) /
-            self.Z_isotensional(0, approach=approach)
-        )
+        if approach == 'asymptotic':
+            k_isotensional = (
+                self.Z_isotensional(p, transition_state=True) /
+                self.Z_isotensional(p)
+            ) / (
+                self.Z_isotensional(0, transition_state=True) /
+                self.Z_isotensional(0)
+            )
+        elif approach == 'monte carlo':
+            k_isotensional = self.k_isotensional_monte_carlo(p, **kwargs)
+        return k_isotensional
 
     def Z_0_isotensional(self, p, lambda_):
         r"""The nondimensional isotensional partition function
diff --git a/statmechcrack/monte_carlo.py b/statmechcrack/monte_carlo.py
index c205743..c2c9db0 100644
--- a/statmechcrack/monte_carlo.py
+++ b/statmechcrack/monte_carlo.py
@@ -175,9 +175,71 @@ def fun(seed, output):
             return np.mean(process_results)
         return serial_fun(burned_in_config, **kwargs)
 
+    def p_isometric_monte_carlo_serial(self, v, init_config,
+                                       num_samples=1000000, **kwargs):
+        """Serial calculation for :meth:`p_isometric_monte_carlo`.
+
+        Args:
+            v (array_like):
+                The nondimensional end separation.
+            init_config (np.ndarray):
+                The initial configuration, typically already burned-in.
+            num_samples (int, optional, default=1000000):
+                The number of samples to use.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`mh_next_config`.
+
+        Returns:
+            float: The nondimensional end force.
+
+        """
+        p = 0
+        config = init_config
+        for counter in range(1, 1 + num_samples):
+            p_next = self.kappa*(v - (2*config[0] - config[1]))
+            p += (p_next - p)/counter
+            config = self.mh_next_config(config, **kwargs)
+        return p
+
+    def p_isometric_monte_carlo(self, v, **kwargs):
+        r"""The nondimensional end force
+        as a function of the nondimensional end separation
+        in the isometric ensemble, using a
+        Metropolis-Hastings Markov chain Monte Carlo calculation
+        of the ensemble average
+
+        .. math::
+            p(v) = \kappa\langle v - 2s_1 + s_2\rangle.
+
+        Args:
+            v (array_like): The nondimensional end separation.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`parallel_calculation`.
+
+        Returns:
+            float: The nondimensional end force.
+
+        """
+        v = self.np_array(v)
+        p = np.zeros(v.shape)
+        for i, v_i in enumerate(v):
+            self.beta_e = lambda config: self.beta_U(v_i, config)
+
+            def serial_fun(init_config, **kwargs):
+                return self.p_isometric_monte_carlo_serial(
+                        v_i, init_config, **kwargs
+                    )
+
+            p[i] = self.parallel_calculation(
+                serial_fun,
+                self.minimize_beta_U(v_i)[2][:, 0],
+                **kwargs
+            )
+        return p
+
     def v_isotensional_monte_carlo_serial(self, init_config,
                                           num_samples=1000000, **kwargs):
-        """
+        """Serial calculation for :meth:`v_isotensional_monte_carlo`.
 
         Args:
             init_config (np.ndarray):
@@ -200,10 +262,14 @@ def v_isotensional_monte_carlo_serial(self, init_config,
         return v
 
     def v_isotensional_monte_carlo(self, p, **kwargs):
-        """The nondimensional end separation
+        r"""The nondimensional end separation
         as a function of the nondimensional end force
         in the isotensional ensemble, using a
-        Metropolis-Hastings Markov chain Monte Carlo calculation.
+        Metropolis-Hastings Markov chain Monte Carlo calculation
+        of the ensemble average
+
+        .. math::
+            v(p) = \langle s_0\rangle.
 
         Args:
             p (array_like): The nondimensional end force.
@@ -226,9 +292,9 @@ def v_isotensional_monte_carlo(self, p, **kwargs):
             )
         return v
 
-    def p_isometric_monte_carlo_serial(self, v, init_config,
-                                       num_samples=1000000, **kwargs):
-        """
+    def beta_A_isometric_monte_carlo_serial(self, v, init_config,
+                                            num_samples=1000000, **kwargs):
+        """Serial calculation for :meth:`beta_A_isometric_monte_carlo`.
 
         Args:
             v (array_like):
@@ -241,22 +307,23 @@ def p_isometric_monte_carlo_serial(self, v, init_config,
                 Passed to :meth:`mh_next_config`.
 
         Returns:
-            float: The nondimensional end force.
+            numpy.ndarray: The relative nondimensional Helmholtz free energy.
 
         """
-        p = 0
-        config = init_config
-        for counter in range(1, 1 + num_samples):
-            p_next = self.kappa*(v - (2*config[0] - config[1]))
-            p += (p_next - p)/counter
-            config = self.mh_next_config(config, **kwargs)
-        return p
+        pass
 
-    def p_isometric_monte_carlo(self, v, **kwargs):
-        """The nondimensional end force
+    def beta_A_isometric_monte_carlo(self, v, **kwargs):
+        r"""The relative nondimensional Helmholtz free energy
         as a function of the nondimensional end separation
         in the isometric ensemble, using a
-        Metropolis-Hastings Markov chain Monte Carlo calculation.
+        Metropolis-Hastings Markov chain Monte Carlo calculation
+        of the ensemble average
+
+        .. math::
+            \beta\Delta A(v) =
+            -\ln\left\langle e^{-\beta\Delta U}\right\rangle_{v=1},
+
+        where :math:`\Delta U\equiv U(v,s)-U(1,s)`.
 
         Args:
             v (array_like): The nondimensional end separation.
@@ -264,22 +331,186 @@ def p_isometric_monte_carlo(self, v, **kwargs):
                 Passed to :meth:`parallel_calculation`.
 
         Returns:
-            float: The nondimensional end force.
+            numpy.ndarray: The relative nondimensional Helmholtz free energy.
 
         """
         v = self.np_array(v)
-        p = np.zeros(v.shape)
+        beta_A = np.zeros(v.shape)
         for i, v_i in enumerate(v):
             self.beta_e = lambda config: self.beta_U(v_i, config)
 
             def serial_fun(init_config, **kwargs):
-                return self.p_isometric_monte_carlo_serial(
+                return self.beta_A_isometric_monte_carlo_serial(
                         v_i, init_config, **kwargs
                     )
 
-            p[i] = self.parallel_calculation(
+            beta_A[i] = self.parallel_calculation(
                 serial_fun,
                 self.minimize_beta_U(v_i)[2][:, 0],
                 **kwargs
             )
-        return p
+        return beta_A
+
+    def beta_G_isotensional_monte_carlo_serial(self, init_config,
+                                               num_samples=1000000, **kwargs):
+        """Serial calculation for :meth:`beta_G_isotensional_monte_carlo`.
+
+        Args:
+            init_config (np.ndarray):
+                The initial configuration, typically already burned-in.
+            num_samples (int, optional, default=1000000):
+                The number of samples to use.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`mh_next_config`.
+
+        Returns:
+            numpy.ndarray: The relative nondimensional Gibbs free energy.
+
+        """
+        pass
+
+    def beta_G_isotensional_monte_carlo(self, p, **kwargs):
+        r"""The relative nondimensional Gibbs free energy
+        as a function of the nondimensional end force
+        in the isometric ensemble, using a
+        Metropolis-Hastings Markov chain Monte Carlo calculation
+        of the ensemble average
+
+        .. math::
+            \beta\Delta G(v) =
+            -\ln\left\langle e^{-\beta\Delta\Pi}\right\rangle_{p=0},
+
+        where :math:`\Delta\Pi\equiv\Pi(p,s)-\Pi(0,s)`.
+
+        Args:
+            p (array_like): The nondimensional end force.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`parallel_calculation`.
+
+        Returns:
+            numpy.ndarray: The relative nondimensional Gibbs free energy.
+
+        """
+        p = self.np_array(p)
+        beta_G_isotensional = np.zeros(p.shape)
+        for i, p_i in enumerate(p):
+            self.beta_e = lambda vs: self.beta_Pi(p_i, vs[0], vs[1:])
+            v_guess, s_guess = self.minimize_beta_Pi(p_i)[1:]
+            beta_G_isotensional[i] = self.parallel_calculation(
+                self.beta_G_isotensional_monte_carlo_serial,
+                np.concatenate((v_guess, s_guess[:, 0])),
+                **kwargs
+            )
+        return beta_G_isotensional
+
+    def k_isometric_monte_carlo_serial(self, v, init_config,
+                                       num_samples=1000000, **kwargs):
+        """Serial calculation for :meth:`k_isometric_monte_carlo`.
+
+        Args:
+            v (array_like): The nondimensional end separation.
+            init_config (np.ndarray):
+                The initial configuration, typically already burned-in.
+            num_samples (int, optional, default=1000000):
+                The number of samples to use.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`mh_next_config`.
+
+        Returns:
+            numpy.ndarray: The nondimensional forward reaction rate.
+
+        """
+        pass
+
+    def k_isometric_monte_carlo(self, v, **kwargs):
+        r"""The nondimensional forward reaction rate coefficient
+        as a function of the nondimensional end force
+        in the isometric ensemble, using a
+        Metropolis-Hastings Markov chain Monte Carlo calculation
+        of the ensemble average
+
+        .. math::
+            k = \frac{
+                    \left\langle e^{-\beta U}\right\rangle_{v=1}^\ddagger
+                }{
+                    \left\langle e^{-\beta U}\right\rangle_{v=1}
+                }.
+
+        Args:
+            v (array_like): The nondimensional end separation.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`parallel_calculation`.
+
+        Returns:
+            numpy.ndarray: The nondimensional forward reaction rate.
+
+        """
+        v = self.np_array(v)
+        k_isometric = np.zeros(v.shape)
+        for i, v_i in enumerate(v):
+            self.beta_e = lambda config: self.beta_U(v_i, config)
+
+            def serial_fun(init_config, **kwargs):
+                return self.k_isometric_monte_carlo_serial(
+                        v_i, init_config, **kwargs
+                    )
+
+            k_isometric[i] = self.parallel_calculation(
+                serial_fun,
+                self.minimize_beta_U(v_i)[2][:, 0],
+                **kwargs
+            )
+        return k_isometric
+
+    def k_isotensional_monte_carlo_serial(self, init_config,
+                                          num_samples=1000000, **kwargs):
+        """Serial calculation for :meth:`k_isotensional_monte_carlo`.
+
+        Args:
+            init_config (np.ndarray):
+                The initial configuration, typically already burned-in.
+            num_samples (int, optional, default=1000000):
+                The number of samples to use.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`mh_next_config`.
+
+        Returns:
+            numpy.ndarray: The nondimensional forward reaction rate.
+
+        """
+        pass
+
+    def k_isotensional_monte_carlo(self, p, **kwargs):
+        r"""The nondimensional forward reaction rate coefficient
+        as a function of the nondimensional end force
+        in the isotensional ensemble, using a
+        Metropolis-Hastings Markov chain Monte Carlo calculation
+        of the ensemble average
+
+        .. math::
+            k = \frac{
+                    \left\langle e^{-\beta\Pi}\right\rangle_{p=0}^\ddagger
+                }{
+                    \left\langle e^{-\beta\Pi}\right\rangle_{p=0}
+                }.
+
+        Args:
+            p (array_like): The nondimensional end force.
+            **kwargs: Arbitrary keyword arguments.
+                Passed to :meth:`parallel_calculation`.
+
+        Returns:
+            numpy.ndarray: The nondimensional forward reaction rate.
+
+        """
+        p = self.np_array(p)
+        k_isotensional = np.zeros(p.shape)
+        for i, p_i in enumerate(p):
+            self.beta_e = lambda vs: self.beta_Pi(p_i, vs[0], vs[1:])
+            v_guess, s_guess = self.minimize_beta_Pi(p_i)[1:]
+            k_isotensional[i] = self.parallel_calculation(
+                self.k_isotensional_monte_carlo_serial,
+                np.concatenate((v_guess, s_guess[:, 0])),
+                **kwargs
+            )
+        return k_isotensional