From 4a3a62ef0da7afd482ec215c51fc379cf0aef0f9 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Mon, 26 Feb 2024 21:15:19 +0000 Subject: [PATCH] anisotropic stishovite model --- contrib/anisotropic_stishovite/Readme.md | 36 + ...ndrault_2003_post_stishovite_unit_cell.dat | 14 + .../Andrault_2003_stishovite_unit_cell.dat | 28 + .../data/Fischer_2018_stishovite.dat | 138 ++++ .../Ito_1974_stishovite_room_pressure.dat | 14 + .../data/Wang_2012_stv.dat | 62 ++ .../data/Zhang_2021_stishovite_angles.dat | 26 + .../Zhang_2021_stishovite_elastic_tensor.dat | 31 + .../data/Zhang_2021_stishovite_unit_cell.dat | 39 + .../anisotropic_stishovite/stishovite_data.py | 271 +++++++ .../stishovite_fit_eos.py | 547 ++++++++++++++ .../stishovite_model.py | 681 ++++++++++++++++++ .../stishovite_model_plots.py | 534 ++++++++++++++ .../stishovite_parameters.py | 92 +++ 14 files changed, 2513 insertions(+) create mode 100644 contrib/anisotropic_stishovite/Readme.md create mode 100644 contrib/anisotropic_stishovite/data/Andrault_2003_post_stishovite_unit_cell.dat create mode 100644 contrib/anisotropic_stishovite/data/Andrault_2003_stishovite_unit_cell.dat create mode 100644 contrib/anisotropic_stishovite/data/Fischer_2018_stishovite.dat create mode 100644 contrib/anisotropic_stishovite/data/Ito_1974_stishovite_room_pressure.dat create mode 100644 contrib/anisotropic_stishovite/data/Wang_2012_stv.dat create mode 100644 contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_angles.dat create mode 100644 contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_elastic_tensor.dat create mode 100644 contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_unit_cell.dat create mode 100644 contrib/anisotropic_stishovite/stishovite_data.py create mode 100644 contrib/anisotropic_stishovite/stishovite_fit_eos.py create mode 100644 contrib/anisotropic_stishovite/stishovite_model.py create mode 100644 contrib/anisotropic_stishovite/stishovite_model_plots.py create mode 100644 contrib/anisotropic_stishovite/stishovite_parameters.py diff --git a/contrib/anisotropic_stishovite/Readme.md b/contrib/anisotropic_stishovite/Readme.md new file mode 100644 index 00000000..dc6de24e --- /dev/null +++ b/contrib/anisotropic_stishovite/Readme.md @@ -0,0 +1,36 @@ +## 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 +## R. Myhill + +The python scripts contained in this directory accompany the paper +submitted by Myhill (2024). +There are comments in each of the scripts, so it is hoped that by reading +these in tandem with the original paper, the reader will get a feel +for how to use the code for their own purposes. + +The scripts contained in this directory are as follows: + +stishovite_data.py +------------------ +This script packages all of the data in the data directory into convenient +dictionaries. + +stishovite_parameters.py +------------------------ +This script provides the final fitted parameters for the model, +along with the bounds used during fitting. + +stishovite_model.py +------------------- +This file makes the burnman solution objects from the parameters. +The function get_models() also modifies the Zhang elasticity data +based on the model parameters, as described in the paper. + +stishovite_fit_eos.py +--------------------- +This script allows the user to refine the model parameters +based on the experimental data. + +stishovite_model_plots.py +------------------------- +This script creates all of the plots presented in the paper. diff --git a/contrib/anisotropic_stishovite/data/Andrault_2003_post_stishovite_unit_cell.dat b/contrib/anisotropic_stishovite/data/Andrault_2003_post_stishovite_unit_cell.dat new file mode 100644 index 00000000..5e910f8e --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Andrault_2003_post_stishovite_unit_cell.dat @@ -0,0 +1,14 @@ +# P P_err a a_err b b_err c c_err V V_err +66.31 0.19893 3.8971 0.0007 4.0008 0.0008 2.5658 0.0002 40.005 0.013 +69.58 0.20874 3.8829 0.0007 3.9955 0.0007 2.5608 0.0003 39.729 0.013 +76.42 0.22926 3.8467 0.0011 4.0114 0.0012 2.5509 0.0004 39.363 0.021 +81.72 0.24516 3.8325 0.0009 3.9848 0.001 2.5438 0.0002 38.848 0.016 +89.39 0.26817 3.817 0.001 3.9787 0.0011 2.5339 0.0002 38.481 0.017 +94.64 0.28392 3.8014 0.001 3.975 0.0011 2.5279 0.0002 38.198 0.017 +99.41 0.29823 3.7882 0.0009 3.9709 0.001 2.5221 0.0002 37.939 0.016 +103.7 0.3111 3.7773 0.001 3.9664 0.0011 2.5169 0.0003 37.709 0.017 +106.9 0.3207 3.7675 0.0006 3.9576 0.0007 2.513 0.0002 37.469 0.011 +114.6 0.3438 3.7505 0.0007 3.9533 0.0008 2.5065 0.0003 37.164 0.013 +120.4 0.3612 3.7364 0.0007 3.9482 0.0008 2.4998 0.0002 36.876 0.013 +124.1 0.3723 3.727 0.0008 3.9449 0.0009 2.496 0.0003 36.699 0.013 +128 0.384 3.7201 0.0008 3.9422 0.0009 2.4913 0.0003 36.535 0.014 diff --git a/contrib/anisotropic_stishovite/data/Andrault_2003_stishovite_unit_cell.dat b/contrib/anisotropic_stishovite/data/Andrault_2003_stishovite_unit_cell.dat new file mode 100644 index 00000000..614af997 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Andrault_2003_stishovite_unit_cell.dat @@ -0,0 +1,28 @@ +# P P_err a a_err c c_err V V_err +0.0001 0.00001 4.17755 0.00016 2.66518 0.00034 46.5126 0.0061 +1.168 0.006 4.17119 0.00027 2.66356 0.0003 46.3429 0.0053 +2.299 0.008 4.16504 0.00011 2.6618 0.00024 46.1756 0.0043 +3.137 0.008 4.16051 0.00014 2.66062 0.00029 46.055 0.0051 +4.252 0.008 4.15488 0.00013 2.65868 0.00026 45.8969 0.0045 +5.037 0.011 4.15084 0.00015 2.65767 0.0003 45.7902 0.0053 +5.851 0.009 4.14669 0.00012 2.65612 0.00022 45.6721 0.0038 +6.613 0.008 4.14323 0.00014 2.6547 0.00027 45.5715 0.005 +7.504 0.009 4.13888 0.00012 2.6534 0.00022 45.4536 0.0041 +8.264 0.011 4.13555 0.00015 2.65226 0.0003 45.3609 0.0056 +9.635 0.012 4.12962 0.00011 2.64976 0.00023 45.1885 0.0042 +11.69 0.03507 4.1217 0.0001 2.6458 0.0001 44.947 0.002 +17.67 0.05301 4.0962 0.0001 2.6381 0.0001 44.264 0.002 +22.38 0.06714 4.0781 0.0001 2.6321 0.0002 43.776 0.003 +29.38 0.08814 4.0552 0.0003 2.6192 0.0005 43.073 0.009 +37.71 0.11313 4.0287 0.0003 2.6048 0.0004 42.278 0.008 +46.03 0.13809 4.0035 0.0004 2.5919 0.001 41.544 0.017 +52.73 0.15819 3.9859 0.0004 2.5806 0.0004 40.999 0.009 +26.32 0.07896 4.0576 0.0002 2.6217 0.0003 43.164 0.006 +30.98 0.09294 4.044 0.0002 2.6155 0.0003 42.772 0.005 +34.21 0.10263 4.0308 0.0001 2.6101 0.0002 42.407 0.003 +38.45 0.11535 4.0207 0.0002 2.6037 0.0002 42.093 0.004 +43.37 0.13011 4.003 0.0002 2.5966 0.0002 41.61 0.004 +47.49 0.14247 3.9907 0.0003 2.5921 0.0004 41.28 0.007 +50.47 0.15141 3.9836 0.0002 2.5837 0.0004 41 0.008 +55.91 0.16773 3.9755 0.0003 2.5767 0.0008 40.724 0.013 +57.56 0.17268 3.9725 0.0003 2.5788 0.0009 40.695 0.015 diff --git a/contrib/anisotropic_stishovite/data/Fischer_2018_stishovite.dat b/contrib/anisotropic_stishovite/data/Fischer_2018_stishovite.dat new file mode 100644 index 00000000..1d971040 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Fischer_2018_stishovite.dat @@ -0,0 +1,138 @@ +# American Mineralogist: May 2018 Deposit AM-18-56267 +# Fischer et al.: Equations of state and phase boundary for SiO2 +# Table S2: Pressures (P), temperatures (T), SiO2 phase identification, and lattice parameters (a, b, c) from X-ray diffraction experiments. +# +# File phase line PKBr err PPt err Tsample err TKBr err aKBr err aPt err aSiO2 err bSiO2 err cSiO2 err +R133_1021 stishovite 13-ID-D 22.0 1.4 18.8 0.8 1672 109 1360 355 3.492 0.003 3.886 0.000 4.121 0.003 4.121 0.003 2.648 0.005 +R133_1026 stishovite 13-ID-D 21.8 1.6 18.0 0.9 1947 113 1572 426 3.500 0.003 3.898 0.001 4.124 0.001 4.124 0.001 2.648 0.002 +R133_1027 stishovite 13-ID-D 22.1 1.5 18.1 0.9 1850 114 1498 401 3.494 0.002 3.895 0.001 4.125 0.002 4.125 0.002 2.648 0.002 +R133_1028 stishovite 13-ID-D 21.9 1.4 19.5 1.5 1757 111 1425 377 3.495 0.002 3.886 0.005 4.120 0.001 4.120 0.001 2.647 0.001 +R133_1030 stishovite 13-ID-D 21.5 1.3 18.9 0.9 1603 109 1306 337 3.498 0.003 3.884 0.002 4.116 0.001 4.116 0.001 2.646 0.001 +R133_1031 stishovite 13-ID-D 21.3 1.2 18.8 0.9 1516 108 1239 315 3.498 0.002 3.881 0.001 4.114 0.001 4.114 0.001 2.646 0.001 +R133_1033 stishovite 13-ID-D 21.2 1.1 18.9 0.9 1404 110 1152 286 3.498 0.002 3.878 0.001 4.111 0.001 4.111 0.001 2.644 0.001 +R133_1034 stishovite 13-ID-D 21.3 1.0 19.2 0.9 1298 107 1071 259 3.495 0.001 3.874 0.001 4.110 0.001 4.110 0.001 2.643 0.001 +R133_1035 stishovite 13-ID-D 21.3 0.9 19.0 1.0 1226 107 1015 240 3.493 0.002 3.872 0.002 4.109 0.001 4.109 0.001 2.643 0.001 +R164_160 stishovite 16-ID-B 34.2 1.4 34.7 1.0 1614 121 1315 340 3.374 0.003 3.828 0.001 4.068 0.001 4.068 0.001 2.627 0.001 +R164_162 stishovite 16-ID-B 34.2 1.2 33.8 0.8 1478 107 1210 305 3.373 0.002 3.827 0.000 4.064 0.001 4.064 0.001 2.629 0.001 +R133_1055 stishovite 13-ID-D 36.3 1.5 0.0 0.0 1895 113 1532 412 3.362 0.002 0.000 0.000 4.061 0.001 4.061 0.001 2.629 0.001 +R133_1060 stishovite 13-ID-D 36.0 1.6 0.0 0.0 2009 113 1620 442 3.364 0.002 0.000 0.000 4.066 0.001 4.066 0.001 2.628 0.002 +R133_1061 stishovite 13-ID-D 36.2 1.7 0.0 0.0 2086 120 1680 462 3.364 0.002 0.000 0.000 4.064 0.001 4.064 0.001 2.628 0.001 +R133_1063 stishovite 13-ID-D 35.6 1.5 0.0 0.0 1784 112 1447 384 3.365 0.002 0.000 0.000 4.062 0.001 4.062 0.001 2.625 0.001 +R133_1064 stishovite 13-ID-D 35.6 1.4 0.0 0.0 1681 109 1366 357 3.364 0.002 0.000 0.000 4.060 0.001 4.060 0.001 2.626 0.001 +R133_1066 stishovite 13-ID-D 35.7 1.2 0.0 0.0 1539 108 1257 321 3.362 0.001 0.000 0.000 4.059 0.001 4.059 0.001 2.626 0.001 +R133_1068 stishovite 13-ID-D 35.7 1.1 0.0 0.0 1412 106 1159 288 3.360 0.001 0.000 0.000 4.056 0.001 4.056 0.001 2.624 0.001 +R133_1071 stishovite 13-ID-D 35.9 0.9 0.0 0.0 1258 105 1039 248 3.358 0.001 0.000 0.000 4.052 0.001 4.052 0.001 2.623 0.001 +R130_090 stishovite 13-ID-D 40.2 2.0 0.0 0.0 2501 125 2001 569 3.350 0.002 0.000 0.000 4.056 0.005 4.056 0.005 2.620 0.007 +R130_091 stishovite 13-ID-D 40.1 2.7 0.0 0.0 3001 228 2387 697 3.344 0.004 0.000 0.000 4.056 0.006 4.056 0.006 2.621 0.007 +R130_092 stishovite 13-ID-D 39.7 2.1 0.0 0.0 2325 120 1865 523 3.339 0.003 0.000 0.000 4.053 0.008 4.053 0.008 2.617 0.008 +R130_093 stishovite 13-ID-D 40.4 2.1 0.0 0.0 2401 136 1924 543 3.337 0.003 0.000 0.000 4.055 0.002 4.055 0.002 2.622 0.002 +R130_094 stishovite 13-ID-D 40.3 2.3 0.0 0.0 2645 253 2112 606 3.339 0.003 0.000 0.000 4.056 0.006 4.056 0.006 2.625 0.006 +R130_099 stishovite 13-ID-D 40.2 2.0 42.1 2.4 2271 129 1823 509 3.336 0.003 3.820 0.007 4.051 0.002 4.051 0.002 2.617 0.003 +R130_100 stishovite 13-ID-D 40.1 2.0 43.3 1.1 2226 120 1788 498 3.337 0.003 3.816 0.002 4.051 0.003 4.051 0.003 2.617 0.003 +R130_101 stishovite 13-ID-D 40.2 1.9 41.3 1.6 2081 119 1676 460 3.335 0.003 3.818 0.004 4.050 0.003 4.050 0.003 2.616 0.003 +R130_102 stishovite 13-ID-D 40.0 1.8 40.9 1.4 1913 115 1546 417 3.334 0.003 3.815 0.003 4.049 0.003 4.049 0.003 2.615 0.004 +R130_103 stishovite 13-ID-D 40.2 1.7 40.1 1.3 1847 114 1495 400 3.333 0.003 3.817 0.003 4.050 0.002 4.050 0.002 2.614 0.003 +R130_104 stishovite 13-ID-D 40.0 1.7 40.6 1.1 1783 111 1446 384 3.333 0.003 3.814 0.002 4.049 0.002 4.049 0.002 2.615 0.003 +R130_105 stishovite 13-ID-D 39.6 1.5 39.5 1.3 1663 109 1353 353 3.335 0.003 3.814 0.003 4.048 0.002 4.048 0.002 2.615 0.003 +R130_106 stishovite 13-ID-D 39.6 1.4 38.5 1.5 1575 109 1285 330 3.334 0.003 3.815 0.004 4.048 0.002 4.048 0.002 2.614 0.003 +R130_107 stishovite 13-ID-D 39.7 1.4 39.5 0.7 1467 108 1201 302 3.334 0.003 3.812 0.002 4.047 0.003 4.047 0.003 2.614 0.004 +R130_108 stishovite 13-ID-D 39.7 1.4 37.8 1.1 1420 108 1165 290 3.332 0.003 3.813 0.002 4.044 0.003 4.044 0.003 2.609 0.004 +R130_190 stishovite 13-ID-D 47.2 1.7 0.0 0.0 1955 119 1579 428 3.289 0.002 0.0 0.0 4.027 0.002 4.027 0.002 2.613 0.003 +R130_193 stishovite 13-ID-D 48.1 2.0 0.0 0.0 2187 124 1758 488 3.285 0.003 0.0 0.0 4.030 0.001 4.030 0.001 2.613 0.002 +R130_195 stishovite 13-ID-D 48.3 2.1 49.3 1.9 2355 130 1888 531 3.286 0.003 3.801 0.004 4.031 0.001 4.031 0.001 2.615 0.001 +R130_196 stishovite 13-ID-D 48.3 2.2 50.1 1.6 2518 137 2014 573 3.287 0.002 3.802 0.003 4.033 0.002 4.033 0.002 2.616 0.003 +R130_197 stishovite 13-ID-D 48.4 2.4 49.9 1.6 2703 130 2157 621 3.288 0.003 3.807 0.004 4.035 0.001 4.035 0.001 2.618 0.002 +R130_198 stishovite 13-ID-D 48.4 2.5 48.8 1.1 2833 141 2257 654 3.289 0.003 3.813 0.001 4.036 0.002 4.036 0.002 2.620 0.003 +R130_202 stishovite 13-ID-D 48.1 2.1 49.9 1.2 2280 132 1830 512 3.286 0.003 3.798 0.002 4.028 0.001 4.028 0.001 2.616 0.002 +R130_205 stishovite 13-ID-D 47.9 1.8 47.8 2.1 2077 133 1673 459 3.285 0.002 3.799 0.005 4.029 0.001 4.029 0.001 2.611 0.002 +R130_207 stishovite 13-ID-D 47.7 1.6 48.4 2.0 1868 120 1512 406 3.285 0.002 3.793 0.005 4.025 0.002 4.025 0.002 2.605 0.002 +R130_208 stishovite 13-ID-D 47.7 1.5 47.8 2.1 1766 115 1432 379 3.284 0.002 3.792 0.005 4.024 0.002 4.024 0.002 2.606 0.002 +R130_211 stishovite 13-ID-D 47.8 1.4 48.0 2.1 1698 110 1380 362 3.283 0.002 3.790 0.005 4.026 0.003 4.026 0.003 2.606 0.004 +R130_212 stishovite 13-ID-D 47.6 1.3 46.4 1.4 1610 108 1312 339 3.284 0.002 3.793 0.003 4.025 0.003 4.025 0.003 2.605 0.003 +R130_213 stishovite 13-ID-D 47.6 1.2 46.3 1.3 1477 107 1209 305 3.282 0.002 3.790 0.003 4.022 0.002 4.022 0.002 2.602 0.003 +R130_214 stishovite 13-ID-D 47.6 1.1 46.3 1.8 1395 106 1146 284 3.282 0.001 3.788 0.005 4.021 0.002 4.021 0.002 2.604 0.002 +R130_215 stishovite 13-ID-D 47.6 1.1 45.8 1.7 1295 105 1068 258 3.281 0.002 3.788 0.004 4.020 0.002 4.020 0.002 2.605 0.003 +R130_216 stishovite 13-ID-D 47.6 1.1 45.0 1.2 1224 104 1013 239 3.281 0.002 3.788 0.002 4.021 0.002 4.021 0.002 2.607 0.002 +R130_217 stishovite 13-ID-D 47.3 0.9 45.6 2.8 1134 104 943 216 3.282 0.001 3.785 0.007 4.019 0.002 4.019 0.002 2.605 0.003 +R130_235 stishovite 13-ID-D 56.5 1.8 0.0 0.0 1928 112 1557 421 3.238 0.002 0.0 0.0 4.006 0.007 4.006 0.007 2.607 0.011 +R130_237 stishovite 13-ID-D 56.5 2.0 0.0 0.0 2076 116 1672 459 3.239 0.003 0.0 0.0 4.005 0.007 4.005 0.007 2.614 0.011 +R130_240 stishovite 13-ID-D 57.1 2.1 0.0 0.0 2281 118 1830 512 3.238 0.002 0.0 0.0 4.008 0.007 4.008 0.007 2.615 0.013 +R130_244 stishovite 13-ID-D 57.6 2.3 0.0 0.0 2659 124 2123 609 3.238 0.002 0.0 0.0 4.011 0.004 4.011 0.004 2.611 0.004 +R130_245 stishovite 13-ID-D 57.6 2.2 0.0 0.0 2551 123 2039 581 3.237 0.002 0.0 0.0 4.007 0.003 4.007 0.003 2.605 0.003 +R130_246 stishovite 13-ID-D 57.5 2.1 0.0 0.0 2385 130 1911 539 3.236 0.002 0.0 0.0 4.007 0.005 4.007 0.005 2.605 0.005 +R130_249 stishovite 13-ID-D 57.4 1.9 0.0 0.0 2180 119 1752 486 3.235 0.002 0.0 0.0 4.003 0.005 4.003 0.005 2.603 0.006 +R130_250 stishovite 13-ID-D 57.3 1.7 0.0 0.0 1987 120 1603 436 3.235 0.001 0.0 0.0 4.000 0.005 4.000 0.005 2.606 0.006 +R130_251 stishovite 13-ID-D 57.3 1.5 0.0 0.0 1796 116 1456 387 3.233 0.001 0.0 0.0 3.995 0.003 3.995 0.003 2.610 0.004 +R130_253 stishovite 13-ID-D 57.2 1.4 0.0 0.0 1717 109 1394 366 3.233 0.001 0.0 0.0 3.998 0.007 3.998 0.007 2.601 0.007 +R130_254 stishovite 13-ID-D 56.6 2.0 0.0 0.0 1658 109 1349 351 3.236 0.004 0.0 0.0 3.992 0.003 3.992 0.003 2.612 0.004 +R130_255 stishovite 13-ID-D 57.2 1.3 0.0 0.0 1564 108 1276 327 3.232 0.001 0.0 0.0 3.989 0.002 3.989 0.002 2.610 0.003 +R130_256 stishovite 13-ID-D 57.2 1.2 0.0 0.0 1502 108 1228 311 3.232 0.001 0.0 0.0 3.992 0.003 3.992 0.003 2.607 0.004 +R130_257 stishovite 13-ID-D 57.0 1.1 0.0 0.0 1374 106 1129 278 3.232 0.001 0.0 0.0 3.992 0.003 3.992 0.003 2.607 0.003 +R130_258 stishovite 13-ID-D 57.0 1.0 0.0 0.0 1276 105 1054 253 3.231 0.001 0.0 0.0 3.993 0.004 3.993 0.004 2.607 0.005 +R130_260 stishovite 13-ID-D 57.1 0.9 0.0 0.0 1160 104 964 223 3.230 0.001 0.0 0.0 3.995 0.001 3.995 0.001 2.591 0.001 +R130_267 stishovite 13-ID-D 57.3 2.5 0.0 0.0 2876 139 2290 665 3.241 0.002 0.0 0.0 4.013 0.002 4.013 0.002 2.610 0.003 +R130_268 stishovite 13-ID-D 57.4 2.3 0.0 0.0 2701 132 2155 620 3.239 0.002 0.0 0.0 4.011 0.003 4.011 0.003 2.618 0.004 +R130_269 stishovite 13-ID-D 58.1 2.7 0.0 0.0 3056 143 2430 712 3.238 0.002 0.0 0.0 4.020 0.003 4.020 0.003 2.613 0.004 +R130_270 stishovite 13-ID-D 58.1 2.7 0.0 0.0 3213 135 2551 752 3.238 0.002 0.0 0.0 4.021 0.004 4.021 0.004 2.610 0.005 +R130_271 stishovite 13-ID-D 58.1 2.9 0.0 0.0 3274 159 2599 768 3.239 0.002 0.0 0.0 4.023 0.005 4.023 0.005 2.612 0.006 +R130_273 stishovite 13-ID-D 57.6 2.7 0.0 0.0 3112 146 2473 726 3.241 0.001 0.0 0.0 4.018 0.007 4.018 0.007 2.613 0.008 +R130_297 stishovite 13-ID-D 61.2 1.7 59.2 0.0 1907 112 1541 415 3.215 0.002 3.765 -1.000 3.990 0.008 3.990 0.008 2.569 0.011 +R130_298 stishovite 13-ID-D 61.0 1.9 60.2 0.0 2023 113 1631 445 3.217 0.002 3.765 -1.000 3.989 0.003 3.989 0.003 2.584 0.004 +R130_300 stishovite 13-ID-D 61.7 2.0 60.7 1.0 2095 115 1687 464 3.214 0.002 3.765 0.001 3.989 0.004 3.989 0.004 2.583 0.006 +R130_302 stishovite 13-ID-D 61.6 2.0 59.4 1.2 2209 116 1775 493 3.215 0.002 3.770 0.002 3.991 0.003 3.991 0.003 2.591 0.005 +R130_303 stishovite 13-ID-D 61.4 2.1 60.3 1.0 2356 120 1889 531 3.217 0.002 3.771 0.001 3.994 0.002 3.994 0.002 2.594 0.003 +R130_307 stishovite 13-ID-D 62.1 2.4 60.9 1.3 2575 130 2058 588 3.215 0.002 3.774 0.002 3.996 0.002 3.996 0.002 2.593 0.004 +R130_309 stishovite 13-ID-D 62.8 2.5 62.0 1.2 2785 127 2221 642 3.214 0.002 3.775 0.002 3.998 0.003 3.998 0.003 2.592 0.004 +R130_312 stishovite 13-ID-D 62.7 2.8 0.0 0.0 3130 134 2487 731 3.216 0.002 0.0 0.0 3.996 0.003 3.996 0.003 2.592 0.004 +R130_313 stishovite 13-ID-D 62.9 3.0 0.0 0.0 3228 136 2563 756 3.216 0.003 0.0 0.0 3.999 0.003 3.999 0.003 2.590 0.003 +R130_323 stishovite 13-ID-D 61.6 1.9 59.4 1.7 2160 134 1737 481 3.215 0.002 3.769 0.003 3.989 0.001 3.989 0.001 2.589 0.003 +R130_325 stishovite 13-ID-D 61.2 1.7 60.7 0.9 2021 118 1629 445 3.216 0.001 3.763 0.000 3.987 0.004 3.987 0.004 2.592 0.008 +R130_326 stishovite 13-ID-D 61.2 1.5 59.9 1.2 1860 152 1505 403 3.215 0.001 3.762 0.000 3.981 0.006 3.981 0.006 2.595 0.008 +R130_329 stishovite 13-ID-D 61.5 1.6 59.7 0.0 1821 122 1475 393 3.213 0.002 3.762 -1.000 3.985 0.003 3.985 0.003 2.580 0.003 +R130_331 stishovite 13-ID-D 61.4 1.6 58.2 1.0 1633 109 1330 345 3.213 0.002 3.762 0.001 3.984 0.004 3.984 0.004 2.579 0.006 +R130_333 stishovite 13-ID-D 61.0 1.2 58.5 1.0 1540 108 1258 321 3.214 0.001 3.759 0.001 3.982 0.005 3.982 0.005 2.581 0.006 +R130_334 stishovite 13-ID-D 60.8 1.1 57.8 0.8 1405 106 1153 286 3.214 0.001 3.758 0.000 3.985 0.005 3.985 0.005 2.581 0.010 +R130_335 stishovite 13-ID-D 60.9 1.0 59.0 1.7 1300 105 1072 259 3.213 0.001 3.753 0.004 3.983 0.007 3.983 0.007 2.573 0.011 +R133_1107 stishovite 13-ID-D 65.8 1.8 0.0 0.0 2258 117 1813 506 3.197 0.001 0.0 0.0 3.979 0.003 3.979 0.003 2.584 0.003 +R133_1108 stishovite 13-ID-D 65.8 1.9 0.0 0.0 2324 119 1864 523 3.197 0.001 0.0 0.0 3.979 0.003 3.979 0.003 2.586 0.004 +R133_1111 stishovite 13-ID-D 65.7 1.8 0.0 0.0 2174 117 1748 484 3.197 0.001 0.0 0.0 3.977 0.003 3.977 0.003 2.586 0.004 +R133_1113 stishovite 13-ID-D 65.7 1.7 0.0 0.0 2052 115 1654 453 3.196 0.001 0.0 0.0 3.975 0.003 3.975 0.003 2.582 0.004 +R133_1115 stishovite 13-ID-D 65.8 1.6 0.0 0.0 1958 113 1581 429 3.195 0.001 0.0 0.0 3.977 0.001 3.977 0.001 2.581 0.002 +R133_1117 stishovite 13-ID-D 65.8 1.4 0.0 0.0 1820 111 1474 393 3.194 0.001 0.0 0.0 3.974 0.002 3.974 0.002 2.580 0.002 +R133_1120 stishovite 13-ID-D 65.9 1.3 0.0 0.0 1702 111 1383 363 3.193 0.000 0.0 0.0 3.972 0.002 3.972 0.002 2.581 0.003 +R133_1122 stishovite 13-ID-D 66.1 1.2 0.0 0.0 1570 110 1281 329 3.192 0.000 0.0 0.0 3.972 0.002 3.972 0.002 2.578 0.003 +R133_1123 stishovite 13-ID-D 66.0 1.1 0.0 0.0 1480 108 1211 305 3.192 0.000 0.0 0.0 3.969 0.002 3.969 0.002 2.578 0.003 +R133_1124 stishovite 13-ID-D 66.0 1.0 0.0 0.0 1388 107 1140 282 3.191 0.000 0.0 0.0 3.969 0.002 3.969 0.002 2.579 0.002 +R164_198 stishovite 16-ID-B 74.1 2.5 0.0 0.0 2613 138 2088 598 3.165 0.002 0.0 0.0 3.963 0.003 3.963 0.003 2.569 0.004 +R164_199 stishovite 16-ID-B 74.4 2.1 0.0 0.0 2449 179 1961 555 3.163 0.001 0.0 0.0 3.963 0.001 3.963 0.001 2.575 0.002 +R164_206 stishovite 16-ID-B 74.7 2.0 0.0 0.0 2351 123 1885 530 3.161 0.001 0.0 0.0 3.961 0.002 3.961 0.002 2.576 0.002 +R164_208 stishovite 16-ID-B 74.5 1.8 73.2 1.0 2242 125 1800 502 3.161 0.001 3.738 0.001 3.960 0.002 3.960 0.002 2.577 0.002 +R164_213 stishovite 16-ID-B 74.5 1.7 73.6 1.0 2126 118 1711 472 3.161 0.001 3.735 0.001 3.959 0.002 3.959 0.002 2.575 0.002 +R164_214 stishovite 16-ID-B 74.5 1.7 72.8 1.1 2055 117 1656 454 3.160 0.001 3.735 0.001 3.957 0.001 3.957 0.001 2.575 0.001 +R130_339 CaCl2-type 13-ID-D 60.6 0.8 0.0 0.0 1054 104 882 196 3.213 0.001 0.0 0.0 3.986 0.021 3.961 0.017 2.580 0.010 +R133_1125 CaCl2-type 13-ID-D 65.9 1.0 0.0 0.0 1313 105 1082 262 3.191 0.000 0.0 0.0 3.975 0.004 3.958 0.004 2.574 0.003 +R133_1127 CaCl2-type 13-ID-D 66.0 0.8 0.0 0.0 1204 106 998 234 3.190 0.000 0.0 0.0 3.978 0.004 3.958 0.004 2.573 0.003 +R133_1128 CaCl2-type 13-ID-D 66.0 0.8 0.0 0.0 1133 105 943 216 3.190 0.000 0.0 0.0 3.981 0.005 3.955 0.006 2.574 0.006 +R164_216 CaCl2-type 16-ID-B 74.6 1.6 72.7 1.0 1949 113 1574 426 3.159 0.001 3.735 0.001 3.974 0.003 3.945 0.003 2.573 0.002 +R164_219 CaCl2-type 16-ID-B 74.7 1.5 72.4 1.0 1856 114 1502 402 3.159 0.001 3.733 0.001 3.975 0.004 3.945 0.004 2.574 0.002 +R164_223 CaCl2-type 16-ID-B 74.6 1.4 71.6 1.0 1733 111 1407 371 3.158 0.001 3.732 0.000 3.968 0.004 3.948 0.004 2.572 0.002 +R164_226 CaCl2-type 16-ID-B 74.5 1.3 71.1 1.0 1604 110 1307 337 3.158 0.001 3.731 0.000 3.966 0.002 3.945 0.002 2.573 0.001 +R164_229 CaCl2-type 16-ID-B 74.2 1.1 70.3 1.0 1488 110 1218 308 3.159 0.001 3.731 0.001 3.962 0.004 3.940 0.005 2.575 0.003 +R164_237 CaCl2-type 16-ID-B 74.3 1.1 70.8 0.9 1367 106 1124 276 3.157 0.001 3.727 0.001 3.969 0.002 3.937 0.002 2.572 0.001 +R164_251 CaCl2-type 16-ID-B 74.5 1.0 71.0 1.0 1298 105 1071 259 3.156 0.001 3.726 0.001 3.969 0.004 3.930 0.004 2.573 0.002 +R164_259 CaCl2-type 16-ID-B 84.6 1.6 84.0 0.9 1255 104 1037 247 3.121 0.002 3.698 0.001 3.944 0.007 3.898 0.007 2.549 0.004 +R164_261 CaCl2-type 16-ID-B 84.4 1.7 84.3 0.9 1314 105 1083 263 3.121 0.002 3.698 0.001 3.934 0.009 3.906 0.009 2.550 0.005 +R164_262 CaCl2-type 16-ID-B 84.6 1.9 85.1 0.9 1450 106 1188 298 3.121 0.003 3.699 0.001 3.943 0.007 3.910 0.007 2.548 0.005 +R164_264 CaCl2-type 16-ID-B 84.9 1.8 85.7 0.9 1569 108 1280 328 3.121 0.002 3.700 0.001 3.929 0.007 3.912 0.007 2.550 0.004 +R164_265 CaCl2-type 16-ID-B 85.2 1.8 86.1 1.0 1676 116 1363 356 3.120 0.002 3.701 0.000 3.935 0.006 3.908 0.006 2.551 0.004 +R164_268 CaCl2-type 16-ID-B 85.5 1.9 87.1 1.0 1771 122 1436 380 3.120 0.002 3.700 0.000 3.950 0.009 3.905 0.009 2.548 0.005 +R164_269 CaCl2-type 16-ID-B 85.6 1.9 87.6 1.0 1862 117 1506 404 3.120 0.002 3.701 0.000 3.953 0.008 3.901 0.007 2.550 0.005 +R164_270 CaCl2-type 16-ID-B 86.2 2.0 88.9 1.0 2051 124 1653 453 3.119 0.002 3.701 0.001 3.952 0.008 3.906 0.007 2.550 0.005 +R164_271 CaCl2-type 16-ID-B 86.1 2.2 88.6 1.1 2136 118 1719 475 3.120 0.002 3.703 0.001 3.954 0.008 3.905 0.007 2.550 0.005 +R164_274 CaCl2-type 16-ID-B 87.1 2.1 88.9 1.0 2387 119 1912 539 3.117 0.001 3.704 0.001 3.964 0.007 3.900 0.006 2.555 0.003 +R164_276 CaCl2-type 16-ID-B 86.7 2.1 89.5 1.1 2257 118 1812 506 3.118 0.002 3.705 0.001 3.960 0.005 3.907 0.005 2.554 0.003 +R164_278 CaCl2-type 16-ID-B 88.5 2.2 89.3 2.9 2551 126 2040 582 3.114 0.001 3.707 0.001 3.959 0.005 3.916 0.005 2.553 0.003 +R164_279 CaCl2-type 16-ID-B 88.6 2.4 90.0 1.0 2835 237 2259 655 3.115 0.001 3.708 0.000 3.957 0.007 3.913 0.007 2.556 0.004 +R164_281 CaCl2-type 16-ID-B 89.0 2.3 89.7 1.2 2735 127 2181 629 3.113 0.001 3.711 0.001 3.961 0.006 3.910 0.006 2.556 0.003 +R164_282 CaCl2-type 16-ID-B 88.8 2.2 0.0 0.0 2645 225 2112 606 3.113 0.001 0.0 0.0 3.952 0.007 3.913 0.007 2.557 0.004 +# -1 implies only one diffraction peak was clearly resolved, so no uncertainty can be assigned diff --git a/contrib/anisotropic_stishovite/data/Ito_1974_stishovite_room_pressure.dat b/contrib/anisotropic_stishovite/data/Ito_1974_stishovite_room_pressure.dat new file mode 100644 index 00000000..47454137 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Ito_1974_stishovite_room_pressure.dat @@ -0,0 +1,14 @@ +# T(K) a a_err b b_err c c_err V V_err +291.15 4.1789 0.00208945 4.1789 0.00208945 2.6664 0.00208945 46.56389037 0.049147819 +413.15 4.1834 0.0020917 4.1834 0.0020917 2.666 0.0020917 46.6572276 0.049279652 +473.15 4.1852 0.0020926 4.1852 0.0020926 2.6667 0.0020926 46.70964797 0.049339583 +523.15 4.1864 0.0020932 4.1864 0.0020932 2.6674 0.0020932 46.74870559 0.0493815 +573.15 4.1881 0.00209405 4.1881 0.00209405 2.6672 0.00209405 46.78317239 0.049431031 +773.15 4.195 0.0020975 4.195 0.0020975 2.6699 0.0020975 46.98496695 0.049661643 +291.15 4.1795 0.00208975 4.1795 0.00208975 2.6649 0.00208975 46.55106014 0.049153415 + +663.15 4.1922 0.0020961 4.1922 0.0020961 2.6686 0.0020961 46.89941969 0.049566275 +723.15 4.1937 0.00209685 4.1937 0.00209685 2.6698 0.00209685 46.95409215 0.04962154 +773.15 4.1951 0.00209755 4.1951 0.00209755 2.67 0.00209755 46.98896691 0.049665497 +823.15 4.1967 0.00209835 4.1967 0.00209835 2.6713 0.00209835 47.04771265 0.049724693 +873.15 4.1987 0.00209935 4.1987 0.00209935 2.6713 0.00209935 47.09256592 0.049785202 diff --git a/contrib/anisotropic_stishovite/data/Wang_2012_stv.dat b/contrib/anisotropic_stishovite/data/Wang_2012_stv.dat new file mode 100644 index 00000000..37ea7d9a --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Wang_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 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/Zhang_2021_stishovite_angles.dat b/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_angles.dat new file mode 100644 index 00000000..943f8267 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_angles.dat @@ -0,0 +1,26 @@ +# P Perr V Verr D Derr thetavar thetavarerr angle angleerr +0 0.00001 7.351 0.002 0.01278 0.00003 27.186 0.007 0 0.001 +2.8 0.1 7.319 0.004 0.01196 0.00007 27.75 0.02 0 0.001 +7.8 0.1 7.224 0.008 0.00923 0.0001 27.35 0.03 0 0.001 +13 0.2 7.128 0.008 0.00745 0.00008 26.7 0.03 0 0.001 +16 0.1 7.075 0.008 0.00685 0.00008 27.52 0.03 0 0.001 +16.9 0.3 7.057 0.009 0.00711 0.00011 27.16 0.03 0 0.001 +19.7 0.1 6.988 0.008 0.00537 0.00006 25.96 0.03 0 0.001 +21.4 0.1 6.994 0.009 0.00582 0.00007 27.92 0.03 0 0.001 +26.8 0.2 6.895 0.008 0.00662 0.00008 28.91 0.03 0 0.001 +28.5 0.2 6.828 0.009 0.00412 0.00005 27.16 0.03 0 0.001 +33.8 0.2 6.738 0.005 0.00325 0.00003 27.32 0.02 0 0.001 +40.4 0.2 6.659 0.005 0.00276 0.00002 28.01 0.02 0 0.001 +48.7 0.2 6.542 0.008 0.00075 0.00001 26.66 0.03 0 0.001 +49.8 0.2 6.522 0.01 0.00053 0.00001 25.85 0.04 0 0.2 +52.4 0.2 6.53 0.013 0.00036 0.00004 26.85 0.05 2.1 0.2 +54.2 0.2 6.511 0.013 0.00225 0.00001 25.52 0.05 2.6 0.2 +55.6 0.2 6.473 0.008 0.00076 0.00001 26.71 0.03 2.5 0.1 +58.6 0.3 6.466 0.013 0.00204 0.00004 24.52 0.05 3.3 0.2 +62 0.3 6.435 0.016 0.0027 0.00007 23.63 0.06 3.8 0.2 +64.4 0.3 6.418 0.018 0.00274 0.00011 24.92 0.07 4 0.2 +65.8 0.3 6.405 0.014 0.00221 0.00005 25.31 0.05 4.2 0.2 +68 0.3 6.386 0.013 0.00366 0.00007 23.79 0.05 4.6 0.2 +71 0.3 6.356 0.009 0.00304 0.00004 24.26 0.03 4.9 0.1 +73.8 0.3 6.341 0.013 0.00396 0.00008 23.68 0.05 5.1 0.2 +75.3 0.3 6.3 0.008 0.00259 0.00003 24.29 0.03 5.4 0.1 diff --git a/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_elastic_tensor.dat b/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_elastic_tensor.dat new file mode 100644 index 00000000..3f31e827 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_elastic_tensor.dat @@ -0,0 +1,31 @@ +# P  err C11  err C12  err C13  err C22  err C23  err C33  err C44  err C55  err C66  err ρ  err +# (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (GPa)  (g/cm3)  (g/cm3)  +# Platelet P1  +0.0001 0.00001 454.2 2.4 196.1 3.1 193 3.3 454.2 2.4 193 3.3 760.2 6.9 261.6 1.6 261.6 1.6 319.7 2 4.286 0.001 +26.4 0.8 558.2 4 371.6 5 233.5 2.8 558.2 4 233.5 2.8 847.5 6.2 302.6 1.9 302.6 1.9 396.2 3.2 4.598 0.001 +# Platelet P2  +8.2 0.3 491.3 1.9 248.7 2.3 205.7 6.8 491.3 1.9 205.7 6.8 787.8 33.6 275.8 2.5 275.8 2.5 344.8 1.6 4.392 0.001 +10.9 0.1 502.2 0.5 266.2 0.7 210.1 1.9 502.2 0.5 210.1 1.9 800.8 38 280 0.4 280 0.4 352.8 0.4 4.425 0.001 +14.6 0.1 516.9 1.3 290.8 1.6 215.3 4.4 516.9 1.3 215.3 4.4 812.6 42 286 1.4 286 1.4 363.6 0.9 4.469 0.001 +23.8 0.1 548.9 1.4 353.8 1.7 228.6 4.9 548.9 1.4 228.6 4.9 843.2 45 299.4 1.4 299.4 1.4 389.4 1.2 4.571 0.001 +32.6 0.1 574.2 1.2 418 1.4 240.4 4 574.2 1.2 240.4 4 867.1 49 311.8 1.4 311.8 1.4 413.2 1 4.662 0.001 +38.7 0.4 587.9 1.8 464.5 2.1 248.1 6.4 587.9 1.8 248.1 6.4 881.8 52 319.7 1.9 319.7 1.9 429.1 1.5 4.721 0.001 +42.3 0.2 594.4 1.8 494.1 2.1 252.2 5.9 594.4 1.8 252.2 5.9 889.8 55 324.5 1.6 324.5 1.6 438.8 1.7 4.756 0.001 +50.1 0.2 603.5 5.3 559.9 6 262.1 16.3 603.5 5.3 262.1 16.3 905 58 334.1 5.3 334.1 5.3 458.5 4.8 4.826 0.001 +# Platelet P3  +17.9 0.1 531.2 3.5 310.7 7.8 220.8 3.2 531.2 3.5 220.8 3.2 823.4 3 290.6 1.5 290.6 1.5 371.4 3.6 4.507 0.001 +20.8 0.1 539.2 4.8 331.9 10.3 224.1 4.1 539.2 4.8 224.1 4.1 832.6 4 294.5 2.1 294.5 2.1 382.2 5 4.538 0.001 +28.5 0.5 564.2 4.8 386.7 10.9 230.2 3.6 564.2 4.8 230.2 3.6 855.7 3.4 306.6 1.7 306.6 1.7 403.7 4.5 4.62 0.001 +35.3 0.3 583.9 4.5 440.1 9.9 240.3 3.9 583.9 4.5 240.3 3.9 873.2 3.4 315.9 1.8 315.9 1.8 416.4 4.4 4.688 0.001 +41 0.3 592.7 3.2 482.6 8.2 251.4 2.2 592.7 3.2 251.4 2.2 886.6 2 322.8 1 322.8 1 434.4 5.1 4.743 0.001 +42.3 0.1 594.7 3.9 492.7 7.7 250.2 3.3 594.7 3.9 250.2 3.3 890.9 2.9 324.3 1.5 324.3 1.5 439.7 4.1 4.755 0.001 +44.1 0.1 597.4 5.8 505.5 12 253.5 5.2 597.4 5.8 253.5 5.2 894.4 4.5 327.2 2.3 327.2 2.3 444.6 6.1 4.772 0.001 +45.8 0.2 599.7 3.8 521.9 8.7 256.3 3.3 599.7 3.8 256.3 3.3 898.2 2.9 329.1 1.6 329.1 1.6 447.4 4.4 4.788 0.001 +47.9 0.1 600.6 3.7 540 6.7 258.6 3 600.6 3.7 258.6 3 903.7 2.9 331.4 1.4 331.4 1.4 454.9 3.5 4.807 0.001 +49.1 0.2 603 3.8 549.4 7.6 260.4 3.2 603 3.8 260.4 3.2 905.7 3 333.1 1.4 333.1 1.4 456.6 3.9 4.818 0.001 +52.1 0.2 605.5 4 578.6 8.3 264.2 3.5 605.5 4 264.2 3.5 912.2 3.1 336.4 1.5 336.4 1.5 462.4 4.1 4.844 0.001 +55 0.5 605.9 5 605.5 9.4 267.6 4.2 605.9 5 267.6 4.2 918.5 3.8 339.8 1.8 339.8 1.8 470.2 5 4.869 0.002 +57.2 0.1 593.8 36.8 594.9 5.8 237.4 9.6 673.8 3.3 298.4 2.4 920.1 1.6 346 1.2 338.8 1.5 475.9 2.7 4.899 0.002 +60.1 0.2 609 92 585.6 20.6 226.1 29.8 723.5 10.2 310.8 7.4 922.2 4.8 351.8 3.5 340.4 4.3 483.3 8.8 4.929 0.002 +64.2 0.2 648.5 96.8 576 17.4 219.3 25.8 779.6 8.9 322 6.4 926.3 4.2 357.9 3.4 343 4.1 494.5 6.3 4.971 0.002 +69.7 0.3 695.2 99.8 572.7 19.9 215.5 26.4 838.7 9.4 331.4 6.3 932.3 4.1 365.5 3.7 347 4.6 507.2 6.5 5.026 0.002 diff --git a/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_unit_cell.dat b/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_unit_cell.dat new file mode 100644 index 00000000..26174d7e --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_unit_cell.dat @@ -0,0 +1,39 @@ +# P (GPa)  err a (A)  err b (A)  err c (A)  err V (Å)  err ρ (g/cm3)  err +0.0001  0.000001 4.1771 0.001 4.1771 0.001 2.669 0.0008 46.57 0.03 4.286 0.003 +4.3 0.1 4.1539 0.0006 4.1539 0.0006 2.6611 0.0008 45.92 0.01 4.347 0.001 +5.9 0.1 4.1466 0.0006 4.1466 0.0006 2.6596 0.0008 45.73 0.01 4.365 0.001 +7.4 0.1 4.1387 0.0006 4.1387 0.0006 2.6561 0.0008 45.5 0.01 4.388 0.001 +9 0.1 4.1319 0.0006 4.1319 0.0006 2.6541 0.0008 45.31 0.01 4.405 0.001 +10.4 0.1 4.1255 0.0006 4.1255 0.0006 2.6518 0.0008 45.13 0.01 4.423 0.001 +12.8 0.1 4.1156 0.0006 4.1156 0.0006 2.6479 0.0008 44.85 0.01 4.451 0.001 +14.8 0.1 4.1074 0.0006 4.1074 0.0006 2.6446 0.0008 44.62 0.01 4.474 0.001 +16.1 0.1 4.1023 0.0006 4.1023 0.0006 2.6424 0.0008 44.47 0.01 4.489 0.001 +17.1 0.1 4.0992 0.0006 4.0992 0.0006 2.6416 0.0008 44.39 0.01 4.497 0.001 +18.3 0.1 4.0949 0.0006 4.0949 0.0006 2.6395 0.0008 44.26 0.01 4.51 0.001 +19.6 0.1 4.09 0.0006 4.09 0.0006 2.6375 0.0008 44.12 0.01 4.524 0.001 +19.8 0.1 4.0889 0.0006 4.0889 0.0006 2.6371 0.0008 44.09 0.01 4.527 0.001 +21.9 0.1 4.082 0.0006 4.082 0.0006 2.6322 0.0008 43.86 0.01 4.551 0.001 +23 0.2 4.0785 0.0006 4.0785 0.0006 2.631 0.0008 43.76 0.01 4.561 0.001 +25.3 0.1 4.0705 0.0006 4.0705 0.0006 2.6278 0.0008 43.54 0.01 4.585 0.001 +26.8 0.1 4.0654 0.0006 4.0654 0.0006 2.6261 0.0008 43.4 0.01 4.599 0.001 +27.6 0.2 4.0612 0.0006 4.0612 0.0006 2.6238 0.0008 43.27 0.01 4.613 0.001 +29.2 0.2 4.0566 0.0006 4.0566 0.0006 2.6221 0.0008 43.15 0.01 4.626 0.001 +30.8 0.2 4.0528 0.0006 4.0528 0.0006 2.6203 0.0008 43.04 0.01 4.638 0.001 +32.1 0.2 4.0485 0.0006 4.0485 0.0006 2.6174 0.0008 42.9 0.01 4.653 0.001 +33.8 0.2 4.0422 0.0005 4.0422 0.0005 2.6166 0.0008 42.75 0.01 4.669 0.001 +35.2 0.2 4.0393 0.0005 4.0393 0.0005 2.6133 0.0007 42.64 0.01 4.682 0.001 +37.3 0.2 4.0328 0.0006 4.0328 0.0006 2.6105 0.0008 42.45 0.01 4.702 0.001 +38.6 0.2 4.0293 0.0006 4.0293 0.0006 2.6096 0.0008 42.37 0.01 4.711 0.001 +40.3 0.2 4.024 0.0006 4.024 0.0006 2.6062 0.0008 42.2 0.01 4.73 0.001 +42.4 0.2 4.016 0.0006 4.016 0.0006 2.6039 0.0008 42 0.01 4.753 0.001 +43.9 0.2 4.011 0.0006 4.011 0.0006 2.6015 0.0007 41.85 0.01 4.769 0.001 +45.2 0.2 4.008 0.0006 4.008 0.0006 2.5992 0.0008 41.75 0.01 4.781 0.001 +47.6 0.2 4.001 0.0006 4.001 0.0006 2.5965 0.0008 41.56 0.01 4.802 0.001 +49.4 0.2 3.997 0.0006 3.997 0.0006 2.5929 0.0008 41.42 0.01 4.819 0.001 +50.8 0.2 3.994 0.0006 3.994 0.0006 2.5895 0.0007 41.31 0.01 4.832 0.001 +52.7 0.2 3.96 0.0112 4.0131 0.0009 2.5886 0.0008 41.14 0.11 4.852 0.013 +54.8 0.2 3.9487 0.0015 4.0106 0.0006 2.5844 0.0007 40.93 0.02 4.877 0.002 +56.9 0.2 3.9385 0.0015 4.0082 0.0006 2.58 0.0007 40.73 0.02 4.901 0.002 +58.4 0.2 3.9319 0.0015 4.0065 0.0006 2.5783 0.0007 40.62 0.02 4.915 0.002 +59.9 0.2 3.9242 0.0015 4.0048 0.0006 2.5762 0.0007 40.49 0.02 4.931 0.002 +61.5 0.2 3.9168 0.0015 4.0029 0.0006 2.5741 0.0007 40.36 0.02 4.946 0.002 diff --git a/contrib/anisotropic_stishovite/stishovite_data.py b/contrib/anisotropic_stishovite/stishovite_data.py new file mode 100644 index 00000000..c0303742 --- /dev/null +++ b/contrib/anisotropic_stishovite/stishovite_data.py @@ -0,0 +1,271 @@ +import numpy as np +from burnman.utils.unitcell import molar_volume_from_unit_cell_volume + + +def add_to_data(data, publication, compilation): + if publication not in data: + data[publication] = {} + for name, header, dataset in compilation: + data[publication][name] = { + header[i]: dataset[:, i] for i in range(len(dataset[0])) + } + + +def get_data(): + data = {} + + # Zhang data + CT_data = np.loadtxt("data/Zhang_2021_stishovite_elastic_tensor.dat") + unit_cell_data = np.loadtxt("data/Zhang_2021_stishovite_unit_cell.dat") + CT_header = [ + "P", + "P_err", + "C11", + "C11_err", + "C12", + "C12_err", + "C13", + "C13_err", + "C22", + "C22_err", + "C23", + "C23_err", + "C33", + "C33_err", + "C44", + "C44_err", + "C55", + "C55_err", + "C66", + "C66_err", + "rho", + "rho_err", + ] + unit_cell_header = [ + "P", + "P_err", + "a", + "a_err", + "b", + "b_err", + "c", + "c_err", + "V", + "V_err", + "rho", + "rho_err", + ] + + isstv = unit_cell_data[:, 2] == unit_cell_data[:, 4] + + compilation = [ + ["CT", CT_header, CT_data], + ["stv", unit_cell_header, unit_cell_data[isstv]], + ["poststv", unit_cell_header, unit_cell_data[~isstv]], + ] + add_to_data(data, "Zhang_2021", compilation) + for d in [ + data["Zhang_2021"]["CT"], + data["Zhang_2021"]["stv"], + data["Zhang_2021"]["poststv"], + ]: + d["T"] = 0.0 * d["P"] + 298.15 + d["T_err"] = 0.0 * d["P"] + 5.0 + + # Andrault data + stv_data = np.loadtxt("data/Andrault_2003_stishovite_unit_cell.dat") + poststv_data = np.loadtxt("data/Andrault_2003_post_stishovite_unit_cell.dat") + stv_header = ["P", "P_err", "a", "a_err", "c", "c_err", "V", "V_err"] + poststv_header = [ + "P", + "P_err", + "a", + "a_err", + "b", + "b_err", + "c", + "c_err", + "V", + "V_err", + ] + compilation = [ + ["stv", stv_header, stv_data], + ["poststv", poststv_header, poststv_data], + ] + add_to_data(data, "Andrault_2003", compilation) + d = data["Andrault_2003"]["stv"] + d["b"] = d["a"] + d["b_err"] = d["a_err"] + for d in [data["Andrault_2003"]["stv"], data["Andrault_2003"]["poststv"]]: + d["T"] = 0.0 * d["P"] + 298.15 + d["T_err"] = 0.0 * d["P"] + 5.0 + + # Fischer data (swap a and b) + F_data = np.genfromtxt("data/Fischer_2018_stishovite.dat", dtype=str) + file, phase, beamline = F_data.T[:3] + mask = phase == "stishovite" + + F_data = np.genfromtxt("data/Fischer_2018_stishovite.dat")[:, 3:] + header = [ + "P", + "P_err", + "P_Pt", + "P_Pt_err", + "T", + "T_err", + "T_KBr", + "T_KBr_err", + "a_KBr", + "a_KBr_err", + "a_Pt", + "a_Pt_err", + "b", + "b_err", + "a", + "a_err", + "c", + "c_err", + ] + compilation = [["stv", header, F_data[mask]], ["poststv", header, F_data[~mask]]] + add_to_data(data, "Fischer_2018", compilation) + + for d in [data["Fischer_2018"]["stv"], data["Fischer_2018"]["poststv"]]: + d["V"] = d["a"] * d["b"] * d["c"] + d["V_err"] = d["V"] * np.sqrt( + ( + np.power(d["a_err"] / d["a"], 2.0) + + np.power(d["b_err"] / d["b"], 2.0) + + np.power(d["c_err"] / d["c"], 2.0) + ) + ) + + # Ito data (only stishovite) + I_data = np.loadtxt("data/Ito_1974_stishovite_room_pressure.dat") + header = ["T", "a", "a_err", "b", "b_err", "c", "c_err", "V", "V_err"] + compilation = [["stv", header, I_data]] + add_to_data(data, "Ito_1974", compilation) + + d = data["Ito_1974"]["stv"] + d["T_err"] = d["T"] * 0.0 + 5.0 + d["P"] = d["T"] * 0.0 + 0.0001 + d["P_err"] = d["T"] * 0.0 + 0.000001 + + # Convert to SI units + d = data["Zhang_2021"]["CT"] + d["P"] = d["P"] * 1.0e9 + d["P_err"] = d["P_err"] * 1.0e9 + + f_ln = np.cbrt(molar_volume_from_unit_cell_volume(1.0, 2.0)) + + for d in [ + data["Andrault_2003"]["stv"], + data["Andrault_2003"]["poststv"], + data["Ito_1974"]["stv"], + data["Fischer_2018"]["stv"], + data["Fischer_2018"]["poststv"], + data["Zhang_2021"]["stv"], + data["Zhang_2021"]["poststv"], + ]: + 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) + d["P"] = d["P"] * 1.0e9 + d["P_err"] = d["P_err"] * 1.0e9 + + for ln in ["a", "b", "c"]: + d[f"{ln}"] = d[f"{ln}"] * f_ln + d[f"{ln}_err"] = d[f"{ln}_err"] * f_ln + + # Molar mass (M) is apparently rounded by Zhang to get rho. + # As rho is a derived value, we correct rho here and + # return the derived volumes + M_Zhang = 0.0601 + M = 0.06008 + f = M / M_Zhang + data["Zhang_2021"]["CT"]["rho"] = 1.0e3 * f * data["Zhang_2021"]["CT"]["rho"] + data["Zhang_2021"]["CT"]["rho_err"] = ( + 1.0e3 * f * data["Zhang_2021"]["CT"]["rho_err"] + ) + data["Zhang_2021"]["CT"]["V"] = M / data["Zhang_2021"]["CT"]["rho"] + data["Zhang_2021"]["CT"]["V_err"] = M * ( + data["Zhang_2021"]["CT"]["rho_err"] + / np.power(data["Zhang_2021"]["CT"]["rho"], 2.0) + ) + for phase in ["stv", "poststv"]: + data["Zhang_2021"][phase]["rho"] = 1.0e3 * f * data["Zhang_2021"][phase]["rho"] + data["Zhang_2021"][phase]["rho_err"] = ( + 1.0e3 * f * data["Zhang_2021"][phase]["rho_err"] + ) + data["Zhang_2021"][phase]["V2"] = M / data["Zhang_2021"][phase]["rho"] + data["Zhang_2021"][phase]["V2_err"] = M * ( + data["Zhang_2021"][phase]["rho_err"] + / np.power(data["Zhang_2021"][phase]["rho"], 2.0) + ) + + return data + + +def common_data(): + # First, let's read in the PVT equation of state data + data = get_data() + + d = {} + d["PTV"] = {"stv": {}, "poststv": {}} + d["PTV_err"] = {"stv": {}, "poststv": {}} + d["abc"] = {"stv": {}, "poststv": {}} + d["abc_err"] = {"stv": {}, "poststv": {}} + + for pub in ["Ito_1974", "Andrault_2003", "Fischer_2018", "Zhang_2021"]: + for phase in ["stv", "poststv"]: + try: + Z_data = data[pub][phase] + d["PTV"][phase][pub] = np.array( + [Z_data["P"], Z_data["T"], Z_data["V"]] + ).T + d["PTV_err"][phase][pub] = np.array( + [Z_data["P_err"], Z_data["T_err"], Z_data["V_err"]] + ).T + except: + pass + + try: + Z_data = data[pub][phase] + d["abc"][phase][pub] = np.array( + [Z_data["a"], Z_data["b"], Z_data["c"]] + ).T + d["abc_err"][phase][pub] = np.array( + [Z_data["a_err"], Z_data["b_err"], Z_data["c_err"]] + ).T + except: + pass + return d + + +def get_stiffnesses(CN): + """ """ + n = len(CN["C11"]) + Cs = np.zeros((n, 6, 6)) + Cs_err = np.zeros((n, 6, 6)) + for ij in ["11", "12", "13", "22", "23", "33", "44", "55", "66"]: + i, j = (int(ij[0]) - 1, int(ij[1]) - 1) + Cs[:, i, j] = CN[f"C{ij}"] + Cs[:, j, i] = CN[f"C{ij}"] + Cs_err[:, i, j] = CN[f"C{ij}_err"] + Cs_err[:, j, i] = CN[f"C{ij}_err"] + K_NR = 1.0 / np.sum(np.linalg.inv(Cs)[:, :3, :3], axis=(1, 2)) + return Cs, Cs_err, K_NR + + +CN_data = get_data()["Zhang_2021"]["CT"] +CN_GPa, CN_err_GPa, KNR_GPa = get_stiffnesses(CN_data) +P_for_CN = CN_data["P"] +T_for_CN = CN_data["T"] +V_for_CN = CN_data["V"] +SN_invGPa = np.linalg.inv(CN_GPa) + +# approximate error in K_NR +KNR_err_GPa = np.sqrt(np.sum(np.power(CN_err_GPa[:, :3, :3], 2.0), axis=(1, 2))) / 9.0 + +beta_NR_invGPa = np.sum(SN_invGPa[:, :3, :3], axis=(1, 2)) +SNoverbetaNR_obs = np.einsum("ijk, i->ijk", SN_invGPa, 1.0 / beta_NR_invGPa) +betaNRoverSN_err_obs = np.ones_like(SNoverbetaNR_obs) * 0.2 +lnV_for_CN = np.log(V_for_CN) diff --git a/contrib/anisotropic_stishovite/stishovite_fit_eos.py b/contrib/anisotropic_stishovite/stishovite_fit_eos.py new file mode 100644 index 00000000..137f7a41 --- /dev/null +++ b/contrib/anisotropic_stishovite/stishovite_fit_eos.py @@ -0,0 +1,547 @@ +from __future__ import absolute_import +from __future__ import print_function + +import numpy as np +from scipy.optimize import minimize, differential_evolution + +from stishovite_data import ( + common_data, + P_for_CN, + V_for_CN, + T_for_CN, + KNR_GPa, + KNR_err_GPa, + CN_GPa, + CN_err_GPa, +) +from scipy.optimize import root_scalar +from stishovite_model import make_models, make_scalar_model, modify_Zhang_elasticity +from copy import deepcopy +from burnman import RelaxedSolution +import matplotlib.pyplot as plt +from stishovite_parameters import scalar_bounds, cell_bounds, elastic_bounds +from stishovite_parameters import scalar_args, cell_args, elastic_args +from stishovite_parameters import scalar_and_cell_args, all_args + +data = common_data() + +delta_ijk = np.einsum("ij, lk, l->ijk", np.eye(6), np.eye(6), np.ones(6)) + +min_misfit_scalar = [1.0e7] +min_misfit_cell = [1.0e7] +min_misfit_elastic = [1.0e7] +min_misfit_combined = [1.0e7] + +plot = False + + +def transition_pressure(solution, T): + def diff_mu(P_GPa): + solution.set_state(P_GPa * 1.0e9, T) + solution.set_composition([0.499, 0.501]) + mu = solution.partial_gibbs + return mu[1] - mu[0] + + sol = root_scalar(diff_mu, x0=80.0, x1=60.0) + assert sol.converged + return sol.root * 1.0e9 + + +# We want to fit the Zhang room temperature data (a, b, c) +# and also the symmetry breaking pressure from Fischer at 3000 K +# (approximately 88 GPa) + + +def misfit_scalar(args): + ( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + dDebye_0, + P_tr_GPa, + fP_Zhang, + fP_Andrault, + fP2_Zhang, + fP2_Andrault, + ) = args + + scalar_prms = make_scalar_model( + dVQ0, dKQ0, dKpQ0, dgrQ0, dqQ0, V0Q1overV0Q0, dDebye_0, P_tr_GPa + ) + stishovite_Q0, stishovite_Q1, scalar_stv, ESV_interactions = scalar_prms + + if plot: + fig = plt.figure() + ax = [fig.add_subplot(1, 3, i) for i in range(1, 4)] + + relaxed_stv = RelaxedSolution( + scalar_stv, + relaxation_vectors=np.array([[1.0, -1.0]]), + unrelaxed_vectors=np.array([[0.5, 0.5]]), + ) + relaxed_stv.set_composition([1.0]) + + # Fit the HT transition pressure from Fischer + 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)) + + # Fit volumes for stishovite from Zhang, Andrault + for phase, publication in [ + ("stv", "Zhang_2021"), + ("poststv", "Zhang_2021"), + ("stv", "Andrault_2003"), + ("poststv", "Andrault_2003"), + ]: + if publication == "Zhang_2021": + fP = fP_Zhang + fP2 = fP2_Zhang + elif publication == "Andrault_2003": + fP = fP_Andrault + fP2 = fP2_Andrault + else: + fP = 1.0 + fP2 = 0.0 + + PTV = deepcopy(data["PTV"][phase][publication]) + PTV_err = deepcopy(data["PTV_err"][phase][publication]) + PTV_err[:, 0] = np.max([1.0e8 + 0.01 * PTV[:, 0], PTV_err[:, 0]]) + P_model = relaxed_stv.evaluate_with_volumes(["pressure"], PTV[:, 2], PTV[:, 1])[ + 0 + ] + P_actual = PTV[:, 0] * (fP + fP2 * PTV[:, 0]) + chisqr += np.sum(np.power((P_model - P_actual) / PTV_err[:, 0], 2.0)) + + if plot: + sort = np.argsort(P_model) + ax[0].plot(P_model[sort], PTV[:, 2][sort]) + ax[0].scatter(P_actual[sort], PTV[:, 2][sort]) + + # Fit relative volumes from Ito + PTV = deepcopy(data["PTV"]["stv"]["Ito_1974"]) + PTV_err = deepcopy(data["PTV_err"]["stv"]["Ito_1974"]) + scalar_stv.set_composition([0.5, 0.5]) + 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)) + + if plot: + sort = np.argsort(PTV[:, 1]) + ax[1].plot(PTV[:, 1][sort], V_model[sort]) + ax[1].scatter(PTV[:, 1][sort], PTV[:, 2][sort]) + + # Fit isentropic bulk moduli from Zhang + 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)) + + if plot: + sort = np.argsort(V_for_CN) + ax[2].plot(V_for_CN[sort], KNR_obs[sort] / 1.0e9) + ax[2].scatter(V_for_CN[sort], KNR_GPa[sort]) + plt.show() + + if chisqr < min_misfit_scalar[0]: + min_misfit_scalar[0] = chisqr + print(repr(args)) + print(chisqr) + return chisqr + + +def misfit_cell(args, scalar_args): + ( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + dDebye_0, + P_tr_GPa, + fP_Zhang, + fP_Andrault, + fP2_Zhang, + fP2_Andrault, + ) = scalar_args + (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = args + ( + a11, + a22, + a33, + a44, + a55, + a66, + b11, + b22, + b33, + b44, + b55, + b66, + c44, + c55, + c66, + b112, + b222, + b332, + ) = np.zeros(18) + + models = make_models( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + a0Q1, + b0Q1, + dDebye_0, + P_tr_GPa, + PsiI_33_a, + PsiI_33_b, + PsiI_33_c, + PsiI_33_b2, + PsiI_33_c2, + a11, + a22, + a33, + a44, + a55, + a66, + b11, + b22, + b33, + b44, + b55, + b66, + c44, + c55, + c66, + b112, + b222, + b332, + ) + _, _, _, stishovite_relaxed = models + stishovite_relaxed.set_composition([1.0]) + + chisqr = 0.0 + + if plot: + fig = plt.figure() + ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)] + + # fit a, b and c for stishovite from Zhang and Andrault + for phase, publication in [ + ("stv", "Zhang_2021"), + ("poststv", "Zhang_2021"), + # ("stv", "Andrault_2003"), + ("poststv", "Andrault_2003"), + ]: + if publication == "Zhang_2021": + fP = fP_Zhang + fP2 = fP2_Zhang + elif publication == "Andrault_2003": + fP = fP_Andrault + fP2 = fP2_Andrault + else: + fP = 1.0 + fP2 = 0.0 + + PTV = data["PTV"][phase][publication] + abc_obs = data["abc"][phase][publication] + abc_err_obs = np.max( + [ + data["abc"][phase][publication] * 0.0005, + data["abc_err"][phase][publication], + ], + axis=0, + ) + + P_actual = PTV[:, 0] * (fP + fP2 * PTV[:, 0]) + cell_parameters_model = stishovite_relaxed.evaluate( + ["cell_parameters"], P_actual, PTV[:, 1] + )[0] + abc_model = cell_parameters_model[:, :3] + + chisqr += np.sum(np.power((abc_obs - abc_model) / abc_err_obs, 2.0)) + if plot: + for i in range(3): + j = 1 + if i < 2: + j = 0 + sort = np.argsort(P_actual) + ax[j].errorbar( + P_actual[sort], + abc_obs[sort, i], + yerr=abc_err_obs[sort, i], + ls="none", + ) + + if plot: + pressures = np.linspace(1.0e5, 120.0e9, 101) + temperatures = pressures * 0.0 + 298.15 + cell_parameters_model = stishovite_relaxed.evaluate( + ["cell_parameters"], pressures, temperatures + )[0] + abc_model = cell_parameters_model[:, :3] + for i in range(3): + j = 1 + if i < 2: + j = 0 + ax[j].plot(pressures, abc_model[:, i]) + plt.show() + + if chisqr < min_misfit_cell[0]: + min_misfit_cell[0] = chisqr + print(repr(args)) + print(chisqr) + return chisqr + + +def misfit_elastic(args, scalar_args, cell_args): + ( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + dDebye_0, + P_tr_GPa, + fP_Zhang, + fP_Andrault, + fP2_Zhang, + fP2_Andrault, + ) = 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) = ( + args + ) + + b55i = b44i + b22i = b11i + b222i = b112i + c55 = c44 + + frel = -0.14 + b11 = b11i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b22 = b22i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b33 = b33i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b112 = b112i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b222 = b222i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b332 = b332i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b44 = b44i / (c44 * np.exp(c44 * frel) - 1.0) + b55 = b55i / (c55 * np.exp(c55 * frel) - 1.0) + b66 = b66i / (c66 * np.exp(c66 * frel) - 1.0) + + models = make_models( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + a0Q1, + b0Q1, + dDebye_0, + P_tr_GPa, + PsiI_33_a, + PsiI_33_b, + PsiI_33_c, + PsiI_33_b2, + PsiI_33_c2, + a11, + a22, + a33, + a44, + a55, + a66, + b11, + b22, + b33, + b44, + b55, + b66, + c44, + c55, + c66, + b112, + b222, + b332, + ) + _, _, _, stishovite_relaxed = models + stishovite_relaxed.set_composition([1.0]) + + chisqr = 0.0 + + if plot: + fig = plt.figure() + ax = [fig.add_subplot(1, 3, i) for i in range(1, 4)] + + # Fit compliance data + fP = fP_Zhang + fP2 = fP2_Zhang + P_actual = P_for_CN * (fP + fP2 * P_for_CN) + CN_model = stishovite_relaxed.evaluate( + ["isentropic_stiffness_tensor"], P_actual, T_for_CN + )[0] + sort = np.argsort(P_actual) + for i, j in [ + (0, 0), + (1, 1), + (2, 2), + (0, 1), + (0, 2), + (1, 2), + (3, 3), + (4, 4), + (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)) + + if plot: + axi = 2 + if i < 2 and j < 2: + axi = 0 + elif i < 3: + axi = 1 + + lc = ax[axi].plot( + P_actual[sort], CN_model[sort, i, j] / 1.0e9, label=f"{i+1}{j+1}" + ) + ax[axi].errorbar( + P_actual[sort], + CN_GPa[sort, i, j], + yerr=CN_err_GPa[sort, i, j], + ls="none", + c=lc[0].get_color(), + ) + + if plot: + for axi in range(3): + ax[axi].legend() + plt.show() + + if chisqr < min_misfit_elastic[0]: + min_misfit_elastic[0] = chisqr + print(repr(args)) + print(chisqr) + return chisqr + + +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) + + if chisqr < min_misfit_combined[0]: + min_misfit_combined[0] = chisqr + print(repr(args)) + print(chisqr) + return chisqr + + +def misfit_cell_and_elastic(args, scalar_args): + cell_args = args[:7] + elastic_args = args[7:] + chisqr = misfit_cell(cell_args, scalar_args) + modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) + chisqr += misfit_elastic(elastic_args, scalar_args, cell_args) + + if chisqr < min_misfit_combined[0]: + min_misfit_combined[0] = chisqr + print(repr(args)) + print(chisqr) + return chisqr + + +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) + modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) + chisqr += misfit_elastic(elastic_args, scalar_args, cell_args) + + 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(chisqr) + return chisqr + + +if __name__ == "__main__": + if False: + sol = minimize( + misfit_scalar, scalar_args, bounds=scalar_bounds, method="Nelder-Mead" + ) + + if False: + + def misfit_de(cell_args): + return misfit_cell(cell_args, scalar_args) + + sol = differential_evolution(func=misfit_de, bounds=cell_bounds, x0=cell_args) + + sol = minimize( + misfit_cell, + cell_args, + bounds=cell_bounds, + args=(scalar_args), + method="Nelder-Mead", + ) + + if False: + sol = minimize( + misfit_scalar_and_cell, + scalar_and_cell_args, + bounds=scalar_bounds + cell_bounds, + method="Nelder-Mead", + ) + + if True: + sol = minimize( + misfit_all, + all_args, + bounds=scalar_bounds + cell_bounds + elastic_bounds, + method="Nelder-Mead", + ) + + if False: + cell_and_elastic_args = np.concatenate((cell_args, elastic_args)) + modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) + sol = minimize( + misfit_cell_and_elastic, + cell_and_elastic_args, + bounds=cell_bounds + elastic_bounds, + args=(scalar_args), + method="Nelder-Mead", + ) + + if False: + sol = minimize( + misfit_elastic, + elastic_args, + bounds=elastic_bounds, + args=(scalar_args, cell_args), + method="Powell", + ) + + if False: + sol = minimize( + misfit_cell, + cell_args, + bounds=cell_bounds, + args=(scalar_args), + method="Nelder-Mead", + ) diff --git a/contrib/anisotropic_stishovite/stishovite_model.py b/contrib/anisotropic_stishovite/stishovite_model.py new file mode 100644 index 00000000..0a43c094 --- /dev/null +++ b/contrib/anisotropic_stishovite/stishovite_model.py @@ -0,0 +1,681 @@ +import burnman +from burnman.classes.solutionmodel import PolynomialSolution +from burnman import AnisotropicMineral, AnisotropicSolution +from burnman import RelaxedAnisotropicSolution +import numpy as np +from copy import copy, deepcopy +from stishovite_parameters import all_args +from stishovite_data import P_for_CN, T_for_CN, SN_invGPa, CN_GPa +from tabulate import tabulate + +stv_SLB = burnman.minerals.SLB_2022.st() +stv_SLB.property_modifiers = [] + + +def modify_Zhang_elasticity(scalar_args, cell_args, elastic_args): + ( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + dDebye_0, + P_tr_GPa, + fP_Zhang, + fP_Andrault, + fP2_Zhang, + fP2_Andrault, + ) = 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) = ( + elastic_args + ) + + b55i = b44i + b22i = b11i + b222i = b112i + c55 = c44 + + frel = -0.14 + b11 = b11i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b22 = b22i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b33 = b33i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b112 = b112i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b222 = b222i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b332 = b332i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b44 = b44i / (c44 * np.exp(c44 * frel) - 1.0) + b55 = b55i / (c55 * np.exp(c55 * frel) - 1.0) + b66 = b66i / (c66 * np.exp(c66 * frel) - 1.0) + + models = make_models( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + a0Q1, + b0Q1, + dDebye_0, + P_tr_GPa, + PsiI_33_a, + PsiI_33_b, + PsiI_33_c, + PsiI_33_b2, + PsiI_33_c2, + a11, + a22, + a33, + a44, + a55, + a66, + b11, + b22, + b33, + b44, + b55, + b66, + c44, + c55, + c66, + b112, + b222, + b332, + ) + _, _, _, stishovite_relaxed = models + stishovite_relaxed.set_composition([1.0]) + + # Fit compliance data + fP = fP_Zhang + fP2 = fP2_Zhang + P_actual = P_for_CN * (fP + fP2 * P_for_CN) + + beta_T_model = stishovite_relaxed.evaluate( + ["isentropic_compressibility_tensor"], P_actual, T_for_CN + )[0] + beta_ii_model = np.einsum("ijj->ij", beta_T_model) + SN = deepcopy(SN_invGPa) / 1.0e9 + beta_ii_obs = np.einsum("ijk->ij", SN[:, :3, :3]) + beta_RT_obs = np.einsum("ijk->i", SN[:, :3, :3]) + + if False: + # Modify in an "L" shape (C13, C23, C33, C32, C31) + f = beta_ii_model[:, 2] / beta_ii_obs[:, 2] + g_obs = 2.0 * np.sum(SN[:, :3, 2], axis=1) - SN[:, 2, 2] + SN[:, 2, :] = np.einsum("ij, i->ij", SN[:, 2, :], f) + SN[:, :, 2] = np.einsum("ij, i->ij", SN[:, :, 2], f) + SN[:, 2, 2] = SN[:, 2, 2] / f + + f2 = (beta_RT_obs - g_obs * f) / (beta_RT_obs - g_obs) + SN[:, :2, :2] = np.einsum("ijk, i->ijk", SN[:, :2, :2], f2) + elif False: + # Modify only the diagonals + for i in range(3): + SN[:, i, i] = SN[:, i, i] + beta_ii_model[:, i] - beta_ii_obs[:, i] + else: + # Modify only the off-diagonals + dbeta = beta_ii_model - beta_ii_obs + SN[:, 0, 1] = SN[:, 0, 1] + (dbeta[:, 0] + dbeta[:, 1] - dbeta[:, 2]) / 2.0 + SN[:, 0, 2] = SN[:, 0, 2] + (dbeta[:, 0] + dbeta[:, 2] - dbeta[:, 1]) / 2.0 + SN[:, 1, 2] = SN[:, 1, 2] + (dbeta[:, 1] + dbeta[:, 2] - dbeta[:, 0]) / 2.0 + SN[:, 1, 0] = SN[:, 0, 1] + SN[:, 2, 0] = SN[:, 0, 2] + SN[:, 2, 1] = SN[:, 1, 2] + + # modify in place + # SN is used to calculate, so overwriting repeatedly is ok + CN_GPa[:, :, :] = np.linalg.inv(SN) / 1.0e9 + + +def make_scalar_model(dVQ0, dKQ0, dKpQ0, dgrQ0, dqQ0, V0Q1overV0Q0, dDebye_0, P_tr_GPa): + # Make Q0 + stishovite_Q0 = deepcopy(stv_SLB) + for prm, dv in [ + ("V_0", dVQ0 * 1.0e-6), # cm^3 + ("K_0", dKQ0 * 1.0e9), # GPa + ("Kprime_0", dKpQ0), + ("grueneisen_0", dgrQ0), + ("q_0", dqQ0), + ]: + stishovite_Q0.params[prm] = stv_SLB.params[prm] + dv + + # Make Q1 + stishovite_Q1 = deepcopy(stv_SLB) + Q_scale_factor = 6.25e10 # user-chosen variable (scales Q) + V0Q1 = V0Q1overV0Q0 * stishovite_Q0.params["V_0"] + dV = V0Q1 - stishovite_Q0.params["V_0"] + dF = -dV * Q_scale_factor + for prm, dv in [ + ("V_0", dV), + ("K_0", 0.0), + ("Kprime_0", 0.0), + ("grueneisen_0", 0.0), + ("q_0", 0.0), + ("F_0", dF), + ("Debye_0", dDebye_0), + ]: + stishovite_Q1.params[prm] = stishovite_Q0.params[prm] + dv + + """ + # (1^6 - Q^6) = 12*(pa^5*pb + pb^5*pa) + 40*(pa^3*pb^3) + # (1^4 - Q^4) = 8*(pa^3*pb + pb^3*pa) + # (1^2 - Q^2) = 4*(pa*pb) + + # G_xs = (a-b)*(1 - Q^2) + b*(1 - Q^n) + + # At Q = 0 + # G_xs = a + # Let a = GQ0 - GQ1 + + # At Q = 1 + # G_xs = 0 + + # At Q = 0 + # d2G/dQdQ = -2*(a - b) - c*b*Q^(n-2) + # a = b + + # Thus at Q0, the symmetry breaking occurs when b = (GQ0-GQ1) + """ + + T_tr = 298.15 # user-chosen variable + P_tr = P_tr_GPa * 1.0e9 + stishovite_Q0.set_state(P_tr, T_tr) + stishovite_Q1.set_state(P_tr, T_tr) + b_xs = stishovite_Q0.gibbs - stishovite_Q1.gibbs + + # Q^4 + ESV_interactions = [ + [-4.0 * b_xs, 0.0, 0.0, 0, 1, 1, 1], + [8.0 * b_xs, 0.0, 0.0, 0, 1, 1, 3], + [8.0 * b_xs, 0.0, 0.0, 0, 3, 1, 1], + ] + + """ + # Q^6 + ESV_interactions = [[-4.*b_xs, 0., 0., 0, 1, 1, 1], + [12.*b_xs, 0., 0., 0, 5, 1, 1], + [40.*b_xs, 0., 0., 0, 3, 1, 3], + [12.*b_xs, 0., 0., 0, 1, 1, 5]] + """ + + class stv_poststv(burnman.Solution): + def __init__(self, molar_fractions=None): + self.name = "stishovite-post-stishovite" + self.solution_model = PolynomialSolution( + endmembers=[[stishovite_Q1, "[Si]O2"], [stishovite_Q1, "[Si]O2"]], + ESV_interactions=ESV_interactions, + interaction_endmembers=[stishovite_Q0, stishovite_Q1], + endmember_coefficients_and_interactions=[[4.0, -4.0, 0, 1, 1, 1]], + ) + self.b_xs = b_xs # for easy access later + self.endmember_coefficients_and_interactions = [[4.0, -4.0, 0, 1, 1, 1]] + burnman.Solution.__init__(self, molar_fractions=molar_fractions) + + return stishovite_Q0, stishovite_Q1, stv_poststv(), ESV_interactions + + +def make_models( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + a0Q1, + b0Q1, + dDebye_0, + P_tr_GPa, + PsiI_33_a, + PsiI_33_b, + PsiI_33_c, + PsiI_33_b2, + PsiI_33_c2, + a11, + a22, + a33, + a44, + a55, + a66, + b11, + b22, + b33, + b44, + b55, + b66, + c44, + c55, + c66, + b112, + b222, + b332, +): + + scalar_prms = make_scalar_model( + dVQ0, dKQ0, dKpQ0, dgrQ0, dqQ0, V0Q1overV0Q0, dDebye_0, P_tr_GPa + ) + stishovite_Q0, stishovite_Q1, scalar_stv, ESV_interactions = scalar_prms + + c0Q1 = stishovite_Q1.params["V_0"] / (a0Q1 * b0Q1) + + """ + ANISOTROPIC MODELS + """ + + def psi_func(f, Pth, params): + b1 = params["b1"] + b2 = params["b2"] + dPsidf = ( + params["a"] + + b1 * params["c1"] * (np.exp(params["c1"] * f) - 1.0) + + b2 * params["c2"] * (np.exp(params["c2"] * f) - 1.0) + ) + Psi = ( + 0.0 + + (params["a"] - b1 * params["c1"] - b2 * params["c2"]) * f + + b1 * (np.exp(params["c1"] * f) - 1.0) + + b2 * (np.exp(params["c2"] * f) - 1.0) + ) + dPsidPth = np.zeros((6, 6)) + return (Psi, dPsidf, dPsidPth) + + # Make the ordered model + # Note that c is not dependent on Q, so take the value of + # c from the standard state volume of the Q1 structure + # in normal Q0 stishovite (at room temperature): + # c_Q1(V0, T0) = c_Q0(V0(Q1), T0) + + # Note also that we assume that the Q-related strain is + # not a function of V or T (note not in agreement with + # Fischer data), so we can use the splitting + # for the Q1-structure at any V or T where it is stable + # to calculate the relative splitting of a and b at any + # other state where Q is known. + + cell_parameters = np.array([a0Q1, b0Q1, c0Q1, 90.0, 90.0, 90.0]) # m, degrees + anisotropic_parameters = { + "a": np.zeros((6, 6)), + "b1": np.zeros((6, 6)), + "c1": np.zeros((6, 6)), + "b2": np.zeros((6, 6)), + "c2": np.zeros((6, 6)), + } + + """ + beta_T/beta_RT = d(Psi*I)/df + F_ij = exp(Psi I) + + M=FM0 + + for orthotropic systems: + ln(M_ii) = ln(F_ii) + ln(M0_ii) + ln(M_ii) - ln(M0_ii) = ln(F_ii) = (Psi I)_ii + ln(a) - ln(a0) = (Psi I)_11 + ln(b) - ln(b0) = (Psi I)_22 + ln(c) - ln(c0) = (Psi I)_33 + ln(c/c0Q1) = aQ1 * fQ1 + bQ1 * (np.exp(cQ1 * fQ1) - 1.) + + note that bQ1 needs to be scaled relative to b fit using Q0, + because of the change in volume: + Let f(Q1) be the constant ln(V0Q1/V0Q0) + fQ1 is ln(VQ1/V0Q1) + ln(c) - ln(c0Q1) = a * (f - f(Q1)) + b * (np.exp(c * f) - np.exp(c * f(Q1))) + ln(c) - ln(c0Q1) = a * fQ1 + b * (np.exp(c)*(np.exp(f) - np.exp(f(Q1)))) + ln(c) - ln(c0Q1) = a * fQ1 + b*np.exp(f(Q1)) * (np.exp(c)*(np.exp(f - f(Q1)) - 1)) + ln(c) - ln(c0Q1) = a * fQ1 + b*np.exp(f(Q1)) * (np.exp(c * fQ1) - 1) + + Thus + aQ1 = a + bQ1 = b*np.exp(f(Q1)) = b*V0Q1/V0Q0 + cQ1 = c + + # lna_Q1 = 0.5*(lnV - lnc_Q1) - lna_minus_lnb/2. + # ln(a/a0)_Q1 = 0.5*(f - ln(c/c0)_Q1) + """ + + # the following give the a, b1, c1 + PsiI_33 = { + "a": PsiI_33_a, + "b1": PsiI_33_b, + "c1": PsiI_33_c, + "b2": PsiI_33_b2, + "c2": PsiI_33_c2, + } + # The following assumes that the compression of the a and + # b axes is the same + PsiI_11 = { + "a": 0.5 * (1.0 - PsiI_33_a), + "b1": -0.5 * PsiI_33_b, + "c1": PsiI_33_c, + "b2": -0.5 * PsiI_33_b2, + "c2": PsiI_33_c2, + } + PsiI_22 = { + "a": 0.5 * (1.0 - PsiI_33_a), + "b1": -0.5 * PsiI_33_b, + "c1": PsiI_33_c, + "b2": -0.5 * PsiI_33_b2, + "c2": PsiI_33_c2, + } + + anisotropic_parameters["a"][0, 0] = a11 + anisotropic_parameters["b1"][0, 0] = b11 + anisotropic_parameters["c1"][0, 0] = PsiI_33_c + anisotropic_parameters["b2"][0, 0] = b112 + anisotropic_parameters["c2"][0, 0] = PsiI_33_c2 + anisotropic_parameters["a"][1, 1] = a22 + anisotropic_parameters["b1"][1, 1] = b22 + anisotropic_parameters["c1"][1, 1] = PsiI_33_c + anisotropic_parameters["b2"][1, 1] = b222 + anisotropic_parameters["c2"][1, 1] = PsiI_33_c2 + anisotropic_parameters["a"][2, 2] = a33 + anisotropic_parameters["b1"][2, 2] = b33 + anisotropic_parameters["c1"][2, 2] = PsiI_33_c + anisotropic_parameters["b2"][2, 2] = b332 + anisotropic_parameters["c2"][2, 2] = PsiI_33_c2 + anisotropic_parameters["a"][3, 3] = a44 + anisotropic_parameters["b1"][3, 3] = b44 + anisotropic_parameters["c1"][3, 3] = c44 + anisotropic_parameters["a"][4, 4] = a55 + anisotropic_parameters["b1"][4, 4] = b55 + anisotropic_parameters["c1"][4, 4] = c55 + anisotropic_parameters["a"][5, 5] = a66 + anisotropic_parameters["b1"][5, 5] = b66 + anisotropic_parameters["c1"][5, 5] = c66 + + # Fill the rest + for p in ["a", "b1", "b2"]: + + anisotropic_parameters[p][0, 1] = 0.5 * ( + (PsiI_11[p] + PsiI_22[p] - PsiI_33[p]) + - ( + anisotropic_parameters[p][0, 0] + + anisotropic_parameters[p][1, 1] + - anisotropic_parameters[p][2, 2] + ) + ) + anisotropic_parameters[p][0, 2] = 0.5 * ( + (PsiI_11[p] + PsiI_33[p] - PsiI_22[p]) + - ( + anisotropic_parameters[p][0, 0] + + anisotropic_parameters[p][2, 2] + - anisotropic_parameters[p][1, 1] + ) + ) + anisotropic_parameters[p][1, 2] = 0.5 * ( + (PsiI_22[p] + PsiI_33[p] - PsiI_11[p]) + - ( + anisotropic_parameters[p][1, 1] + + anisotropic_parameters[p][2, 2] + - anisotropic_parameters[p][0, 0] + ) + ) + + anisotropic_parameters[p][1, 0] = anisotropic_parameters[p][0, 1] + anisotropic_parameters[p][2, 0] = anisotropic_parameters[p][0, 2] + anisotropic_parameters[p][2, 1] = anisotropic_parameters[p][1, 2] + + for i, j in [[0, 1], [0, 2], [1, 2]]: + anisotropic_parameters["c1"][i, j] = PsiI_33_c + anisotropic_parameters["c1"][j, i] = PsiI_33_c + anisotropic_parameters["c2"][i, j] = PsiI_33_c2 + anisotropic_parameters["c2"][j, i] = PsiI_33_c2 + + # Make the antiordered model + # Remember, Voigt form to standard notation is: + # 11->1, 22->2, 33->3, 23->4, 13->5, 12->6 + # So ordered to antiordered transition through + # a tetragonal phase involves the + # following change of rows and columns + # 1 <-> 2 (i.e. indices 0 and 1) + # 4 <-> 5 (i.e. indices 3 and 4) + + anti_cell_parameters = deepcopy(cell_parameters) + anti_cell_parameters[0] = cell_parameters[1] + anti_cell_parameters[1] = cell_parameters[0] + anti_anisotropic_parameters = deepcopy(anisotropic_parameters) + for p in ["a", "b1", "c1", "b2", "c2"]: + for i, j in [(0, 1), (3, 4)]: + l1 = copy(anti_anisotropic_parameters[p][j, :]) + l2 = copy(anti_anisotropic_parameters[p][i, :]) + anti_anisotropic_parameters[p][i, :] = l1 + anti_anisotropic_parameters[p][j, :] = l2 + l1 = copy(anti_anisotropic_parameters[p][:, j]) + l2 = copy(anti_anisotropic_parameters[p][:, i]) + anti_anisotropic_parameters[p][:, i] = l1 + anti_anisotropic_parameters[p][:, j] = l2 + + # print(np.sum(anti_anisotropic_parameters["a"][:3], axis=0)[:3]) + # print(np.sum(anti_anisotropic_parameters["b1"][:3], axis=0)[:3]) + # print(np.sum(anti_anisotropic_parameters["c1"][:3], axis=0)[:3]) + # exit() + anisotropic_stv_Q1 = AnisotropicMineral( + stishovite_Q1, + cell_parameters, + anisotropic_parameters, + psi_function=psi_func, + orthotropic=True, + ) + anisotropic_anti_stv_Q1 = AnisotropicMineral( + stishovite_Q1, + anti_cell_parameters, + anti_anisotropic_parameters, + psi_function=psi_func, + orthotropic=True, + ) + + # An ideal mixing model + def fn(lnV, Pth, X, params): + z = np.zeros((6, 6)) + f = np.zeros((3, 3, 2)) + return (z, z, z, f) + + prm = {} + + stishovite_anisotropic = AnisotropicSolution( + name="stishovite", + solution_model=PolynomialSolution( + endmembers=[ + [anisotropic_stv_Q1, "[Si]O2"], + [anisotropic_anti_stv_Q1, "[Si]O2"], + ], + ESV_interactions=ESV_interactions, + interaction_endmembers=[stishovite_Q0, stishovite_Q1], + endmember_coefficients_and_interactions=[[4.0, -4.0, 0, 1, 1, 1]], + ), + psi_excess_function=fn, + anisotropic_parameters=prm, + ) + + stishovite_relaxed = RelaxedAnisotropicSolution( + stishovite_anisotropic, + relaxation_vectors=np.array([[1.0, -1.0]]), + unrelaxed_vectors=np.array([[0.5, 0.5]]), + ) + + return stishovite_Q1, scalar_stv, stishovite_anisotropic, stishovite_relaxed + + +def get_models(): + + scalar_args = all_args[:12] + cell_args = all_args[12:19] + elastic_args = all_args[19:] + + modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) + + ( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + dDebye_0, + P_tr_GPa, + fP_Zhang, + fP_Andrault, + fP2_Zhang, + fP2_Andrault, + ) = 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) = ( + elastic_args + ) + + b55i = b44i + b22i = b11i + b222i = b112i + c55 = c44 + + frel = -0.14 + b11 = b11i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b22 = b22i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b33 = b33i / (PsiI_33_c * np.exp(PsiI_33_c * frel) - 1.0) + b112 = b112i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b222 = b222i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b332 = b332i / (PsiI_33_c2 * np.exp(PsiI_33_c2 * frel) - 1.0) + b44 = b44i / (c44 * np.exp(c44 * frel) - 1.0) + b55 = b55i / (c55 * np.exp(c55 * frel) - 1.0) + b66 = b66i / (c66 * np.exp(c66 * frel) - 1.0) + + models = make_models( + dVQ0, + dKQ0, + dKpQ0, + dgrQ0, + dqQ0, + V0Q1overV0Q0, + a0Q1, + b0Q1, + dDebye_0, + P_tr_GPa, + PsiI_33_a, + PsiI_33_b, + PsiI_33_c, + PsiI_33_b2, + PsiI_33_c2, + a11, + a22, + a33, + a44, + a55, + a66, + b11, + b22, + b33, + b44, + b55, + b66, + c44, + c55, + c66, + b112, + b222, + b332, + ) + + # The following assumes that the compression of the a and + # b axes is the same + PsiI_11_a = 0.5 * (1.0 - PsiI_33_a) + PsiI_11_b = -0.5 * PsiI_33_b + PsiI_11_b2 = -0.5 * PsiI_33_b2 + + sm = models[1] + Q0, Q1 = sm.solution_model.interaction_endmembers + properties = ["F_0", "V_0", "K_0", "Kprime_0", "Debye_0", "grueneisen_0", "q_0"] + pps = [ + "$\\mathcal{F}_0$", + "$V_0$", + "$K_0$", + "$K'_0$", + "Debye $T_0$", + "$\\gamma_0$", + "$q_0$", + ] + 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(["$a_0$", "-", a0Q1]) + table.append(["$b_0$", "-", b0Q1]) + + print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".6e")) + + print("b_xs", sm.b_xs) + + headers = [""] + headers.extend([f"$\\Psi_{{{i+1}}}$" for i in range(3)]) + headers.extend([f"$\\Psi_{{{i+1}{i+1}}}$" for i in range(6)]) + table = [headers] + table.extend( + [ + ["$a$", PsiI_11_a, "as above", PsiI_33_a, a11, a22, a33, a44, a55, a66], + [ + "$b_1$", + PsiI_11_b, + "as above", + PsiI_33_b, + b11, + "as above", + b33, + b44, + "as above", + b66, + ], + [ + "$c_1$", + PsiI_33_c, + "as above", + "as above", + "as above", + "as above", + "as above", + c44, + "as above", + c66, + ], + [ + "$b_2$", + PsiI_11_b2, + "as above", + PsiI_33_b2, + b112, + "as above", + b332, + "-", + "-", + "-", + ], + [ + "$c_2$", + PsiI_33_c2, + "as above", + "as above", + "as above", + "as above", + "as above", + "-", + "-", + "-", + ], + ] + ) + + table = [ + [f"{item:.6e}" if type(item) is np.float64 else item for item in row] + for row in table + ] + print( + tabulate( + list(map(list, zip(*table))), + headers="firstrow", + tablefmt="latex_raw", + floatfmt=".6e", + ) + ) + + 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 diff --git a/contrib/anisotropic_stishovite/stishovite_model_plots.py b/contrib/anisotropic_stishovite/stishovite_model_plots.py new file mode 100644 index 00000000..145c7f12 --- /dev/null +++ b/contrib/anisotropic_stishovite/stishovite_model_plots.py @@ -0,0 +1,534 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib +import matplotlib.cm as cm +from matplotlib.gridspec import GridSpec +from copy import deepcopy +from scipy.optimize import minimize +from stishovite_data import common_data +from stishovite_model import get_models +from burnman.tools.eos import check_anisotropic_eos_consistency +from stishovite_data import CN_GPa, P_for_CN, KNR_GPa, CN_err_GPa +from stishovite_fit_eos import transition_pressure +from burnman.tools.plot import plot_projected_elastic_properties +from burnman.utils.math import is_positive_definite +import burnman + +CN_GPa_orig = deepcopy(CN_GPa) + +models = get_models() +_, stishovite_scalar, stishovite_unrelaxed, stishovite_relaxed = models[0] +fP_Zhang, fP_Andrault, fP2_Zhang, fP2_Andrault = 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_relaxed = True +plot_Q_V_abc = True +plot_G = True +plot_lnabc = True + +if plot_isotropic_velocities: + pressures = np.linspace(9.0e9, 128e9, 501) + seismic_model = burnman.seismic.PREM() + depths = seismic_model.depth(pressures) + temperatures = burnman.geotherm.brown_shankland(depths) + + v = stishovite_relaxed.evaluate( + [ + "isentropic_shear_modulus_voigt", + "isentropic_shear_modulus_reuss", + "isentropic_shear_modulus_vrh", + "isentropic_bulk_modulus_voigt", + "isentropic_bulk_modulus_reuss", + "isentropic_bulk_modulus_vrh", + "density", + ], + pressures, + temperatures, + ) + + G_V, G_R, G_VRH, K_V, K_R, K_VRH, rho = v + + VS_V = np.sqrt(G_V / rho) + VS_R = np.sqrt(G_R / rho) + VS_VRH = np.sqrt(G_VRH / rho) + + Vphi_V = np.sqrt(K_V / rho) + Vphi_R = np.sqrt(K_R / rho) + Vphi_VRH = np.sqrt(K_VRH / rho) + + VP_V = np.sqrt((K_V + 4.0 / 3.0 * G_V) / rho) + VP_R = np.sqrt((K_R + 4.0 / 3.0 * G_R) / rho) + VP_VRH = np.sqrt((K_VRH + 4.0 / 3.0 * G_VRH) / rho) + + fig = plt.figure() + ax = [fig.add_subplot(1, 1, 1)] + + ax[0].fill_between( + depths / 1.0e3, VP_R / 1.0e3, VP_V / 1.0e3, alpha=0.3, color="blue" + ) + ax[0].plot(depths / 1.0e3, VP_VRH / 1.0e3, color="blue", label="$V_P$") + + ax[0].fill_between( + depths / 1.0e3, Vphi_R / 1.0e3, Vphi_V / 1.0e3, alpha=0.3, color="purple" + ) + ax[0].plot(depths / 1.0e3, Vphi_VRH / 1.0e3, color="purple", label="$V_{\\Phi}$") + + ax[0].fill_between( + depths / 1.0e3, VS_R / 1.0e3, VS_V / 1.0e3, alpha=0.3, color="red" + ) + ax[0].plot(depths / 1.0e3, VS_VRH / 1.0e3, color="red", label="$V_S$") + + ax[0].set_ylim( + 0.0, + ) + ax[0].set_xlim(np.min(depths) / 1.0e3, np.max(depths) / 1.0e3) + + ax[0].set_xlabel("Depths (km)") + ax[0].set_ylabel("Velocities (km/s)") + ax[0].legend() + fig.set_tight_layout(True) + fig.savefig("figures/stv_isotropic_velocities.pdf") + 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]: + print( + f"Plotting seismic properties at {delta_P/1.e9:.2f} GPa above the transition" + ) + P = transition_pressure(stishovite_scalar, T) + delta_P + + stishovite_relaxed.set_state(P, T) + stishovite_unrelaxed.set_state(P, T) + + print( + f"Is elastic tensor positive definite?: {is_positive_definite(stishovite_relaxed.isentropic_compliance_tensor)}" + ) + print( + f"Isotropic Poisson ratio: {stishovite_relaxed.isentropic_isotropic_poisson_ratio:.2f}" + ) + + plot_types = [ + "vp", + "vs1", + "vs2", + "linear compressibility", + "minimum poisson ratio", + "maximum poisson ratio", + ] + + fig = plt.figure(figsize=(12, 7)) + ax = [fig.add_subplot(2, 3, i, projection="polar") for i in range(1, 7)] + + if low_res_seismic_properties: + contour_sets, ticks, lines = plot_projected_elastic_properties( + stishovite_relaxed, plot_types, ax + ) + else: + contour_sets, ticks, lines = plot_projected_elastic_properties( + stishovite_relaxed, plot_types, ax, 181, 721, 100, 361 + ) + + for i in range(len(contour_sets)): + cbar = fig.colorbar(contour_sets[i], ax=ax[i], ticks=ticks[i], pad=0.1) + cbar.add_lines(lines[i]) + + fig.set_tight_layout(True) + fig.savefig( + f"figures/stishovite_seismic_properties_{P/1.e9:.2f}_GPa_{int(T)}_K.pdf" + ) + plt.show() + +data = common_data() + +PTV = np.concatenate( + ( + data["PTV"]["stv"]["Ito_1974"], + data["PTV"]["stv"]["Zhang_2021"], + data["PTV"]["poststv"]["Zhang_2021"], + data["PTV"]["stv"]["Andrault_2003"], + data["PTV"]["poststv"]["Andrault_2003"], + data["PTV"]["stv"]["Fischer_2018"], + data["PTV"]["poststv"]["Fischer_2018"], + ) +) +PTV_err = np.concatenate( + ( + data["PTV_err"]["stv"]["Ito_1974"], + data["PTV_err"]["stv"]["Zhang_2021"], + data["PTV_err"]["poststv"]["Zhang_2021"], + data["PTV_err"]["stv"]["Andrault_2003"], + data["PTV_err"]["poststv"]["Andrault_2003"], + data["PTV_err"]["stv"]["Fischer_2018"], + data["PTV_err"]["poststv"]["Fischer_2018"], + ) +) +abc = np.concatenate( + ( + data["abc"]["stv"]["Ito_1974"], + data["abc"]["stv"]["Zhang_2021"], + data["abc"]["poststv"]["Zhang_2021"], + data["abc"]["stv"]["Andrault_2003"], + data["abc"]["poststv"]["Andrault_2003"], + data["abc"]["stv"]["Fischer_2018"], + data["abc"]["poststv"]["Fischer_2018"], + ) +) +abc_err = np.concatenate( + ( + data["abc_err"]["stv"]["Ito_1974"], + data["abc_err"]["stv"]["Zhang_2021"], + data["abc_err"]["poststv"]["Zhang_2021"], + data["abc_err"]["stv"]["Andrault_2003"], + data["abc_err"]["poststv"]["Andrault_2003"], + data["abc_err"]["stv"]["Fischer_2018"], + data["abc_err"]["poststv"]["Fischer_2018"], + ) +) + +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, + ) +) + +Ito_mask = id == 1 +Zhang_mask = id == 2 +Andrault_mask = id == 3 +Fischer_mask = id == 4 + +PTV[Zhang_mask, 0] *= fP_Zhang + PTV[Zhang_mask, 0] * fP2_Zhang +PTV[Andrault_mask, 0] *= fP_Andrault + PTV[Andrault_mask, 0] * fP2_Andrault +P_for_CN *= fP_Zhang + P_for_CN * fP2_Zhang + +print( + f"Transition pressure at 3000 K: {transition_pressure(stishovite_scalar, 3000.0) / 1.0e9} GPa" +) + +stishovite_relaxed.set_composition([1]) +print(f"Consistent?: {check_anisotropic_eos_consistency(stishovite_relaxed)}") + +if plot_lnabc: + fig_lnabc = plt.figure(figsize=(8, 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"] + labels = [ + "Ito et al. (1974)", + "Zhang et al. (2021)", + "Andrault et al. (2003)", + "Fischer et al. (2018)", + ] + for i, mask in enumerate([Ito_mask, Zhang_mask, Andrault_mask, Fischer_mask]): + c = colors[i] + + lnV = deepcopy(np.log(PTV[mask, 2])) + if i == 0: + stishovite_scalar.set_composition([0.5, 0.5]) + stishovite_scalar.set_state(1.0e5, 298.15) + lnV += np.log(stishovite_scalar.V / PTV[mask, 2][0]) + + ax_lnabc[0].scatter(lnV, np.log(abc[mask, 0]), color=c) + ax_lnabc[0].errorbar( + lnV, + np.log(abc[mask, 0]), + xerr=PTV_err[mask, 2] / PTV[mask, 2], + yerr=abc_err[mask, 0] / abc[mask, 0], + ls="none", + color=c, + ) + ax_lnabc[0].scatter(lnV, np.log(abc[mask, 1]), color=c) + ax_lnabc[0].errorbar( + lnV, + np.log(abc[mask, 1]), + xerr=PTV_err[mask, 2] / PTV[mask, 2], + yerr=abc_err[mask, 1] / abc[mask, 1], + ls="none", + color=c, + ) + + ax_lnabc[0].scatter( + lnV, 0.5 * np.log(abc[mask, 0] * abc[mask, 1]), s=5, color=c + ) + ax_lnabc[1].scatter(lnV, np.log(abc[mask, 2]), color=c, label=labels[i]) + ax_lnabc[1].errorbar( + lnV, + np.log(abc[mask, 2]), + xerr=PTV_err[mask, 2] / PTV[mask, 2], + yerr=abc_err[mask, 2] / abc[mask, 2], + ls="none", + color=c, + ) + + for i in range(2): + ax_lnabc[i].set_xlabel("$\\ln(V)$") + + 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() + fig_lnabc.set_tight_layout(True) + fig_lnabc.savefig("figures/stv_lnabc.pdf") + plt.show() + + +if plot_G: + fig_G = plt.figure(figsize=(5, 4)) + ax_G = [fig_G.add_subplot(1, 1, 1)] + pressures = np.linspace(40.0e9, 100.0e9, 13) + ps = np.linspace(-0.25, 1.25, 201) + pressure_grid, p_grid = np.meshgrid(pressures, ps) + molar_fractions = np.moveaxis(np.array([p_grid, 1.0 - p_grid]), 0, 2) + T_grid = pressure_grid * 0.0 + 298.15 + G_xs = stishovite_scalar.evaluate( + ["excess_gibbs"], pressure_grid, T_grid, molar_fractions + )[0] + 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" + ) + else: + ax_G[0].plot(p_grid[:, i], G_xs[:, i] / 1000.0, linestyle=":") + 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].set_xlim(-0.25, 1.25) + ax_G[0].set_xlabel("$p_{Q1}$") + ax_G[0].set_ylabel("$\\mathcal{G}_{xs}$ (kJ/mol)") + ax_G[0].legend() + + # instantiate a second axes that shares the same y-axis + ax_G2 = ax_G[0].twiny() + ax_G2.set_xlim(-1.5, 1.5) + ax_G2.tick_params(axis="x") + ax_G2.set_xlabel("$Q$") + fig_G.set_tight_layout(True) + fig_G.savefig("figures/stv_G.pdf") + plt.show() + + +if plot_Q_V_abc: + fig_Q = plt.figure(figsize=(5, 4)) + gs = GridSpec(1, 2, width_ratios=[20, 1]) # 1 row, 2 columns + ax_Q = [fig_Q.add_subplot(gs[i]) for i in range(2)] + + fig_V = plt.figure(figsize=(8, 4)) + ax_V = [fig_V.add_subplot(1, 2, i) for i in range(1, 3)] + + fig_abc = plt.figure(figsize=(8, 4)) + ax_abc = [fig_abc.add_subplot(1, 2, i) for i in range(1, 3)] + + Tmin = np.min(PTV[:, 1]) + Tmax = 4000.0 + + cmap = "rainbow" + norm = matplotlib.colors.Normalize(vmin=Tmin, vmax=Tmax, clip=True) + mapper = cm.ScalarMappable(norm=norm, cmap=cmap) + + temperatures = np.array([4000.0, 3000.0, 2000.0, 1000.0, 298.15]) + stishovite_relaxed.set_composition([1.0]) + + for i, T in enumerate(temperatures): + print(f"Plotting Qs and volumes at {T} K") + color = mapper.to_rgba(T) + Pmin = np.max([1.0e5, (T - 300.0) * 1.0e7 - 10.0e9]) + pressures = np.linspace(Pmin, 120.0e9, 151) + Ts = pressures * 0.0 + T + (volumes, cell_params, molar_fractions) = stishovite_relaxed.evaluate( + ["V", "cell_parameters", "molar_fractions"], pressures, Ts + ) + stishovite_unrelaxed.set_composition([0.5, 0.5]) + (volumes_Q0, cell_params_Q0) = stishovite_unrelaxed.evaluate( + ["V", "cell_parameters"], pressures, Ts + ) + ax_Q[0].plot( + pressures / 1.0e9, molar_fractions[:, 0] - molar_fractions[:, 1], c=color + ) + a = cell_params[:, 0] + b = cell_params[:, 1] + c = cell_params[:, 2] + ax_V[0].plot(pressures / 1.0e9, volumes * 1.0e6, c=color) + ax_V[1].plot(pressures / 1.0e9, (volumes - volumes_Q0) * 1.0e6, c=color) + ax_abc[0].plot(pressures / 1.0e9, cell_params[:, 0] * 1.0e2, c=color) + ax_abc[0].plot(pressures / 1.0e9, cell_params[:, 1] * 1.0e2, c=color) + ax_abc[1].plot(pressures / 1.0e9, cell_params[:, 2] * 1.0e2, c=color) + + scatter = ax_V[0].scatter( + PTV[:, 0] / 1.0e9, PTV[:, 2] * 1.0e6, c=PTV[:, 1], cmap=cmap, norm=norm + ) + ax_abc[0].scatter( + PTV[:, 0] / 1.0e9, abc[:, 0] * 1.0e2, c=PTV[:, 1], cmap=cmap, norm=norm + ) + ax_abc[0].scatter( + PTV[:, 0] / 1.0e9, abc[:, 1] * 1.0e2, c=PTV[:, 1], cmap=cmap, norm=norm + ) + ax_abc[1].scatter( + PTV[:, 0] / 1.0e9, abc[:, 2] * 1.0e2, c=PTV[:, 1], cmap=cmap, norm=norm + ) + + ax_Q[0].set_ylim(0.0, 1.5) + + # plot rotation angles on Q plot + d = np.loadtxt("data/Zhang_2021_stishovite_angles.dat", unpack=True) + pressures_GPa = d[0] # *(fP_Zhang + d[0]*1.e9*fP2_Zhang) + Ts = pressures_GPa * 0.0 + 298.15 + angles, angles_err = [d[-2], d[-1]] + molar_fractions = stishovite_relaxed.evaluate( + ["molar_fractions"], pressures_GPa * 1.0e9, Ts + )[0] + Q_model = molar_fractions[:, 0] - molar_fractions[:, 1] + + def Q_misfit(args): + return np.sum(np.power((angles * args[0] - Q_model) / angles_err, 2.0)) + + sol = minimize(Q_misfit, [0.0]) + + # instantiate a second axes that shares the same x-axis + color = "tab:blue" + ax_Q2 = ax_Q[0].twinx() + ax_Q2.scatter(pressures_GPa, angles, color=color) + ax_Q2.errorbar(pressures_GPa, angles, yerr=angles_err, ls="none", color=color) + ax_Q2.set_ylim(0.0, 1.5 / sol.x) + ax_Q2.tick_params(axis="y", labelcolor=color) + ax_Q2.set_ylabel("SiO$_6$ rotation angle ($^{\\circ}$)", color=color) + + fig_Q.colorbar( + scatter, cax=ax_Q[1], orientation="vertical", label="Temperature (K)" + ) + + fig_V.colorbar(scatter, ax=ax_V[1], label="Temperature (K)") + fig_abc.colorbar(scatter, ax=ax_abc[1], label="Temperature (K)") + + for i in range(2): + ax_V[i].set_xlabel("Pressure (GPa)") + ax_abc[i].set_xlabel("Pressure (GPa)") + + ax_Q[0].set_xlabel("Pressure (GPa)") + ax_Q[0].set_ylabel("Q") + ax_V[0].set_ylabel("Volume (cm$^{3}$/mol)") + ax_V[1].set_ylabel("Excess volume (cm$^{3}$/mol)") + ax_abc[0].set_ylabel("$a$, $b$ length (cm/mol$^{\\frac{1}{3}}$)") + ax_abc[1].set_ylabel("$c$ length (cm/mol$^{\\frac{1}{3}}$)") + + fig_Q.set_tight_layout(True) + fig_V.set_tight_layout(True) + fig_abc.set_tight_layout(True) + + fig_Q.savefig("figures/stv_Q.pdf") + fig_V.savefig("figures/stv_V.pdf") + fig_abc.savefig("figures/stv_abc.pdf") + +if plot_relaxed: + fig_relaxed = plt.figure(figsize=(12, 4)) + fig_Q0 = plt.figure(figsize=(12, 4)) + + ax_relaxed = [fig_relaxed.add_subplot(1, 3, i) for i in range(1, 4)] + ax_Q0 = [fig_Q0.add_subplot(1, 3, i) for i in range(1, 4)] + + T = 298.15 + pressures = np.linspace(1.0e5, 70.0e9, 101) + Ts = pressures * 0.0 + T + + (molar_fractions, C_N, K_NR) = stishovite_relaxed.evaluate( + [ + "molar_fractions", + "isentropic_stiffness_tensor", + "isentropic_bulk_modulus_reuss", + ], + pressures, + Ts, + ) + (C_N_u, K_NR_u) = stishovite_relaxed.unrelaxed.evaluate( + ["isentropic_stiffness_tensor", "isentropic_bulk_modulus_reuss"], + pressures, + Ts, + molar_fractions, + ) + + stishovite_unrelaxed.set_composition([0.5, 0.5]) + (C_N_Q0, K_NR_Q0) = stishovite_unrelaxed.evaluate( + ["isentropic_stiffness_tensor", "isentropic_bulk_modulus_reuss"], pressures, Ts + ) + stishovite_unrelaxed.set_composition([1.0, 0.0]) + (C_N_Q1, K_NR_Q1) = stishovite_unrelaxed.evaluate( + ["isentropic_stiffness_tensor", "isentropic_bulk_modulus_reuss"], pressures, Ts + ) + + c = ax_relaxed[1].plot(pressures / 1.0e9, K_NR / 1.0e9, label="$K_{NR}$") + 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_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)$") + + for axi, i, j in ( + (0, 0, 0), + (0, 0, 1), + (0, 1, 1), + (1, 0, 2), + (1, 1, 2), + (1, 2, 2), + (2, 3, 3), + (2, 4, 4), + (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}}}$" + ) + color = c[0].get_color() + ax_relaxed[axi].plot( + pressures / 1.0e9, C_N_Q0[:, i, j] / 1.0e9, linestyle=":", c=color + ) + ax_relaxed[axi].scatter( + P_for_CN / 1.0e9, CN_GPa[:, i, j], label=f"{i+1}{j+1}", color=color + ) + ax_relaxed[axi].errorbar( + P_for_CN / 1.0e9, + CN_GPa_orig[:, i, j], + yerr=CN_err_GPa[:, i, j], + alpha=0.5, + color=color, + ls="none", + ) + ax_relaxed[axi].scatter( + P_for_CN / 1.0e9, CN_GPa_orig[:, i, j], s=5, alpha=0.5, color=color + ) + + ax_Q0[axi].plot( + pressures / 1.0e9, C_N_Q0[:, i, j] / 1.0e9, label=f"{i+1}{j+1} (Q0)" + ) + ax_Q0[axi].plot( + pressures / 1.0e9, C_N_Q1[:, i, j] / 1.0e9, label=f"{i+1}{j+1} (Q1)" + ) + + for i in range(3): + ax_relaxed[i].legend() + ax_Q0[i].legend() + + ax_relaxed[i].set_xlabel("Pressure (GPa)") + ax_Q0[i].set_xlabel("Pressure (GPa)") + + ax_relaxed[i].set_ylabel("Modulus (GPa)") + ax_Q0[i].set_ylabel("Modulus (GPa)") + + fig_relaxed.set_tight_layout(True) + fig_Q0.set_tight_layout(True) + + fig_relaxed.savefig("figures/stv_relaxed.pdf") + + plt.show() diff --git a/contrib/anisotropic_stishovite/stishovite_parameters.py b/contrib/anisotropic_stishovite/stishovite_parameters.py new file mode 100644 index 00000000..316c8c06 --- /dev/null +++ b/contrib/anisotropic_stishovite/stishovite_parameters.py @@ -0,0 +1,92 @@ +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_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, +] + + +# a0Q1, b0Q1, +# PsiI_33_a, PsiI_33_b, PsiI_33_c, +# PsiI_33_b2, PsiI_33_c2 +# 1151 +cell_args = [ + 2.72829954e-02, + 2.85720276e-02, + 2.10496560e-01, + -6.02250364e-01, + 1.11223359e00, + 8.20161603e-05, + -1.04689736e01, +] + + +# a11, a22, a33, a44, a55, a66 +# b11, b33, b44, b66 +# c44, c66 +# b112, b332 + +# let the b parameters be equal to +# frel = -0.14 +# b1_input = b1 * c1 * exp(c1 * frel) - 1) +# 178 +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, +] + +scalar_and_cell_args = np.concatenate((scalar_args, cell_args)) +all_args = np.concatenate((scalar_args, cell_args, elastic_args)) + +scalar_bounds = ( + (-8.0e-3, -4.0e-3), # dVQ0 + (-10.0, 0.0), # dKQ0 + (-0.5, 2.0), # dKpQ0 4.0292, + (-0.6, 1.0), # dgrQ0 1.55674, + (-2.0, 3.0), # dqQ0 2.2096, + (0.99, 0.999), # V0Q1overV0Q0 + (0.0, 80.0), # dDebye_0 + (46, 53), # P_tr_GPa, + (0.9, 1.1), # s + (0.9, 1.1), # s + (-0.1, 0.1), # s + (-0.1, 0.1), +) # s + +cell_bounds = ( + (2.6e-2, 2.8e-2), # a0Q1 + (2.7e-2, 3.0e-2), # b0Q1 + (0.0, 1.0), # PsiI_33_a + (-1.0, -0.1), # PsiI_33_b + (0.5, 2), # PsiI_33_c + (-1.0e-2, 1.0e-2), # PsiI_33_b2 + (-20, -2), +) # PsiI_33_c2 + +elastic_bounds = tuple((None, None) for i in range(14))