From c9746ba9628b333ff07036ba9f0886599b450fee Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Tue, 5 Mar 2024 10:16:45 +0000 Subject: [PATCH] assert Phi=0 at standard state --- burnman/classes/anisotropicmineral.py | 6 ++++++ contrib/anisotropic_eos/cubic_fitting_second_form.py | 8 ++++---- .../anisotropic_eos/cubic_fitting_second_form_atherm.py | 4 ++-- .../cubic_fitting_second_form_only_therm.py | 4 ++-- .../anisotropic_eos/orthorhombic_fitting_second_form.py | 8 ++++---- 5 files changed, 18 insertions(+), 12 deletions(-) diff --git a/burnman/classes/anisotropicmineral.py b/burnman/classes/anisotropicmineral.py index 6c1d1c9a0..694ced167 100644 --- a/burnman/classes/anisotropicmineral.py +++ b/burnman/classes/anisotropicmineral.py @@ -204,6 +204,12 @@ def __init__( self.anisotropic_params = anisotropic_parameters self.psi_function = psi_function + Psi_Voigt0 = self.psi_function(0., 0., self.anisotropic_params)[0] + if not np.all(Psi_Voigt0 == 0.): + raise ValueError("All elements of Psi should evaluate to zero at " + "standard state. The current array evaluates to: " + f"{Psi_Voigt0}") + # cell_vectors is the transpose of the cell tensor M self.cell_vectors_0 = cell_parameters_to_vectors(cell_parameters) diff --git a/contrib/anisotropic_eos/cubic_fitting_second_form.py b/contrib/anisotropic_eos/cubic_fitting_second_form.py index d55a0dfbd..955e5e2f1 100644 --- a/contrib/anisotropic_eos/cubic_fitting_second_form.py +++ b/contrib/anisotropic_eos/cubic_fitting_second_form.py @@ -61,8 +61,8 @@ def psi_func(f, Pth, params): Psi = ( 0.0 + params["a"] * f - + params["b_1"] * np.exp(params["c_1"] * f) - + params["b_2"] * np.exp(params["c_2"] * f) + + params["b_1"] * (np.exp(params["c_1"] * f) - 1.) + + params["b_2"] * (np.exp(params["c_2"] * f) - 1.) ) Psi = Psi + Pth / 1.0e9 * ( params["b_3"] * np.exp(params["c_3"] * f) @@ -379,8 +379,8 @@ def cubic_misfit(x, imin): f'the isotropic model: $V_0$: {m.params["V_0"]*1.e6:.5f} cm$^3$/mol, ' f'$K_0$: {m.params["K_0"]/1.e9:.5f} GPa, ' f'$K\'_0$: {m.params["Kprime_0"]:.5f}, ' - f'$\Theta_0$: {m.params["Debye_0"]:.5f} K, ' - f'$\gamma_0$: {m.params["grueneisen_0"]:.5f}, ' + f'$\\Theta_0$: {m.params["Debye_0"]:.5f} K, ' + f'$\\gamma_0$: {m.params["grueneisen_0"]:.5f}, ' f'and $q_0$: {m.params["q_0"]:.5f}.' ) diff --git a/contrib/anisotropic_eos/cubic_fitting_second_form_atherm.py b/contrib/anisotropic_eos/cubic_fitting_second_form_atherm.py index 6a48fb046..e6534283c 100644 --- a/contrib/anisotropic_eos/cubic_fitting_second_form_atherm.py +++ b/contrib/anisotropic_eos/cubic_fitting_second_form_atherm.py @@ -53,8 +53,8 @@ def psi_func(f, Pth, params): Psi = ( 0.0 + params["a"] * f - + params["b_1"] * np.exp(params["c_1"] * f) - + params["b_2"] * np.exp(params["c_2"] * f) + + params["b_1"] * (np.exp(params["c_1"] * f) - 1.) + + params["b_2"] * (np.exp(params["c_2"] * f) - 1.) ) dPsidPth = 0.0 * params["b_1"] return (Psi, dPsidf, dPsidPth) diff --git a/contrib/anisotropic_eos/cubic_fitting_second_form_only_therm.py b/contrib/anisotropic_eos/cubic_fitting_second_form_only_therm.py index 2e8bc606a..d808eec0f 100644 --- a/contrib/anisotropic_eos/cubic_fitting_second_form_only_therm.py +++ b/contrib/anisotropic_eos/cubic_fitting_second_form_only_therm.py @@ -114,8 +114,8 @@ def psi_func(f, Pth, params): Psi = ( 0.0 + params["a"] * f - + params["b_1"] * np.exp(params["c_1"] * f) - + params["b_2"] * np.exp(params["c_2"] * f) + + params["b_1"] * (np.exp(params["c_1"] * f) - 1.) + + params["b_2"] * (np.exp(params["c_2"] * f) - 1.) ) Psi = Psi + Pth / 1.0e9 * ( params["b_3"] * np.exp(params["c_3"] * f) diff --git a/contrib/anisotropic_eos/orthorhombic_fitting_second_form.py b/contrib/anisotropic_eos/orthorhombic_fitting_second_form.py index 5615ffcf5..c23d9efaa 100644 --- a/contrib/anisotropic_eos/orthorhombic_fitting_second_form.py +++ b/contrib/anisotropic_eos/orthorhombic_fitting_second_form.py @@ -100,8 +100,8 @@ def psi_func(f, Pth, params): Psi = ( 0.0 + params["a"] * f - + params["b_1"] * np.exp(params["c_1"] * f) - + params["b_2"] * np.exp(params["c_2"] * f) + + params["b_1"] * (np.exp(params["c_1"] * f) - 1.) + + params["b_2"] * (np.exp(params["c_2"] * f) - 1.) + params["d"] * Pth / 1.0e9 ) dPsidPth = params["d"] / 1.0e9 @@ -393,8 +393,8 @@ def orthorhombic_misfit(x, imin): f'the isotropic model: $V_0$: {m.params["V_0"]*1.e6:.5f} cm$^3$/mol, ' f'$K_0$: {m.params["K_0"]/1.e9:.5f} GPa, ' f'$K\'_0$: {m.params["Kprime_0"]:.5f}, ' - f'$\Theta_0$: {m.params["Debye_0"]:.5f} K, ' - f'$\gamma_0$: {m.params["grueneisen_0"]:.5f}, ' + f'$\\Theta_0$: {m.params["Debye_0"]:.5f} K, ' + f'$\\gamma_0$: {m.params["grueneisen_0"]:.5f}, ' f'and $q_0$: {m.params["q_0"]:.5f}.' )