From b06d353596413ecf74481e24c7018207509adac6 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Sun, 6 Oct 2024 16:58:04 +0100 Subject: [PATCH] update stishovite model files --- burnman/classes/material.py | 3 +- contrib/anisotropic_stishovite/Readme.md | 17 +- .../data/Nishihara_et_al_2005_stv.dat | 47 +++ .../data/Wang_2012_stv.dat | 62 ---- .../data/Wang_et_al_2012_stv.dat | 62 ++++ .../model_output/covariance_matrix.dat | 37 ++ .../anisotropic_stishovite/stishovite_data.py | 47 ++- .../stishovite_fit_eos.py | 144 +++++-- .../stishovite_model.py | 43 ++- .../stishovite_model_Carpenter_2000.py | 303 +++++++++++++++ .../stishovite_model_Zhang_2021.py | 350 ++++++++++++++++++ .../stishovite_model_covariance_matrix.py | 98 +++++ .../stishovite_model_plots.py | 141 +++++-- .../stishovite_parameters.py | 136 ++++--- 14 files changed, 1293 insertions(+), 197 deletions(-) create mode 100644 contrib/anisotropic_stishovite/data/Nishihara_et_al_2005_stv.dat delete mode 100644 contrib/anisotropic_stishovite/data/Wang_2012_stv.dat create mode 100644 contrib/anisotropic_stishovite/data/Wang_et_al_2012_stv.dat create mode 100644 contrib/anisotropic_stishovite/model_output/covariance_matrix.dat create mode 100644 contrib/anisotropic_stishovite/stishovite_model_Carpenter_2000.py create mode 100644 contrib/anisotropic_stishovite/stishovite_model_Zhang_2021.py create mode 100644 contrib/anisotropic_stishovite/stishovite_model_covariance_matrix.py diff --git a/burnman/classes/material.py b/burnman/classes/material.py index 3e7c36ac..462ea7de 100644 --- a/burnman/classes/material.py +++ b/burnman/classes/material.py @@ -293,9 +293,10 @@ def evaluate(self, vars_list, pressures, temperatures, molar_fractions=None): assert pressures.shape == temperatures.shape + first_index = list(np.ndenumerate(temperatures))[0][0] if molar_fractions is not None: molar_fractions = np.array(molar_fractions) - self.set_composition(molar_fractions[0]) + self.set_composition(molar_fractions[first_index]) assert temperatures.shape == molar_fractions.shape[:-1] # First, check the output types of all the requested variables: diff --git a/contrib/anisotropic_stishovite/Readme.md b/contrib/anisotropic_stishovite/Readme.md index dc6de24e..4423d3a1 100644 --- a/contrib/anisotropic_stishovite/Readme.md +++ b/contrib/anisotropic_stishovite/Readme.md @@ -1,5 +1,5 @@ ## Companion code to -# A self-consistent Landau formalism to model instantaneous and time-dependent elastic softening and second order phase transitions in anisotropic phases, with application to stishovite +# A model of elastic softening and second order phase transitions in anisotropic phases, with application to stishovite and post-stishovite ## R. Myhill The python scripts contained in this directory accompany the paper @@ -34,3 +34,18 @@ based on the experimental data. stishovite_model_plots.py ------------------------- This script creates all of the plots presented in the paper. + +stishovite_model_covariance_matrix.py +------------------------------------- +This script outputs the model uncertainties and correlation matrix +from the covariance matrix. + +stishovite_model_Carpenter_2000.py +---------------------------------- +Plots the Gibbs energy, volumes, Reuss bulk modulus and elastic moduli +from the model of Carpenter et al. (2000) + +stishovite_model_Zhang_2021.py +------------------------------ +Plots the Gibbs energy, volumes, Reuss bulk modulus and elastic moduli +from the model of Zhang et al. (2021) diff --git a/contrib/anisotropic_stishovite/data/Nishihara_et_al_2005_stv.dat b/contrib/anisotropic_stishovite/data/Nishihara_et_al_2005_stv.dat new file mode 100644 index 00000000..fdfaf8e8 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Nishihara_et_al_2005_stv.dat @@ -0,0 +1,47 @@ +# This paper was solely studying stv. +# P T a a_err c c_err V V_err +0 300 4.178 0.001 2.668 0.001 46.56 0.02 +4.76 300 4.152 0.001 2.658 0.001 45.82 0.02 +7.51 300 4.139 0.001 2.653 0.001 45.46 0.01 +10.7 300 4.124 0.001 2.647 0.001 45.03 0.02 +13.55 300 4.111 0.001 2.641 0.001 44.64 0.01 +14.63 300 4.106 0.001 2.64 0.001 44.51 0.02 +16.19 300 4.1 0.001 2.637 0.001 44.32 0.02 +17.87 300 4.093 0.001 2.634 0.001 44.11 0.02 +19.53 300 4.088 0.001 2.629 0.001 43.92 0.02 +21.23 300 4.08 0.001 2.626 0.001 43.73 0.02 +5.03 473 4.155 0.001 2.658 0.001 45.88 0.02 +7.79 473 4.141 0.001 2.655 0.001 45.53 0.02 +10.95 473 4.127 0.001 2.647 0.001 45.08 0.02 +13.82 473 4.113 0.001 2.642 0.001 44.69 0.01 +14.91 473 4.107 0.001 2.641 0.001 44.55 0.02 +16.39 473 4.101 0.001 2.638 0.001 44.37 0.02 +18.1 473 4.095 0.001 2.634 0.001 44.17 0.02 +19.77 473 4.088 0.001 2.63 0.001 43.95 0.02 +21.47 473 4.082 0.001 2.627 0.001 43.78 0.03 +5.49 673 4.158 0.001 2.659 0.001 45.97 0.01 +8.23 673 4.144 0.001 2.655 0.001 45.59 0.01 +11.39 673 4.129 0.001 2.648 0.001 45.16 0.02 +14.19 673 4.115 0.001 2.643 0.001 44.76 0.02 +15.29 673 4.11 0.001 2.642 0.001 44.61 0.02 +16.81 673 4.104 0.001 2.639 0.001 44.44 0.02 +18.44 673 4.097 0.001 2.635 0.001 44.23 0.02 +20.07 673 4.091 0.001 2.631 0.001 44.02 0.02 +21.75 673 4.084 0.001 2.628 0.001 43.82 0.02 +6.03 873 4.161 0.001 2.661 0.001 46.08 0.01 +8.76 873 4.147 0.001 2.656 0.001 45.69 0.01 +11.9 873 4.132 0.001 2.65 0.001 45.25 0.02 +14.66 873 4.118 0.001 2.644 0.001 44.85 0.01 +15.73 873 4.112 0.001 2.643 0.001 44.69 0.02 +17.25 873 4.106 0.001 2.64 0.001 44.51 0.02 +18.83 873 4.1 0.001 2.636 0.001 44.31 0.01 +20.46 873 4.093 0.001 2.632 0.001 44.09 0.01 +22.09 873 4.086 0.001 2.63 0.001 43.91 0.02 +9.31 1073 4.151 0.001 2.658 0.001 45.8 0.01 +12.37 1073 4.135 0.001 2.651 0.001 45.33 0.02 +15.01 1073 4.121 0.001 2.645 0.001 44.92 0.01 +16.14 1073 4.115 0.001 2.644 0.001 44.77 0.01 +17.64 1073 4.109 0.001 2.641 0.001 44.59 0.03 +19.24 1073 4.103 0.001 2.637 0.001 44.4 0.03 +20.87 1073 4.096 0.001 2.632 0.001 44.15 0.02 +22.47 1073 4.088 0.001 2.631 0.001 43.97 0.02 diff --git a/contrib/anisotropic_stishovite/data/Wang_2012_stv.dat b/contrib/anisotropic_stishovite/data/Wang_2012_stv.dat deleted file mode 100644 index 37ea7d9a..00000000 --- a/contrib/anisotropic_stishovite/data/Wang_2012_stv.dat +++ /dev/null @@ -1,62 +0,0 @@ - -# This paper was almost solely studying stv. -# Just a few post-stv data points, not in this table -# (see their Fig 2) -# Data Number T (K) P (GPa) Stishovite Au V (Å3) -a (Å) c (Å) V (Å3) -M1034007 1300 19.32(10) 4.1145(3) 2.6464(4) 44.80(2) 63.84(2) -M1034008 1100 18.9(2) 4.1105(3) 2.6443(2) 44.68(2) 63.57(3) -M1034009 900 18.5(2) 4.1062(4) 2.6422(3) 44.55(2) 63.30(4) -M1034010 700 17.9(2) 4.1031(2) 2.6390(2) 44.43(3) 63.08(3) -M1034011 500 17.5(2) 4.1002(4) 2.6377(4) 44.34(3) 62.84(2) -M1034012 300 16.85(12) 4.0987(3) 2.6359(2) 44.28(2) 62.66(2) -M1036003 1300 25.6(2) 4.0866(2) 2.6345(3) 44.00(2) 62.24(2) -M1036004 1100 25.02(2) 4.0828(3) 2.6337(4) 43.90(2) 62.04(2) -M1036005 900 24.6(3) 4.0792(6) 2.6304(9) 43.77(4) 61.82(4) -M1036006 700 24.1(2) 4.07704(5) 2.62788(8) 43.681(4) 61.64(2) -M1036007 500 23.4(2) 4.0751(4) 2.6235(9) 43.60(5) 61.48(3) -M1036008 300 23.05(10) 4.0743(4) 2.6237(5) 43.55(4) 61.28(1) -M1036012 1300 30.9(2) 4.0630(2) 2.6260(4) 43.35(4) 61.06(3) -M1036013 1100 30.3(2) 4.0600(4) 2.6232(3) 43.24(2) 60.90(2) -M1036014 900 29.98(2) 4.0591(4) 2.6197(3) 43.16(3) 60.70(3) -M1036015 700 29.5(2) 4.0545(3) 2.6192(2) 43.06(3) 60.5(2) -M1036016 500 29.1(3) 4.0538(4) 2.6181(5) 43.02(2) 60.36(3) -M1036017 300 28.5(2) 4.0513(7) 2.6167(5) 42.95(4) 60.22(2) -M914003 1700 38.5(4) 4.0462(4) 2.6211(6) 42.91(3) 60.08(4) -M914005 1500 37.7(3) 4.0425(4) 2.6193(7) 42.80(3) 59.98(3) -M914006 1300 37.4(4) 4.0392(5) 2.6172(8) 42.70(4) 59.80(4) -M914007 1100 36.7(4) 4.0365(6) 2.6155(10) 42.62(5) 59.68(5) -M914008 900 36.4(4) 4.0340(7) 2.6135(10) 42.53(5) 59.50(4) -M914009 700 35.9(3) 4.0322(5) 2.6121(9) 42.47(4) 59.36(3) -M914010 500 35.1(2) 4.0308(4) 2.6109(7) 42.42(3) 59.28(3) -M914011 300 34.2(3) 4.0302(4) 2.6105(7) 42.40(3) 59.21(3) -M942004 1700 40.4(4) 4.0406(4) 2.6177(7) 42.74(3) 59.74(4) -M942005 1500 40.1(4) 4.0362(3) 2.6153(5) 42.60(2) 59.54(4) -M942006 1300 39.3(3) 4.0339(4) 2.6125(7) 42.51(3) 59.46(3) -M942007 1100 39.0(2) 4.0310(3) 2.6106(6) 42.42(3) 59.27(2) -M942008 900 38.5(2) 4.0290(5) 2.6085(9) 42.34(4) 59.13(2) -M942009 700 38.0(2) 4.0263(4) 2.6070(7) 42.26(3) 59.00(2) -M942010 500 37.4(2) 4.0247(2) 2.6062(4) 42.21(2) 58.89(3) -M942011 300 36.8(3) 4.0234(4) 2.6058(8) 42.18(4) 58.77(3) -M915004 1700 46.1(5) 4.0201(5) 2.6081(8) 42.15(4) 58.75(5) -M915005 1500 45.9(2) 4.0168(4) 2.6056(7) 42.04(4) 58.56(2) -M915006 1300 45.6(4) 4.0143(3) 2.6039(5) 41.96(2) 58.40(4) -M915007 1100 44.88(11) 4.0117(3) 2.6024(6) 41.88(3) 58.31(1) -M915008 900 44.5(2) 4.0095(3) 2.6010(5) 41.81(2) 58.17(2) -M915009 700 43.9(2) 4.0079(2) 2.5993(4) 41.75(2) 58.07(2) -M915010 500 43.3(2) 4.0064(4) 2.5988(8) 41.71(4) 57.97(4) -M915011 300 42.5(4) 4.0058(3) 2.5977(5) 41.68(2) 57.89(5) -M941009 1700 48.1(3) 4.01369(10) 2.60407(10) 41.951(14) 58.41(2) -M941010 1500 47.8(4) 4.01017(10) 2.60158(10) 41.837(14) 58.26(3) -M941011 1300 47.4(2) 4.0063(3) 2.5999(6) 41.73(3) 58.12(4) -M941012 1100 47.3(2) 4.0036(2) 2.5977(4) 41.64(2) 57.94(2) -M941013 900 47.10(13) 4.0004(3) 2.5954(5) 41.54(2) 57.78(2) -M941014 700 47.1(2) 3.9972(4) 2.5939(6) 41.45(3) 57.60(3) -M941015 500 46.6(2) 3.9957(4) 2.5924(6) 41.39(3) 57.48(2) -M941016 300 46.16(12) 3.9948(3) 2.5916(6) 41.36(3) 57.373(10) -M922006 1300 54.4(3) 3.9863(4) 2.5901(7) 41.16(3) 57.10(4) -M922007 1100 54.5(3) 3.9833(4) 2.5877(7) 41.06(3) 56.91(3) -M922008 900 54.4(2) 3.9801(7) 2.5857(12) 40.96(6) 56.76(3) -M922009 700 54.0(3) 3.9781(4) 2.5839(6) 40.89(3) 56.64(2) -M922010 500 53.7(4) 3.9758(7) 2.5822(13) 40.82(6) 56.53(3) -M922011 300 53.4(4) 3.9750(7) 2.5819(12) 40.80(5) 56.41(3) diff --git a/contrib/anisotropic_stishovite/data/Wang_et_al_2012_stv.dat b/contrib/anisotropic_stishovite/data/Wang_et_al_2012_stv.dat new file mode 100644 index 00000000..df9b8a95 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Wang_et_al_2012_stv.dat @@ -0,0 +1,62 @@ + +# This paper was almost solely studying stv. +# Just a few post-stv data points, not in this table +# (see their Fig 2) +# Data Number +# Analysis T P P_err a a_err c c_err V V_err Vau Vau_err +1034007 1300 19.32 0.1 4.1145 0.0003 2.6464 0.0004 44.8 0.02 63.84 0.02 +1034008 1100 18.9 0.2 4.1105 0.0003 2.6443 0.0002 44.68 0.02 63.57 0.03 +1034009 900 18.5 0.2 4.1062 0.0004 2.6422 0.0003 44.55 0.02 63.3 0.04 +1034010 700 17.9 0.2 4.1031 0.0002 2.639 0.0002 44.43 0.03 63.08 0.03 +1034011 500 17.5 0.2 4.1002 0.0004 2.6377 0.0004 44.34 0.03 62.84 0.02 +1034012 300 16.85 0.12 4.0987 0.0003 2.6359 0.0002 44.28 0.02 62.66 0.02 +1036003 1300 25.6 0.2 4.0866 0.0002 2.6345 0.0003 44 0.02 62.24 0.02 +1036004 1100 25.02 0.2 4.0828 0.0003 2.6337 0.0004 43.9 0.02 62.04 0.02 +1036005 900 24.6 0.3 4.0792 0.0006 2.6304 0.0009 43.77 0.04 61.82 0.04 +1036006 700 24.1 0.2 4.07704 0.0005 2.62788 0.0008 43.681 0.04 61.64 0.02 +1036007 500 23.4 0.2 4.0751 0.0004 2.6235 0.0009 43.6 0.05 61.48 0.03 +1036008 300 23.05 0.1 4.0743 0.0004 2.6237 0.0005 43.55 0.04 61.28 0.01 +1036012 1300 30.9 0.2 4.063 0.0002 2.626 0.0004 43.35 0.04 61.06 0.03 +1036013 1100 30.3 0.2 4.06 0.0004 2.6232 0.0003 43.24 0.02 60.9 0.02 +1036014 900 29.98 0.2 4.0591 0.0004 2.6197 0.0003 43.16 0.03 60.7 0.03 +1036015 700 29.5 0.2 4.0545 0.0003 2.6192 0.0002 43.06 0.03 60.5 0.02 +1036016 500 29.1 0.3 4.0538 0.0004 2.6181 0.0005 43.02 0.02 60.36 0.03 +1036017 300 28.5 0.2 4.0513 0.0007 2.6167 0.0005 42.95 0.04 60.22 0.02 +914003 1700 38.5 0.4 4.0462 0.0004 2.6211 0.0006 42.91 0.03 60.08 0.04 +914005 1500 37.7 0.3 4.0425 0.0004 2.6193 0.0007 42.8 0.03 59.98 0.03 +914006 1300 37.4 0.4 4.0392 0.0005 2.6172 0.0008 42.7 0.04 59.8 0.04 +914007 1100 36.7 0.4 4.0365 0.0006 2.6155 0.001 42.62 0.05 59.68 0.05 +914008 900 36.4 0.4 4.034 0.0007 2.6135 0.001 42.53 0.05 59.5 0.04 +914009 700 35.9 0.3 4.0322 0.0005 2.6121 0.0009 42.47 0.04 59.36 0.03 +914010 500 35.1 0.2 4.0308 0.0004 2.6109 0.0007 42.42 0.03 59.28 0.03 +914011 300 34.2 0.3 4.0302 0.0004 2.6105 0.0007 42.4 0.03 59.21 0.03 +942004 1700 40.4 0.4 4.0406 0.0004 2.6177 0.0007 42.74 0.03 59.74 0.04 +942005 1500 40.1 0.4 4.0362 0.0003 2.6153 0.0005 42.6 0.02 59.54 0.04 +942006 1300 39.3 0.3 4.0339 0.0004 2.6125 0.0007 42.51 0.03 59.46 0.03 +942007 1100 39 0.2 4.031 0.0003 2.6106 0.0006 42.42 0.03 59.27 0.02 +942008 900 38.5 0.2 4.029 0.0005 2.6085 0.0009 42.34 0.04 59.13 0.02 +942009 700 38 0.2 4.0263 0.0004 2.607 0.0007 42.26 0.03 59 0.02 +942010 500 37.4 0.2 4.0247 0.0002 2.6062 0.0004 42.21 0.02 58.89 0.03 +942011 300 36.8 0.3 4.0234 0.0004 2.6058 0.0008 42.18 0.04 58.77 0.03 +915004 1700 46.1 0.5 4.0201 0.0005 2.6081 0.0008 42.15 0.04 58.75 0.05 +915005 1500 45.9 0.2 4.0168 0.0004 2.6056 0.0007 42.04 0.04 58.56 0.02 +915006 1300 45.6 0.4 4.0143 0.0003 2.6039 0.0005 41.96 0.02 58.4 0.04 +915007 1100 44.88 0.11 4.0117 0.0003 2.6024 0.0006 41.88 0.03 58.31 0.01 +915008 900 44.5 0.2 4.0095 0.0003 2.601 0.0005 41.81 0.02 58.17 0.02 +915009 700 43.9 0.2 4.0079 0.0002 2.5993 0.0004 41.75 0.02 58.07 0.02 +915010 500 43.3 0.2 4.0064 0.0004 2.5988 0.0008 41.71 0.04 57.97 0.04 +915011 300 42.5 0.4 4.0058 0.0003 2.5977 0.0005 41.68 0.02 57.89 0.05 +941009 1700 48.1 0.3 4.01369 0.0001 2.60407 0.0001 41.951 0.014 58.41 0.02 +941010 1500 47.8 0.4 4.01017 0.0001 2.60158 0.0001 41.837 0.014 58.26 0.03 +941011 1300 47.4 0.2 4.0063 0.0003 2.5999 0.0006 41.73 0.03 58.12 0.04 +941012 1100 47.3 0.2 4.0036 0.0002 2.5977 0.0004 41.64 0.02 57.94 0.02 +941013 900 47.1 0.13 4.0004 0.0003 2.5954 0.0005 41.54 0.02 57.78 0.02 +941014 700 47.1 0.2 3.9972 0.0004 2.5939 0.0006 41.45 0.03 57.6 0.03 +941015 500 46.6 0.2 3.9957 0.0004 2.5924 0.0006 41.39 0.03 57.48 0.02 +941016 300 46.16 0.12 3.9948 0.0003 2.5916 0.0006 41.36 0.03 57.373 0.01 +922006 1300 54.4 0.3 3.9863 0.0004 2.5901 0.0007 41.16 0.03 57.1 0.04 +922007 1100 54.5 0.3 3.9833 0.0004 2.5877 0.0007 41.06 0.03 56.91 0.03 +922008 900 54.4 0.2 3.9801 0.0007 2.5857 0.0012 40.96 0.06 56.76 0.03 +922009 700 54 0.3 3.9781 0.0004 2.5839 0.0006 40.89 0.03 56.64 0.02 +922010 500 53.7 0.4 3.9758 0.0007 2.5822 0.0013 40.82 0.06 56.53 0.03 +922011 300 53.4 0.4 3.975 0.0007 2.5819 0.0012 40.8 0.05 56.41 0.03 diff --git a/contrib/anisotropic_stishovite/model_output/covariance_matrix.dat b/contrib/anisotropic_stishovite/model_output/covariance_matrix.dat new file mode 100644 index 00000000..2ceb515f --- /dev/null +++ b/contrib/anisotropic_stishovite/model_output/covariance_matrix.dat @@ -0,0 +1,37 @@ +2.969072323847245294e-05 -2.414319325043906995e-03 9.735316303185922192e-07 -2.034079978431028023e-05 1.597751501944453352e-04 5.032599281155131899e-07 -2.088564195265813286e-03 1.141858121336443432e-03 3.949357292634677550e-05 1.850758822971476200e-05 3.299016538863077687e-05 9.077296660207598603e-05 -5.137345609827280432e-16 -8.143086031774118376e-17 -4.589020455106308953e-16 -3.098661009272525452e-15 2.582109938010748615e-08 2.624612085038632830e-08 3.608749591974849653e-06 -6.937855671113378833e-05 -1.371705693477332281e-04 -2.071765376120491374e-07 -2.516517250145258906e-03 -3.928954203215390470e-05 1.279490468007478382e-06 9.068426317865086325e-08 -1.010411452316174498e-05 -1.002915587787911641e-05 -7.096411239012776800e-06 -6.166408270828315896e-05 7.954518997506867072e-06 -1.457372693470361177e-03 -1.665472713324716442e-03 5.402267720757495664e-04 1.943321742220476944e-03 -3.333863796457967564e-04 -5.468710926687100291e-05 +-2.414319325043906995e-03 4.249826663546418248e-01 -3.094884841608068937e-04 5.573670640395235942e-04 1.149432319847632812e-02 -4.490399119970610087e-05 1.665240607391086369e-01 -6.504618095173082815e-02 -2.282360652002456136e-03 -7.476555221169896813e-04 -2.364247942489634824e-03 -7.079281708105459453e-03 3.475237159535944607e-14 5.352818153943201318e-15 3.945412998590303000e-14 2.598336304972241831e-13 -2.151349431467241060e-06 -2.167558592309546433e-06 -2.808918726555411629e-04 1.185104206621878306e-02 1.657883033879654217e-02 1.543514452236958691e-05 2.136226422370665345e-01 2.137533177873366264e-03 2.287660214046896914e-03 2.979745879017003073e-04 1.676699508613893001e-03 1.732916157635252218e-03 1.298896908940160239e-03 9.282782630378589445e-03 -6.763338260323719427e-04 9.451541906012153604e-02 1.143881756135775030e-01 -3.500456106541598861e-02 -1.335569570692768659e-01 2.491333162046089483e-02 4.670416326687474098e-03 +9.735316303185922192e-07 -3.094884841608068937e-04 7.262712937098192892e-06 -9.849114812544193482e-06 -8.091894307406850438e-05 -8.384368839662062259e-08 2.757286574523949750e-04 -3.302108663012920522e-05 -7.940751942138309900e-07 4.870195465031972041e-07 -5.201307034047055083e-08 8.923622864770794612e-07 2.587147620834972909e-17 -8.072235707204832032e-18 4.295548432640169746e-18 -3.705586445323490784e-17 -3.382128028653160845e-10 -3.237635627080133576e-10 2.880580686211860440e-07 -5.630122258543856047e-05 -5.550398399081168347e-05 -1.368692423501269688e-08 -5.511266268928062019e-04 5.445097707514176965e-06 1.881736127791208387e-06 1.074217982496532016e-06 -9.573321138379633936e-07 -8.955180902632183681e-07 -8.519837545226716281e-07 -3.229239994559254898e-05 3.022377409523564049e-07 1.005144316480618646e-04 6.291873262212429902e-05 -3.646627150018448262e-05 -7.162025230758143908e-05 -2.353051652724181582e-05 -4.818973351892917084e-06 +-2.034079978431028023e-05 5.573670640395235942e-04 -9.849114812544193482e-06 3.362871704584756774e-03 7.747114374398873272e-03 -1.609643254583494532e-06 6.691987623949255913e-03 -3.052808866339202370e-03 6.414389714812255498e-05 -1.154995575737634931e-05 3.264042564485330100e-04 5.367620605078604222e-04 -1.764022828001739892e-15 -8.026163323787587436e-17 -4.187898242504952461e-15 -1.643837790116418740e-14 -2.415063062091929880e-08 -3.005924930827426758e-08 2.969477106024997500e-07 -2.599807913875601048e-03 -2.502578179730502710e-03 -3.448229919827389252e-08 9.894412705418731673e-04 1.614015226545799700e-04 1.096321509004083357e-05 7.773663567167847920e-07 6.743863193988200014e-06 6.281162398372516905e-06 4.516749168953107673e-06 -1.923508018533116436e-03 1.331924023050925137e-05 -5.849646589347031259e-03 -5.178846973454946830e-03 2.153268188959647257e-03 6.003639016901835684e-03 1.387439896431741414e-03 8.672656714053429282e-05 +1.597751501944453352e-04 1.149432319847632812e-02 -8.091894307406850438e-05 7.747114374398873272e-03 1.331325052127259045e-01 -2.029319015312089713e-05 5.925313582364208354e-02 -3.080991714587031846e-02 1.977215830163714296e-03 2.741585182335783783e-04 2.896860335651843684e-04 1.694321705139406218e-03 -4.032963359981584820e-14 -3.038866244634352617e-15 -1.406216877566180905e-14 -7.731203689154381191e-14 -3.333970888846127249e-10 -8.709526033092263374e-08 5.768850215532460621e-05 -1.808405967871253364e-02 -1.896972119286157929e-02 -3.448716959724578149e-06 3.799553542857860011e-02 1.185191496548233239e-03 1.182857261486174993e-03 5.547822027831921226e-05 1.003402092902857769e-04 1.077696970302018857e-04 8.883519413119686178e-05 -1.588843399738775614e-02 3.661169141778293729e-05 -1.294007493555653365e-01 -1.216020000364273906e-01 4.772874293545138064e-02 1.412309683815566597e-01 1.866643744552765397e-02 1.065831122547543450e-03 +5.032599281155131899e-07 -4.490399119970610087e-05 -8.384368839662062259e-08 -1.609643254583494532e-06 -2.029319015312089713e-05 1.534491722260848744e-07 -5.010148782737223420e-04 1.389084317312687150e-04 8.027025917434763237e-07 -3.086014703537749928e-07 5.402179638462372378e-07 1.348672238299724075e-06 -1.243086105691209076e-17 1.406908605920547081e-17 -5.669532384363459182e-18 -4.265439752377840368e-17 2.258055017134934322e-09 1.973339258111572039e-09 -2.312332839147157345e-07 -2.916663369361920704e-05 -2.519716699320902141e-05 6.390334241489579055e-09 -2.181936752914516461e-05 -2.885741513929382642e-06 -9.584663865549323316e-06 -5.294172336870067481e-09 -5.146313275530546443e-07 -5.539881244564768786e-07 -2.523146163460764631e-07 -1.594898464460605111e-05 8.444660670564603750e-08 -2.640464914727091680e-05 -3.561425506840147387e-05 9.788159208161594391e-06 4.150511446075701752e-05 -4.390083089488345825e-05 4.229752159922561635e-07 +-2.088564195265813286e-03 1.665240607391086369e-01 2.757286574523949750e-04 6.691987623949255913e-03 5.925313582364208354e-02 -5.010148782737223420e-04 2.104896430272050534e+00 -4.981610783074122972e-01 -3.365689227855006831e-03 6.453674787356451911e-04 -2.134808163348684597e-03 -5.585079353706507398e-03 4.763216963682577079e-14 -4.384607976308321193e-14 2.390016003294783340e-14 1.802000319076651299e-13 -7.553474453448608289e-06 -6.879771852291717019e-06 6.747976667041734062e-04 1.029989985721349172e-01 9.065041706401480437e-02 -2.040043730075414874e-05 4.725134464369972614e-02 1.027426075386230972e-02 3.051814272148788490e-02 1.825216372649425331e-05 1.755644760324465278e-03 1.870870624454052007e-03 8.708104448430653436e-04 5.720248449212836372e-02 -3.219229529387673773e-04 1.057903669450520889e-01 1.403403673864646317e-01 -3.924802079810447292e-02 -1.636564359123428292e-01 1.451439494535149144e-01 -1.738338909140980243e-03 +1.141858121336443432e-03 -6.504618095173082815e-02 -3.302108663012920522e-05 -3.052808866339202370e-03 -3.080991714587031846e-02 1.389084317312687150e-04 -4.981610783074122972e-01 2.170061849111732388e-01 1.828695319578236243e-03 3.001159231045504765e-04 1.312031061506022458e-03 3.225448143816585494e-03 -1.887583244373389328e-14 9.531646394130322699e-15 -1.395286465047703843e-14 -1.010094527683821402e-13 2.222821203875844844e-06 2.537007182150678717e-06 -6.392355604498717206e-05 -3.473801794131257253e-02 -3.257967242021848203e-02 3.079233917323632859e-06 -3.340170502487335713e-02 -5.494265105720975290e-03 -5.915567986028016229e-03 3.383003929474542214e-05 -5.605255102331018402e-04 -5.695318433638057315e-04 -2.820483943677134196e-04 -1.955720914959569248e-02 2.528818449588239310e-04 -4.322324541853273666e-02 -6.023684334790740946e-02 1.611303878009478704e-02 7.044389440078994979e-02 -4.483966495114385248e-02 -8.233919291224370057e-05 +3.949357292634677550e-05 -2.282360652002456136e-03 -7.940751942138309900e-07 6.414389714812255498e-05 1.977215830163714296e-03 8.027025917434763237e-07 -3.365689227855006831e-03 1.828695319578236243e-03 1.016967765522213269e-04 2.704638970645951757e-05 4.532050448181657753e-05 1.357918390547996671e-04 -1.592926120230101003e-15 -9.821249544289911861e-17 -7.439606753073941868e-16 -4.845003724100701487e-15 3.813363256363073580e-08 3.857585521230416037e-08 3.200685394482180146e-06 -2.112929067660789938e-04 -2.858688383419348751e-04 -2.408182879736639907e-07 -2.713414985469090809e-03 -2.657470988356218109e-05 -2.803719763090280236e-05 2.580429643797891438e-07 -1.010510257238622547e-05 -9.686357516880338036e-06 -5.814282625649671241e-06 -1.795567582710501570e-04 1.072466086869132534e-05 -4.802813362752815160e-03 -4.963194473939220967e-03 1.775867457136094804e-03 5.777840977318744395e-03 -3.405953240386222296e-04 -6.056307379456799997e-05 +1.850758822971476200e-05 -7.476555221169896813e-04 4.870195465031972041e-07 -1.154995575737634931e-05 2.741585182335783783e-04 -3.086014703537749928e-07 6.453674787356451911e-04 3.001159231045504765e-04 2.704638970645951757e-05 8.634980512679053723e-05 2.160796390805643907e-05 5.824445473519590950e-05 -3.264625985174926743e-16 -7.923909448595848997e-16 -2.869213732277017375e-16 -1.943758295631656079e-15 8.266752272149200309e-09 9.778994440663878912e-09 3.524928150674835131e-06 9.244866045590156400e-05 2.724159914307262731e-05 -1.754112212171557053e-07 -1.639613353212979532e-03 -1.838827713483549386e-05 4.881136721097782943e-05 1.113170450201948483e-06 -2.085181637611447406e-06 -1.662825435048065551e-06 -1.529648396480678433e-06 3.754504312323681382e-05 5.070866914767413610e-06 -9.756958387685016505e-04 -1.072946360077371318e-03 3.617995149859560026e-04 1.251930415222073795e-03 -5.528075146326992567e-05 -3.912271590726792503e-05 +3.299016538863077687e-05 -2.364247942489634824e-03 -5.201307034047055083e-08 3.264042564485330100e-04 2.896860335651843684e-04 5.402179638462372378e-07 -2.134808163348684597e-03 1.312031061506022458e-03 4.532050448181657753e-05 2.160796390805643907e-05 3.741988925585817429e-04 1.659927024235176589e-04 -5.879810778647083869e-16 -9.391283520432811530e-17 -7.571720499217147661e-15 -5.255400233547673271e-15 2.856267415224482285e-08 2.904112999090483503e-08 4.142496374942100873e-06 -2.570082280690285368e-04 -3.251801223006570694e-04 -2.402523089244516050e-07 -3.216676545935040574e-03 -3.961294658191956365e-05 2.107728167369185196e-06 6.525664614377412270e-07 -9.965538983237963089e-06 -9.808018743190203645e-06 -6.833063562445468623e-06 -1.870407675260386544e-04 1.107990032257984879e-05 -1.679299878847527938e-03 -1.899267739233110111e-03 6.223804919963958409e-04 2.215661331552164909e-03 -3.730348743630668195e-04 -6.364872634283788350e-05 +9.077296660207598603e-05 -7.079281708105459453e-03 8.923622864770794612e-07 5.367620605078604222e-04 1.694321705139406218e-03 1.348672238299724075e-06 -5.585079353706507398e-03 3.225448143816585494e-03 1.357918390547996671e-04 5.824445473519590950e-05 1.659927024235176589e-04 9.861110209393142770e-03 -1.891966634045982456e-15 -2.711700981994465871e-16 -2.204199540613405901e-15 -5.253591335895288997e-13 7.794795518439994591e-08 7.846052278910172042e-08 1.148661919795845201e-05 -6.454285445130538692e-04 -8.440052609354726587e-04 -6.640530256914299064e-07 -7.926725382679215490e-03 -1.011002018751854773e-04 9.225664148879019028e-06 1.051179776515564498e-06 -2.918386730194719129e-05 -2.889439043922104893e-05 -2.028183930935848183e-05 -5.055290379226907082e-04 2.773567462221827859e-05 -5.498729953080172886e-03 -6.050663216414568953e-03 2.036131827085406725e-03 7.054155462547743484e-03 -8.578412202787384210e-04 -1.606945507961290737e-04 +-5.137345609827280432e-16 3.475237159535944607e-14 2.587147620834972909e-17 -1.764022828001739892e-15 -4.032963359981584820e-14 -1.243086105691209076e-17 4.763216963682577079e-14 -1.887583244373389328e-14 -1.592926120230101003e-15 -3.264625985174926743e-16 -5.879810778647083869e-16 -1.891966634045982456e-15 2.851148476213537420e-26 1.058426796312021171e-27 1.107481097947859209e-26 7.026483273915552917e-26 -5.379740661803572684e-19 -5.402147542475392543e-19 -1.533973283924579781e-17 5.186818943221322513e-15 5.756938989258702090e-15 2.461311477771505841e-18 3.091641941796989494e-14 -1.759372847258843210e-16 1.070459744475160279e-15 9.401517695428871141e-18 1.520685182370616798e-16 1.525755060186724085e-16 8.834165471558264173e-17 4.235084795229275426e-15 -1.343201925832196203e-16 8.757869065918129794e-14 8.734305031308171818e-14 -3.234243315335766641e-14 -1.015755225261617099e-13 2.162259177191088287e-15 6.313068106737503019e-16 +-8.143086031774118376e-17 5.352818153943201318e-15 -8.072235707204832032e-18 -8.026163323787587436e-17 -3.038866244634352617e-15 1.406908605920547081e-17 -4.384607976308321193e-14 9.531646394130322699e-15 -9.821249544289911861e-17 -7.923909448595848997e-16 -9.391283520432811530e-17 -2.711700981994465871e-16 1.058426796312021171e-27 8.811082299008530012e-27 1.487801391128000374e-27 9.563873128644333457e-27 1.229211106917345516e-19 9.086668192364899074e-20 -4.005920918774881884e-17 -2.847060601679682921e-15 -2.127964203851167272e-15 1.628865061860011811e-18 8.722981696312964399e-15 -1.245990632231858332e-16 -1.035329050434101653e-15 -1.089157915915628154e-18 -1.201790912326496301e-17 -1.672240764147307434e-17 2.748078604662010185e-18 -1.462373406033409279e-15 -2.706459684488763294e-17 3.979391298362956564e-15 3.914280900809954267e-15 -1.474664675128080572e-15 -4.571334409002748931e-15 -3.184500434619874255e-15 2.915127137224032614e-16 +-4.589020455106308953e-16 3.945412998590303000e-14 4.295548432640169746e-18 -4.187898242504952461e-15 -1.406216877566180905e-14 -5.669532384363459182e-18 2.390016003294783340e-14 -1.395286465047703843e-14 -7.439606753073941868e-16 -2.869213732277017375e-16 -7.571720499217147661e-15 -2.204199540613405901e-15 1.107481097947859209e-26 1.487801391128000374e-27 1.670781716139887458e-25 7.303558544180336786e-26 -3.863256970748645692e-19 -3.853063443221715582e-19 -6.025552387684062132e-17 4.625905237017853671e-15 5.623605591899965859e-15 3.465867171663132686e-18 3.897465960433185356e-14 4.136114348206219717e-16 -5.577700348424316488e-17 -2.850597210977476503e-19 1.578558198322722499e-16 1.574596966400089160e-16 1.116888499524688689e-16 3.637569907625921999e-15 -1.478423120060089212e-16 3.269968615819430543e-14 3.510560327678775616e-14 -1.209817091222388069e-14 -4.090315881827249918e-14 3.317222791750073899e-15 7.690375928819282975e-16 +-3.098661009272525452e-15 2.598336304972241831e-13 -3.705586445323490784e-17 -1.643837790116418740e-14 -7.731203689154381191e-14 -4.265439752377840368e-17 1.802000319076651299e-13 -1.010094527683821402e-13 -4.845003724100701487e-15 -1.943758295631656079e-15 -5.255400233547673271e-15 -5.253591335895288997e-13 7.026483273915552917e-26 9.563873128644333457e-27 7.303558544180336786e-26 2.936862964043296133e-23 -2.642810303210343532e-18 -2.643764778721222251e-18 -3.963086126803719396e-16 2.376937110121502144e-14 3.064589526524977282e-14 2.284806361195078194e-17 2.579873988473421512e-13 3.168273772135320469e-15 -3.091775265829969024e-16 -1.241280744886867784e-17 1.054145316287755313e-15 1.049161398503431286e-15 7.417407914940684902e-16 1.908816841781816436e-14 -9.268810345429511462e-16 2.059572385345710425e-13 2.235959995215321721e-13 -7.622820185235721193e-14 -2.605932125881225939e-13 2.553960731936533514e-14 5.242908346781543878e-15 +2.582109938010748615e-08 -2.151349431467241060e-06 -3.382128028653160845e-10 -2.415063062091929880e-08 -3.333970888846127249e-10 2.258055017134934322e-09 -7.553474453448608289e-06 2.222821203875844844e-06 3.813363256363073580e-08 8.266752272149200309e-09 2.856267415224482285e-08 7.794795518439994591e-08 -5.379740661803572684e-19 1.229211106917345516e-19 -3.863256970748645692e-19 -2.642810303210343532e-18 7.403427778043011779e-11 3.792255765611489853e-11 -1.752501317773718225e-08 -3.892781849899091207e-07 -2.772317334686791469e-07 -5.393331797473890970e-11 -2.529673605135954261e-06 -4.363570597273100514e-08 -1.095915722966923104e-07 9.977143842030775020e-10 -1.309301747407404281e-08 -1.353873145121602425e-08 -7.794890802066248781e-09 -1.475428397776662541e-07 6.079910763427108562e-09 -1.429439820800710253e-06 -1.674430426264918292e-06 5.296222786552825868e-07 1.952425956793901815e-06 -7.209465997261584918e-07 -4.236010055304580633e-08 +2.624612085038632830e-08 -2.167558592309546433e-06 -3.237635627080133576e-10 -3.005924930827426758e-08 -8.709526033092263374e-08 1.973339258111572039e-09 -6.879771852291717019e-06 2.537007182150678717e-06 3.857585521230416037e-08 9.778994440663878912e-09 2.904112999090483503e-08 7.846052278910172042e-08 -5.402147542475392543e-19 9.086668192364899074e-20 -3.853063443221715582e-19 -2.643764778721222251e-18 3.792255765611489853e-11 7.115769517967051503e-11 -1.816743239538513489e-08 -4.007424093499594748e-07 -2.822264494255077849e-07 -5.807776819282940936e-11 -2.854449519605377725e-06 -7.525430127036434062e-08 -1.081976407863369931e-07 1.566013896545607172e-09 -1.251205534443805253e-08 -1.287411686617860229e-08 -7.601134826475198951e-09 -1.479189936788522961e-07 6.683421504262534207e-09 -1.455166275855798905e-06 -1.690749426185897373e-06 5.391347495603226767e-07 1.971511674726411526e-06 -8.267564687813066038e-07 -4.700992181711387086e-08 +3.608749591974849653e-06 -2.808918726555411629e-04 2.880580686211860440e-07 2.969477106024997500e-07 5.768850215532460621e-05 -2.312332839147157345e-07 6.747976667041734062e-04 -6.392355604498717206e-05 3.200685394482180146e-06 3.524928150674835131e-06 4.142496374942100873e-06 1.148661919795845201e-05 -1.533973283924579781e-17 -4.005920918774881884e-17 -6.025552387684062132e-17 -3.963086126803719396e-16 -1.752501317773718225e-08 -1.816743239538513489e-08 3.264114120969943217e-05 -2.548394869425140520e-05 -2.779137482701451935e-04 -7.915047909730105148e-08 1.886790910630985283e-03 1.090754866631791889e-05 7.801050689041739985e-06 -3.432197661963476500e-06 -5.301808378832312896e-07 -4.152882865697312327e-07 -6.031727788188544236e-07 -2.128547701336814147e-04 1.234274341509008116e-06 -4.493057470519621835e-05 -7.206594708028664050e-05 1.697945549754080645e-05 8.504953299522758909e-05 1.777752907677872599e-04 2.585254330549545762e-05 +-6.937855671113378833e-05 1.185104206621878306e-02 -5.630122258543856047e-05 -2.599807913875601048e-03 -1.808405967871253364e-02 -2.916663369361920704e-05 1.029989985721349172e-01 -3.473801794131257253e-02 -2.112929067660789938e-04 9.244866045590156400e-05 -2.570082280690285368e-04 -6.454285445130538692e-04 5.186818943221322513e-15 -2.847060601679682921e-15 4.625905237017853671e-15 2.376937110121502144e-14 -3.892781849899091207e-07 -4.007424093499594748e-07 -2.548394869425140520e-05 5.406641666172327942e-02 5.174314523077627043e-02 -2.119053630567342706e-08 4.233491535612741996e-02 -1.462953360525826117e-03 3.121355558429767408e-03 -4.192856834587190711e-04 1.080013838424443047e-04 1.186614724903431369e-04 5.628102500579559636e-05 3.656786376748527578e-02 -2.546513518647200729e-06 1.417633226812336188e-02 1.486108430352409832e-02 -5.230857934383274961e-03 -1.726575396760205230e-02 2.391948208788366448e-03 -1.790482763215092395e-03 +-1.371705693477332281e-04 1.657883033879654217e-02 -5.550398399081168347e-05 -2.502578179730502710e-03 -1.896972119286157929e-02 -2.519716699320902141e-05 9.065041706401480437e-02 -3.257967242021848203e-02 -2.858688383419348751e-04 2.724159914307262731e-05 -3.251801223006570694e-04 -8.440052609354726587e-04 5.756938989258702090e-15 -2.127964203851167272e-15 5.623605591899965859e-15 3.064589526524977282e-14 -2.772317334686791469e-07 -2.822264494255077849e-07 -2.779137482701451935e-04 5.174314523077627043e-02 5.185024329073833343e-02 1.255540786621290009e-06 1.302520162131516959e-02 -1.618341401727035021e-03 2.953214562481817212e-03 -3.416134893996327412e-04 1.167556037089365365e-04 1.252744296441264085e-04 6.603921468351778485e-05 3.689859664766347025e-02 -2.563449933879626929e-05 1.590642365676888520e-02 1.704477207252333079e-02 -5.876392370890174723e-03 -1.982649448940359757e-02 -3.922392057424447923e-05 -1.976172892504302666e-03 +-2.071765376120491374e-07 1.543514452236958691e-05 -1.368692423501269688e-08 -3.448229919827389252e-08 -3.448716959724578149e-06 6.390334241489579055e-09 -2.040043730075414874e-05 3.079233917323632859e-06 -2.408182879736639907e-07 -1.754112212171557053e-07 -2.402523089244516050e-07 -6.640530256914299064e-07 2.461311477771505841e-18 1.628865061860011811e-18 3.465867171663132686e-18 2.284806361195078194e-17 -5.393331797473890970e-11 -5.807776819282940936e-11 -7.915047909730105148e-08 -2.119053630567342706e-08 1.255540786621290009e-06 2.923778826857790215e-09 8.101263196662729328e-06 -1.850550483830866121e-07 -4.340263833203487453e-07 1.100083824610656931e-08 4.306354748750698602e-08 3.889510608936747429e-08 3.734877925264959754e-08 8.633797122601068312e-07 -6.417850997642396732e-08 7.220878837043200342e-06 8.583828081722996397e-06 -2.684398257591137100e-06 -1.003887185171714319e-05 -8.037912468433751875e-07 3.038103882007981565e-07 +-2.516517250145258906e-03 2.136226422370665345e-01 -5.511266268928062019e-04 9.894412705418731673e-04 3.799553542857860011e-02 -2.181936752914516461e-05 4.725134464369972614e-02 -3.340170502487335713e-02 -2.713414985469090809e-03 -1.639613353212979532e-03 -3.216676545935040574e-03 -7.926725382679215490e-03 3.091641941796989494e-14 8.722981696312964399e-15 3.897465960433185356e-14 2.579873988473421512e-13 -2.529673605135954261e-06 -2.854449519605377725e-06 1.886790910630985283e-03 4.233491535612741996e-02 1.302520162131516959e-02 8.101263196662729328e-06 2.561886360437845678e+00 1.712688632339961700e-02 -1.495969548298849670e-02 -4.941174031154322611e-03 8.451753728352438908e-04 8.312588628538810937e-04 6.207798751111286498e-04 -1.486990425054944553e-02 -6.548694347379695625e-04 8.400416410730938976e-02 1.040558310541357112e-01 -3.123910806711079069e-02 -1.216717023960072891e-01 1.668395658113263902e-01 2.568679781949267243e-02 +-3.928954203215390470e-05 2.137533177873366264e-03 5.445097707514176965e-06 1.614015226545799700e-04 1.185191496548233239e-03 -2.885741513929382642e-06 1.027426075386230972e-02 -5.494265105720975290e-03 -2.657470988356218109e-05 -1.838827713483549386e-05 -3.961294658191956365e-05 -1.011002018751854773e-04 -1.759372847258843210e-16 -1.245990632231858332e-16 4.136114348206219717e-16 3.168273772135320469e-15 -4.363570597273100514e-08 -7.525430127036434062e-08 1.090754866631791889e-05 -1.462953360525826117e-03 -1.618341401727035021e-03 -1.850550483830866121e-07 1.712688632339961700e-02 1.114660806185970585e-03 -7.730335285552777892e-04 -4.452308713267074861e-06 1.481836476482716150e-05 1.451846216162652682e-05 8.452527659270013381e-06 -1.359203485518305915e-03 -1.601010014264435757e-06 -1.002409728823324721e-03 -3.485597911598559660e-04 3.639571529618867399e-04 3.878510497305322842e-04 1.604008248187792666e-03 1.834410652820026033e-04 +1.279490468007478382e-06 2.287660214046896914e-03 1.881736127791208387e-06 1.096321509004083357e-05 1.182857261486174993e-03 -9.584663865549323316e-06 3.051814272148788490e-02 -5.915567986028016229e-03 -2.803719763090280236e-05 4.881136721097782943e-05 2.107728167369185196e-06 9.225664148879019028e-06 1.070459744475160279e-15 -1.035329050434101653e-15 -5.577700348424316488e-17 -3.091775265829969024e-16 -1.095915722966923104e-07 -1.081976407863369931e-07 7.801050689041739985e-06 3.121355558429767408e-03 2.953214562481817212e-03 -4.340263833203487453e-07 -1.495969548298849670e-02 -7.730335285552777892e-04 1.688998948499471605e-03 1.159026578533110570e-05 2.989596973369903929e-05 3.373840088359871794e-05 1.465319842208428671e-05 1.862022016073073931e-03 -5.124694057679437498e-06 2.798040993920551233e-03 2.876093839470280755e-03 -1.029387046857675172e-03 -3.333290952544525499e-03 3.469864107285604033e-03 -1.855656174130421361e-04 +9.068426317865086325e-08 2.979745879017003073e-04 1.074217982496532016e-06 7.773663567167847920e-07 5.547822027831921226e-05 -5.294172336870067481e-09 1.825216372649425331e-05 3.383003929474542214e-05 2.580429643797891438e-07 1.113170450201948483e-06 6.525664614377412270e-07 1.051179776515564498e-06 9.401517695428871141e-18 -1.089157915915628154e-18 -2.850597210977476503e-19 -1.241280744886867784e-17 9.977143842030775020e-10 1.566013896545607172e-09 -3.432197661963476500e-06 -4.192856834587190711e-04 -3.416134893996327412e-04 1.100083824610656931e-08 -4.941174031154322611e-03 -4.452308713267074861e-06 1.159026578533110570e-05 2.103990647271761701e-05 1.114955092000448355e-06 1.206700317689306529e-06 9.328667613641622459e-07 -2.012489084703237706e-04 -1.281890104442271536e-06 3.024086810582569051e-05 2.773027802293759217e-05 -1.096665998980926882e-05 -3.201338373551938023e-05 -2.511310858789703626e-04 -4.049192332247160327e-05 +-1.010411452316174498e-05 1.676699508613893001e-03 -9.573321138379633936e-07 6.743863193988200014e-06 1.003402092902857769e-04 -5.146313275530546443e-07 1.755644760324465278e-03 -5.605255102331018402e-04 -1.010510257238622547e-05 -2.085181637611447406e-06 -9.965538983237963089e-06 -2.918386730194719129e-05 1.520685182370616798e-16 -1.201790912326496301e-17 1.578558198322722499e-16 1.054145316287755313e-15 -1.309301747407404281e-08 -1.251205534443805253e-08 -5.301808378832312896e-07 1.080013838424443047e-04 1.167556037089365365e-04 4.306354748750698602e-08 8.451753728352438908e-04 1.481836476482716150e-05 2.989596973369903929e-05 1.114955092000448355e-06 2.884322043042034708e-04 -2.538126012564094220e-04 5.387122904190521982e-06 6.894907238685718954e-05 -2.703157629667282969e-06 2.554210195701856097e-03 4.889853800134472996e-04 -9.511869428156902347e-04 -5.708735812038205929e-04 1.938946049089669489e-04 1.649637945841375946e-05 +-1.002915587787911641e-05 1.732916157635252218e-03 -8.955180902632183681e-07 6.281162398372516905e-06 1.077696970302018857e-04 -5.539881244564768786e-07 1.870870624454052007e-03 -5.695318433638057315e-04 -9.686357516880338036e-06 -1.662825435048065551e-06 -9.808018743190203645e-06 -2.889439043922104893e-05 1.525755060186724085e-16 -1.672240764147307434e-17 1.574596966400089160e-16 1.049161398503431286e-15 -1.353873145121602425e-08 -1.287411686617860229e-08 -4.152882865697312327e-07 1.186614724903431369e-04 1.252744296441264085e-04 3.889510608936747429e-08 8.312588628538810937e-04 1.451846216162652682e-05 3.373840088359871794e-05 1.206700317689306529e-06 -2.538126012564094220e-04 3.102166925762899776e-04 5.611063091380308048e-06 7.541673326503943767e-05 -2.661931171835422640e-06 1.470998787372507915e-02 4.955205309483954376e-04 -5.437021090677873411e-03 -5.783272781411675492e-04 1.993184776291691035e-04 1.588876821035091380e-05 +-7.096411239012776800e-06 1.298896908940160239e-03 -8.519837545226716281e-07 4.516749168953107673e-06 8.883519413119686178e-05 -2.523146163460764631e-07 8.708104448430653436e-04 -2.820483943677134196e-04 -5.814282625649671241e-06 -1.529648396480678433e-06 -6.833063562445468623e-06 -2.028183930935848183e-05 8.834165471558264173e-17 2.748078604662010185e-18 1.116888499524688689e-16 7.417407914940684902e-16 -7.794890802066248781e-09 -7.601134826475198951e-09 -6.031727788188544236e-07 5.628102500579559636e-05 6.603921468351778485e-05 3.734877925264959754e-08 6.207798751111286498e-04 8.452527659270013381e-06 1.465319842208428671e-05 9.328667613641622459e-07 5.387122904190521982e-06 5.611063091380308048e-06 1.988825337167201494e-05 3.812443652156480859e-05 -1.944066131252356749e-06 2.255075602087547727e-04 9.370667484578088652e-03 -8.357229684017524840e-05 -1.091244912444888637e-02 1.067404421836058503e-04 1.284101720431881188e-05 +-6.166408270828315896e-05 9.282782630378589445e-03 -3.229239994559254898e-05 -1.923508018533116436e-03 -1.588843399738775614e-02 -1.594898464460605111e-05 5.720248449212836372e-02 -1.955720914959569248e-02 -1.795567582710501570e-04 3.754504312323681382e-05 -1.870407675260386544e-04 -5.055290379226907082e-04 4.235084795229275426e-15 -1.462373406033409279e-15 3.637569907625921999e-15 1.908816841781816436e-14 -1.475428397776662541e-07 -1.479189936788522961e-07 -2.128547701336814147e-04 3.656786376748527578e-02 3.689859664766347025e-02 8.633797122601068312e-07 -1.486990425054944553e-02 -1.359203485518305915e-03 1.862022016073073931e-03 -2.012489084703237706e-04 6.894907238685718954e-05 7.541673326503943767e-05 3.812443652156480859e-05 2.724309422537568928e-02 -7.792580236021764118e-06 1.214317040165824486e-02 1.235390285211402685e-02 -4.479740771749052969e-03 -1.435268803652801212e-02 -6.216048243951985408e-03 -1.671552675289353493e-03 +7.954518997506867072e-06 -6.763338260323719427e-04 3.022377409523564049e-07 1.331924023050925137e-05 3.661169141778293729e-05 8.444660670564603750e-08 -3.219229529387673773e-04 2.528818449588239310e-04 1.072466086869132534e-05 5.070866914767413610e-06 1.107990032257984879e-05 2.773567462221827859e-05 -1.343201925832196203e-16 -2.706459684488763294e-17 -1.478423120060089212e-16 -9.268810345429511462e-16 6.079910763427108562e-09 6.683421504262534207e-09 1.234274341509008116e-06 -2.546513518647200729e-06 -2.563449933879626929e-05 -6.417850997642396732e-08 -6.548694347379695625e-04 -1.601010014264435757e-06 -5.124694057679437498e-06 -1.281890104442271536e-06 -2.703157629667282969e-06 -2.661931171835422640e-06 -1.944066131252356749e-06 -7.792580236021764118e-06 3.938618265045845181e-06 -3.810800813135566405e-04 -4.397826624549711397e-04 1.413248991076312738e-04 5.133429895970311298e-04 -8.965737950838059862e-05 -2.129537338315999433e-05 +-1.457372693470361177e-03 9.451541906012153604e-02 1.005144316480618646e-04 -5.849646589347031259e-03 -1.294007493555653365e-01 -2.640464914727091680e-05 1.057903669450520889e-01 -4.322324541853273666e-02 -4.802813362752815160e-03 -9.756958387685016505e-04 -1.679299878847527938e-03 -5.498729953080172886e-03 8.757869065918129794e-14 3.979391298362956564e-15 3.269968615819430543e-14 2.059572385345710425e-13 -1.429439820800710253e-06 -1.455166275855798905e-06 -4.493057470519621835e-05 1.417633226812336188e-02 1.590642365676888520e-02 7.220878837043200342e-06 8.400416410730938976e-02 -1.002409728823324721e-03 2.798040993920551233e-03 3.024086810582569051e-05 2.554210195701856097e-03 1.470998787372507915e-02 2.255075602087547727e-04 1.214317040165824486e-02 -3.810800813135566405e-04 7.124318197334406300e+00 2.676044558680029106e-01 -2.632743994364230122e+00 -3.111554357622127176e-01 1.764057209661643380e-03 1.756704866057536658e-03 +-1.665472713324716442e-03 1.143881756135775030e-01 6.291873262212429902e-05 -5.178846973454946830e-03 -1.216020000364273906e-01 -3.561425506840147387e-05 1.403403673864646317e-01 -6.023684334790740946e-02 -4.963194473939220967e-03 -1.072946360077371318e-03 -1.899267739233110111e-03 -6.050663216414568953e-03 8.734305031308171818e-14 3.914280900809954267e-15 3.510560327678775616e-14 2.235959995215321721e-13 -1.674430426264918292e-06 -1.690749426185897373e-06 -7.206594708028664050e-05 1.486108430352409832e-02 1.704477207252333079e-02 8.583828081722996397e-06 1.040558310541357112e-01 -3.485597911598559660e-04 2.876093839470280755e-03 2.773027802293759217e-05 4.889853800134472996e-04 4.955205309483954376e-04 9.370667484578088652e-03 1.235390285211402685e-02 -4.397826624549711397e-04 2.676044558680029106e-01 7.686030370760271957e+00 -9.884182661136756143e-02 -8.940093905890782011e+00 7.354504365195378701e-03 2.178830498571339251e-03 +5.402267720757495664e-04 -3.500456106541598861e-02 -3.646627150018448262e-05 2.153268188959647257e-03 4.772874293545138064e-02 9.788159208161594391e-06 -3.924802079810447292e-02 1.611303878009478704e-02 1.775867457136094804e-03 3.617995149859560026e-04 6.223804919963958409e-04 2.036131827085406725e-03 -3.234243315335766641e-14 -1.474664675128080572e-15 -1.209817091222388069e-14 -7.622820185235721193e-14 5.296222786552825868e-07 5.391347495603226767e-07 1.697945549754080645e-05 -5.230857934383274961e-03 -5.876392370890174723e-03 -2.684398257591137100e-06 -3.123910806711079069e-02 3.639571529618867399e-04 -1.029387046857675172e-03 -1.096665998980926882e-05 -9.511869428156902347e-04 -5.437021090677873411e-03 -8.357229684017524840e-05 -4.479740771749052969e-03 1.413248991076312738e-04 -2.632743994364230122e+00 -9.884182661136756143e-02 9.729171339957245479e-01 1.149288680405661928e-01 -6.978562179520457964e-04 -6.534298927597297166e-04 +1.943321742220476944e-03 -1.335569570692768659e-01 -7.162025230758143908e-05 6.003639016901835684e-03 1.412309683815566597e-01 4.150511446075701752e-05 -1.636564359123428292e-01 7.044389440078994979e-02 5.777840977318744395e-03 1.251930415222073795e-03 2.215661331552164909e-03 7.054155462547743484e-03 -1.015755225261617099e-13 -4.571334409002748931e-15 -4.090315881827249918e-14 -2.605932125881225939e-13 1.952425956793901815e-06 1.971511674726411526e-06 8.504953299522758909e-05 -1.726575396760205230e-02 -1.982649448940359757e-02 -1.003887185171714319e-05 -1.216717023960072891e-01 3.878510497305322842e-04 -3.333290952544525499e-03 -3.201338373551938023e-05 -5.708735812038205929e-04 -5.783272781411675492e-04 -1.091244912444888637e-02 -1.435268803652801212e-02 5.133429895970311298e-04 -3.111554357622127176e-01 -8.940093905890782011e+00 1.149288680405661928e-01 1.039881525384822325e+01 -8.678148783162666774e-03 -2.548835131946574056e-03 +-3.333863796457967564e-04 2.491333162046089483e-02 -2.353051652724181582e-05 1.387439896431741414e-03 1.866643744552765397e-02 -4.390083089488345825e-05 1.451439494535149144e-01 -4.483966495114385248e-02 -3.405953240386222296e-04 -5.528075146326992567e-05 -3.730348743630668195e-04 -8.578412202787384210e-04 2.162259177191088287e-15 -3.184500434619874255e-15 3.317222791750073899e-15 2.553960731936533514e-14 -7.209465997261584918e-07 -8.267564687813066038e-07 1.777752907677872599e-04 2.391948208788366448e-03 -3.922392057424447923e-05 -8.037912468433751875e-07 1.668395658113263902e-01 1.604008248187792666e-03 3.469864107285604033e-03 -2.511310858789703626e-04 1.938946049089669489e-04 1.993184776291691035e-04 1.067404421836058503e-04 -6.216048243951985408e-03 -8.965737950838059862e-05 1.764057209661643380e-03 7.354504365195378701e-03 -6.978562179520457964e-04 -8.678148783162666774e-03 4.943448466562492399e-02 1.944024799375950992e-03 +-5.468710926687100291e-05 4.670416326687474098e-03 -4.818973351892917084e-06 8.672656714053429282e-05 1.065831122547543450e-03 4.229752159922561635e-07 -1.738338909140980243e-03 -8.233919291224370057e-05 -6.056307379456799997e-05 -3.912271590726792503e-05 -6.364872634283788350e-05 -1.606945507961290737e-04 6.313068106737503019e-16 2.915127137224032614e-16 7.690375928819282975e-16 5.242908346781543878e-15 -4.236010055304580633e-08 -4.700992181711387086e-08 2.585254330549545762e-05 -1.790482763215092395e-03 -1.976172892504302666e-03 3.038103882007981565e-07 2.568679781949267243e-02 1.834410652820026033e-04 -1.855656174130421361e-04 -4.049192332247160327e-05 1.649637945841375946e-05 1.588876821035091380e-05 1.284101720431881188e-05 -1.671552675289353493e-03 -2.129537338315999433e-05 1.756704866057536658e-03 2.178830498571339251e-03 -6.534298927597297166e-04 -2.548835131946574056e-03 1.944024799375950992e-03 4.524222784337420585e-04 diff --git a/contrib/anisotropic_stishovite/stishovite_data.py b/contrib/anisotropic_stishovite/stishovite_data.py index c0303742..f48ba4a7 100644 --- a/contrib/anisotropic_stishovite/stishovite_data.py +++ b/contrib/anisotropic_stishovite/stishovite_data.py @@ -150,6 +150,42 @@ def get_data(): d["P"] = d["T"] * 0.0 + 0.0001 d["P_err"] = d["T"] * 0.0 + 0.000001 + # Wang data (only stishovite) + W_data = np.loadtxt("data/Wang_et_al_2012_stv.dat") + header = [ + "id", + "T", + "P", + "P_err", + "a", + "a_err", + "c", + "c_err", + "V", + "V_err", + "VAu", + "VAu_err", + ] + compilation = [["stv", header, W_data]] + add_to_data(data, "Wang_2012", compilation) + + d = data["Wang_2012"]["stv"] + d["T_err"] = d["T"] / 100.0 + d["b"] = d["a"] + d["b_err"] = d["a_err"] + + # Nishihara data (only stishovite) + N_data = np.loadtxt("data/Nishihara_et_al_2005_stv.dat") + header = ["P", "T", "a", "a_err", "c", "c_err", "V", "V_err"] + compilation = [["stv", header, N_data]] + add_to_data(data, "Nishihara_2005", compilation) + + d = data["Nishihara_2005"]["stv"] + d["T_err"] = d["T"] / 100.0 + d["P_err"] = d["P"] * 0.05 + 0.2 # conservative estimate of the pressure error + d["b"] = d["a"] + d["b_err"] = d["a_err"] + # Convert to SI units d = data["Zhang_2021"]["CT"] d["P"] = d["P"] * 1.0e9 @@ -165,6 +201,8 @@ def get_data(): data["Fischer_2018"]["poststv"], data["Zhang_2021"]["stv"], data["Zhang_2021"]["poststv"], + data["Wang_2012"]["stv"], + data["Nishihara_2005"]["stv"], ]: d["V"] = molar_volume_from_unit_cell_volume(d["V"], 2.0) d["V_err"] = molar_volume_from_unit_cell_volume(d["V_err"], 2.0) @@ -214,7 +252,14 @@ def common_data(): d["abc"] = {"stv": {}, "poststv": {}} d["abc_err"] = {"stv": {}, "poststv": {}} - for pub in ["Ito_1974", "Andrault_2003", "Fischer_2018", "Zhang_2021"]: + for pub in [ + "Ito_1974", + "Andrault_2003", + "Nishihara_2005", + "Wang_2012", + "Fischer_2018", + "Zhang_2021", + ]: for phase in ["stv", "poststv"]: try: Z_data = data[pub][phase] diff --git a/contrib/anisotropic_stishovite/stishovite_fit_eos.py b/contrib/anisotropic_stishovite/stishovite_fit_eos.py index 137f7a41..0fcc55c8 100644 --- a/contrib/anisotropic_stishovite/stishovite_fit_eos.py +++ b/contrib/anisotropic_stishovite/stishovite_fit_eos.py @@ -2,7 +2,7 @@ from __future__ import print_function import numpy as np -from scipy.optimize import minimize, differential_evolution +from scipy.optimize import minimize, differential_evolution, leastsq from stishovite_data import ( common_data, @@ -64,8 +64,12 @@ def misfit_scalar(args): P_tr_GPa, fP_Zhang, fP_Andrault, + fP_Wang, + fP_Nishihara, fP2_Zhang, fP2_Andrault, + fP2_Wang, + fP2_Nishihara, ) = args scalar_prms = make_scalar_model( @@ -88,10 +92,12 @@ def misfit_scalar(args): P_tr_3000K_model = transition_pressure(scalar_stv, 3000.0) / 1.0e9 P_tr_3000K_obs = 90.0 P_tr_3000K_err = 1.0 - chisqr = np.sum(np.power((P_tr_3000K_model - P_tr_3000K_obs) / P_tr_3000K_err, 2.0)) + misfits = [(P_tr_3000K_model - P_tr_3000K_obs) / P_tr_3000K_err] # Fit volumes for stishovite from Zhang, Andrault for phase, publication in [ + ("stv", "Nishihara_2005"), + ("stv", "Wang_2012"), ("stv", "Zhang_2021"), ("poststv", "Zhang_2021"), ("stv", "Andrault_2003"), @@ -103,6 +109,12 @@ def misfit_scalar(args): elif publication == "Andrault_2003": fP = fP_Andrault fP2 = fP2_Andrault + elif publication == "Wang_2012": + fP = fP_Wang + fP2 = fP2_Wang + elif publication == "Nishihara_2005": + fP = fP_Nishihara + fP2 = fP2_Nishihara else: fP = 1.0 fP2 = 0.0 @@ -114,7 +126,7 @@ def misfit_scalar(args): 0 ] P_actual = PTV[:, 0] * (fP + fP2 * PTV[:, 0]) - chisqr += np.sum(np.power((P_model - P_actual) / PTV_err[:, 0], 2.0)) + misfits.extend(list((P_model - P_actual) / PTV_err[:, 0])) if plot: sort = np.argsort(P_model) @@ -128,7 +140,7 @@ def misfit_scalar(args): scalar_stv.set_state(1.0e5, 298.15) PTV[:, 2] = PTV[:, 2] / PTV[0, 2] * scalar_stv.V V_model = relaxed_stv.evaluate(["V"], PTV[:, 0], PTV[:, 1])[0] - chisqr += np.sum(np.power((V_model - PTV[:, 2]) / PTV_err[:, 2], 2.0)) + misfits.extend(list((V_model - PTV[:, 2]) / PTV_err[:, 2])) if plot: sort = np.argsort(PTV[:, 1]) @@ -139,7 +151,7 @@ def misfit_scalar(args): KNR_obs = relaxed_stv.evaluate_with_volumes( ["isentropic_bulk_modulus_reuss"], V_for_CN, T_for_CN )[0] - chisqr += np.sum(np.power((KNR_GPa - KNR_obs / 1.0e9) / KNR_err_GPa, 2.0)) + misfits.extend(list((KNR_GPa - KNR_obs / 1.0e9) / KNR_err_GPa)) if plot: sort = np.argsort(V_for_CN) @@ -147,11 +159,14 @@ def misfit_scalar(args): ax[2].scatter(V_for_CN[sort], KNR_GPa[sort]) plt.show() + misfits = np.array(misfits) + chisqr = np.sum(np.power(misfits, 2.0)) + if chisqr < min_misfit_scalar[0]: min_misfit_scalar[0] = chisqr print(repr(args)) print(chisqr) - return chisqr + return misfits def misfit_cell(args, scalar_args): @@ -166,8 +181,12 @@ def misfit_cell(args, scalar_args): P_tr_GPa, fP_Zhang, fP_Andrault, + fP_Wang, + fP_Nishihara, fP2_Zhang, fP2_Andrault, + fP2_Wang, + fP2_Nishihara, ) = scalar_args (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = args ( @@ -229,7 +248,7 @@ def misfit_cell(args, scalar_args): _, _, _, stishovite_relaxed = models stishovite_relaxed.set_composition([1.0]) - chisqr = 0.0 + misfits = [] if plot: fig = plt.figure() @@ -248,6 +267,12 @@ def misfit_cell(args, scalar_args): elif publication == "Andrault_2003": fP = fP_Andrault fP2 = fP2_Andrault + elif publication == "Wang_2012": + fP = fP_Wang + fP2 = fP2_Wang + elif publication == "Nishihara_2005": + fP = fP_Nishihara + fP2 = fP2_Nishihara else: fP = 1.0 fP2 = 0.0 @@ -268,7 +293,7 @@ def misfit_cell(args, scalar_args): )[0] abc_model = cell_parameters_model[:, :3] - chisqr += np.sum(np.power((abc_obs - abc_model) / abc_err_obs, 2.0)) + misfits.extend(list(np.ravel((abc_obs - abc_model) / abc_err_obs))) if plot: for i in range(3): j = 1 @@ -296,11 +321,14 @@ def misfit_cell(args, scalar_args): ax[j].plot(pressures, abc_model[:, i]) plt.show() + misfits = np.array(misfits) + chisqr = np.sum(np.power(misfits, 2.0)) + if chisqr < min_misfit_cell[0]: min_misfit_cell[0] = chisqr print(repr(args)) print(chisqr) - return chisqr + return misfits def misfit_elastic(args, scalar_args, cell_args): @@ -315,8 +343,12 @@ def misfit_elastic(args, scalar_args, cell_args): P_tr_GPa, fP_Zhang, fP_Andrault, + fP_Wang, + fP_Nishihara, fP2_Zhang, fP2_Andrault, + fP2_Wang, + fP2_Nishihara, ) = scalar_args (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = cell_args (a11, a22, a33, a44, a55, a66, b11i, b33i, b44i, b66i, c44, c66, b112i, b332i) = ( @@ -377,7 +409,7 @@ def misfit_elastic(args, scalar_args, cell_args): _, _, _, stishovite_relaxed = models stishovite_relaxed.set_composition([1.0]) - chisqr = 0.0 + misfits = [] if plot: fig = plt.figure() @@ -403,7 +435,7 @@ def misfit_elastic(args, scalar_args, cell_args): (5, 5), ]: chis = (CN_GPa[:, i, j] - CN_model[:, i, j] / 1.0e9) / CN_err_GPa[:, i, j] - chisqr += np.sum(np.power(chis, 2.0)) + misfits.extend(list(chis)) if plot: axi = 2 @@ -428,63 +460,74 @@ def misfit_elastic(args, scalar_args, cell_args): ax[axi].legend() plt.show() + misfits = np.array(misfits) + chisqr = np.sum(np.power(misfits, 2.0)) + if chisqr < min_misfit_elastic[0]: min_misfit_elastic[0] = chisqr print(repr(args)) print(chisqr) - return chisqr + return misfits def misfit_scalar_and_cell(args): - scalar_args = args[:12] - cell_args = args[12:] - chisqr = misfit_scalar(scalar_args) - chisqr += misfit_cell(cell_args, scalar_args) + scalar_args = args[:16] + cell_args = args[16:] + misfits = np.concatenate( + (misfit_scalar(scalar_args), misfit_cell(cell_args, scalar_args)) + ) + chisqr = np.sum(np.power(misfits, 2.0)) if chisqr < min_misfit_combined[0]: min_misfit_combined[0] = chisqr print(repr(args)) print(chisqr) - return chisqr + return misfits def misfit_cell_and_elastic(args, scalar_args): cell_args = args[:7] elastic_args = args[7:] - chisqr = misfit_cell(cell_args, scalar_args) + misfits_1 = misfit_cell(cell_args, scalar_args) modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) - chisqr += misfit_elastic(elastic_args, scalar_args, cell_args) + misfits_2 = misfit_elastic(elastic_args, scalar_args, cell_args) + + misfits = np.concatenate((misfits_1, misfits_2)) + chisqr = np.sum(np.power(misfits, 2.0)) if chisqr < min_misfit_combined[0]: min_misfit_combined[0] = chisqr print(repr(args)) print(chisqr) - return chisqr + return misfits def misfit_all(args): - scalar_args = args[:12] - cell_args = args[12:19] - elastic_args = args[19:] - chisqr = misfit_scalar(scalar_args) - chisqr += misfit_cell(cell_args, scalar_args) + scalar_args = args[:16] + cell_args = args[16:23] + elastic_args = args[23:] + misfits_1 = misfit_scalar(scalar_args) + misfits_2 = misfit_cell(cell_args, scalar_args) modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) - chisqr += misfit_elastic(elastic_args, scalar_args, cell_args) + misfits_3 = misfit_elastic(elastic_args, scalar_args, cell_args) + + misfits = np.concatenate((misfits_1, misfits_2, misfits_3)) + chisqr = np.sum(np.power(misfits, 2.0)) if chisqr < min_misfit_combined[0]: min_misfit_combined[0] = chisqr - print(repr(args[:12])) - print(repr(args[12:19])) - print(repr(args[19:])) + print(repr(args[:16])) + print(repr(args[16:23])) + print(repr(args[23:])) print(chisqr) - return chisqr + return misfits if __name__ == "__main__": if False: - sol = minimize( - misfit_scalar, scalar_args, bounds=scalar_bounds, method="Nelder-Mead" - ) + print("Starting misfit_scalar") + sol = leastsq(misfit_scalar, scalar_args, full_output=True) + print(sol) if False: @@ -502,14 +545,37 @@ def misfit_de(cell_args): ) if False: - sol = minimize( - misfit_scalar_and_cell, - scalar_and_cell_args, - bounds=scalar_bounds + cell_bounds, - method="Nelder-Mead", - ) + print("Starting misfit_scalar_and_cell") + sol = leastsq( + misfit_scalar_and_cell, scalar_and_cell_args, factor=0.01, full_output=True + ) # factor = 0.01 definitely helps + print(sol) if True: + print("Starting misfit_all") + sol = leastsq( + misfit_all, all_args, factor=0.01, full_output=True + ) # factor = 0.01 definitely helps + + # Estimating covariance. See here: + # https://stackoverflow.com/q/14854339/6272561 + x = sol[0] + residuals = misfit_all(x) + n_data = len(residuals) + n_params = len(x) + residual_variance = np.sum(np.power(residuals, 2.0)) / (n_data - n_params) + + sig = np.sqrt(np.diag(sol[1] * residual_variance)) + np.savetxt("model_output/covariance_matrix.dat", sol[1] * residual_variance) + print(repr(sig)) + for i in range(len(x)): + print(f"{x[i]:.4e} +/- {sig[i]:.4e}") + + if False: + + def lsq(args): + return np.sum(np.power(misfit_all(args), 2.0)) + sol = minimize( misfit_all, all_args, diff --git a/contrib/anisotropic_stishovite/stishovite_model.py b/contrib/anisotropic_stishovite/stishovite_model.py index 0a43c094..9a289208 100644 --- a/contrib/anisotropic_stishovite/stishovite_model.py +++ b/contrib/anisotropic_stishovite/stishovite_model.py @@ -24,8 +24,12 @@ def modify_Zhang_elasticity(scalar_args, cell_args, elastic_args): P_tr_GPa, fP_Zhang, fP_Andrault, + fP_Wang, + fP_Nishihara, fP2_Zhang, fP2_Andrault, + fP2_Wang, + fP2_Nishihara, ) = scalar_args (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = cell_args (a11, a22, a33, a44, a55, a66, b11i, b33i, b44i, b66i, c44, c66, b112i, b332i) = ( @@ -498,9 +502,9 @@ def fn(lnV, Pth, X, params): def get_models(): - scalar_args = all_args[:12] - cell_args = all_args[12:19] - elastic_args = all_args[19:] + scalar_args = all_args[:16] + cell_args = all_args[16:23] + elastic_args = all_args[23:] modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) @@ -515,8 +519,12 @@ def get_models(): P_tr_GPa, fP_Zhang, fP_Andrault, + fP_Wang, + fP_Nishihara, fP2_Zhang, fP2_Andrault, + fP2_Wang, + fP2_Nishihara, ) = scalar_args (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = cell_args (a11, a22, a33, a44, a55, a66, b11i, b33i, b44i, b66i, c44, c66, b112i, b332i) = ( @@ -583,6 +591,8 @@ def get_models(): sm = models[1] Q0, Q1 = sm.solution_model.interaction_endmembers + Q1.params["F_0"] -= Q0.params["F_0"] + Q0.params["F_0"] = 0 properties = ["F_0", "V_0", "K_0", "Kprime_0", "Debye_0", "grueneisen_0", "q_0"] pps = [ "$\\mathcal{F}_0$", @@ -595,12 +605,13 @@ def get_models(): ] table = [["", "$Q = 0$", "$Q = 1$"]] for i, p in enumerate(properties): - table.append([pps[i], f"{Q0.params[p]:.6e}", f"{Q1.params[p]:.6e}"]) + table.append([pps[i], f"{Q0.params[p]:.4e}", f"{Q1.params[p]:.4e}"]) + table[1][1] = "-" table.append(["$a_0$", "-", a0Q1]) table.append(["$b_0$", "-", b0Q1]) - print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".6e")) + print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".4e")) print("b_xs", sm.b_xs) @@ -663,7 +674,7 @@ def get_models(): ) table = [ - [f"{item:.6e}" if type(item) is np.float64 else item for item in row] + [f"{item:.4e}" if type(item) is np.float64 else item for item in row] for row in table ] print( @@ -671,11 +682,23 @@ def get_models(): list(map(list, zip(*table))), headers="firstrow", tablefmt="latex_raw", - floatfmt=".6e", + floatfmt=".4e", ) ) print("Corrections to pressure:") - print(f" Zhang: {fP_Zhang}, {fP2_Zhang}") - print(f" Andrault: {fP_Andrault}, {fP2_Andrault}") - return models, fP_Zhang, fP_Andrault, fP2_Zhang, fP2_Andrault + print(f"${fP_Zhang:.4f}{fP2_Zhang*1e9:+.4f}P$ (Zhang)") + print(f"${fP_Andrault:.4f}{fP2_Andrault*1e9:+.4f}P$ (Andrault)") + print(f"${fP_Wang:.4f}{fP2_Wang*1e9:+.4f}P$ (Wang)") + print(f"${fP_Nishihara:.4f}{fP2_Nishihara*1e9:+.4f}P$ (Nishihara)") + return ( + models, + fP_Zhang, + fP_Andrault, + fP_Wang, + fP_Nishihara, + fP2_Zhang, + fP2_Andrault, + fP2_Wang, + fP2_Nishihara, + ) diff --git a/contrib/anisotropic_stishovite/stishovite_model_Carpenter_2000.py b/contrib/anisotropic_stishovite/stishovite_model_Carpenter_2000.py new file mode 100644 index 00000000..1a559142 --- /dev/null +++ b/contrib/anisotropic_stishovite/stishovite_model_Carpenter_2000.py @@ -0,0 +1,303 @@ +import numpy as np +from scipy.integrate import cumulative_trapezoid +import matplotlib.pyplot as plt +from burnman.minerals.HP_2011_ds62 import stv +from string import ascii_lowercase + +# Expressions for the bare elastic constant values +# are taken from Table 2 +# All units are GPa + + +def C11(P_GPa): + return 578 + 5.38 * P_GPa + + +def C33(P_GPa): + return 776 + 4.94 * P_GPa + + +def C12(P_GPa): + return 86 + 5.38 * P_GPa + + +def C13(P_GPa): + return 191 + 2.72 * P_GPa + + +def C44(P_GPa): + return 252 + 1.88 * P_GPa + + +def C66(P_GPa): + return 323 + 3.10 * P_GPa + + +# The following parameters are also taken from Table 2 +Pastc = 49 +Pc = 49.0 + 50.7 + +l1 = -8.0 +l2 = 24.62 +l3 = 17.0 +a = -0.04856 +b = 10.94 + +# l4 and l6 are taken from the text +l4 = 20.0 +l6 = 20.0 + + +def gibbs_over_tetragonal_volume_1(P_GPa): + """ + The Gibbs energy divided by the + volume of the (bare) tetragonal phase, + as defined in Equation 1. + """ + tet_mask = np.array([P < Pastc for P in P_GPa]) + gibbs_values = np.empty_like(P_GPa) + gibbs_values[tet_mask] = 0.0 + + P = P_GPa[~tet_mask] + Q2 = Qsqr(P) + Q = np.sqrt(Q2) + Q4 = Q2 * Q2 + + e1se2 = -l2 / (0.5 * (C11(P) - C12(P))) * Q + e1pe2 = -l1 / (0.5 * (C11(P) + C12(P))) * Q2 + e1 = (e1se2 + e1pe2) / 2.0 + e2 = (-e1se2 + e1pe2) / 2.0 + e3 = -l3 / C33(P) * Q2 + + gibbs_values[~tet_mask] = ( + 0.5 * a * (P - Pc) * Q2 + + 0.25 * b * Q4 + + l1 * (e1 + e2) * Q2 + + l2 * (e1 - e2) * Q + + l3 * e3 * Q2 + + 0.25 * (C11(P) + C12(P)) * e1pe2 * e1pe2 + + 0.25 * (C11(P) - C12(P)) * e1se2 * e1se2 + + C13(P) * e1pe2 * e3 + + 0.5 * C33(P) * e3 * e3 + ) + return gibbs_values + + +def gibbs_over_tetragonal_volume_4(P_GPa): + """ + The Gibbs energy divided by the + volume of the (bare) tetragonal phase, + as defined in Equation 4. + """ + tet_mask = np.array([P < Pastc for P in P_GPa]) + P = P_GPa[~tet_mask] + gibbs_values = np.empty_like(P_GPa) + gibbs_values[tet_mask] = 0.0 + gibbs_values[~tet_mask] = 0.5 * a * (P - Pastc) * Qsqr(P) + 0.25 * bast(P) * Qsqr( + P + ) * Qsqr(P) + return gibbs_values + + +def Qsqr(P_GPa): + """ + The value of Q squared, as defined in Equation 3 + and used in both Equations 1 and 4. + """ + tet_mask = np.array([P < Pastc for P in P_GPa]) + Q2 = np.empty_like(P_GPa) + Q2[tet_mask] = 0 + P = P_GPa[~tet_mask] + Q2[~tet_mask] = a / bast(P) * (Pastc - P) + return Q2 + + +def bast(P_GPa): + """ + The asterisked value of b, as given in + Equation 5 and used in Equation 4. + """ + return b - 2.0 * ( + ( + l3 * l3 * (C11(P_GPa) + C12(P_GPa)) + + 2.0 * l1 * l1 * C33(P_GPa) + - 4.0 * l1 * l3 * C13(P_GPa) + ) + / ((C11(P_GPa) + C12(P_GPa)) * C33(P_GPa) - 2.0 * C13(P_GPa) * C13(P_GPa)) + ) + + +def chi(P_GPa): + """ + The value of chi, as defined in + Equations 9 and 10 and used to define the Cijs. + """ + tet_mask = np.array([P < Pastc for P in P_GPa]) + chi_values = np.empty_like(P_GPa) + chi_values[tet_mask] = 1.0 / (a * (P_GPa[tet_mask] - Pc)) + chi_values[~tet_mask] = 1.0 / ( + 2.0 * a * b / bast(P_GPa[~tet_mask]) * (Pastc - P_GPa[~tet_mask]) + + a * (Pastc - Pc) + ) + return chi_values + + +def CT(P_GPa, mode="eqm"): + # Expressions in Table 1 + C = np.zeros((6, 6, len(P_GPa))) + + if mode == "eqm": + X = chi(P_GPa) + Q = np.sqrt(Qsqr(P_GPa)) + else: + X = np.ones_like(P_GPa) + Q = np.ones_like(P_GPa) + + tet_mask = np.array([(P < Pastc and mode == "eqm") or mode == "tet" for P in P_GPa]) + C[0, 0, tet_mask] = C11(P_GPa[tet_mask]) - l2 * l2 * X[tet_mask] + C[1, 1, tet_mask] = C[0, 0, tet_mask] + C[2, 2, tet_mask] = C33(P_GPa[tet_mask]) + C[0, 1, tet_mask] = C12(P_GPa[tet_mask]) + l2 * l2 * X[tet_mask] + C[0, 2, tet_mask] = C13(P_GPa[tet_mask]) + C[1, 2, tet_mask] = C[0, 2, tet_mask] + C[3, 3, tet_mask] = C44(P_GPa[tet_mask]) + C[4, 4, tet_mask] = C44(P_GPa[tet_mask]) + C[5, 5, tet_mask] = C66(P_GPa[tet_mask]) + + Q = Q[~tet_mask] + X = X[~tet_mask] + C[0, 0, ~tet_mask] = ( + C11(P_GPa[~tet_mask]) + - (4.0 * l1 * l1 * Q * Q + l2 * l2 + 4.0 * l1 * l2 * Q) * X + ) + C[1, 1, ~tet_mask] = ( + C11(P_GPa[~tet_mask]) + - (4.0 * l1 * l1 * Q * Q + l2 * l2 - 4.0 * l1 * l2 * Q) * X + ) + C[2, 2, ~tet_mask] = C33(P_GPa[~tet_mask]) - 4.0 * l3 * l3 * Q * Q * X + C[0, 1, ~tet_mask] = C12(P_GPa[~tet_mask]) - (4.0 * l1 * l1 * Q * Q - l2 * l2) * X + C[0, 2, ~tet_mask] = ( + C13(P_GPa[~tet_mask]) - (4 * l1 * l3 * Q * Q + 2.0 * l2 * l3 * Q) * X + ) + C[1, 2, ~tet_mask] = ( + C13(P_GPa[~tet_mask]) - (4 * l1 * l3 * Q * Q - 2.0 * l2 * l3 * Q) * X + ) + C[3, 3, ~tet_mask] = C44(P_GPa[~tet_mask]) + 2.0 * l4 * Q + C[4, 4, ~tet_mask] = C44(P_GPa[~tet_mask]) - 2.0 * l4 * Q + C[5, 5, ~tet_mask] = C66(P_GPa[~tet_mask]) + 2.0 * l6 * Q * Q + + C[1, 0] = C[0, 1] + C[2, 0] = C[0, 2] + C[2, 1] = C[1, 2] + + return np.moveaxis(C, 2, 0) + + +def strains(P_GPa): + tet_mask = np.array([P < Pastc for P in P_GPa]) + epsilons = np.zeros((3, len(P_GPa))) + P = P_GPa[~tet_mask] + Q2 = Qsqr(P) + Q = np.sqrt(Q2) + + e1se2 = -l2 / (0.5 * (C11(P) - C12(P))) * Q + e1pe2 = -l1 / (0.5 * (C11(P) + C12(P))) * Q2 + epsilons[0, ~tet_mask] = (e1se2 + e1pe2) / 2.0 + epsilons[1, ~tet_mask] = (-e1se2 + e1pe2) / 2.0 + epsilons[2, ~tet_mask] = -l3 / C33(P) * Q2 + return epsilons + + +pressures = np.linspace(0.0e9, 120.0e9, 1001) +temperatures = 300.0 + pressures * 0.0 + +Ceqm = CT(pressures / 1.0e9, mode="eqm") * 1.0e9 +Ctet = CT(pressures / 1.0e9, mode="tet") * 1.0e9 +Cort = CT(pressures / 1.0e9, mode="ort") * 1.0e9 +ST = np.linalg.inv(Ceqm) +KTR = 1.0 / np.sum(ST[:, :3, :3], axis=(1, 2)) +ST_tet = np.linalg.inv(Ctet) +KTR_tet = 1.0 / np.sum(ST_tet[:, :3, :3], axis=(1, 2)) + +eps = strains(pressures / 1.0e9) +plt.plot(pressures / 1.0e9, np.power(eps[0] - eps[1], 2.0)) +plt.show() + +for i, j in [(1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (1, 2), (1, 3), (2, 3)]: + (ln,) = plt.plot( + pressures / 1.0e9, Ceqm[:, i - 1, j - 1] / 1.0e9, label=f"$C_{{{i}{j}}}$" + ) + # plt.plot(pressures, Cort[:, i, j], + # label=f'$C_{{{i}{j}}}$', linestyle=':', c=ln.get_color()) + +plt.plot(pressures / 1.0e9, KTR / 1.0e9, label="$K_{{TR}}$") +plt.legend() +plt.show() + +V_0 = 1.401e-05 +# V_0 = 1.401e-05 / 6.0223e23 * 1.e30 * 2. +lnV = cumulative_trapezoid(-1.0 / KTR, pressures, initial=0) + np.log(V_0) +V_orig = np.exp(lnV) + +fig = plt.figure(figsize=(12, 4)) +ax = [fig.add_subplot(1, 3, i) for i in range(1, 4)] + +G = gibbs_over_tetragonal_volume_1(pressures / 1.0e9) * V_orig * 1.0e9 +ax[0].plot(pressures / 1.0e9, G, label="$\\mathcal{G}_{{eq 1}} \\cdot V_{{tet}}$") +G = gibbs_over_tetragonal_volume_4(pressures / 1.0e9) * V_orig * 1.0e9 +ax[0].plot(pressures / 1.0e9, G, label="$\\mathcal{G}_{{eq 4}} \\cdot V_{{tet}}$") + +ax[0].set_ylabel("Excess Gibbs energy ($\\mathcal{{G}}_{{xs}}$; J/mol)") + +ax[1].plot( + pressures / 1.0e9, + V_orig * 1.0e6, + linestyle=":", + linewidth=3.0, + label="tetragonal: from $V_0$ and $C_{{Tij0}}$s", +) + +V_diff = np.gradient(G, pressures, edge_order=2) +V = V_orig + V_diff + +ax[1].plot( + pressures / 1.0e9, + V * 1.0e6, + linestyle="--", + label="eqm: from $V_0 + C_{{Tij0}}$s + $\\partial \\mathcal{{G}}_{{xs}}/\\partial P$", +) + +V_HP = stv().evaluate(["V"], pressures, temperatures)[0] +# ax[1].plot(pressures / 1.0e9, V_HP * 1.0e6, linestyle=":", label="eqm: V(HP)") + +ax[1].set_ylabel("Volume (cm$^3$/mol)") + +KTR2 = -np.gradient(pressures, np.log(V)) + +ax[2].plot( + pressures / 1.0e9, + KTR_tet / 1.0e9, + linestyle=":", + linewidth=3.0, + label="tetragonal: from $C_{{Tij0}}$s", +) +ax[2].plot( + pressures / 1.0e9, KTR / 1.0e9, label="eqm: from eqm $C_{{Tij}}$s", c="black" +) +ax[2].plot( + pressures / 1.0e9, + KTR2 / 1.0e9, + linestyle="--", + label="eqm: from $-V (\\partial P / \\partial V)_T$", +) + +ax[2].set_ylabel("Isothermal Reuss bulk modulus (GPa)") + +for i in range(3): + ax[i].set_title(f"({ascii_lowercase[i]})") + ax[i].set_xlabel("Pressure (GPa)") + ax[i].legend() + +fig.set_tight_layout(True) +fig.savefig("figures/Carpenter_2000_model.pdf") +plt.show() diff --git a/contrib/anisotropic_stishovite/stishovite_model_Zhang_2021.py b/contrib/anisotropic_stishovite/stishovite_model_Zhang_2021.py new file mode 100644 index 00000000..12c50ec6 --- /dev/null +++ b/contrib/anisotropic_stishovite/stishovite_model_Zhang_2021.py @@ -0,0 +1,350 @@ +import numpy as np +from scipy.integrate import cumulative_trapezoid +import matplotlib.pyplot as plt +from burnman.minerals.HP_2011_ds62 import stv +from burnman.minerals.SLB_2011 import periclase +from string import ascii_lowercase + + +stishovite = periclase() +stishovite.params["V_0"] = 46.57e-6 # only relative matters +stishovite.params["K_0"] = 317.2e9 +stishovite.params["Kprime_0"] = 4.8 + + +V_0p = stishovite.params["V_0"] +K_0p = stishovite.params["K_0"] / 1.0e9 +Kprime_0p = stishovite.params["Kprime_0"] + +dijdkl = np.einsum("ij, kl->ijkl", np.eye(3), np.eye(3)) +kron = -(dijdkl + np.einsum("ijkl->iljk", dijdkl) + np.einsum("ijkl->jlik", dijdkl)) + + +def Cij(P_GPa, Cij0, Cijprime0, dijkl): + if len(P_GPa) > 0: + V = stishovite.evaluate(["V"], P_GPa * 1.0e9, P_GPa * 0.0 + 300.0)[0] + fE = (np.power(V_0p / V, 2.0 / 3.0) - 1) / 2.0 + Cij = np.power(1.0 + 2 * fE, 2.5) * ( + Cij0 + + (3.0 * K_0p * Cijprime0 - 5.0 * Cij0) * fE + + ( + 6.0 * K_0p * Cijprime0 + - 14 * Cij0 + - 3.0 / 2.0 * K_0p * dijkl * (3.0 * Kprime_0p - 16.0) + ) + * fE + * fE + ) + return Cij + else: + return [] + + +# Expressions for the bare elastic constant values +# are taken from Table SVI +# All units are GPa + +Pastc = 55.0 +Pc = 55.2 + 55 + +l1 = 11.03 +l2 = 27.61 +l3 = 16.79 +l4 = 18.94 +l6 = 15.15 +a = -0.0501 +b = 11.0 + + +def C11(P_GPa): + return Cij(P_GPa, 592.3, 10.8, kron[0, 0, 0, 0]) + + +def C12(P_GPa): + return Cij(P_GPa, 57.9, 8.81, kron[0, 0, 1, 1]) + + +def C13(P_GPa): + return Cij(P_GPa, 193, 2.91, kron[0, 0, 2, 2]) + + +def C33(P_GPa): + return Cij(P_GPa, 760.2, 7.07, kron[2, 2, 2, 2]) + + +def C44(P_GPa): + return Cij(P_GPa, 261.6, 3.18, kron[1, 2, 1, 2]) + + +def C66(P_GPa): + return Cij(P_GPa, 319.7, 5.60, kron[0, 1, 0, 1]) + + +def gibbs_over_tetragonal_volume_1(P_GPa): + """ + The Gibbs energy divided by the + volume of the (bare) tetragonal phase, + as defined in Equation 1. + """ + tet_mask = np.array([P < Pastc for P in P_GPa]) + gibbs_values = np.empty_like(P_GPa) + gibbs_values[tet_mask] = 0.0 + + P = P_GPa[~tet_mask] + Q2 = Qsqr(P) + Q = np.sqrt(Q2) + Q4 = Q2 * Q2 + + e1se2 = -l2 / (0.5 * (C11(P) - C12(P))) * Q + e1pe2 = -l1 / (0.5 * (C11(P) + C12(P))) * Q2 + e1 = (e1se2 + e1pe2) / 2.0 + e2 = (-e1se2 + e1pe2) / 2.0 + e3 = -l3 / C33(P) * Q2 + + gibbs_values[~tet_mask] = ( + 0.5 * a * (P - Pc) * Q2 + + 0.25 * b * Q4 + + l1 * (e1 + e2) * Q2 + + l2 * (e1 - e2) * Q + + l3 * e3 * Q2 + + 0.25 * (C11(P) + C12(P)) * e1pe2 * e1pe2 + + 0.25 * (C11(P) - C12(P)) * e1se2 * e1se2 + + C13(P) * e1pe2 * e3 + + 0.5 * C33(P) * e3 * e3 + ) + return gibbs_values + + +def gibbs_over_tetragonal_volume_4(P_GPa): + """ + The Gibbs energy divided by the + volume of the (bare) tetragonal phase, + as defined in Equation 4. + """ + tet_mask = np.array([P < Pastc for P in P_GPa]) + P = P_GPa[~tet_mask] + gibbs_values = np.empty_like(P_GPa) + gibbs_values[tet_mask] = 0.0 + gibbs_values[~tet_mask] = 0.5 * a * (P - Pastc) * Qsqr(P) + 0.25 * bast(P) * Qsqr( + P + ) * Qsqr(P) + return gibbs_values + + +def Qsqr(P_GPa): + """ + The value of Q squared, as defined in Equation 3 + and used in both Equations 1 and 4. + """ + tet_mask = np.array([P < Pastc for P in P_GPa]) + Q2 = np.empty_like(P_GPa) + Q2[tet_mask] = 0 + P = P_GPa[~tet_mask] + Q2[~tet_mask] = a / bast(P) * (Pastc - P) + return Q2 + + +def bast(P_GPa): + """ + The asterisked value of b, as given in + Equation 5 and used in Equation 4. + """ + return b - 2.0 * ( + ( + l3 * l3 * (C11(P_GPa) + C12(P_GPa)) + + 2.0 * l1 * l1 * C33(P_GPa) + - 4.0 * l1 * l3 * C13(P_GPa) + ) + / ((C11(P_GPa) + C12(P_GPa)) * C33(P_GPa) - 2.0 * C13(P_GPa) * C13(P_GPa)) + ) + + +def chi(P_GPa): + """ + The value of chi, as defined in + Equations 9 and 10 and used to define the Cijs. + """ + tet_mask = np.array([P < Pastc for P in P_GPa]) + chi_values = np.empty_like(P_GPa) + chi_values[tet_mask] = 1.0 / (a * (P_GPa[tet_mask] - Pc)) + chi_values[~tet_mask] = 1.0 / ( + 2.0 * a * b / bast(P_GPa[~tet_mask]) * (Pastc - P_GPa[~tet_mask]) + + a * (Pastc - Pc) + ) + return chi_values + + +def CT(P_GPa, mode="eqm"): + # Expressions in Table 1 + C = np.zeros((6, 6, len(P_GPa))) + + if mode == "eqm": + X = chi(P_GPa) + Q = np.sqrt(Qsqr(P_GPa)) + else: + X = np.ones_like(P_GPa) + Q = np.ones_like(P_GPa) + + tet_mask = np.array([(P < Pastc and mode == "eqm") or mode == "tet" for P in P_GPa]) + C[0, 0, tet_mask] = C11(P_GPa[tet_mask]) - l2 * l2 * X[tet_mask] + C[1, 1, tet_mask] = C[0, 0, tet_mask] + C[2, 2, tet_mask] = C33(P_GPa[tet_mask]) + C[0, 1, tet_mask] = C12(P_GPa[tet_mask]) + l2 * l2 * X[tet_mask] + C[0, 2, tet_mask] = C13(P_GPa[tet_mask]) + C[1, 2, tet_mask] = C[0, 2, tet_mask] + C[3, 3, tet_mask] = C44(P_GPa[tet_mask]) + C[4, 4, tet_mask] = C44(P_GPa[tet_mask]) + C[5, 5, tet_mask] = C66(P_GPa[tet_mask]) + + Q = Q[~tet_mask] + X = X[~tet_mask] + C[0, 0, ~tet_mask] = ( + C11(P_GPa[~tet_mask]) + - (4.0 * l1 * l1 * Q * Q + l2 * l2 + 4.0 * l1 * l2 * Q) * X + ) + C[1, 1, ~tet_mask] = ( + C11(P_GPa[~tet_mask]) + - (4.0 * l1 * l1 * Q * Q + l2 * l2 - 4.0 * l1 * l2 * Q) * X + ) + C[2, 2, ~tet_mask] = C33(P_GPa[~tet_mask]) - 4.0 * l3 * l3 * Q * Q * X + C[0, 1, ~tet_mask] = C12(P_GPa[~tet_mask]) - (4.0 * l1 * l1 * Q * Q - l2 * l2) * X + C[0, 2, ~tet_mask] = ( + C13(P_GPa[~tet_mask]) - (4 * l1 * l3 * Q * Q + 2.0 * l2 * l3 * Q) * X + ) + C[1, 2, ~tet_mask] = ( + C13(P_GPa[~tet_mask]) - (4 * l1 * l3 * Q * Q - 2.0 * l2 * l3 * Q) * X + ) + C[3, 3, ~tet_mask] = C44(P_GPa[~tet_mask]) + 2.0 * l4 * Q + C[4, 4, ~tet_mask] = C44(P_GPa[~tet_mask]) - 2.0 * l4 * Q + C[5, 5, ~tet_mask] = C66(P_GPa[~tet_mask]) + 2.0 * l6 * Q * Q + + C[1, 0] = C[0, 1] + C[2, 0] = C[0, 2] + C[2, 1] = C[1, 2] + + return np.moveaxis(C, 2, 0) + + +def strains(P_GPa): + tet_mask = np.array([P < Pastc for P in P_GPa]) + epsilons = np.zeros((3, len(P_GPa))) + P = P_GPa[~tet_mask] + Q2 = Qsqr(P) + Q = np.sqrt(Q2) + + e1se2 = -l2 / (0.5 * (C11(P) - C12(P))) * Q + e1pe2 = -l1 / (0.5 * (C11(P) + C12(P))) * Q2 + epsilons[0, ~tet_mask] = (e1se2 + e1pe2) / 2.0 + epsilons[1, ~tet_mask] = (-e1se2 + e1pe2) / 2.0 + epsilons[2, ~tet_mask] = -l3 / C33(P) * Q2 + return epsilons + + +pressures = np.linspace(0.0e9, 120.0e9, 1001) +temperatures = 300.0 + pressures * 0.0 + +Ceqm = CT(pressures / 1.0e9, mode="eqm") * 1.0e9 +Ctet = CT(pressures / 1.0e9, mode="tet") * 1.0e9 +Cort = CT(pressures / 1.0e9, mode="ort") * 1.0e9 +ST = np.linalg.inv(Ceqm) +KTR = 1.0 / np.sum(ST[:, :3, :3], axis=(1, 2)) +ST_tet = np.linalg.inv(Ctet) +KTR_tet = 1.0 / np.sum(ST_tet[:, :3, :3], axis=(1, 2)) + +eps = strains(pressures / 1.0e9) +plt.plot(pressures / 1.0e9, np.power(eps[0] - eps[1], 2.0)) +plt.show() + +for i, j in [(1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (1, 2), (1, 3), (2, 3)]: + (ln,) = plt.plot( + pressures / 1.0e9, Ceqm[:, i - 1, j - 1] / 1.0e9, label=f"$C_{{{i}{j}}}$" + ) + # plt.plot(pressures, Cort[:, i, j], + # label=f'$C_{{{i}{j}}}$', linestyle=':', c=ln.get_color()) + +plt.plot(pressures / 1.0e9, KTR / 1.0e9, label="$K_{{TR}}$") + + +d = np.loadtxt("data/Zhang_2021_stishovite_elastic_tensor.dat") +# C11 C12 C13 C22 C23 C33 C44 C55 C66 +plt.scatter(d[:, 0], d[:, 2]) +plt.scatter(d[:, 0], d[:, 8]) +plt.scatter(d[:, 0], d[:, 12]) +plt.scatter(d[:, 0], d[:, 14]) +plt.scatter(d[:, 0], d[:, 16]) +plt.scatter(d[:, 0], d[:, 18]) +plt.scatter(d[:, 0], d[:, 4]) +plt.scatter(d[:, 0], d[:, 6]) +plt.scatter(d[:, 0], d[:, 10]) + + +plt.legend() +plt.show() + +V_0 = 1.401e-05 +# V_0 = 1.401e-05 / 6.0223e23 * 1.e30 * 2. +lnV = cumulative_trapezoid(-1.0 / KTR, pressures, initial=0) + np.log(V_0) +V_orig = np.exp(lnV) + +fig = plt.figure(figsize=(12, 4)) +ax = [fig.add_subplot(1, 3, i) for i in range(1, 4)] + +G = gibbs_over_tetragonal_volume_1(pressures / 1.0e9) * V_orig * 1.0e9 +ax[0].plot(pressures / 1.0e9, G, label="$\\mathcal{G}_{{eq 1}} \\cdot V_{{tet}}$") +G = gibbs_over_tetragonal_volume_4(pressures / 1.0e9) * V_orig * 1.0e9 +ax[0].plot(pressures / 1.0e9, G, label="$\\mathcal{G}_{{eq 4}} \\cdot V_{{tet}}$") + +ax[0].set_ylabel("Excess Gibbs energy ($\\mathcal{{G}}_{{xs}}$; J/mol)") + +ax[1].plot( + pressures / 1.0e9, + V_orig * 1.0e6, + linestyle=":", + linewidth=3.0, + label="tetragonal: from $V_0$ and $C_{{Tij0}}$s", +) + +V_diff = np.gradient(G, pressures, edge_order=2) +V = V_orig + V_diff + +ax[1].plot( + pressures / 1.0e9, + V * 1.0e6, + linestyle="--", + label="eqm: from $V_0 + C_{{Tij0}}$s + $\\partial \\mathcal{{G}}_{{xs}}/\\partial P$", +) + +V_HP = stv().evaluate(["V"], pressures, temperatures)[0] +# ax[1].plot(pressures / 1.0e9, V_HP * 1.0e6, linestyle=":", label="eqm: V(HP)") + +ax[1].set_ylabel("Volume (cm$^3$/mol)") + +KTR2 = -np.gradient(pressures, np.log(V)) + +ax[2].plot( + pressures / 1.0e9, + KTR_tet / 1.0e9, + linestyle=":", + linewidth=3.0, + label="tetragonal: from $C_{{Tij0}}$s", +) +ax[2].plot( + pressures / 1.0e9, KTR / 1.0e9, label="eqm: from eqm $C_{{Tij}}$s", c="black" +) +ax[2].plot( + pressures / 1.0e9, + KTR2 / 1.0e9, + linestyle="--", + label="eqm: from $-V (\\partial P / \\partial V)_T$", +) + +ax[2].set_ylabel("Isothermal Reuss bulk modulus (GPa)") + +for i in range(3): + ax[i].set_title(f"({ascii_lowercase[i]})") + ax[i].set_xlabel("Pressure (GPa)") + ax[i].legend() + +fig.set_tight_layout(True) +fig.savefig("figures/Zhang_2021_model.pdf") +plt.show() diff --git a/contrib/anisotropic_stishovite/stishovite_model_covariance_matrix.py b/contrib/anisotropic_stishovite/stishovite_model_covariance_matrix.py new file mode 100644 index 00000000..42f67123 --- /dev/null +++ b/contrib/anisotropic_stishovite/stishovite_model_covariance_matrix.py @@ -0,0 +1,98 @@ +import numpy as np +import pandas as pd +import seaborn as sns +import matplotlib.pyplot as plt +from stishovite_parameters import ( + all_args, + scalar_param_names, + cell_param_names, + all_param_names, +) +import burnman +from tabulate import tabulate + +popt = all_args +pcov = np.loadtxt("model_output/covariance_matrix.dat") + +stv_SLB = burnman.minerals.SLB_2022.st() +stv_SLB.property_modifiers = [] + +popt[0] += stv_SLB.params["V_0"] * 1.0e6 +popt[1] += stv_SLB.params["K_0"] / 1.0e9 +popt[2] += stv_SLB.params["Kprime_0"] +popt[3] += stv_SLB.params["grueneisen_0"] +popt[4] += stv_SLB.params["q_0"] + +sigma = np.sqrt(np.diag(pcov)) +rel_error = [ + "{:g}".format(float("{:.1g}".format(i))) for i in sigma / np.abs(popt) * 100.0 +] + +table = [all_param_names, popt, sigma, rel_error] +table = list(zip(*table)) +print( + tabulate( + table, + headers=["Parameter", "Value", "Sigma", "Relative error (\\%, 1 s.f.)"], + tablefmt="latex_raw", + floatfmt=["s", ".6e", ".6e"], + ) +) + + +def correlation_from_covariance(covariance): + v = np.sqrt(np.diag(covariance)) + outer_v = np.outer(v, v) + correlation = covariance / outer_v + correlation[covariance == 0] = 0 + return correlation + + +corr = np.around(correlation_from_covariance(pcov) * 100.0, 0) + +d = pd.DataFrame(data=corr, columns=all_param_names, index=all_param_names) + + +sns.set_theme(style="white") + +# Set up the matplotlib figure +fig, ax = plt.subplots(figsize=(11, 9)) + +# Generate a custom diverging colormap +cmap = sns.diverging_palette(230, 20, as_cmap=True) + +# Draw the heatmap with the mask and correct aspect ratio +sns.heatmap( + d, + cmap=cmap, + center=0, + vmin=-100, + annot=True, + fmt="g", + square=True, + linewidths=0.5, + annot_kws={"size": 6}, + cbar_kws={"shrink": 0.5, "label": "Correlation coefficient (%)"}, +) + +ax.hlines( + [ + len(scalar_param_names) - 8, + len(scalar_param_names), + len(scalar_param_names + cell_param_names), + ], + *ax.get_xlim() +) +ax.vlines( + [ + len(scalar_param_names) - 8, + len(scalar_param_names), + len(scalar_param_names + cell_param_names), + ], + *ax.get_ylim() +) + + +fig.set_tight_layout(True) +fig.savefig("figures/correlation_matrix.pdf") +plt.show() diff --git a/contrib/anisotropic_stishovite/stishovite_model_plots.py b/contrib/anisotropic_stishovite/stishovite_model_plots.py index 145c7f12..eddb1984 100644 --- a/contrib/anisotropic_stishovite/stishovite_model_plots.py +++ b/contrib/anisotropic_stishovite/stishovite_model_plots.py @@ -18,18 +18,29 @@ models = get_models() _, stishovite_scalar, stishovite_unrelaxed, stishovite_relaxed = models[0] -fP_Zhang, fP_Andrault, fP2_Zhang, fP2_Andrault = models[1:] +( + fP_Zhang, + fP_Andrault, + fP_Wang, + fP_Nishihara, + fP2_Zhang, + fP2_Andrault, + fP2_Wang, + fP2_Nishihara, +) = models[1:] stishovite_relaxed.set_composition([1.0]) stishovite_unrelaxed.set_composition([0.5, 0.5]) -plot_isotropic_velocities = True -plot_seismic_properties = True -low_res_seismic_properties = True +plot_isotropic_velocities = False +plot_seismic_properties = False +low_res_seismic_properties = False plot_relaxed = True -plot_Q_V_abc = True -plot_G = True -plot_lnabc = True +plot_Q_V_abc = False +plot_G = False +plot_lnabc = False + +show_plots_one_by_one = False if plot_isotropic_velocities: pressures = np.linspace(9.0e9, 128e9, 501) @@ -93,12 +104,15 @@ ax[0].legend() fig.set_tight_layout(True) fig.savefig("figures/stv_isotropic_velocities.pdf") - plt.show() + if show_plots_one_by_one: + plt.show() if plot_seismic_properties: T = 2200.0 # for delta_P in [0.1e9]: - for delta_P in [-40.0e9, -1.0e9, -0.1e9, 0.1e9, 1.0e9, 40.0e9]: + # for delta_P in [-40.0e9, -1.0e9, -0.1e9, 0.1e9, 1.0e9, 40.0e9]: + # for delta_P in [-40.0e9, -30.e9, -20.e9, -10.e9, -1.0e9, -0.1e9, 0.1e9, 1.0e9, 10.e9, 20.e9, 30.e9, 40.0e9]: + for delta_P in [-10.0e9]: print( f"Plotting seismic properties at {delta_P/1.e9:.2f} GPa above the transition" ) @@ -143,13 +157,19 @@ fig.savefig( f"figures/stishovite_seismic_properties_{P/1.e9:.2f}_GPa_{int(T)}_K.pdf" ) - plt.show() + fig.savefig( + f"figures/stishovite_seismic_properties_{P/1.e9:.2f}_GPa_{int(T)}_K.jpg" + ) + if show_plots_one_by_one: + plt.show() data = common_data() PTV = np.concatenate( ( data["PTV"]["stv"]["Ito_1974"], + data["PTV"]["stv"]["Wang_2012"], + data["PTV"]["stv"]["Nishihara_2005"], data["PTV"]["stv"]["Zhang_2021"], data["PTV"]["poststv"]["Zhang_2021"], data["PTV"]["stv"]["Andrault_2003"], @@ -161,6 +181,8 @@ PTV_err = np.concatenate( ( data["PTV_err"]["stv"]["Ito_1974"], + data["PTV_err"]["stv"]["Wang_2012"], + data["PTV_err"]["stv"]["Nishihara_2005"], data["PTV_err"]["stv"]["Zhang_2021"], data["PTV_err"]["poststv"]["Zhang_2021"], data["PTV_err"]["stv"]["Andrault_2003"], @@ -172,6 +194,8 @@ abc = np.concatenate( ( data["abc"]["stv"]["Ito_1974"], + data["abc"]["stv"]["Wang_2012"], + data["abc"]["stv"]["Nishihara_2005"], data["abc"]["stv"]["Zhang_2021"], data["abc"]["poststv"]["Zhang_2021"], data["abc"]["stv"]["Andrault_2003"], @@ -183,6 +207,8 @@ abc_err = np.concatenate( ( data["abc_err"]["stv"]["Ito_1974"], + data["abc_err"]["stv"]["Wang_2012"], + data["abc_err"]["stv"]["Nishihara_2005"], data["abc_err"]["stv"]["Zhang_2021"], data["abc_err"]["poststv"]["Zhang_2021"], data["abc_err"]["stv"]["Andrault_2003"], @@ -195,22 +221,29 @@ id = np.concatenate( ( np.ones(len(data["abc_err"]["stv"]["Ito_1974"])) * 1, - np.ones(len(data["abc_err"]["stv"]["Zhang_2021"])) * 2, - np.ones(len(data["abc_err"]["poststv"]["Zhang_2021"])) * 2, - np.ones(len(data["abc_err"]["stv"]["Andrault_2003"])) * 3, - np.ones(len(data["abc_err"]["poststv"]["Andrault_2003"])) * 3, - np.ones(len(data["abc_err"]["stv"]["Fischer_2018"])) * 4, - np.ones(len(data["abc_err"]["poststv"]["Fischer_2018"])) * 4, + np.ones(len(data["abc_err"]["stv"]["Wang_2012"])) * 2, + np.ones(len(data["abc_err"]["stv"]["Nishihara_2005"])) * 3, + np.ones(len(data["abc_err"]["stv"]["Zhang_2021"])) * 4, + np.ones(len(data["abc_err"]["poststv"]["Zhang_2021"])) * 4, + np.ones(len(data["abc_err"]["stv"]["Andrault_2003"])) * 5, + np.ones(len(data["abc_err"]["poststv"]["Andrault_2003"])) * 5, + np.ones(len(data["abc_err"]["stv"]["Fischer_2018"])) * 6, + np.ones(len(data["abc_err"]["poststv"]["Fischer_2018"])) * 6, ) ) Ito_mask = id == 1 -Zhang_mask = id == 2 -Andrault_mask = id == 3 -Fischer_mask = id == 4 +Wang_mask = id == 2 +Nishihara_mask = id == 3 +Zhang_mask = id == 4 +Andrault_mask = id == 5 +Fischer_mask = id == 6 PTV[Zhang_mask, 0] *= fP_Zhang + PTV[Zhang_mask, 0] * fP2_Zhang PTV[Andrault_mask, 0] *= fP_Andrault + PTV[Andrault_mask, 0] * fP2_Andrault +PTV[Wang_mask, 0] *= fP_Wang + PTV[Wang_mask, 0] * fP2_Wang +PTV[Nishihara_mask, 0] *= fP_Nishihara + PTV[Nishihara_mask, 0] * fP2_Nishihara +P_for_CN_orig = deepcopy(P_for_CN) P_for_CN *= fP_Zhang + P_for_CN * fP2_Zhang print( @@ -221,17 +254,28 @@ print(f"Consistent?: {check_anisotropic_eos_consistency(stishovite_relaxed)}") if plot_lnabc: - fig_lnabc = plt.figure(figsize=(8, 4)) + fig_lnabc = plt.figure(figsize=(10, 4)) ax_lnabc = [fig_lnabc.add_subplot(1, 2, i) for i in range(1, 3)] - colors = ["tab:blue", "tab:red", "tab:purple", "tab:orange"] + colors = [ + "tab:olive", + "tab:cyan", + "tab:orange", + "tab:blue", + "tab:purple", + "tab:red", + ] labels = [ - "Ito et al. (1974)", "Zhang et al. (2021)", - "Andrault et al. (2003)", "Fischer et al. (2018)", + "Wang et al. (2012)", + "Nishihara et al. (2005)", + "Andrault et al. (2003)", + "Ito et al. (1974)", ] - for i, mask in enumerate([Ito_mask, Zhang_mask, Andrault_mask, Fischer_mask]): + for i, mask in enumerate( + [Zhang_mask, Fischer_mask, Wang_mask, Nishihara_mask, Andrault_mask, Ito_mask] + ): c = colors[i] lnV = deepcopy(np.log(PTV[mask, 2])) @@ -277,10 +321,12 @@ ax_lnabc[0].set_ylabel("$\\ln(a)$, $\\ln(b)$ (cm/mol$^{\\frac{1}{3}}$)") ax_lnabc[1].set_ylabel("$\\ln(c)$ (cm/mol$^{\\frac{1}{3}}$)") - ax_lnabc[1].legend() + handles, labels = ax_lnabc[1].get_legend_handles_labels() + ax_lnabc[1].legend(handles[::-1], labels[::-1]) fig_lnabc.set_tight_layout(True) fig_lnabc.savefig("figures/stv_lnabc.pdf") - plt.show() + if show_plots_one_by_one: + plt.show() if plot_G: @@ -294,17 +340,29 @@ G_xs = stishovite_scalar.evaluate( ["excess_gibbs"], pressure_grid, T_grid, molar_fractions )[0] + + n_lines = len(pressures) + cmap = matplotlib.colormaps["rainbow"] + + # Take colors at regular intervals spanning the colormap. + colors = cmap(np.linspace(0, 1, n_lines)) + for i in range(len(pressures)): if i % 2 == 0: ax_G[0].plot( - p_grid[:, i], G_xs[:, i] / 1000.0, label=f"{pressures[i]/1.e9} GPa" + p_grid[:, i], + G_xs[:, i] / 1000.0, + label=f"{pressures[i]/1.e9} GPa", + color=colors[i], ) else: - ax_G[0].plot(p_grid[:, i], G_xs[:, i] / 1000.0, linestyle=":") + ax_G[0].plot( + p_grid[:, i], G_xs[:, i] / 1000.0, linestyle=":", color=colors[i] + ) imin = np.argmin(G_xs[:, i]) G_min = G_xs[imin, i] / 1.0e3 p_min = p_grid[imin, i] - ax_G[0].scatter([p_min, 1.0 - p_min], [G_min, G_min]) + ax_G[0].scatter([p_min, 1.0 - p_min], [G_min, G_min], color=colors[i]) ax_G[0].set_xlim(-0.25, 1.25) ax_G[0].set_xlabel("$p_{Q1}$") ax_G[0].set_ylabel("$\\mathcal{G}_{xs}$ (kJ/mol)") @@ -317,7 +375,8 @@ ax_G2.set_xlabel("$Q$") fig_G.set_tight_layout(True) fig_G.savefig("figures/stv_G.pdf") - plt.show() + if show_plots_one_by_one: + plt.show() if plot_Q_V_abc: @@ -430,6 +489,8 @@ def Q_misfit(args): fig_Q.savefig("figures/stv_Q.pdf") fig_V.savefig("figures/stv_V.pdf") fig_abc.savefig("figures/stv_abc.pdf") + if show_plots_one_by_one: + plt.show() if plot_relaxed: fig_relaxed = plt.figure(figsize=(12, 4)) @@ -467,14 +528,14 @@ def Q_misfit(args): ["isentropic_stiffness_tensor", "isentropic_bulk_modulus_reuss"], pressures, Ts ) - c = ax_relaxed[1].plot(pressures / 1.0e9, K_NR / 1.0e9, label="$K_{NR}$") + c = ax_relaxed[1].plot(pressures / 1.0e9, K_NR / 1.0e9, label="$K_{SR}$") ax_relaxed[1].plot( pressures / 1.0e9, K_NR_u / 1.0e9, linestyle=":", c=c[0].get_color() ) - ax_relaxed[1].scatter(P_for_CN / 1.0e9, KNR_GPa, label="$K_{NR}$") + ax_relaxed[1].scatter(P_for_CN / 1.0e9, KNR_GPa, label="$K_{SR}$") - ax_Q0[1].plot(pressures / 1.0e9, K_NR_Q0 / 1.0e9, label="$K_{NR} (Q0)$") - ax_Q0[1].plot(pressures / 1.0e9, K_NR_Q1 / 1.0e9, label="$K_{NR} (Q1)$") + ax_Q0[1].plot(pressures / 1.0e9, K_NR_Q0 / 1.0e9, label="$K_{SR} (Q0)$") + ax_Q0[1].plot(pressures / 1.0e9, K_NR_Q1 / 1.0e9, label="$K_{SR} (Q1)$") for axi, i, j in ( (0, 0, 0), @@ -488,7 +549,7 @@ def Q_misfit(args): (2, 5, 5), ): c = ax_relaxed[axi].plot( - pressures / 1.0e9, C_N[:, i, j] / 1.0e9, label=f"$C_{{N{i+1}{j+1}}}$" + pressures / 1.0e9, C_N[:, i, j] / 1.0e9, label=f"$C_{{S{i+1}{j+1}}}$" ) color = c[0].get_color() ax_relaxed[axi].plot( @@ -497,8 +558,10 @@ def Q_misfit(args): ax_relaxed[axi].scatter( P_for_CN / 1.0e9, CN_GPa[:, i, j], label=f"{i+1}{j+1}", color=color ) + + # Plot original data ax_relaxed[axi].errorbar( - P_for_CN / 1.0e9, + P_for_CN_orig / 1.0e9, CN_GPa_orig[:, i, j], yerr=CN_err_GPa[:, i, j], alpha=0.5, @@ -506,7 +569,7 @@ def Q_misfit(args): ls="none", ) ax_relaxed[axi].scatter( - P_for_CN / 1.0e9, CN_GPa_orig[:, i, j], s=5, alpha=0.5, color=color + P_for_CN_orig / 1.0e9, CN_GPa_orig[:, i, j], s=5, alpha=0.5, color=color ) ax_Q0[axi].plot( @@ -531,4 +594,8 @@ def Q_misfit(args): fig_relaxed.savefig("figures/stv_relaxed.pdf") + if show_plots_one_by_one: + plt.show() + +if not show_plots_one_by_one: plt.show() diff --git a/contrib/anisotropic_stishovite/stishovite_parameters.py b/contrib/anisotropic_stishovite/stishovite_parameters.py index 316c8c06..9d93e02a 100644 --- a/contrib/anisotropic_stishovite/stishovite_parameters.py +++ b/contrib/anisotropic_stishovite/stishovite_parameters.py @@ -1,68 +1,108 @@ import numpy as np -# dVQ0, dKQ0, dKpQ0, dgrQ0, dqQ0, -# V0Q1overV0Q0, dDebye_0, P_tr_GPa, -# fP_Zhang, fP_Andrault, fP2_Zhang, fP2_Andrault -# 74 +scalar_param_names = [ + "$V_0 (Q=0)$", + "$K_0 (Q=0)$", + "$K'_0 (Q=0)$", + "$\\gamma_0 (Q=0)$", + "$q_0 (Q=0)$", + "$V_0 (Q=1) / V_0 (Q=0)$", + "$\\Theta_0 (Q=1) - \\Theta_0 (Q=0)$", + "$P_{{tr}} (T = T_{{ref}})$", + "$f_{{P}}$ (Zhang)", + "$f_{{P}}$ (Andrault)", + "$f_{{P}}$ (Wang)", + "$f_{{P}}$ (Nishihara)", + "$f_{{P2}}$ (Zhang)", + "$f_{{P2}}$ (Andrault)", + "$f_{{P2}}$ (Wang)", + "$f_{{P2}}$ (Nishihara)", +] + +# 76 scalar_args = [ - -5.28259962e-03, - -2.78232366e00, - 1.86241455e-03, - -2.16473596e-01, - -8.76004122e-02, - 9.93918688e-01, - 1.76286306e01, - 4.97878625e01, - 9.30217075e-01, - 9.93494742e-01, - -4.54240420e-13, - -9.20410207e-13, + -1.15714715e-02, + -2.49554064e00, + -7.69118482e-04, + -2.08731542e-01, + -6.08850674e-01, + 9.93834483e-01, + 1.76227286e01, + 4.99610469e01, + 8.91475407e-01, + 9.89762819e-01, + 9.87845187e-01, + 9.95720410e-01, + 2.92917650e-13, + -9.08667999e-13, + -2.20268875e-13, + 1.55926921e-12, ] -# a0Q1, b0Q1, -# PsiI_33_a, PsiI_33_b, PsiI_33_c, -# PsiI_33_b2, PsiI_33_c2 -# 1151 +cell_param_names = [ + "$a_0 (Q=1)$", + "$b_0 (Q=1)$", + "$\\Psi_{{3a}}$", + "$\\Psi_{{3b}}$", + "$\\Psi_{{3c}}$", + "$\\Psi_{{3b2}}$", + "$\\Psi_{{3c2}}$", +] + +# 1117 cell_args = [ - 2.72829954e-02, - 2.85720276e-02, - 2.10496560e-01, - -6.02250364e-01, - 1.11223359e00, - 8.20161603e-05, - -1.04689736e01, + 2.72779892e-02, + 2.85698440e-02, + 2.10385858e-01, + -6.02813854e-01, + 1.10903361e00, + 8.86818639e-05, + -1.09619314e01, ] -# a11, a22, a33, a44, a55, a66 -# b11, b33, b44, b66 -# c44, c66 -# b112, b332 +elastic_param_names = [ + "$\\Psi_{{11a}}$", + "$\\Psi_{{22a}}$", + "$\\Psi_{{33a}}$", + "$\\Psi_{{44a}}$", + "$\\Psi_{{55a}}$", + "$\\Psi_{{66a}}$", + "$\\Psi_{{11b}}$", + "$\\Psi_{{33b}}$", + "$\\Psi_{{44b}}$", + "$\\Psi_{{66b}}$", + "$\\Psi_{{44c}}$", + "$\\Psi_{{66c}}$", + "$\\Psi_{{11b2}}$", + "$\\Psi_{{33b2}}$", +] # let the b parameters be equal to # frel = -0.14 # b1_input = b1 * c1 * exp(c1 * frel) - 1) -# 178 +# 197 elastic_args = [ - 3.21306666e-01, - 8.34212845e-01, - 4.93022578e-01, - 1.12011984e00, - 1.20957700e00, - 9.49364288e-01, - -4.47565581e-02, - 7.06681926e-04, - 3.85553457e-01, - 1.13867050e-01, - 9.74629124e-01, - 9.98028989e-01, - 5.64360733e-01, - 2.12770903e-01, + 3.24807928e-01, + 8.37663793e-01, + 4.93000185e-01, + 1.11959128e00, + 1.20914946e00, + 9.49237642e-01, + -4.29973046e-02, + 7.10161257e-04, + 3.85475084e-01, + 1.13312302e-01, + 9.74598492e-01, + 9.97390953e-01, + 5.85741559e-01, + 2.13037695e-01, ] scalar_and_cell_args = np.concatenate((scalar_args, cell_args)) all_args = np.concatenate((scalar_args, cell_args, elastic_args)) +all_param_names = scalar_param_names + cell_param_names + elastic_param_names scalar_bounds = ( (-8.0e-3, -4.0e-3), # dVQ0 @@ -75,6 +115,10 @@ (46, 53), # P_tr_GPa, (0.9, 1.1), # s (0.9, 1.1), # s + (0.9, 1.1), # s + (0.9, 1.1), # s + (-0.1, 0.1), # s + (-0.1, 0.1), (-0.1, 0.1), # s (-0.1, 0.1), ) # s