diff --git a/src/polykin/thermo/acm/base.py b/src/polykin/thermo/acm/base.py
index bd8fcc3..7749128 100644
--- a/src/polykin/thermo/acm/base.py
+++ b/src/polykin/thermo/acm/base.py
@@ -4,12 +4,12 @@
 
 from abc import ABC, abstractmethod
 
-import numpy as np
-from numpy import log, sqrt
+from numpy import dot, log, sqrt
 from scipy.constants import R
 
 from polykin.utils.math import eps
 from polykin.utils.types import FloatVector
+from polykin.math import derivative_centered
 
 
 class ActivityCoefficientModel(ABC):
@@ -31,12 +31,31 @@ def Dgmix(self, T: float, x: FloatVector) -> float:
         float
             Molar Gibbs energy of mixing. Unit = J/mol.
         """
-        return self.Dhmix(T, x) - T*self.Dsmix(T, x)
+        return self.gE(T, x) - T*self._Dsmix_ideal(T, x)
+
+    def Dhmix(self, T: float, x: FloatVector) -> float:
+        r"""Molar enthalpy of mixing, $\Delta h_{mix}$.
+
+        $$ \Delta h_{mix} = h^{E} $$
+
+        Parameters
+        ----------
+        T : float
+            Temperature. Unit = K.
+        x : FloatVector
+            Mole fractions of all components. Unit = mol/mol.
+
+        Returns
+        -------
+        float
+            Molar enthalpy of mixing. Unit = J/mol.
+        """
+        return self.hE(T, x)
 
     def Dsmix(self, T: float, x: FloatVector) -> float:
         r"""Molar entropy of mixing, $\Delta s_{mix}$.
 
-        $$ \Delta s_{mix} = s^{E} + R \sum_i {x_i \ln{x_i}} $$
+        $$ \Delta s_{mix} = s^{E} - R \sum_i {x_i \ln{x_i}} $$
 
         Parameters
         ----------
@@ -50,12 +69,12 @@ def Dsmix(self, T: float, x: FloatVector) -> float:
         float
             Molar entropy of mixing. Unit = J/(mol·K).
         """
-        return self.sE(T, x) + R*np.dot(x, log(x))
+        return self.sE(T, x) + self._Dsmix_ideal(T, x)
 
-    def Dhmix(self, T: float, x: FloatVector) -> float:
-        r"""Molar enthalpy of mixing, $\Delta h_{mix}$.
+    def _Dsmix_ideal(self, T: float, x: FloatVector) -> float:
+        r"""Molar entropy of mixing of ideal solution, $\Delta s_{mix}^{ideal}$.
 
-        $$ \Delta h_{mix} = h^{E} $$
+        $$ \Delta s_{mix}^{ideal} = - R \sum_i {x_i \ln{x_i}} $$
 
         Parameters
         ----------
@@ -67,9 +86,10 @@ def Dhmix(self, T: float, x: FloatVector) -> float:
         Returns
         -------
         float
-            Molar enthalpy of mixing. Unit = J/mol.
+            Molar entropy of mixing. Unit = J/(mol·K).
         """
-        return self.hE(T, x)
+        xp = x[x > 0]
+        return -R*dot(xp, log(xp))
 
     def hE(self, T: float, x: FloatVector) -> float:
         r"""Molar excess enthalpy, $h^{E}$.
@@ -107,10 +127,7 @@ def sE(self, T: float, x: FloatVector) -> float:
         float
             Molar excess entropy. Unit = J/(mol·K).
         """
-        dT = 2*sqrt(eps)*T
-        gE_plus = self.gE(T + dT, x)
-        gE_minus = self.gE(T - dT, x)
-        return (gE_minus - gE_plus)/(2*dT)
+        return derivative_centered(lambda t: self.gE(t, x), T)[0]
 
     def a(self, T: float, x: FloatVector) -> FloatVector:
         r"""Activities, $a_i$.