scale anharmonic entropy

Author: bobmyhill

burnman/eos/__init__.py

@@ -18,12 +18,6 @@
 from .mie_grueneisen_debye import MGD2, MGD3
 from .slb import SLB2, SLB3, SLB3Conductive
 from .modular_mie_grueneisen_debye import ModularMGD
-from .debye_temperature_models import (
-    DebyeTemperatureModelBase,
-    SLB,
-    PowerLawGammaSimple,
-    PowerLawGamma,
-)
 from .modified_tait import MT
 from .hp import HP_TMT
 from .hp import HP_TMTL

burnman/eos/anharmonic_debye_pade.py

@@ -153,7 +153,7 @@ def _nondimensional_helmholtz_energy(T, debye_T):
     :rtype: float
     """
     t = T / debye_T
-    return _helmholtz_pade(t)
+    return _helmholtz_pade(t) * debye_T
 
 
 @jit(nopython=True)
@@ -170,7 +170,7 @@ def _nondimensional_entropy(T, debye_T):
     :rtype: float
     """
     t = T / debye_T
-    return -_dhelmholtzdt_pade(t) / debye_T
+    return -_dhelmholtzdt_pade(t)
 
 
 @jit(nopython=True)
@@ -187,7 +187,7 @@ def _nondimensional_heat_capacity(T, debye_T):
     :rtype: float
     """
     t = T / debye_T
-    return -t * _d2helmholtzdt2_pade(t) / debye_T
+    return -t * _d2helmholtzdt2_pade(t)
 
 
 @jit(nopython=True)
@@ -205,7 +205,7 @@ def _nondimensional_dhelmholtz_dTheta(T, debye_T):
     :rtype: float
     """
     t = T / debye_T
-    return -_dhelmholtzdt_pade(t) * t / debye_T
+    return _helmholtz_pade(t) - _dhelmholtzdt_pade(t) * t
 
 
 @jit(nopython=True)
@@ -223,7 +223,7 @@ def _nondimensional_d2helmholtz_dTheta2(T, debye_T):
     :rtype: float
     """
     t = T / debye_T
-    return t * (t * _d2helmholtzdt2_pade(t) + 2.0 * _dhelmholtzdt_pade(t)) / debye_T**2
+    return t * t * _d2helmholtzdt2_pade(t) / debye_T
 
 
 @jit(nopython=True)
@@ -241,7 +241,7 @@ def _nondimensional_dentropy_dTheta(T, debye_T):
     :rtype: float
     """
     t = T / debye_T
-    return (_d2helmholtzdt2_pade(t) * t + _dhelmholtzdt_pade(t)) / debye_T**2
+    return _d2helmholtzdt2_pade(t) * t / debye_T
 
 
 class AnharmonicDebyePade:
@@ -256,17 +256,22 @@ class AnharmonicDebyePade:
     :math:`A = a_{anh} * (V/V_0)^{m_{anh}}`, with both :math:`a_{anh}`
     and :math:`m_{anh}` being parameters of the model.
     The term :math:`F_a` is calculated using the 3-5 Pade approximant to
-    the function: :math:`\\int_0^x (E_{D}/3nR) dt / x^{4}`, then
+    the function: :math:`\\int_0^x E_{D}(T)/3nR dT / x^{4}`, then
     post-multiplied by :math:`x^{4}`.
 
     The :math:`E_{D}` term inside the integral is the thermal energy of a
-    Debye solid per mole of atoms. This expression is chosen because it
-    matches the behaviour of the anharmonic contribution to the entropy
-    at low and high temperatures - i.e., it is equal to zero at low temperature
-    (with all derivatives also equal to zero) and linear at high temperature.
+    Debye solid per mole of atoms where the Debye temperature is 1 K.
+    This expression is chosen because it matches the behaviour of the
+    anharmonic contribution to the entropy at low and high temperatures;
+    i.e., it is equal to zero at low temperature (with all derivatives
+    also equal to zero) and linear at high temperature.
     See Figure 3 in Oganov and Dorogokupets
     (2004; dx.doi.org/10.1088/0953-8984/16/8/018).
 
+    The model is such that the entropy at the Debye temperature is equal to
+    -A * debye.thermal_energy(T=1., debye_T=1., n=1.) / (3R), which is
+    roughly -0.6744 J/K/mol when A = 1. Note the negative sign.
+
     :return: _description_
     :rtype: _type_
     """

docs/eos.rst

@@ -1,9 +1,10 @@
-Equations of state
-==================
+Isothermal Equations of state
+=============================
 
 Base class
 ----------
-.. autoclass:: burnman.eos.EquationOfState
+
+.. autoclass:: burnman.eos.IsothermalEquationOfState
 
 Murnaghan
 ---------
@@ -49,6 +50,19 @@ Reciprocal K-prime
 
 .. autoclass:: burnman.eos.RKprime
 
+SPOCK
+-----
+
+.. autoclass:: burnman.eos.SPOCK
+
+Thermal Equations of state
+==========================
+
+Base class
+^^^^^^^^^^
+
+.. autoclass:: burnman.eos.EquationOfState
+
 Stixrude and Lithgow-Bertelloni Formulation
 -------------------------------------------
 
@@ -68,7 +82,7 @@ SLB3
 .. autoclass:: burnman.eos.SLB3
 
 Mie-Grüneisen-Debye
--------------------------------------------
+-------------------
 
 Base class
 ^^^^^^^^^^
@@ -85,6 +99,23 @@ MGD3
 
 .. autoclass:: burnman.eos.MGD3
 
+Modular Mie-Grüneisen-Debye
+---------------------------
+
+.. autoclass:: burnman.eos.ModularMGD
+
+Debye Temperature Models
+^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autoclass:: burnman.eos.debye_temperature_models.DebyeTemperatureModelBase
+
+.. autoclass:: burnman.eos.debye_temperature_models.SLB
+
+.. autoclass:: burnman.eos.debye_temperature_models.PowerLawGammaSimple
+
+.. autoclass:: burnman.eos.debye_temperature_models.PowerLawGamma
+
+
 Modified Tait
 -------------
 

tests/test_anharmonic.py

@@ -1,12 +1,14 @@
 import unittest
 from util import BurnManTest
 
+import burnman
 from burnman.eos.anharmonic_debye_pade import AnharmonicDebyePade
 from burnman.eos.debye_temperature_models import SLB as theta_SLB
-
+from burnman.eos import debye
+from burnman.constants import gas_constant
 
 params = {
-    "a_anh": 1500.0,
+    "a_anh": 1.5,
     "V_0": 2.0e-6,
     "m_anh": 3.0,
     "Debye_0": 1000.0,
@@ -67,6 +69,16 @@ def test_derivatives(self):
         self.assertAlmostEqual(alphaK_T / 1.0e3, dSdV / 1.0e3, places=6)
         self.assertAlmostEqual(alphaK_T / 1.0e3, dPdT / 1.0e3, places=6)
 
+    def test_standard_entropy(self):
+        V = params["V_0"]
+        T = params["debye_temperature_model"](Vrel=1.0, params=params)
+        self.assertAlmostEqual(T, params["Debye_0"], places=6)
+
+        S = AnharmonicDebyePade.entropy(T, V, params) / params["a_anh"]
+        S1 = -debye.thermal_energy(1.0, 1.0, 1.0) / gas_constant / 3.0
+        self.assertAlmostEqual(S, S1, places=4)
+        self.assertAlmostEqual(S, -0.6744, places=4)
+
 
 if __name__ == "__main__":
     unittest.main()