diff --git a/contrib/anisotropic_stishovite/data/Chen_2024_stishovite_velocities.dat b/contrib/anisotropic_stishovite/data/Chen_2024_stishovite_velocities.dat new file mode 100644 index 00000000..61331049 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Chen_2024_stishovite_velocities.dat @@ -0,0 +1,31 @@ +# Pressure(GPa) Uncertainty(GPa) Temperature(°C) Uncertainty(°C) VP(km/s) Uncertainty(km/s) VS(km/s) Uncertainty(km/s) Density(g/cm³) Uncertainty(g/cm³) KS(GPa) Uncertainty(GPa) G(GPa) Uncertainty(GPa) +14.5 0.2 800 5 12.277 0.061 7.223 0.036 4.417 0.002 353.9 0.53 227.7 0.34 +13.8 0.2 800 5 12.238 0.061 7.214 0.036 4.407 0.002 349.7 0.52 226.7 0.34 +13.4 0.2 800 5 12.222 0.061 7.201 0.036 4.402 0.002 348.7 0.52 225.6 0.34 +12.1 0.2 800 5 12.163 0.061 7.200 0.036 4.386 0.002 341.4 0.51 224.7 0.34 +10.9 0.2 800 5 12.099 0.060 7.196 0.036 4.370 0.002 333.9 0.50 223.6 0.34 +14.3 0.2 600 5 12.316 0.062 7.247 0.036 4.431 0.002 357.2 0.54 230.0 0.35 +13.5 0.2 600 5 12.299 0.061 7.256 0.036 4.421 0.002 353.8 0.53 230.0 0.35 +13.4 0.2 600 5 12.264 0.061 7.244 0.036 4.420 0.002 351.0 0.53 229.2 0.34 +11.7 0.2 600 5 12.210 0.061 7.234 0.036 4.399 0.002 344.5 0.52 227.5 0.34 +10.6 0.2 600 5 12.162 0.061 7.241 0.036 4.384 0.002 337.7 0.51 227.1 0.34 +13.9 0.2 400 5 12.378 0.062 7.300 0.036 4.443 0.002 360.4 0.54 233.9 0.35 +13.0 0.2 400 5 12.349 0.062 7.301 0.037 4.432 0.002 356.2 0.53 233.4 0.35 +12.9 0.2 400 5 12.325 0.062 7.290 0.036 4.430 0.002 354.4 0.53 232.7 0.35 +11.6 0.2 400 5 12.267 0.061 7.293 0.036 4.414 0.002 346.8 0.52 232.0 0.35 +10.2 0.2 400 5 12.207 0.061 7.284 0.036 4.396 0.002 339.8 0.51 230.5 0.35 +13.5 0.2 200 5 12.424 0.062 7.336 0.037 4.452 0.002 363.0 0.54 236.7 0.36 +12.7 0.2 200 5 12.398 0.062 7.339 0.037 4.442 0.002 359.1 0.54 236.4 0.35 +12.6 0.2 200 5 12.383 0.062 7.323 0.037 4.440 0.002 358.7 0.54 235.3 0.35 +11.3 0.2 200 5 12.337 0.062 7.332 0.037 4.424 0.002 351.7 0.53 235.0 0.35 +9.9 0.2 200 5 12.275 0.061 7.329 0.037 4.407 0.002 344.1 0.52 233.9 0.35 +13.1 0.2 27 5 12.445 0.062 7.362 0.037 4.458 0.002 363.5 0.55 238.7 0.36 +12.2 0.2 27 5 12.425 0.062 7.357 0.037 4.447 0.002 360.9 0.54 237.8 0.36 +12.2 0.2 27 5 12.424 0.062 7.356 0.037 4.447 0.002 360.9 0.54 237.7 0.36 +10.7 0.2 27 5 12.362 0.062 7.347 0.037 4.428 0.002 353.4 0.53 236.2 0.35 +9.8 0.2 27 5 12.320 0.062 7.352 0.037 4.417 0.002 347.6 0.52 235.9 0.35 +9.6 0.2 27 5 12.317 0.062 7.355 0.037 4.415 0.002 346.9 0.52 235.9 0.35 +7.6 0.2 27 5 12.242 0.061 7.347 0.037 4.389 0.002 337.6 0.51 234.0 0.35 +6.3 0.2 27 5 12.196 0.061 7.342 0.037 4.372 0.002 331.9 0.50 232.8 0.35 +4.7 0.2 27 5 12.163 0.061 7.331 0.037 4.352 0.002 327.9 0.49 231.1 0.35 +3.5 0.2 27 5 12.124 0.061 7.326 0.037 4.336 0.002 323.1 0.48 229.9 0.34 diff --git a/contrib/anisotropic_stishovite/data/Grocholski_2013_stishovite_post_stishovite_unit_cell.dat b/contrib/anisotropic_stishovite/data/Grocholski_2013_stishovite_post_stishovite_unit_cell.dat new file mode 100644 index 00000000..9e44c778 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Grocholski_2013_stishovite_post_stishovite_unit_cell.dat @@ -0,0 +1,14 @@ +# P_Gpa P_Uncertainty a_p (Å) a_p_Uncertainty a a_Uncertainty b b_Uncertainty c c_Uncertainty V V_Uncertainty +0.0001 0.00001 4.028 0.003 4.182 0.002 4.182 0.002 2.666 0.002 14.04 0.01 +7 0.3 4.01 0.003 4.153 0.002 4.153 0.002 2.652 0.002 13.77 0.01 +18 0.4 3.966 0.003 4.105 0.003 4.105 0.003 2.641 0.004 13.39 0.02 +31.2 0.6 3.909 0.003 4.048 0.004 4.048 0.004 2.632 0.005 12.98 0.03 +39.2 0.8 3.88 0.003 4.031 0.004 4.031 0.004 2.605 0.003 12.74 0.02 +51 1.2 3.842 0.003 4.024 0.009 3.952 0.009 2.599 0.005 12.44 0.02 +60.5 1.4 3.815 0.003 4.023 0.009 3.899 0.009 2.567 0.005 12.13 0.02 +68.4 1.6 3.795 0.003 4.023 0.009 3.865 0.008 2.558 0.005 11.98 0.02 +72.6 1.8 3.785 0.003 4.015 0.009 3.848 0.004 2.551 0.005 11.86 0.02 +79.7 2.1 3.769 0.004 4.005 0.009 3.831 0.002 2.54 0.005 11.75 0.02 +84.6 2.2 3.758 0.003 4 0.004 3.82 0.003 2.534 0.005 11.66 0.02 +89.2 2.4 3.748 0.002 3.986 0.004 3.812 0.005 2.528 0.003 11.56 0.02 +95 2.6 3.737 0.002 3.981 0.002 3.801 0.004 2.523 0.002 11.5 0.01 \ No newline at end of file diff --git a/contrib/anisotropic_stishovite/data/Jiang_2009_elastic.dat b/contrib/anisotropic_stishovite/data/Jiang_2009_elastic.dat new file mode 100644 index 00000000..eb3dd0a8 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Jiang_2009_elastic.dat @@ -0,0 +1,11 @@ +# P C11 C33 C13 C12 C44 C66 +0.00 453.96 761.13 192.32 199.49 258.16 320.53 +0.80 455.95 780.80 202.15 205.89 261.31 322.27 +1.60 455.94 766.16 196.05 209.20 256.49 326.25 +3.00 464.22 770.07 203.66 219.58 262.38 327.28 +4.70 465.89 780.01 199.85 225.82 265.21 324.86 +6.40 483.06 803.88 215.77 227.54 265.89 333.94 +8.00 475.39 818.42 212.65 237.11 279.44 331.23 +9.40 507.55 814.30 199.46 266.85 276.75 340.35 +10.80 512.27 818.92 201.81 281.35 281.18 347.03 +12.00 517.94 824.27 208.02 283.46 277.17 350.44 \ No newline at end of file diff --git a/contrib/anisotropic_stishovite/data/Murakami_2003_post_stishovite_unit_cell.dat b/contrib/anisotropic_stishovite/data/Murakami_2003_post_stishovite_unit_cell.dat new file mode 100644 index 00000000..78c37247 --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Murakami_2003_post_stishovite_unit_cell.dat @@ -0,0 +1,5 @@ +# Run_No. P(GPa) a(Å) a_uncertainty(Å) b(Å) b_uncertainty(Å) c(Å) c_uncertainty(Å) V(cm3/mol) V_uncertainty(cm3/mol) +1 87 3.982 0.006 3.809 0.006 2.525 0.004 11.53 0.03 +1 93 3.973 0.004 3.802 0.005 2.508 0.002 11.41 0.02 +3 124 3.929 0.002 3.745 0.005 2.465 0.001 10.92 0.02 +3 136 3.927 0.004 3.718 0.009 2.452 0.002 10.78 0.03 diff --git a/contrib/anisotropic_stishovite/data/Sun_et_al_2019_post_stishovite_volumes.dat b/contrib/anisotropic_stishovite/data/Sun_et_al_2019_post_stishovite_volumes.dat new file mode 100644 index 00000000..614e5fbf --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Sun_et_al_2019_post_stishovite_volumes.dat @@ -0,0 +1,42 @@ +# P(GPa) P_uncertainty(GPa) T(K) T_uncertainty(K) V_SiO2(ų) V_SiO2_uncertainty(ų) V_NaCl(ų) V_NaCl_uncertainty(ų) +64.8 0.8 300 5 40.28 0.04 23.06 0.07 +68.1 0.6 300 5 40.00 0.05 22.77 0.05 +79.5 0.6 300 5 38.91 0.6 21.88 0.04 +81.7 0.7 300 5 38.73 0.4 21.76 0.04 +87.7 1.2 300 5 38.49 0.4 21.32 0.09 +89.8 1.8 300 5 38.22 0.3 21.19 0.12 +91.1 1.8 300 5 38.17 0.6 21.10 0.13 +96.2 1.4 300 5 37.29 0.4 20.47 0.06 +103.1 1.4 300 5 37.16 0.7 20.41 0.8 +109.5 1.6 300 5 36.87 0.5 20.08 0.12 +119.3 1.3 300 5 36.33 0.4 19.59 0.10 +125.4 2.6 300 5 36.05 0.7 19.36 0.07 +83.4 1.6 1150 120 38.42 0.5 21.88 0.9 +116.0 2.0 1150 120 36.82 0.5 19.95 0.12 +64.4 1.4 1180 120 40.46 0.4 23.46 0.9 +65.1 1.8 1580 160 40.56 0.3 23.56 0.6 +85.6 9.9 1450 150 38.88 0.4 21.87 0.9 +92.7 9.9 1550 160 38.21 0.5 21.36 0.6 +109.9 9.9 1500 150 37.19 0.3 20.35 0.6 +127.3 3.2 1500 150 36.28 0.6 19.50 0.12 +141.2 2.8 1600 160 35.73 0.7 18.94 0.15 +65.7 1.4 1700 170 40.65 0.3 23.57 0.9 +71.8 1.4 1750 180 40.16 0.3 23.02 0.10 +92.8 1.8 1800 180 38.45 0.4 21.44 0.11 +96.9 1.8 1700 170 38.21 0.5 21.15 0.12 +117.8 2.4 1700 170 36.89 0.5 19.95 0.15 +128.3 2.6 1700 170 36.28 0.8 19.50 0.16 +142.0 4.8 1700 170 35.73 1.0 18.93 0.17 +86.8 1.8 1900 190 38.95 0.3 21.88 0.10 +97.5 2.0 1930 190 38.21 0.5 21.17 0.12 +117.9 2.4 1930 190 37.02 0.5 20.04 0.13 +128.9 2.6 1900 190 36.30 0.6 19.52 0.15 +73.9 1.4 2100 210 40.30 0.3 22.98 0.8 +90.6 1.8 2120 210 38.89 0.3 21.68 0.11 +96.6 2.0 2100 210 38.31 0.4 21.29 0.12 +97.7 2.0 2250 230 38.29 0.5 21.25 0.12 +98.3 2.0 2130 210 38.25 0.5 21.18 0.13 +114.7 2.4 2100 210 37.33 0.7 20.25 0.14 +119.5 2.8 2100 210 37.00 0.7 20.01 0.14 +128.5 3.6 2200 220 36.38 0.4 19.61 0.14 +131.1 4.5 2200 220 36.34 0.9 19.49 0.14 diff --git a/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_elastic_tensor.dat b/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_elastic_tensor.dat index 3f31e827..6df7ee6d 100644 --- a/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_elastic_tensor.dat +++ b/contrib/anisotropic_stishovite/data/Zhang_2021_stishovite_elastic_tensor.dat @@ -1,31 +1,30 @@ -# 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)  +# P  err C11  err C12  err C13  err C22  err C23  err C33  err C44  err C55  err C66  err ρ  err # 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 +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 +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 +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.299999 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_2023_stishovite_unit_cell.dat b/contrib/anisotropic_stishovite/data/Zhang_2023_stishovite_unit_cell.dat new file mode 100644 index 00000000..b069d39f --- /dev/null +++ b/contrib/anisotropic_stishovite/data/Zhang_2023_stishovite_unit_cell.dat @@ -0,0 +1,27 @@ +# P σ_P a σ_a b σ_b c σ_c V σ_V Unique_refl Rint_percent R1_percent +0 0.1 4.1752 0.0001 4.1752 0.0001 2.6642 0.0001 46.443 0.003 2788 4.63 1.29 +2.8 0.1 4.166 0.0003 4.166 0.0003 2.664 0.0003 46.24 0.01 55 0.61 3.48 +7.8 0.2 4.1416 0.0005 4.1416 0.0005 2.6564 0.0003 45.57 0.01 48 1.52 5.93 +13 0.1 4.12 0.0004 4.12 0.0004 2.6458 0.0003 44.91 0.01 37 0.78 6.28 +16 0.3 4.1066 0.0004 4.1066 0.0004 2.6433 0.0004 44.58 0.01 50 5.89 6.14 +16.9 0.1 4.1045 0.0004 4.1045 0.0004 2.6393 0.0003 44.464 0.008 48 0.35 6.11 +19.7 0.1 4.0891 0.0004 4.0891 0.0004 2.6298 0.0017 43.97 0.03 50 1.65 3.78 +21.4 0.2 4.0875 0.0004 4.0875 0.0004 2.6366 0.0003 44.05 0.01 41 12.74 6 +26.8 0.2 4.0681 0.0005 4.0681 0.0005 2.6259 0.0005 43.46 0.01 43 4.37 7.84 +28.5 0.2 4.052 0.0004 4.052 0.0004 2.6162 0.0018 42.95 0.03 52 1.7 7.73 +33.8 0.2 4.0318 0.0006 4.0318 0.0006 2.607 0.0008 42.38 0.02 48 5.8 6.17 +40.4 0.2 4.0133 0.0007 4.0133 0.0007 2.6 0.003 41.88 0.05 34 8.91 5.61 +48.7 0.2 3.9875 0.0007 3.9875 0.0007 2.584 0.003 41.09 0.05 46 1.06 5.65 +49.8 0.2 3.9819 0.001 3.9819 0.001 2.581 0.003 40.92 0.1 36 1.01 5.35 +52.4 0.2 3.944 0.003 4.015 0.0019 2.5851 0.0013 40.93 0.04 47 1.25 5.66 +54.2 0.2 3.932 0.003 4.0128 0.0018 2.5817 0.0012 40.73 0.04 56 0.87 5.38 +55.6 0.2 3.93 0.004 4.0097 0.001 2.575 0.0005 40.58 0.04 54 4.82 6.22 +58.6 0.3 3.914 0.005 4.0118 0.0012 2.5771 0.0006 40.38 0.05 31 14.89 4.97 +62 0.3 3.901 0.004 4.0089 0.0012 2.5656 0.0008 40.12 0.04 41 2.94 6.88 +64.4 0.3 3.888 0.004 4.008 0.002 2.5681 0.0014 40.01 0.05 46 4.8 6.66 +65.8 0.3 3.882 0.004 4.007 0.001 2.5667 0.0008 39.93 0.04 41 12.9 6.62 +68 0.3 3.871 0.004 4.0051 0.001 2.5661 0.0007 39.71 0.04 41 1.4 6.85 +71 0.3 3.858 0.003 4.004 0.0009 2.558 0.0007 39.51 0.03 46 1.91 4.57 +73.8 0.3 3.85 0.003 4.0016 0.0008 2.5557 0.0006 39.37 0.03 45 2.6 6.05 +75.3 0.3 3.838 0.002 3.9974 0.0007 2.5498 0.0005 39.12 0.02 46 0.55 4.42 +# P errors from the 2021 preprint on ESSOAR diff --git a/contrib/anisotropic_stishovite/model_output/covariance_matrix.dat b/contrib/anisotropic_stishovite/model_output/covariance_matrix.dat index 2ceb515f..9895ff89 100644 --- a/contrib/anisotropic_stishovite/model_output/covariance_matrix.dat +++ b/contrib/anisotropic_stishovite/model_output/covariance_matrix.dat @@ -1,37 +1,29 @@ -2.969072323847245294e-05 -2.414319325043906995e-03 9.735316303185922192e-07 -2.034079978431028023e-05 1.597751501944453352e-04 5.032599281155131899e-07 -2.088564195265813286e-03 1.141858121336443432e-03 3.949357292634677550e-05 1.850758822971476200e-05 3.299016538863077687e-05 9.077296660207598603e-05 -5.137345609827280432e-16 -8.143086031774118376e-17 -4.589020455106308953e-16 -3.098661009272525452e-15 2.582109938010748615e-08 2.624612085038632830e-08 3.608749591974849653e-06 -6.937855671113378833e-05 -1.371705693477332281e-04 -2.071765376120491374e-07 -2.516517250145258906e-03 -3.928954203215390470e-05 1.279490468007478382e-06 9.068426317865086325e-08 -1.010411452316174498e-05 -1.002915587787911641e-05 -7.096411239012776800e-06 -6.166408270828315896e-05 7.954518997506867072e-06 -1.457372693470361177e-03 -1.665472713324716442e-03 5.402267720757495664e-04 1.943321742220476944e-03 -3.333863796457967564e-04 -5.468710926687100291e-05 --2.414319325043906995e-03 4.249826663546418248e-01 -3.094884841608068937e-04 5.573670640395235942e-04 1.149432319847632812e-02 -4.490399119970610087e-05 1.665240607391086369e-01 -6.504618095173082815e-02 -2.282360652002456136e-03 -7.476555221169896813e-04 -2.364247942489634824e-03 -7.079281708105459453e-03 3.475237159535944607e-14 5.352818153943201318e-15 3.945412998590303000e-14 2.598336304972241831e-13 -2.151349431467241060e-06 -2.167558592309546433e-06 -2.808918726555411629e-04 1.185104206621878306e-02 1.657883033879654217e-02 1.543514452236958691e-05 2.136226422370665345e-01 2.137533177873366264e-03 2.287660214046896914e-03 2.979745879017003073e-04 1.676699508613893001e-03 1.732916157635252218e-03 1.298896908940160239e-03 9.282782630378589445e-03 -6.763338260323719427e-04 9.451541906012153604e-02 1.143881756135775030e-01 -3.500456106541598861e-02 -1.335569570692768659e-01 2.491333162046089483e-02 4.670416326687474098e-03 -9.735316303185922192e-07 -3.094884841608068937e-04 7.262712937098192892e-06 -9.849114812544193482e-06 -8.091894307406850438e-05 -8.384368839662062259e-08 2.757286574523949750e-04 -3.302108663012920522e-05 -7.940751942138309900e-07 4.870195465031972041e-07 -5.201307034047055083e-08 8.923622864770794612e-07 2.587147620834972909e-17 -8.072235707204832032e-18 4.295548432640169746e-18 -3.705586445323490784e-17 -3.382128028653160845e-10 -3.237635627080133576e-10 2.880580686211860440e-07 -5.630122258543856047e-05 -5.550398399081168347e-05 -1.368692423501269688e-08 -5.511266268928062019e-04 5.445097707514176965e-06 1.881736127791208387e-06 1.074217982496532016e-06 -9.573321138379633936e-07 -8.955180902632183681e-07 -8.519837545226716281e-07 -3.229239994559254898e-05 3.022377409523564049e-07 1.005144316480618646e-04 6.291873262212429902e-05 -3.646627150018448262e-05 -7.162025230758143908e-05 -2.353051652724181582e-05 -4.818973351892917084e-06 --2.034079978431028023e-05 5.573670640395235942e-04 -9.849114812544193482e-06 3.362871704584756774e-03 7.747114374398873272e-03 -1.609643254583494532e-06 6.691987623949255913e-03 -3.052808866339202370e-03 6.414389714812255498e-05 -1.154995575737634931e-05 3.264042564485330100e-04 5.367620605078604222e-04 -1.764022828001739892e-15 -8.026163323787587436e-17 -4.187898242504952461e-15 -1.643837790116418740e-14 -2.415063062091929880e-08 -3.005924930827426758e-08 2.969477106024997500e-07 -2.599807913875601048e-03 -2.502578179730502710e-03 -3.448229919827389252e-08 9.894412705418731673e-04 1.614015226545799700e-04 1.096321509004083357e-05 7.773663567167847920e-07 6.743863193988200014e-06 6.281162398372516905e-06 4.516749168953107673e-06 -1.923508018533116436e-03 1.331924023050925137e-05 -5.849646589347031259e-03 -5.178846973454946830e-03 2.153268188959647257e-03 6.003639016901835684e-03 1.387439896431741414e-03 8.672656714053429282e-05 -1.597751501944453352e-04 1.149432319847632812e-02 -8.091894307406850438e-05 7.747114374398873272e-03 1.331325052127259045e-01 -2.029319015312089713e-05 5.925313582364208354e-02 -3.080991714587031846e-02 1.977215830163714296e-03 2.741585182335783783e-04 2.896860335651843684e-04 1.694321705139406218e-03 -4.032963359981584820e-14 -3.038866244634352617e-15 -1.406216877566180905e-14 -7.731203689154381191e-14 -3.333970888846127249e-10 -8.709526033092263374e-08 5.768850215532460621e-05 -1.808405967871253364e-02 -1.896972119286157929e-02 -3.448716959724578149e-06 3.799553542857860011e-02 1.185191496548233239e-03 1.182857261486174993e-03 5.547822027831921226e-05 1.003402092902857769e-04 1.077696970302018857e-04 8.883519413119686178e-05 -1.588843399738775614e-02 3.661169141778293729e-05 -1.294007493555653365e-01 -1.216020000364273906e-01 4.772874293545138064e-02 1.412309683815566597e-01 1.866643744552765397e-02 1.065831122547543450e-03 -5.032599281155131899e-07 -4.490399119970610087e-05 -8.384368839662062259e-08 -1.609643254583494532e-06 -2.029319015312089713e-05 1.534491722260848744e-07 -5.010148782737223420e-04 1.389084317312687150e-04 8.027025917434763237e-07 -3.086014703537749928e-07 5.402179638462372378e-07 1.348672238299724075e-06 -1.243086105691209076e-17 1.406908605920547081e-17 -5.669532384363459182e-18 -4.265439752377840368e-17 2.258055017134934322e-09 1.973339258111572039e-09 -2.312332839147157345e-07 -2.916663369361920704e-05 -2.519716699320902141e-05 6.390334241489579055e-09 -2.181936752914516461e-05 -2.885741513929382642e-06 -9.584663865549323316e-06 -5.294172336870067481e-09 -5.146313275530546443e-07 -5.539881244564768786e-07 -2.523146163460764631e-07 -1.594898464460605111e-05 8.444660670564603750e-08 -2.640464914727091680e-05 -3.561425506840147387e-05 9.788159208161594391e-06 4.150511446075701752e-05 -4.390083089488345825e-05 4.229752159922561635e-07 --2.088564195265813286e-03 1.665240607391086369e-01 2.757286574523949750e-04 6.691987623949255913e-03 5.925313582364208354e-02 -5.010148782737223420e-04 2.104896430272050534e+00 -4.981610783074122972e-01 -3.365689227855006831e-03 6.453674787356451911e-04 -2.134808163348684597e-03 -5.585079353706507398e-03 4.763216963682577079e-14 -4.384607976308321193e-14 2.390016003294783340e-14 1.802000319076651299e-13 -7.553474453448608289e-06 -6.879771852291717019e-06 6.747976667041734062e-04 1.029989985721349172e-01 9.065041706401480437e-02 -2.040043730075414874e-05 4.725134464369972614e-02 1.027426075386230972e-02 3.051814272148788490e-02 1.825216372649425331e-05 1.755644760324465278e-03 1.870870624454052007e-03 8.708104448430653436e-04 5.720248449212836372e-02 -3.219229529387673773e-04 1.057903669450520889e-01 1.403403673864646317e-01 -3.924802079810447292e-02 -1.636564359123428292e-01 1.451439494535149144e-01 -1.738338909140980243e-03 -1.141858121336443432e-03 -6.504618095173082815e-02 -3.302108663012920522e-05 -3.052808866339202370e-03 -3.080991714587031846e-02 1.389084317312687150e-04 -4.981610783074122972e-01 2.170061849111732388e-01 1.828695319578236243e-03 3.001159231045504765e-04 1.312031061506022458e-03 3.225448143816585494e-03 -1.887583244373389328e-14 9.531646394130322699e-15 -1.395286465047703843e-14 -1.010094527683821402e-13 2.222821203875844844e-06 2.537007182150678717e-06 -6.392355604498717206e-05 -3.473801794131257253e-02 -3.257967242021848203e-02 3.079233917323632859e-06 -3.340170502487335713e-02 -5.494265105720975290e-03 -5.915567986028016229e-03 3.383003929474542214e-05 -5.605255102331018402e-04 -5.695318433638057315e-04 -2.820483943677134196e-04 -1.955720914959569248e-02 2.528818449588239310e-04 -4.322324541853273666e-02 -6.023684334790740946e-02 1.611303878009478704e-02 7.044389440078994979e-02 -4.483966495114385248e-02 -8.233919291224370057e-05 -3.949357292634677550e-05 -2.282360652002456136e-03 -7.940751942138309900e-07 6.414389714812255498e-05 1.977215830163714296e-03 8.027025917434763237e-07 -3.365689227855006831e-03 1.828695319578236243e-03 1.016967765522213269e-04 2.704638970645951757e-05 4.532050448181657753e-05 1.357918390547996671e-04 -1.592926120230101003e-15 -9.821249544289911861e-17 -7.439606753073941868e-16 -4.845003724100701487e-15 3.813363256363073580e-08 3.857585521230416037e-08 3.200685394482180146e-06 -2.112929067660789938e-04 -2.858688383419348751e-04 -2.408182879736639907e-07 -2.713414985469090809e-03 -2.657470988356218109e-05 -2.803719763090280236e-05 2.580429643797891438e-07 -1.010510257238622547e-05 -9.686357516880338036e-06 -5.814282625649671241e-06 -1.795567582710501570e-04 1.072466086869132534e-05 -4.802813362752815160e-03 -4.963194473939220967e-03 1.775867457136094804e-03 5.777840977318744395e-03 -3.405953240386222296e-04 -6.056307379456799997e-05 -1.850758822971476200e-05 -7.476555221169896813e-04 4.870195465031972041e-07 -1.154995575737634931e-05 2.741585182335783783e-04 -3.086014703537749928e-07 6.453674787356451911e-04 3.001159231045504765e-04 2.704638970645951757e-05 8.634980512679053723e-05 2.160796390805643907e-05 5.824445473519590950e-05 -3.264625985174926743e-16 -7.923909448595848997e-16 -2.869213732277017375e-16 -1.943758295631656079e-15 8.266752272149200309e-09 9.778994440663878912e-09 3.524928150674835131e-06 9.244866045590156400e-05 2.724159914307262731e-05 -1.754112212171557053e-07 -1.639613353212979532e-03 -1.838827713483549386e-05 4.881136721097782943e-05 1.113170450201948483e-06 -2.085181637611447406e-06 -1.662825435048065551e-06 -1.529648396480678433e-06 3.754504312323681382e-05 5.070866914767413610e-06 -9.756958387685016505e-04 -1.072946360077371318e-03 3.617995149859560026e-04 1.251930415222073795e-03 -5.528075146326992567e-05 -3.912271590726792503e-05 -3.299016538863077687e-05 -2.364247942489634824e-03 -5.201307034047055083e-08 3.264042564485330100e-04 2.896860335651843684e-04 5.402179638462372378e-07 -2.134808163348684597e-03 1.312031061506022458e-03 4.532050448181657753e-05 2.160796390805643907e-05 3.741988925585817429e-04 1.659927024235176589e-04 -5.879810778647083869e-16 -9.391283520432811530e-17 -7.571720499217147661e-15 -5.255400233547673271e-15 2.856267415224482285e-08 2.904112999090483503e-08 4.142496374942100873e-06 -2.570082280690285368e-04 -3.251801223006570694e-04 -2.402523089244516050e-07 -3.216676545935040574e-03 -3.961294658191956365e-05 2.107728167369185196e-06 6.525664614377412270e-07 -9.965538983237963089e-06 -9.808018743190203645e-06 -6.833063562445468623e-06 -1.870407675260386544e-04 1.107990032257984879e-05 -1.679299878847527938e-03 -1.899267739233110111e-03 6.223804919963958409e-04 2.215661331552164909e-03 -3.730348743630668195e-04 -6.364872634283788350e-05 -9.077296660207598603e-05 -7.079281708105459453e-03 8.923622864770794612e-07 5.367620605078604222e-04 1.694321705139406218e-03 1.348672238299724075e-06 -5.585079353706507398e-03 3.225448143816585494e-03 1.357918390547996671e-04 5.824445473519590950e-05 1.659927024235176589e-04 9.861110209393142770e-03 -1.891966634045982456e-15 -2.711700981994465871e-16 -2.204199540613405901e-15 -5.253591335895288997e-13 7.794795518439994591e-08 7.846052278910172042e-08 1.148661919795845201e-05 -6.454285445130538692e-04 -8.440052609354726587e-04 -6.640530256914299064e-07 -7.926725382679215490e-03 -1.011002018751854773e-04 9.225664148879019028e-06 1.051179776515564498e-06 -2.918386730194719129e-05 -2.889439043922104893e-05 -2.028183930935848183e-05 -5.055290379226907082e-04 2.773567462221827859e-05 -5.498729953080172886e-03 -6.050663216414568953e-03 2.036131827085406725e-03 7.054155462547743484e-03 -8.578412202787384210e-04 -1.606945507961290737e-04 --5.137345609827280432e-16 3.475237159535944607e-14 2.587147620834972909e-17 -1.764022828001739892e-15 -4.032963359981584820e-14 -1.243086105691209076e-17 4.763216963682577079e-14 -1.887583244373389328e-14 -1.592926120230101003e-15 -3.264625985174926743e-16 -5.879810778647083869e-16 -1.891966634045982456e-15 2.851148476213537420e-26 1.058426796312021171e-27 1.107481097947859209e-26 7.026483273915552917e-26 -5.379740661803572684e-19 -5.402147542475392543e-19 -1.533973283924579781e-17 5.186818943221322513e-15 5.756938989258702090e-15 2.461311477771505841e-18 3.091641941796989494e-14 -1.759372847258843210e-16 1.070459744475160279e-15 9.401517695428871141e-18 1.520685182370616798e-16 1.525755060186724085e-16 8.834165471558264173e-17 4.235084795229275426e-15 -1.343201925832196203e-16 8.757869065918129794e-14 8.734305031308171818e-14 -3.234243315335766641e-14 -1.015755225261617099e-13 2.162259177191088287e-15 6.313068106737503019e-16 --8.143086031774118376e-17 5.352818153943201318e-15 -8.072235707204832032e-18 -8.026163323787587436e-17 -3.038866244634352617e-15 1.406908605920547081e-17 -4.384607976308321193e-14 9.531646394130322699e-15 -9.821249544289911861e-17 -7.923909448595848997e-16 -9.391283520432811530e-17 -2.711700981994465871e-16 1.058426796312021171e-27 8.811082299008530012e-27 1.487801391128000374e-27 9.563873128644333457e-27 1.229211106917345516e-19 9.086668192364899074e-20 -4.005920918774881884e-17 -2.847060601679682921e-15 -2.127964203851167272e-15 1.628865061860011811e-18 8.722981696312964399e-15 -1.245990632231858332e-16 -1.035329050434101653e-15 -1.089157915915628154e-18 -1.201790912326496301e-17 -1.672240764147307434e-17 2.748078604662010185e-18 -1.462373406033409279e-15 -2.706459684488763294e-17 3.979391298362956564e-15 3.914280900809954267e-15 -1.474664675128080572e-15 -4.571334409002748931e-15 -3.184500434619874255e-15 2.915127137224032614e-16 --4.589020455106308953e-16 3.945412998590303000e-14 4.295548432640169746e-18 -4.187898242504952461e-15 -1.406216877566180905e-14 -5.669532384363459182e-18 2.390016003294783340e-14 -1.395286465047703843e-14 -7.439606753073941868e-16 -2.869213732277017375e-16 -7.571720499217147661e-15 -2.204199540613405901e-15 1.107481097947859209e-26 1.487801391128000374e-27 1.670781716139887458e-25 7.303558544180336786e-26 -3.863256970748645692e-19 -3.853063443221715582e-19 -6.025552387684062132e-17 4.625905237017853671e-15 5.623605591899965859e-15 3.465867171663132686e-18 3.897465960433185356e-14 4.136114348206219717e-16 -5.577700348424316488e-17 -2.850597210977476503e-19 1.578558198322722499e-16 1.574596966400089160e-16 1.116888499524688689e-16 3.637569907625921999e-15 -1.478423120060089212e-16 3.269968615819430543e-14 3.510560327678775616e-14 -1.209817091222388069e-14 -4.090315881827249918e-14 3.317222791750073899e-15 7.690375928819282975e-16 --3.098661009272525452e-15 2.598336304972241831e-13 -3.705586445323490784e-17 -1.643837790116418740e-14 -7.731203689154381191e-14 -4.265439752377840368e-17 1.802000319076651299e-13 -1.010094527683821402e-13 -4.845003724100701487e-15 -1.943758295631656079e-15 -5.255400233547673271e-15 -5.253591335895288997e-13 7.026483273915552917e-26 9.563873128644333457e-27 7.303558544180336786e-26 2.936862964043296133e-23 -2.642810303210343532e-18 -2.643764778721222251e-18 -3.963086126803719396e-16 2.376937110121502144e-14 3.064589526524977282e-14 2.284806361195078194e-17 2.579873988473421512e-13 3.168273772135320469e-15 -3.091775265829969024e-16 -1.241280744886867784e-17 1.054145316287755313e-15 1.049161398503431286e-15 7.417407914940684902e-16 1.908816841781816436e-14 -9.268810345429511462e-16 2.059572385345710425e-13 2.235959995215321721e-13 -7.622820185235721193e-14 -2.605932125881225939e-13 2.553960731936533514e-14 5.242908346781543878e-15 -2.582109938010748615e-08 -2.151349431467241060e-06 -3.382128028653160845e-10 -2.415063062091929880e-08 -3.333970888846127249e-10 2.258055017134934322e-09 -7.553474453448608289e-06 2.222821203875844844e-06 3.813363256363073580e-08 8.266752272149200309e-09 2.856267415224482285e-08 7.794795518439994591e-08 -5.379740661803572684e-19 1.229211106917345516e-19 -3.863256970748645692e-19 -2.642810303210343532e-18 7.403427778043011779e-11 3.792255765611489853e-11 -1.752501317773718225e-08 -3.892781849899091207e-07 -2.772317334686791469e-07 -5.393331797473890970e-11 -2.529673605135954261e-06 -4.363570597273100514e-08 -1.095915722966923104e-07 9.977143842030775020e-10 -1.309301747407404281e-08 -1.353873145121602425e-08 -7.794890802066248781e-09 -1.475428397776662541e-07 6.079910763427108562e-09 -1.429439820800710253e-06 -1.674430426264918292e-06 5.296222786552825868e-07 1.952425956793901815e-06 -7.209465997261584918e-07 -4.236010055304580633e-08 -2.624612085038632830e-08 -2.167558592309546433e-06 -3.237635627080133576e-10 -3.005924930827426758e-08 -8.709526033092263374e-08 1.973339258111572039e-09 -6.879771852291717019e-06 2.537007182150678717e-06 3.857585521230416037e-08 9.778994440663878912e-09 2.904112999090483503e-08 7.846052278910172042e-08 -5.402147542475392543e-19 9.086668192364899074e-20 -3.853063443221715582e-19 -2.643764778721222251e-18 3.792255765611489853e-11 7.115769517967051503e-11 -1.816743239538513489e-08 -4.007424093499594748e-07 -2.822264494255077849e-07 -5.807776819282940936e-11 -2.854449519605377725e-06 -7.525430127036434062e-08 -1.081976407863369931e-07 1.566013896545607172e-09 -1.251205534443805253e-08 -1.287411686617860229e-08 -7.601134826475198951e-09 -1.479189936788522961e-07 6.683421504262534207e-09 -1.455166275855798905e-06 -1.690749426185897373e-06 5.391347495603226767e-07 1.971511674726411526e-06 -8.267564687813066038e-07 -4.700992181711387086e-08 -3.608749591974849653e-06 -2.808918726555411629e-04 2.880580686211860440e-07 2.969477106024997500e-07 5.768850215532460621e-05 -2.312332839147157345e-07 6.747976667041734062e-04 -6.392355604498717206e-05 3.200685394482180146e-06 3.524928150674835131e-06 4.142496374942100873e-06 1.148661919795845201e-05 -1.533973283924579781e-17 -4.005920918774881884e-17 -6.025552387684062132e-17 -3.963086126803719396e-16 -1.752501317773718225e-08 -1.816743239538513489e-08 3.264114120969943217e-05 -2.548394869425140520e-05 -2.779137482701451935e-04 -7.915047909730105148e-08 1.886790910630985283e-03 1.090754866631791889e-05 7.801050689041739985e-06 -3.432197661963476500e-06 -5.301808378832312896e-07 -4.152882865697312327e-07 -6.031727788188544236e-07 -2.128547701336814147e-04 1.234274341509008116e-06 -4.493057470519621835e-05 -7.206594708028664050e-05 1.697945549754080645e-05 8.504953299522758909e-05 1.777752907677872599e-04 2.585254330549545762e-05 --6.937855671113378833e-05 1.185104206621878306e-02 -5.630122258543856047e-05 -2.599807913875601048e-03 -1.808405967871253364e-02 -2.916663369361920704e-05 1.029989985721349172e-01 -3.473801794131257253e-02 -2.112929067660789938e-04 9.244866045590156400e-05 -2.570082280690285368e-04 -6.454285445130538692e-04 5.186818943221322513e-15 -2.847060601679682921e-15 4.625905237017853671e-15 2.376937110121502144e-14 -3.892781849899091207e-07 -4.007424093499594748e-07 -2.548394869425140520e-05 5.406641666172327942e-02 5.174314523077627043e-02 -2.119053630567342706e-08 4.233491535612741996e-02 -1.462953360525826117e-03 3.121355558429767408e-03 -4.192856834587190711e-04 1.080013838424443047e-04 1.186614724903431369e-04 5.628102500579559636e-05 3.656786376748527578e-02 -2.546513518647200729e-06 1.417633226812336188e-02 1.486108430352409832e-02 -5.230857934383274961e-03 -1.726575396760205230e-02 2.391948208788366448e-03 -1.790482763215092395e-03 --1.371705693477332281e-04 1.657883033879654217e-02 -5.550398399081168347e-05 -2.502578179730502710e-03 -1.896972119286157929e-02 -2.519716699320902141e-05 9.065041706401480437e-02 -3.257967242021848203e-02 -2.858688383419348751e-04 2.724159914307262731e-05 -3.251801223006570694e-04 -8.440052609354726587e-04 5.756938989258702090e-15 -2.127964203851167272e-15 5.623605591899965859e-15 3.064589526524977282e-14 -2.772317334686791469e-07 -2.822264494255077849e-07 -2.779137482701451935e-04 5.174314523077627043e-02 5.185024329073833343e-02 1.255540786621290009e-06 1.302520162131516959e-02 -1.618341401727035021e-03 2.953214562481817212e-03 -3.416134893996327412e-04 1.167556037089365365e-04 1.252744296441264085e-04 6.603921468351778485e-05 3.689859664766347025e-02 -2.563449933879626929e-05 1.590642365676888520e-02 1.704477207252333079e-02 -5.876392370890174723e-03 -1.982649448940359757e-02 -3.922392057424447923e-05 -1.976172892504302666e-03 --2.071765376120491374e-07 1.543514452236958691e-05 -1.368692423501269688e-08 -3.448229919827389252e-08 -3.448716959724578149e-06 6.390334241489579055e-09 -2.040043730075414874e-05 3.079233917323632859e-06 -2.408182879736639907e-07 -1.754112212171557053e-07 -2.402523089244516050e-07 -6.640530256914299064e-07 2.461311477771505841e-18 1.628865061860011811e-18 3.465867171663132686e-18 2.284806361195078194e-17 -5.393331797473890970e-11 -5.807776819282940936e-11 -7.915047909730105148e-08 -2.119053630567342706e-08 1.255540786621290009e-06 2.923778826857790215e-09 8.101263196662729328e-06 -1.850550483830866121e-07 -4.340263833203487453e-07 1.100083824610656931e-08 4.306354748750698602e-08 3.889510608936747429e-08 3.734877925264959754e-08 8.633797122601068312e-07 -6.417850997642396732e-08 7.220878837043200342e-06 8.583828081722996397e-06 -2.684398257591137100e-06 -1.003887185171714319e-05 -8.037912468433751875e-07 3.038103882007981565e-07 --2.516517250145258906e-03 2.136226422370665345e-01 -5.511266268928062019e-04 9.894412705418731673e-04 3.799553542857860011e-02 -2.181936752914516461e-05 4.725134464369972614e-02 -3.340170502487335713e-02 -2.713414985469090809e-03 -1.639613353212979532e-03 -3.216676545935040574e-03 -7.926725382679215490e-03 3.091641941796989494e-14 8.722981696312964399e-15 3.897465960433185356e-14 2.579873988473421512e-13 -2.529673605135954261e-06 -2.854449519605377725e-06 1.886790910630985283e-03 4.233491535612741996e-02 1.302520162131516959e-02 8.101263196662729328e-06 2.561886360437845678e+00 1.712688632339961700e-02 -1.495969548298849670e-02 -4.941174031154322611e-03 8.451753728352438908e-04 8.312588628538810937e-04 6.207798751111286498e-04 -1.486990425054944553e-02 -6.548694347379695625e-04 8.400416410730938976e-02 1.040558310541357112e-01 -3.123910806711079069e-02 -1.216717023960072891e-01 1.668395658113263902e-01 2.568679781949267243e-02 --3.928954203215390470e-05 2.137533177873366264e-03 5.445097707514176965e-06 1.614015226545799700e-04 1.185191496548233239e-03 -2.885741513929382642e-06 1.027426075386230972e-02 -5.494265105720975290e-03 -2.657470988356218109e-05 -1.838827713483549386e-05 -3.961294658191956365e-05 -1.011002018751854773e-04 -1.759372847258843210e-16 -1.245990632231858332e-16 4.136114348206219717e-16 3.168273772135320469e-15 -4.363570597273100514e-08 -7.525430127036434062e-08 1.090754866631791889e-05 -1.462953360525826117e-03 -1.618341401727035021e-03 -1.850550483830866121e-07 1.712688632339961700e-02 1.114660806185970585e-03 -7.730335285552777892e-04 -4.452308713267074861e-06 1.481836476482716150e-05 1.451846216162652682e-05 8.452527659270013381e-06 -1.359203485518305915e-03 -1.601010014264435757e-06 -1.002409728823324721e-03 -3.485597911598559660e-04 3.639571529618867399e-04 3.878510497305322842e-04 1.604008248187792666e-03 1.834410652820026033e-04 -1.279490468007478382e-06 2.287660214046896914e-03 1.881736127791208387e-06 1.096321509004083357e-05 1.182857261486174993e-03 -9.584663865549323316e-06 3.051814272148788490e-02 -5.915567986028016229e-03 -2.803719763090280236e-05 4.881136721097782943e-05 2.107728167369185196e-06 9.225664148879019028e-06 1.070459744475160279e-15 -1.035329050434101653e-15 -5.577700348424316488e-17 -3.091775265829969024e-16 -1.095915722966923104e-07 -1.081976407863369931e-07 7.801050689041739985e-06 3.121355558429767408e-03 2.953214562481817212e-03 -4.340263833203487453e-07 -1.495969548298849670e-02 -7.730335285552777892e-04 1.688998948499471605e-03 1.159026578533110570e-05 2.989596973369903929e-05 3.373840088359871794e-05 1.465319842208428671e-05 1.862022016073073931e-03 -5.124694057679437498e-06 2.798040993920551233e-03 2.876093839470280755e-03 -1.029387046857675172e-03 -3.333290952544525499e-03 3.469864107285604033e-03 -1.855656174130421361e-04 -9.068426317865086325e-08 2.979745879017003073e-04 1.074217982496532016e-06 7.773663567167847920e-07 5.547822027831921226e-05 -5.294172336870067481e-09 1.825216372649425331e-05 3.383003929474542214e-05 2.580429643797891438e-07 1.113170450201948483e-06 6.525664614377412270e-07 1.051179776515564498e-06 9.401517695428871141e-18 -1.089157915915628154e-18 -2.850597210977476503e-19 -1.241280744886867784e-17 9.977143842030775020e-10 1.566013896545607172e-09 -3.432197661963476500e-06 -4.192856834587190711e-04 -3.416134893996327412e-04 1.100083824610656931e-08 -4.941174031154322611e-03 -4.452308713267074861e-06 1.159026578533110570e-05 2.103990647271761701e-05 1.114955092000448355e-06 1.206700317689306529e-06 9.328667613641622459e-07 -2.012489084703237706e-04 -1.281890104442271536e-06 3.024086810582569051e-05 2.773027802293759217e-05 -1.096665998980926882e-05 -3.201338373551938023e-05 -2.511310858789703626e-04 -4.049192332247160327e-05 --1.010411452316174498e-05 1.676699508613893001e-03 -9.573321138379633936e-07 6.743863193988200014e-06 1.003402092902857769e-04 -5.146313275530546443e-07 1.755644760324465278e-03 -5.605255102331018402e-04 -1.010510257238622547e-05 -2.085181637611447406e-06 -9.965538983237963089e-06 -2.918386730194719129e-05 1.520685182370616798e-16 -1.201790912326496301e-17 1.578558198322722499e-16 1.054145316287755313e-15 -1.309301747407404281e-08 -1.251205534443805253e-08 -5.301808378832312896e-07 1.080013838424443047e-04 1.167556037089365365e-04 4.306354748750698602e-08 8.451753728352438908e-04 1.481836476482716150e-05 2.989596973369903929e-05 1.114955092000448355e-06 2.884322043042034708e-04 -2.538126012564094220e-04 5.387122904190521982e-06 6.894907238685718954e-05 -2.703157629667282969e-06 2.554210195701856097e-03 4.889853800134472996e-04 -9.511869428156902347e-04 -5.708735812038205929e-04 1.938946049089669489e-04 1.649637945841375946e-05 --1.002915587787911641e-05 1.732916157635252218e-03 -8.955180902632183681e-07 6.281162398372516905e-06 1.077696970302018857e-04 -5.539881244564768786e-07 1.870870624454052007e-03 -5.695318433638057315e-04 -9.686357516880338036e-06 -1.662825435048065551e-06 -9.808018743190203645e-06 -2.889439043922104893e-05 1.525755060186724085e-16 -1.672240764147307434e-17 1.574596966400089160e-16 1.049161398503431286e-15 -1.353873145121602425e-08 -1.287411686617860229e-08 -4.152882865697312327e-07 1.186614724903431369e-04 1.252744296441264085e-04 3.889510608936747429e-08 8.312588628538810937e-04 1.451846216162652682e-05 3.373840088359871794e-05 1.206700317689306529e-06 -2.538126012564094220e-04 3.102166925762899776e-04 5.611063091380308048e-06 7.541673326503943767e-05 -2.661931171835422640e-06 1.470998787372507915e-02 4.955205309483954376e-04 -5.437021090677873411e-03 -5.783272781411675492e-04 1.993184776291691035e-04 1.588876821035091380e-05 --7.096411239012776800e-06 1.298896908940160239e-03 -8.519837545226716281e-07 4.516749168953107673e-06 8.883519413119686178e-05 -2.523146163460764631e-07 8.708104448430653436e-04 -2.820483943677134196e-04 -5.814282625649671241e-06 -1.529648396480678433e-06 -6.833063562445468623e-06 -2.028183930935848183e-05 8.834165471558264173e-17 2.748078604662010185e-18 1.116888499524688689e-16 7.417407914940684902e-16 -7.794890802066248781e-09 -7.601134826475198951e-09 -6.031727788188544236e-07 5.628102500579559636e-05 6.603921468351778485e-05 3.734877925264959754e-08 6.207798751111286498e-04 8.452527659270013381e-06 1.465319842208428671e-05 9.328667613641622459e-07 5.387122904190521982e-06 5.611063091380308048e-06 1.988825337167201494e-05 3.812443652156480859e-05 -1.944066131252356749e-06 2.255075602087547727e-04 9.370667484578088652e-03 -8.357229684017524840e-05 -1.091244912444888637e-02 1.067404421836058503e-04 1.284101720431881188e-05 --6.166408270828315896e-05 9.282782630378589445e-03 -3.229239994559254898e-05 -1.923508018533116436e-03 -1.588843399738775614e-02 -1.594898464460605111e-05 5.720248449212836372e-02 -1.955720914959569248e-02 -1.795567582710501570e-04 3.754504312323681382e-05 -1.870407675260386544e-04 -5.055290379226907082e-04 4.235084795229275426e-15 -1.462373406033409279e-15 3.637569907625921999e-15 1.908816841781816436e-14 -1.475428397776662541e-07 -1.479189936788522961e-07 -2.128547701336814147e-04 3.656786376748527578e-02 3.689859664766347025e-02 8.633797122601068312e-07 -1.486990425054944553e-02 -1.359203485518305915e-03 1.862022016073073931e-03 -2.012489084703237706e-04 6.894907238685718954e-05 7.541673326503943767e-05 3.812443652156480859e-05 2.724309422537568928e-02 -7.792580236021764118e-06 1.214317040165824486e-02 1.235390285211402685e-02 -4.479740771749052969e-03 -1.435268803652801212e-02 -6.216048243951985408e-03 -1.671552675289353493e-03 -7.954518997506867072e-06 -6.763338260323719427e-04 3.022377409523564049e-07 1.331924023050925137e-05 3.661169141778293729e-05 8.444660670564603750e-08 -3.219229529387673773e-04 2.528818449588239310e-04 1.072466086869132534e-05 5.070866914767413610e-06 1.107990032257984879e-05 2.773567462221827859e-05 -1.343201925832196203e-16 -2.706459684488763294e-17 -1.478423120060089212e-16 -9.268810345429511462e-16 6.079910763427108562e-09 6.683421504262534207e-09 1.234274341509008116e-06 -2.546513518647200729e-06 -2.563449933879626929e-05 -6.417850997642396732e-08 -6.548694347379695625e-04 -1.601010014264435757e-06 -5.124694057679437498e-06 -1.281890104442271536e-06 -2.703157629667282969e-06 -2.661931171835422640e-06 -1.944066131252356749e-06 -7.792580236021764118e-06 3.938618265045845181e-06 -3.810800813135566405e-04 -4.397826624549711397e-04 1.413248991076312738e-04 5.133429895970311298e-04 -8.965737950838059862e-05 -2.129537338315999433e-05 --1.457372693470361177e-03 9.451541906012153604e-02 1.005144316480618646e-04 -5.849646589347031259e-03 -1.294007493555653365e-01 -2.640464914727091680e-05 1.057903669450520889e-01 -4.322324541853273666e-02 -4.802813362752815160e-03 -9.756958387685016505e-04 -1.679299878847527938e-03 -5.498729953080172886e-03 8.757869065918129794e-14 3.979391298362956564e-15 3.269968615819430543e-14 2.059572385345710425e-13 -1.429439820800710253e-06 -1.455166275855798905e-06 -4.493057470519621835e-05 1.417633226812336188e-02 1.590642365676888520e-02 7.220878837043200342e-06 8.400416410730938976e-02 -1.002409728823324721e-03 2.798040993920551233e-03 3.024086810582569051e-05 2.554210195701856097e-03 1.470998787372507915e-02 2.255075602087547727e-04 1.214317040165824486e-02 -3.810800813135566405e-04 7.124318197334406300e+00 2.676044558680029106e-01 -2.632743994364230122e+00 -3.111554357622127176e-01 1.764057209661643380e-03 1.756704866057536658e-03 --1.665472713324716442e-03 1.143881756135775030e-01 6.291873262212429902e-05 -5.178846973454946830e-03 -1.216020000364273906e-01 -3.561425506840147387e-05 1.403403673864646317e-01 -6.023684334790740946e-02 -4.963194473939220967e-03 -1.072946360077371318e-03 -1.899267739233110111e-03 -6.050663216414568953e-03 8.734305031308171818e-14 3.914280900809954267e-15 3.510560327678775616e-14 2.235959995215321721e-13 -1.674430426264918292e-06 -1.690749426185897373e-06 -7.206594708028664050e-05 1.486108430352409832e-02 1.704477207252333079e-02 8.583828081722996397e-06 1.040558310541357112e-01 -3.485597911598559660e-04 2.876093839470280755e-03 2.773027802293759217e-05 4.889853800134472996e-04 4.955205309483954376e-04 9.370667484578088652e-03 1.235390285211402685e-02 -4.397826624549711397e-04 2.676044558680029106e-01 7.686030370760271957e+00 -9.884182661136756143e-02 -8.940093905890782011e+00 7.354504365195378701e-03 2.178830498571339251e-03 -5.402267720757495664e-04 -3.500456106541598861e-02 -3.646627150018448262e-05 2.153268188959647257e-03 4.772874293545138064e-02 9.788159208161594391e-06 -3.924802079810447292e-02 1.611303878009478704e-02 1.775867457136094804e-03 3.617995149859560026e-04 6.223804919963958409e-04 2.036131827085406725e-03 -3.234243315335766641e-14 -1.474664675128080572e-15 -1.209817091222388069e-14 -7.622820185235721193e-14 5.296222786552825868e-07 5.391347495603226767e-07 1.697945549754080645e-05 -5.230857934383274961e-03 -5.876392370890174723e-03 -2.684398257591137100e-06 -3.123910806711079069e-02 3.639571529618867399e-04 -1.029387046857675172e-03 -1.096665998980926882e-05 -9.511869428156902347e-04 -5.437021090677873411e-03 -8.357229684017524840e-05 -4.479740771749052969e-03 1.413248991076312738e-04 -2.632743994364230122e+00 -9.884182661136756143e-02 9.729171339957245479e-01 1.149288680405661928e-01 -6.978562179520457964e-04 -6.534298927597297166e-04 -1.943321742220476944e-03 -1.335569570692768659e-01 -7.162025230758143908e-05 6.003639016901835684e-03 1.412309683815566597e-01 4.150511446075701752e-05 -1.636564359123428292e-01 7.044389440078994979e-02 5.777840977318744395e-03 1.251930415222073795e-03 2.215661331552164909e-03 7.054155462547743484e-03 -1.015755225261617099e-13 -4.571334409002748931e-15 -4.090315881827249918e-14 -2.605932125881225939e-13 1.952425956793901815e-06 1.971511674726411526e-06 8.504953299522758909e-05 -1.726575396760205230e-02 -1.982649448940359757e-02 -1.003887185171714319e-05 -1.216717023960072891e-01 3.878510497305322842e-04 -3.333290952544525499e-03 -3.201338373551938023e-05 -5.708735812038205929e-04 -5.783272781411675492e-04 -1.091244912444888637e-02 -1.435268803652801212e-02 5.133429895970311298e-04 -3.111554357622127176e-01 -8.940093905890782011e+00 1.149288680405661928e-01 1.039881525384822325e+01 -8.678148783162666774e-03 -2.548835131946574056e-03 --3.333863796457967564e-04 2.491333162046089483e-02 -2.353051652724181582e-05 1.387439896431741414e-03 1.866643744552765397e-02 -4.390083089488345825e-05 1.451439494535149144e-01 -4.483966495114385248e-02 -3.405953240386222296e-04 -5.528075146326992567e-05 -3.730348743630668195e-04 -8.578412202787384210e-04 2.162259177191088287e-15 -3.184500434619874255e-15 3.317222791750073899e-15 2.553960731936533514e-14 -7.209465997261584918e-07 -8.267564687813066038e-07 1.777752907677872599e-04 2.391948208788366448e-03 -3.922392057424447923e-05 -8.037912468433751875e-07 1.668395658113263902e-01 1.604008248187792666e-03 3.469864107285604033e-03 -2.511310858789703626e-04 1.938946049089669489e-04 1.993184776291691035e-04 1.067404421836058503e-04 -6.216048243951985408e-03 -8.965737950838059862e-05 1.764057209661643380e-03 7.354504365195378701e-03 -6.978562179520457964e-04 -8.678148783162666774e-03 4.943448466562492399e-02 1.944024799375950992e-03 --5.468710926687100291e-05 4.670416326687474098e-03 -4.818973351892917084e-06 8.672656714053429282e-05 1.065831122547543450e-03 4.229752159922561635e-07 -1.738338909140980243e-03 -8.233919291224370057e-05 -6.056307379456799997e-05 -3.912271590726792503e-05 -6.364872634283788350e-05 -1.606945507961290737e-04 6.313068106737503019e-16 2.915127137224032614e-16 7.690375928819282975e-16 5.242908346781543878e-15 -4.236010055304580633e-08 -4.700992181711387086e-08 2.585254330549545762e-05 -1.790482763215092395e-03 -1.976172892504302666e-03 3.038103882007981565e-07 2.568679781949267243e-02 1.834410652820026033e-04 -1.855656174130421361e-04 -4.049192332247160327e-05 1.649637945841375946e-05 1.588876821035091380e-05 1.284101720431881188e-05 -1.671552675289353493e-03 -2.129537338315999433e-05 1.756704866057536658e-03 2.178830498571339251e-03 -6.534298927597297166e-04 -2.548835131946574056e-03 1.944024799375950992e-03 4.524222784337420585e-04 +2.010435472675886910e-09 -2.675858521978017979e-08 -5.982534266809095175e-09 -2.698030618056596629e-08 -1.989809076874735238e-12 1.269479487876332009e-10 -5.379883140916528167e-07 4.411431339826491038e-08 8.460884350230103591e-11 -8.235579874535305463e-11 5.086303196611991347e-09 3.909799987933974186e-09 -6.270337829424909449e-07 -1.560734009111028709e-09 -1.798068771436684880e-08 9.741659069758279789e-08 1.052901538387951524e-07 -1.403434090874938674e-08 2.955423936136483677e-08 2.952140785935902445e-08 1.085332727849522123e-07 -7.711179591641636079e-08 -1.846199834485905805e-08 -2.095910725117839080e-07 1.212295268261887986e-06 -1.145028511636006787e-06 3.003489585249692202e-05 -5.548432473991033794e-08 1.163159055337209065e-06 +-2.675858521978017979e-08 1.297766063027828726e-03 -3.753669648603687755e-06 2.204948804053052913e-06 -5.760648977086614835e-09 -2.522364024973069784e-07 2.810351414140960452e-04 2.059453949034324586e-04 -9.688275684226148794e-08 1.277079316073990768e-07 -1.700435637905916910e-05 -1.042783194256338995e-05 2.508894046795525257e-04 -4.763427030648742539e-06 2.445678126714358067e-05 -7.768766350936209429e-05 -1.550000746516891001e-04 1.914835953902359177e-06 -1.941649531453191137e-05 -1.919320822886435486e-05 -1.142340243970952425e-04 8.986459785266957064e-05 4.187063139912689347e-06 2.217299814292995278e-04 -1.385430643066194257e-03 1.188125263385224816e-03 -3.449393420040661185e-02 6.806284211416293445e-05 -1.302196147287084044e-03 +-5.982534266809095175e-09 -3.753669648603687755e-06 5.467147788330904279e-06 -1.929446420726715245e-06 3.927915024328045706e-10 9.502851339213074019e-09 -5.731547160736014102e-05 1.495392876565641426e-05 1.071688151308961080e-09 -1.725434567365214470e-10 -1.191868887531181009e-07 -6.836605684076856572e-07 -7.510404330818832544e-05 2.803210685748507891e-07 -1.201099205241394705e-07 -4.782705624218493033e-06 2.068811267157731403e-06 8.597575746850255660e-07 3.463248154055069549e-06 3.532975495996603757e-06 -1.806614884970461895e-05 4.122148316867475904e-07 1.024767440553556745e-06 -5.997920740230977469e-06 -1.845114834708861242e-04 -2.235097988923882958e-05 -4.561999760523118182e-03 2.359450383216521243e-06 -1.862173397419128484e-04 +-2.698030618056596629e-08 2.204948804053052913e-06 -1.929446420726715245e-06 8.877944663226193446e-06 -6.887209628178520836e-10 -1.521223767400015378e-08 1.981431653728706405e-05 -5.714232511650175025e-06 -6.681930418968987725e-09 7.852074010094775080e-09 -7.419927206390406279e-07 -1.003827887550708074e-06 -2.231961288826657386e-05 -3.477150456054405671e-09 1.579290263136871634e-06 -4.648523862594739099e-06 -1.494115171799464481e-05 -1.549645407580595582e-07 1.009957944766114598e-05 1.009326221029824826e-05 -3.588087127642131948e-05 9.272090196410051968e-07 1.021197707334585114e-06 -4.214069500957989611e-05 -3.201145467449226429e-04 -2.307141301817855993e-04 -7.959695123680991802e-03 -5.109658235769883260e-06 -3.475119476432363434e-04 +-1.989809076874735238e-12 -5.760648977086614835e-09 3.927915024328045706e-10 -6.887209628178520836e-10 1.047452371966215647e-12 3.692024030853954120e-12 -2.642804787881784725e-08 4.836244240397750710e-11 -3.225026153927942103e-12 3.372891561747340235e-12 1.687275531821482708e-10 9.001017484816815393e-11 -2.200711982739165082e-08 4.064637270132787020e-11 6.897085277391199193e-10 -3.361049878509963552e-09 -6.557371007649930479e-09 1.801201598159322895e-10 5.036460463942421980e-09 5.048235098404066286e-09 -1.877149415559848179e-08 5.321780894420643688e-09 3.789746867844042224e-10 -2.393386342270401727e-08 -1.599889653366557293e-07 -1.278965889517170458e-07 -3.965513049679856371e-06 -3.785093963734028598e-09 -1.777795032857846728e-07 +1.269479487876332009e-10 -2.522364024973069784e-07 9.502851339213074019e-09 -1.521223767400015378e-08 3.692024030853954120e-12 1.072618923392085248e-08 -2.444250766281064126e-06 -5.556165944175527011e-07 1.079737229591668807e-09 -8.740472329937089759e-10 -5.750858863490159575e-09 4.740987606888703317e-09 2.710816793204982275e-07 3.034175945723352065e-09 -1.999805120433519395e-07 3.764263573571957757e-08 4.181422798691923542e-07 2.166939634446426924e-08 1.434159600630858443e-06 1.436526785418002922e-06 -4.968901227956546984e-06 1.464073563438657625e-06 -8.013697727984940637e-09 -8.021146145905913756e-06 -3.896282249162749667e-05 -4.293640546041120756e-05 -9.640936081544036404e-04 -1.602449575179677763e-06 -4.532299598556670966e-05 +-5.379883140916528167e-07 2.810351414140960452e-04 -5.731547160736014102e-05 1.981431653728706405e-05 -2.642804787881784725e-08 -2.444250766281064126e-06 1.011100176087781528e-02 -3.151729389364358279e-05 -2.210012465133209189e-07 1.430018693056412249e-07 4.727030150893541632e-07 -1.185956556815480327e-05 8.444014896910213971e-04 -8.068582347158779444e-06 3.816302426117755352e-05 -9.509134335701366694e-05 1.290869371424355292e-04 7.999173587881918432e-06 -4.448902829659454115e-04 -4.462615080624281628e-04 5.673066235628386857e-04 -4.368463339528391894e-04 -2.966651419748082100e-06 2.862775861467920131e-03 2.258599654960728037e-03 1.526924499023254306e-02 5.495298627308806139e-02 6.789276878488010973e-04 4.111724444826134889e-03 +4.411431339826491038e-08 2.059453949034324586e-04 1.495392876565641426e-05 -5.714232511650175025e-06 4.836244240397750710e-11 -5.556165944175527011e-07 -3.151729389364358279e-05 8.287635177929500571e-04 -6.743212170134985155e-08 6.485737507432679050e-08 -2.599202544389353039e-06 -3.722761252836061124e-06 -2.202427013627296360e-04 -1.041551004665647798e-07 1.373706124727826798e-05 -2.665527447599296088e-05 -3.483906546259526782e-05 1.181917809097914446e-06 -7.334379604406815261e-05 -7.323818323045451018e-05 2.143945294306314636e-04 -7.481292600847001882e-05 3.559172142580735643e-06 4.659570549436209668e-04 1.488902473850350948e-03 2.523325674687419955e-03 3.682187393506722306e-02 1.071414507704355822e-04 1.857178593453299867e-03 +8.460884350230103591e-11 -9.688275684226148794e-08 1.071688151308961080e-09 -6.681930418968987725e-09 -3.225026153927942103e-12 1.079737229591668807e-09 -2.210012465133209189e-07 -6.743212170134985155e-08 3.152305700313159847e-10 -2.881072600764250739e-10 -5.302398322695695769e-09 -2.817358729753059608e-10 3.819805565475815367e-07 1.770111263217535693e-10 -6.277292817441172622e-08 1.116993027908334727e-07 4.094922280581006152e-07 -8.719975487233701527e-09 9.719937009222355882e-08 9.729488244261514172e-08 -4.779256711338634577e-07 -2.828391647494108599e-07 -1.644794619850250746e-08 -1.706431221052968567e-07 -4.942929899398460255e-06 -9.189580895912764518e-07 -1.229887005227123296e-04 6.618784561824867690e-08 -4.954344793958554478e-06 +-8.235579874535305463e-11 1.277079316073990768e-07 -1.725434567365214470e-10 7.852074010094775080e-09 3.372891561747340235e-12 -8.740472329937089759e-10 1.430018693056412249e-07 6.485737507432679050e-08 -2.881072600764250739e-10 2.851104178984615446e-10 1.571379112124710568e-09 1.686919797292217741e-10 -3.980102346758749090e-07 7.092916960025225752e-11 5.981704433536636329e-08 -1.147320558819750008e-07 -4.081810500505028473e-07 6.279873646548775085e-09 -5.135915139002784187e-09 -5.055020797946611428e-09 2.413037662907759548e-07 3.251551332746034311e-07 1.757670970201573751e-08 -3.658977718676276286e-07 3.245280664429195227e-06 -1.961105585744226006e-06 8.100837084989195424e-05 -1.820466889746353856e-07 2.868363605021105458e-06 +5.086303196611991347e-09 -1.700435637905916910e-05 -1.191868887531181009e-07 -7.419927206390406279e-07 1.687275531821482708e-10 -5.750858863490159575e-09 4.727030150893541632e-07 -2.599202544389353039e-06 -5.302398322695695769e-09 1.571379112124710568e-09 2.100012909662544023e-06 1.067500372544160458e-06 -7.633985400239611885e-05 5.166906854238038822e-08 6.508316636571919618e-07 -9.292632869571105494e-07 -5.437192199874677119e-06 1.135794444260546867e-06 -1.704018825963783847e-05 -1.707420044972905330e-05 4.591104251553039073e-05 6.093201733908097876e-06 1.539966974219876571e-08 8.786098849513336519e-05 3.645306257668421606e-04 4.726146912977818751e-04 9.037471824443470136e-03 1.636191618321751969e-05 4.227424575262562329e-04 +3.909799987933974186e-09 -1.042783194256338995e-05 -6.836605684076856572e-07 -1.003827887550708074e-06 9.001017484816815393e-11 4.740987606888703317e-09 -1.185956556815480327e-05 -3.722761252836061124e-06 -2.817358729753059608e-10 1.686919797292217741e-10 1.067500372544160458e-06 1.815188306379810778e-06 -3.662673368674924463e-05 -5.786456868502590224e-08 2.524888949741236429e-09 1.956461517809685446e-06 -1.610876398027382503e-06 -3.899415111099695314e-07 5.184952793856188186e-06 5.184681251328582246e-06 -9.881016518236210300e-07 6.199989911789667079e-06 -3.264123439185272380e-07 -4.480773228192790525e-05 5.505439686708714516e-05 -2.418200127553242886e-04 1.386268920405101706e-03 -1.331874509802442755e-05 2.193157305120750018e-05 +-6.270337829424909449e-07 2.508894046795525257e-04 -7.510404330818832544e-05 -2.231961288826657386e-05 -2.200711982739165082e-08 2.710816793204982275e-07 8.444014896910213971e-04 -2.202427013627296360e-04 3.819805565475815367e-07 -3.980102346758749090e-07 -7.633985400239611885e-05 -3.662673368674924463e-05 1.131082039421567667e-02 -2.270365598285138139e-05 -8.009363587642599758e-05 3.537792687626881323e-04 5.796572242303914399e-04 -5.093985890945189787e-05 4.678521391953274300e-05 4.595000142320852793e-05 -4.568614174309945848e-04 -8.092956277858830887e-04 -7.586722897012989502e-05 2.927571402887697394e-04 -5.491162368178550100e-03 1.413323656988731396e-03 -1.373143971535892016e-01 2.020372487663332339e-04 -5.125345129111239829e-03 +-1.560734009111028709e-09 -4.763427030648742539e-06 2.803210685748507891e-07 -3.477150456054405671e-09 4.064637270132787020e-11 3.034175945723352065e-09 -8.068582347158779444e-06 -1.041551004665647798e-07 1.770111263217535693e-10 7.092916960025225752e-11 5.166906854238038822e-08 -5.786456868502590224e-08 -2.270365598285138139e-05 1.786518111915964769e-07 -1.795335503922821274e-08 -1.242832365659986865e-06 8.000880296570959303e-07 4.923476882534784764e-08 1.216299141774865806e-06 1.220813387381810254e-06 -1.363990744249371660e-06 1.709234856541364061e-06 1.326762686642541518e-07 -8.824707697377046390e-06 -1.794062539235657335e-06 -4.686681143955169504e-05 -3.945769632847092261e-05 -2.311039907057228941e-06 -8.031008125252086740e-06 +-1.798068771436684880e-08 2.445678126714358067e-05 -1.201099205241394705e-07 1.579290263136871634e-06 6.897085277391199193e-10 -1.999805120433519395e-07 3.816302426117755352e-05 1.373706124727826798e-05 -6.277292817441172622e-08 5.981704433536636329e-08 6.508316636571919618e-07 2.524888949741236429e-09 -8.009363587642599758e-05 -1.795335503922821274e-08 1.363952312190774481e-05 -2.379720837818551153e-05 -8.552223718639338161e-05 1.514011491432705191e-06 -1.102518759229970123e-05 -1.102710239131722966e-05 7.661880580219207948e-05 5.913335059477412340e-05 3.556625851706962321e-06 -1.826666068860154636e-05 8.665422349241999785e-04 -9.741364555841180506e-05 2.158640297574592956e-02 -2.541793048099771652e-05 8.305114282270360023e-04 +9.741659069758279789e-08 -7.768766350936209429e-05 -4.782705624218493033e-06 -4.648523862594739099e-06 -3.361049878509963552e-09 3.764263573571957757e-08 -9.509134335701366694e-05 -2.665527447599296088e-05 1.116993027908334727e-07 -1.147320558819750008e-07 -9.292632869571105494e-07 1.956461517809685446e-06 3.537792687626881323e-04 -1.242832365659986865e-06 -2.379720837818551153e-05 2.686516295581304128e-04 4.185067909142546833e-05 -2.198463577136753164e-06 -1.547470386314452658e-04 -1.550696415335507314e-04 7.235432082383764065e-04 -1.175995195589292014e-04 1.526213108331391269e-06 5.786982774660976185e-04 6.659649989105266565e-03 3.104039763264382096e-03 1.653187072030989913e-01 4.507516723601237617e-05 7.099630861624886893e-03 +1.052901538387951524e-07 -1.550000746516891001e-04 2.068811267157731403e-06 -1.494115171799464481e-05 -6.557371007649930479e-09 4.181422798691923542e-07 1.290869371424355292e-04 -3.483906546259526782e-05 4.094922280581006152e-07 -4.081810500505028473e-07 -5.437192199874677119e-06 -1.610876398027382503e-06 5.796572242303914399e-04 8.000880296570959303e-07 -8.552223718639338161e-05 4.185067909142546833e-05 9.688932706007270150e-04 -2.806879032812250230e-05 -3.183691490652695551e-04 -3.189553665914461079e-04 1.138733396792313262e-03 -5.267418098943810016e-04 -5.091181231887897219e-05 1.763193751609949981e-03 8.995680333825462754e-03 9.471806342709358437e-03 2.227107039551090295e-01 3.609785149330512486e-04 1.044038738942925856e-02 +-1.403434090874938674e-08 1.914835953902359177e-06 8.597575746850255660e-07 -1.549645407580595582e-07 1.801201598159322895e-10 2.166939634446426924e-08 7.999173587881918432e-06 1.181917809097914446e-06 -8.719975487233701527e-09 6.279873646548775085e-09 1.135794444260546867e-06 -3.899415111099695314e-07 -5.093985890945189787e-05 4.923476882534784764e-08 1.514011491432705191e-06 -2.198463577136753164e-06 -2.806879032812250230e-05 1.056362804027514429e-05 -8.481037287228928819e-05 -8.494581159938223056e-05 1.948251320711523059e-04 1.838112797928504610e-05 5.446716556367443957e-06 4.705899113929160418e-04 1.420879765753534909e-03 2.534269964554149244e-03 3.519397139681489789e-02 9.681913184239717069e-05 1.731885439059069332e-03 +2.955423936136483677e-08 -1.941649531453191137e-05 3.463248154055069549e-06 1.009957944766114598e-05 5.036460463942421980e-09 1.434159600630858443e-06 -4.448902829659454115e-04 -7.334379604406815261e-05 9.719937009222355882e-08 -5.135915139002784187e-09 -1.704018825963783847e-05 5.184952793856188186e-06 4.678521391953274300e-05 1.216299141774865806e-06 -1.102518759229970123e-05 -1.547470386314452658e-04 -3.183691490652695551e-04 -8.481037287228928819e-05 6.666883798965164698e-02 6.687990379768993243e-02 1.314713929339253809e-02 3.427719800653128836e-04 6.872189828933031337e-06 -3.352345831071660287e-01 1.067304579002174708e-01 -1.774863683059179964e+00 2.644773702711320063e+00 -5.644851970222779236e-02 1.219901072952075949e-01 +2.952140785935902445e-08 -1.919320822886435486e-05 3.532975495996603757e-06 1.009326221029824826e-05 5.048235098404066286e-09 1.436526785418002922e-06 -4.462615080624281628e-04 -7.323818323045451018e-05 9.729488244261514172e-08 -5.055020797946611428e-09 -1.707420044972905330e-05 5.184681251328582246e-06 4.595000142320852793e-05 1.220813387381810254e-06 -1.102710239131722966e-05 -1.550696415335507314e-04 -3.189553665914461079e-04 -8.494581159938223056e-05 6.687990379768993243e-02 6.717871920148572873e-02 1.317982343604531294e-02 3.437685283381013428e-04 6.900064781338489876e-06 -3.362381769801792708e-01 1.069884849424199424e-01 -1.780018234303760716e+00 2.651156477985307625e+00 -5.653018152053759793e-02 1.222892999845945627e-01 +1.085332727849522123e-07 -1.142340243970952425e-04 -1.806614884970461895e-05 -3.588087127642131948e-05 -1.877149415559848179e-08 -4.968901227956546984e-06 5.673066235628386857e-04 2.143945294306314636e-04 -4.779256711338634577e-07 2.413037662907759548e-07 4.591104251553039073e-05 -9.881016518236210300e-07 -4.568614174309945848e-04 -1.363990744249371660e-06 7.661880580219207948e-05 7.235432082383764065e-04 1.138733396792313262e-03 1.948251320711523059e-04 1.314713929339253809e-02 1.317982343604531294e-02 6.286708702720613662e-01 1.754819185974190818e-03 -1.799179814623763689e-05 -3.846793815786232751e-02 4.883765006164704658e+00 -2.068968073406703134e-01 1.205201353828864939e+02 1.248526762342380431e-03 5.699847013450031419e+00 +-7.711179591641636079e-08 8.986459785266957064e-05 4.122148316867475904e-07 9.272090196410051968e-07 5.321780894420643688e-09 1.464073563438657625e-06 -4.368463339528391894e-04 -7.481292600847001882e-05 -2.828391647494108599e-07 3.251551332746034311e-07 6.093201733908097876e-06 6.199989911789667079e-06 -8.092956277858830887e-04 1.709234856541364061e-06 5.913335059477412340e-05 -1.175995195589292014e-04 -5.267418098943810016e-04 1.838112797928504610e-05 3.427719800653128836e-04 3.437685283381013428e-04 1.754819185974190818e-03 2.044499206723828129e-03 1.720417715470083791e-05 -5.864847084455113331e-03 2.800494224958707812e-02 -3.139570880628861360e-02 7.001535418483920692e-01 -2.229557505069342752e-03 2.307389488616622131e-02 +-1.846199834485905805e-08 4.187063139912689347e-06 1.024767440553556745e-06 1.021197707334585114e-06 3.789746867844042224e-10 -8.013697727984940637e-09 -2.966651419748082100e-06 3.559172142580735643e-06 -1.644794619850250746e-08 1.757670970201573751e-08 1.539966974219876571e-08 -3.264123439185272380e-07 -7.586722897012989502e-05 1.326762686642541518e-07 3.556625851706962321e-06 1.526213108331391269e-06 -5.091181231887897219e-05 5.446716556367443957e-06 6.872189828933031337e-06 6.900064781338489876e-06 -1.799179814623763689e-05 1.720417715470083791e-05 1.736683263133104906e-05 -5.604502271433321135e-05 -8.579744319185587346e-05 -2.982939558778339019e-04 -2.086918431747334171e-03 -1.611499032759180841e-05 -1.368653615970390799e-04 +-2.095910725117839080e-07 2.217299814292995278e-04 -5.997920740230977469e-06 -4.214069500957989611e-05 -2.393386342270401727e-08 -8.021146145905913756e-06 2.862775861467920131e-03 4.659570549436209668e-04 -1.706431221052968567e-07 -3.658977718676276286e-07 8.786098849513336519e-05 -4.480773228192790525e-05 2.927571402887697394e-04 -8.824707697377046390e-06 -1.826666068860154636e-05 5.786982774660976185e-04 1.763193751609949981e-03 4.705899113929160418e-04 -3.352345831071660287e-01 -3.362381769801792708e-01 -3.846793815786232751e-02 -5.864847084455113331e-03 -5.604502271433321135e-05 2.130627813421966721e+00 -1.740538503939484116e-01 1.132457368452565838e+01 -4.251881958073468581e+00 4.922087356490528753e-01 -2.888579202727960538e-01 +1.212295268261887986e-06 -1.385430643066194257e-03 -1.845114834708861242e-04 -3.201145467449226429e-04 -1.599889653366557293e-07 -3.896282249162749667e-05 2.258599654960728037e-03 1.488902473850350948e-03 -4.942929899398460255e-06 3.245280664429195227e-06 3.645306257668421606e-04 5.505439686708714516e-05 -5.491162368178550100e-03 -1.794062539235657335e-06 8.665422349241999785e-04 6.659649989105266565e-03 8.995680333825462754e-03 1.420879765753534909e-03 1.067304579002174708e-01 1.069884849424199424e-01 4.883765006164704658e+00 2.800494224958707812e-02 -8.579744319185587346e-05 -1.740538503939484116e-01 4.403454632902271015e+01 -9.394405695447265447e-01 1.089543930390909054e+03 7.470531431610480388e-02 4.728604203977826614e+01 +-1.145028511636006787e-06 1.188125263385224816e-03 -2.235097988923882958e-05 -2.307141301817855993e-04 -1.278965889517170458e-07 -4.293640546041120756e-05 1.526924499023254306e-02 2.523325674687419955e-03 -9.189580895912764518e-07 -1.961105585744226006e-06 4.726146912977818751e-04 -2.418200127553242886e-04 1.413323656988731396e-03 -4.686681143955169504e-05 -9.741364555841180506e-05 3.104039763264382096e-03 9.471806342709358437e-03 2.534269964554149244e-03 -1.774863683059179964e+00 -1.780018234303760716e+00 -2.068968073406703134e-01 -3.139570880628861360e-02 -2.982939558778339019e-04 1.132457368452565838e+01 -9.394405695447265447e-01 6.020468364136843320e+01 -2.295222916969847304e+01 2.627256703884690570e+00 -1.555259352129249706e+00 +3.003489585249692202e-05 -3.449393420040661185e-02 -4.561999760523118182e-03 -7.959695123680991802e-03 -3.965513049679856371e-06 -9.640936081544036404e-04 5.495298627308806139e-02 3.682187393506722306e-02 -1.229887005227123296e-04 8.100837084989195424e-05 9.037471824443470136e-03 1.386268920405101706e-03 -1.373143971535892016e-01 -3.945769632847092261e-05 2.158640297574592956e-02 1.653187072030989913e-01 2.227107039551090295e-01 3.519397139681489789e-02 2.644773702711320063e+00 2.651156477985307625e+00 1.205201353828864939e+02 7.001535418483920692e-01 -2.086918431747334171e-03 -4.251881958073468581e+00 1.089543930390909054e+03 -2.295222916969847304e+01 2.696033386691932174e+04 1.879690320406714799e+00 1.168352198688235831e+03 +-5.548432473991033794e-08 6.806284211416293445e-05 2.359450383216521243e-06 -5.109658235769883260e-06 -3.785093963734028598e-09 -1.602449575179677763e-06 6.789276878488010973e-04 1.071414507704355822e-04 6.618784561824867690e-08 -1.820466889746353856e-07 1.636191618321751969e-05 -1.331874509802442755e-05 2.020372487663332339e-04 -2.311039907057228941e-06 -2.541793048099771652e-05 4.507516723601237617e-05 3.609785149330512486e-04 9.681913184239717069e-05 -5.644851970222779236e-02 -5.653018152053759793e-02 1.248526762342380431e-03 -2.229557505069342752e-03 -1.611499032759180841e-05 4.922087356490528753e-01 7.470531431610480388e-02 2.627256703884690570e+00 1.879690320406714799e+00 1.454678562368016126e-01 4.338069887155174276e-02 +1.163159055337209065e-06 -1.302196147287084044e-03 -1.862173397419128484e-04 -3.475119476432363434e-04 -1.777795032857846728e-07 -4.532299598556670966e-05 4.111724444826134889e-03 1.857178593453299867e-03 -4.954344793958554478e-06 2.868363605021105458e-06 4.227424575262562329e-04 2.193157305120750018e-05 -5.125345129111239829e-03 -8.031008125252086740e-06 8.305114282270360023e-04 7.099630861624886893e-03 1.044038738942925856e-02 1.731885439059069332e-03 1.219901072952075949e-01 1.222892999845945627e-01 5.699847013450031419e+00 2.307389488616622131e-02 -1.368653615970390799e-04 -2.888579202727960538e-01 4.728604203977826614e+01 -1.555259352129249706e+00 1.168352198688235831e+03 4.338069887155174276e-02 5.316251407340103441e+01 diff --git a/contrib/anisotropic_stishovite/stishovite_data.py b/contrib/anisotropic_stishovite/stishovite_data.py index f48ba4a7..ea4d0dd7 100644 --- a/contrib/anisotropic_stishovite/stishovite_data.py +++ b/contrib/anisotropic_stishovite/stishovite_data.py @@ -1,5 +1,9 @@ import numpy as np from burnman.utils.unitcell import molar_volume_from_unit_cell_volume +from copy import copy +from burnman.calibrants.tools import pressure_to_pressure +from burnman import calibrants +from scipy.integrate import cumulative_simpson def add_to_data(data, publication, compilation): @@ -11,12 +15,73 @@ def add_to_data(data, publication, compilation): } -def get_data(): +Au_Fei = calibrants.Fei_2007.Au() +Pt_Fei = calibrants.Fei_2007.Pt() +Pt_Holmes = calibrants.Holmes_1989.Pt() +Au_Anderson = calibrants.Anderson_1989.Au() +Au_Tsuchiya = calibrants.Tsuchiya_2003.Au() +Pt_Dorogokupets = calibrants.Dorogokupets_2007.Pt() +Pt_Zha = calibrants.Zha_2008.Pt() +Au_Matsui = calibrants.Matsui_2010.Au() +Pt_Matsui = calibrants.Matsui_2009.Pt() +Au_Dorogokupets = calibrants.Dorogokupets_2007.Au() + +calibrant_conversions_none = {} + +calibrant_conversions_Dorogokupets = { + "Andrault_2003": [Pt_Holmes, Pt_Dorogokupets], + "Murakami_2003": [Pt_Holmes, Pt_Dorogokupets], + "Nishihara_2005": [Au_Anderson, Au_Dorogokupets], + "Wang_2012": [Au_Tsuchiya, Au_Dorogokupets], + "Grocholski_2013": [Au_Tsuchiya, Au_Dorogokupets], + "Fischer_2018": [Pt_Dorogokupets, Pt_Dorogokupets], + "Zhang_2021": [Au_Fei, Au_Dorogokupets], + "Zhang_2023": [Au_Fei, Au_Dorogokupets], +} + +calibrant_conversions_Holmes_Tsuchiya = { + "Andrault_2003": [Pt_Holmes, Pt_Holmes], + "Murakami_2003": [Pt_Holmes, Pt_Holmes], + "Nishihara_2005": [Au_Anderson, Au_Tsuchiya], + "Wang_2012": [Au_Tsuchiya, Au_Tsuchiya], + "Grocholski_2013": [Au_Tsuchiya, Au_Tsuchiya], + "Fischer_2018": [Pt_Dorogokupets, Pt_Holmes], + "Zhang_2021": [Au_Fei, Au_Tsuchiya], + "Zhang_2023": [Au_Fei, Au_Tsuchiya], +} + +calibrant_conversions_Matsui = { + "Andrault_2003": [Pt_Holmes, Pt_Matsui], + "Murakami_2003": [Pt_Holmes, Pt_Matsui], + "Nishihara_2005": [Au_Anderson, Au_Matsui], + "Wang_2012": [Au_Tsuchiya, Au_Matsui], + "Grocholski_2013": [Au_Tsuchiya, Au_Matsui], + "Fischer_2018": [Pt_Dorogokupets, Pt_Matsui], + "Zhang_2021": [Au_Fei, Au_Matsui], + "Zhang_2023": [Au_Fei, Au_Matsui], +} + +calibrant_conversions_Fei = { + "Andrault_2003": [Pt_Holmes, Pt_Fei], + "Murakami_2003": [Pt_Holmes, Pt_Fei], + "Nishihara_2005": [Au_Anderson, Au_Fei], + "Wang_2012": [Au_Tsuchiya, Au_Fei], + "Grocholski_2013": [Au_Tsuchiya, Au_Fei], + "Fischer_2018": [Pt_Dorogokupets, Pt_Fei], + "Zhang_2021": [Au_Fei, Au_Fei], + "Zhang_2023": [Au_Fei, Au_Fei], +} + +calibrant_conversions = calibrant_conversions_Fei + + +def get_data(unify_calibrants=True): data = {} - # Zhang data + # Zhang et al., 2021 elastic data CT_data = np.loadtxt("data/Zhang_2021_stishovite_elastic_tensor.dat") - unit_cell_data = np.loadtxt("data/Zhang_2021_stishovite_unit_cell.dat") + idx = np.argsort(CT_data[:, 0]) + CT_data = CT_data[idx] CT_header = [ "P", "P_err", @@ -41,6 +106,11 @@ def get_data(): "rho", "rho_err", ] + + # Zhang et al., 2021 XRD data + unit_cell_data = np.loadtxt("data/Zhang_2021_stishovite_unit_cell.dat") + idx = np.argsort(unit_cell_data[:, 0]) + unit_cell_data = unit_cell_data[idx] unit_cell_header = [ "P", "P_err", @@ -55,7 +125,6 @@ def get_data(): "rho", "rho_err", ] - isstv = unit_cell_data[:, 2] == unit_cell_data[:, 4] compilation = [ @@ -74,7 +143,11 @@ def get_data(): # Andrault data stv_data = np.loadtxt("data/Andrault_2003_stishovite_unit_cell.dat") + idx = np.argsort(stv_data[:, 0]) + stv_data = stv_data[idx] poststv_data = np.loadtxt("data/Andrault_2003_post_stishovite_unit_cell.dat") + idx = np.argsort(poststv_data[:, 0]) + poststv_data = poststv_data[idx] stv_header = ["P", "P_err", "a", "a_err", "c", "c_err", "V", "V_err"] poststv_header = [ "P", @@ -100,12 +173,121 @@ def get_data(): d["T"] = 0.0 * d["P"] + 298.15 d["T_err"] = 0.0 * d["P"] + 5.0 - # Fischer data (swap a and b) + # Murakami data + poststv_data = np.loadtxt("data/Murakami_2003_post_stishovite_unit_cell.dat") + idx = np.argsort(poststv_data[:, 0]) + poststv_data = poststv_data[idx] + poststv_header = [ + "runid", + "P", + "a", + "a_err", + "b", + "b_err", + "c", + "c_err", + "V", + "V_err", + ] + compilation = [["poststv", poststv_header, poststv_data]] + add_to_data(data, "Murakami_2003", compilation) + d = data["Murakami_2003"]["poststv"] + d["T"] = 0.0 * d["P"] + 298.15 + d["T_err"] = 0.0 * d["P"] + 5.0 + + d["V"] = d["V"] * 1.0e-6 / molar_volume_from_unit_cell_volume(1.0, 2.0) + d["V_err"] = d["V_err"] * 1.0e-6 / molar_volume_from_unit_cell_volume(1.0, 2.0) + d["P"] = d["P"] + d["P_err"] = d["P"] * 0.01 + + f_ln = np.cbrt(molar_volume_from_unit_cell_volume(1.0, 2.0)) + for ln in ["a", "b", "c"]: + d[f"{ln}"] = d[f"{ln}"] # * f_ln + d[f"{ln}_err"] = d[f"{ln}_err"] # * f_ln + + # Grocholski data + all_data = np.loadtxt( + "data/Grocholski_2013_stishovite_post_stishovite_unit_cell.dat" + ) + idx = np.argsort(all_data[:, 0]) + all_data = all_data[idx] + stv_header = [ + "P", + "P_err", + "a_p", + "a_p_err", + "a", + "a_err", + "b", + "b_err", + "c", + "c_err", + "V", + "V_err", + ] + stv_mask = all_data[:, 0] < 50.0 + compilation = [ + ["stv", stv_header, copy(all_data[stv_mask])], + ["poststv", stv_header, copy(all_data[~stv_mask])], + ] + add_to_data(data, "Grocholski_2013", compilation) + for d in [data["Grocholski_2013"]["stv"], data["Grocholski_2013"]["poststv"]]: + d["T"] = 0.0 * d["P"] + 298.15 + d["T_err"] = 0.0 * d["P"] + 5.0 + + # Volumes are reported as molar volumes, here we change to unit cell volumes + # We correct back later, with all the other datasets + d["V"] = d["V"] * 1.0e-6 / molar_volume_from_unit_cell_volume(1.0, 2.0) + d["V_err"] = d["V_err"] * 1.0e-6 / molar_volume_from_unit_cell_volume(1.0, 2.0) + + # Zhang data + all_data = np.loadtxt("data/Zhang_2023_stishovite_unit_cell.dat") + idx = np.argsort(all_data[:, 0]) + all_data = all_data[idx] + stv_header = [ + "P", + "P_err", + "a", + "a_err", + "b", + "b_err", + "c", + "c_err", + "V", + "V_err", + "unique_refl", + "Rint", + "R1", + ] + stv_mask = all_data[:, 0] < 50.0 + compilation = [ + ["stv", stv_header, copy(all_data[stv_mask])], + ["poststv", stv_header, copy(all_data[~stv_mask])], + ] + add_to_data(data, "Zhang_2023", compilation) + for d in [data["Zhang_2023"]["stv"], data["Zhang_2023"]["poststv"]]: + d["T"] = 0.0 * d["P"] + 298.15 + d["T_err"] = 0.0 * d["P"] + 5.0 + + # Sun et al., 2019 data + S_data = np.loadtxt("data/Sun_et_al_2019_post_stishovite_volumes.dat") + idx = np.argsort(S_data[:, 0]) + S_data = S_data[idx] + poststv_header = ["P", "P_err", "T", "T_err", "V", "V_err", "V_NaCl", "V_NaCl_err"] + compilation = [ + ["poststv", poststv_header, S_data], + ] + add_to_data(data, "Sun_2019", compilation) + + # Fischer et al., 2018 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:] + idx = np.argsort(F_data[:, 0]) + F_data = F_data[idx] + header = [ "P", "P_err", @@ -150,8 +332,10 @@ def get_data(): d["P"] = d["T"] * 0.0 + 0.0001 d["P_err"] = d["T"] * 0.0 + 0.000001 - # Wang data (only stishovite) + # Wang et al., 2012 data (only stishovite) W_data = np.loadtxt("data/Wang_et_al_2012_stv.dat") + idx = np.argsort(W_data[:, 2]) + W_data = W_data[idx] header = [ "id", "T", @@ -174,8 +358,11 @@ def get_data(): d["b"] = d["a"] d["b_err"] = d["a_err"] - # Nishihara data (only stishovite) + # Nishihara et al., 2005 data (only stishovite) N_data = np.loadtxt("data/Nishihara_et_al_2005_stv.dat") + idx = np.argsort(N_data[:, 0]) + N_data = N_data[idx] + header = ["P", "T", "a", "a_err", "c", "c_err", "V", "V_err"] compilation = [["stv", header, N_data]] add_to_data(data, "Nishihara_2005", compilation) @@ -193,27 +380,83 @@ def get_data(): 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"], - data["Wang_2012"]["stv"], - data["Nishihara_2005"]["stv"], + # Check that P and V are in the correct units + # Unify calibrants if desired + for study in [ + "Andrault_2003", + "Murakami_2003", + "Nishihara_2005", + "Wang_2012", + "Grocholski_2013", + "Fischer_2018", + "Sun_2019", + "Zhang_2021", + "Zhang_2023", ]: + + for phase in ["stv", "poststv"]: + try: + P = data[study][phase]["P"] + T = data[study][phase]["T"] + V = data[study][phase]["V"] + if max(P) > 1.0e6: + print("P too small", study) + exit() + if min(V) < 1.0: + print("V too big", study) + exit() + if unify_calibrants and study in calibrant_conversions: + print(f"Changing calibrant for {study} data ({phase})") + P = [ + pressure_to_pressure( + calibrant_conversions[study][0], + calibrant_conversions[study][1], + P[i] * 1.0e9, + T[i], + ) + for i in range(len(P)) + ] + data[study][phase]["P"] = np.array(P) / 1.0e9 + + except KeyError: + pass + + for study, phase in [ + ("Ito_1974", "stv"), + ("Andrault_2003", "stv"), + ("Andrault_2003", "poststv"), + ("Murakami_2003", "poststv"), + ("Nishihara_2005", "stv"), + ("Wang_2012", "stv"), + ("Grocholski_2013", "stv"), + ("Grocholski_2013", "poststv"), + ("Fischer_2018", "stv"), + ("Fischer_2018", "poststv"), + ("Sun_2019", "poststv"), + ("Zhang_2021", "stv"), + ("Zhang_2021", "poststv"), + ("Zhang_2023", "stv"), + ("Zhang_2023", "poststv"), + ]: + d = data[study][phase] + 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 + try: + for ln in ["a", "b", "c"]: + d[f"{ln}"] = d[f"{ln}"] * f_ln + d[f"{ln}_err"] = d[f"{ln}_err"] * f_ln + except KeyError as e: + if study == "Sun_2019": + pass + else: + print(study, phase) + print(e) - # Molar mass (M) is apparently rounded by Zhang to get rho. + # Molar mass (M) is apparently rounded by Zhang_2021 to get rho. # As rho is a derived value, we correct rho here and # return the derived volumes M_Zhang = 0.0601 @@ -252,6 +495,7 @@ def common_data(): d["abc"] = {"stv": {}, "poststv": {}} d["abc_err"] = {"stv": {}, "poststv": {}} + #### TODO Add Zhang 2023 here? for pub in [ "Ito_1974", "Andrault_2003", @@ -269,7 +513,7 @@ def common_data(): d["PTV_err"][phase][pub] = np.array( [Z_data["P_err"], Z_data["T_err"], Z_data["V_err"]] ).T - except: + except KeyError: pass try: @@ -280,7 +524,47 @@ def common_data(): d["abc_err"][phase][pub] = np.array( [Z_data["a_err"], Z_data["b_err"], Z_data["c_err"]] ).T - except: + except KeyError: + pass + return d + + +def other_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 [ + "Murakami_2003", + "Grocholski_2013", + "Zhang_2023", + ]: + 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 KeyError: + 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 KeyError: pass return d @@ -302,10 +586,14 @@ def get_stiffnesses(CN): CN_data = get_data()["Zhang_2021"]["CT"] CN_GPa, CN_err_GPa, KNR_GPa = get_stiffnesses(CN_data) +CN_err_GPa[CN_err_GPa < 1.0] = 1.0 + +KTR_GPa = KNR_GPa 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) +P_for_CN_new = cumulative_simpson(KTR_GPa * 1.0e9, x=-np.log(V_for_CN), initial=0) # 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 @@ -314,3 +602,48 @@ def get_stiffnesses(CN): 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) + +# Jiang data +J2009_data = np.loadtxt("data/Jiang_2009_elastic.dat", skiprows=1) +J2009_stiffness_matrices = np.zeros((10, 6, 6)) +J2009_pressures = J2009_data[:, 0] * 1.0e9 +J2009_temperatures = J2009_pressures * 0.0 + 298.15 +for i in range(10): + # C11, C33, C13, C12, C44, C66 + J2009_stiffness_matrices[i] = [ + [J2009_data[i, 1], J2009_data[i, 4], J2009_data[i, 3], 0, 0, 0], + [J2009_data[i, 4], J2009_data[i, 1], J2009_data[i, 3], 0, 0, 0], + [J2009_data[i, 3], J2009_data[i, 3], J2009_data[i, 2], 0, 0, 0], + [0, 0, 0, J2009_data[i, 5], 0, 0], + [0, 0, 0, 0, J2009_data[i, 5], 0], + [0, 0, 0, 0, 0, J2009_data[i, 6]], + ] +J2009_stiffness_matrices = J2009_stiffness_matrices * 1.0e9 +J2009_compliance_matrices = np.linalg.inv(J2009_stiffness_matrices) +J2009_KRT = 1.0 / np.sum(J2009_compliance_matrices[:, :3, :3], axis=(1, 2)) +J2009_dPsidfs = np.einsum("i, ijk->ijk", J2009_KRT, J2009_compliance_matrices) + +C2024 = np.loadtxt("data/Chen_2024_stishovite_velocities.dat") +C2024_pressures = C2024[:, 0] * 1.0e9 +C2024_pressures_unc = C2024[:, 1] * 1.0e9 +C2024_temperatures = C2024[:, 2] + 273.15 +C2024_temperatures_unc = C2024[:, 3] + 273.15 +C2024_Vp = C2024[:, 4] * 1.0e3 +C2024_Vp_unc = C2024[:, 5] * 1.0e3 +C2024_Vs = C2024[:, 6] * 1.0e3 +C2024_Vs_unc = C2024[:, 7] * 1.0e3 +C2024_Vphi = ( + np.sqrt(np.power(C2024[:, 4], 2.0) - 4.0 / 3.0 * np.power(C2024[:, 6], 2.0)) * 1.0e3 +) +C2024_Vphi_unc = ( + (1 / C2024_Vphi) + * np.sqrt( + (C2024[:, 4] * C2024[:, 5]) ** 2 + + ((4.0 / 3.0) * C2024[:, 6] * C2024[:, 7]) ** 2 + ) + * 1.0e6 +) +C2024_KS = C2024[:, 10] * 1.0e9 +C2024_KS_unc = C2024[:, 11] * 1.0e9 +C2024_G = C2024[:, 12] * 1.0e9 +C2024_G_unc = C2024[:, 13] * 1.0e9 diff --git a/contrib/anisotropic_stishovite/stishovite_fit_eos.py b/contrib/anisotropic_stishovite/stishovite_fit_eos.py index 0fcc55c8..77cc4b9c 100644 --- a/contrib/anisotropic_stishovite/stishovite_fit_eos.py +++ b/contrib/anisotropic_stishovite/stishovite_fit_eos.py @@ -2,20 +2,29 @@ from __future__ import print_function import numpy as np -from scipy.optimize import minimize, differential_evolution, leastsq +from scipy.optimize import minimize, differential_evolution, leastsq, least_squares from stishovite_data import ( common_data, - P_for_CN, + P_for_CN_new, V_for_CN, T_for_CN, KNR_GPa, KNR_err_GPa, CN_GPa, CN_err_GPa, + J2009_stiffness_matrices, + J2009_pressures, + J2009_temperatures, + C2024_pressures, + C2024_temperatures, + C2024_KS, + C2024_KS_unc, + C2024_G, + C2024_G_unc, ) from scipy.optimize import root_scalar -from stishovite_model import make_models, make_scalar_model, modify_Zhang_elasticity +from stishovite_model import make_models, make_scalar_model from copy import deepcopy from burnman import RelaxedSolution import matplotlib.pyplot as plt @@ -62,14 +71,6 @@ def misfit_scalar(args): V0Q1overV0Q0, dDebye_0, P_tr_GPa, - fP_Zhang, - fP_Andrault, - fP_Wang, - fP_Nishihara, - fP2_Zhang, - fP2_Andrault, - fP2_Wang, - fP2_Nishihara, ) = args scalar_prms = make_scalar_model( @@ -88,36 +89,24 @@ def misfit_scalar(args): ) relaxed_stv.set_composition([1.0]) + misfits = [] # 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 misfits = [(P_tr_3000K_model - P_tr_3000K_obs) / P_tr_3000K_err] - # Fit volumes for stishovite from Zhang, Andrault + # Fit volumes for stishovite from Nishihara, Wang, Andrault for phase, publication in [ ("stv", "Nishihara_2005"), ("stv", "Wang_2012"), - ("stv", "Zhang_2021"), - ("poststv", "Zhang_2021"), + # ("stv", "Zhang_2021"), + # ("poststv", "Zhang_2021"), ("stv", "Andrault_2003"), ("poststv", "Andrault_2003"), + # ("stv", "Fischer_2018"), + # ("poststv", "Fischer_2018"), ]: - if publication == "Zhang_2021": - fP = fP_Zhang - fP2 = fP2_Zhang - elif publication == "Andrault_2003": - fP = fP_Andrault - fP2 = fP2_Andrault - elif publication == "Wang_2012": - fP = fP_Wang - fP2 = fP2_Wang - elif publication == "Nishihara_2005": - fP = fP_Nishihara - fP2 = fP2_Nishihara - else: - fP = 1.0 - fP2 = 0.0 PTV = deepcopy(data["PTV"][phase][publication]) PTV_err = deepcopy(data["PTV_err"][phase][publication]) @@ -125,7 +114,7 @@ def misfit_scalar(args): P_model = relaxed_stv.evaluate_with_volumes(["pressure"], PTV[:, 2], PTV[:, 1])[ 0 ] - P_actual = PTV[:, 0] * (fP + fP2 * PTV[:, 0]) + P_actual = PTV[:, 0] misfits.extend(list((P_model - P_actual) / PTV_err[:, 0])) if plot: @@ -148,14 +137,14 @@ def misfit_scalar(args): ax[1].scatter(PTV[:, 1][sort], PTV[:, 2][sort]) # Fit isentropic bulk moduli from Zhang - KNR_obs = relaxed_stv.evaluate_with_volumes( + KNR_mod = relaxed_stv.evaluate_with_volumes( ["isentropic_bulk_modulus_reuss"], V_for_CN, T_for_CN )[0] - misfits.extend(list((KNR_GPa - KNR_obs / 1.0e9) / KNR_err_GPa)) + misfits.extend(list((KNR_GPa - KNR_mod / 1.0e9) / KNR_err_GPa)) if plot: sort = np.argsort(V_for_CN) - ax[2].plot(V_for_CN[sort], KNR_obs[sort] / 1.0e9) + ax[2].plot(V_for_CN[sort], KNR_mod[sort] / 1.0e9) ax[2].scatter(V_for_CN[sort], KNR_GPa[sort]) plt.show() @@ -179,36 +168,17 @@ def misfit_cell(args, scalar_args): V0Q1overV0Q0, dDebye_0, P_tr_GPa, - fP_Zhang, - fP_Andrault, - fP_Wang, - fP_Nishihara, - fP2_Zhang, - fP2_Andrault, - fP2_Wang, - fP2_Nishihara, ) = scalar_args - (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = args - ( - a11, - a22, - a33, - a44, - a55, - a66, - b11, - b22, - b33, - b44, - b55, - b66, - c44, - c55, - c66, - b112, - b222, - b332, - ) = np.zeros(18) + (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_d, f_PsiI_22) = args + + (a11, a22, a33, a44, a55, a66, b11, b33, b44, b66, d44, d66) = np.zeros(12) + + c44 = 1.0 + c66 = 1.0 + b22 = b11 + b55 = b44 + c55 = c44 + d55 = d44 models = make_models( dVQ0, @@ -224,8 +194,8 @@ def misfit_cell(args, scalar_args): PsiI_33_a, PsiI_33_b, PsiI_33_c, - PsiI_33_b2, - PsiI_33_c2, + PsiI_33_d, + f_PsiI_22, a11, a22, a33, @@ -241,9 +211,9 @@ def misfit_cell(args, scalar_args): c44, c55, c66, - b112, - b222, - b332, + d44, + d55, + d66, ) _, _, _, stishovite_relaxed = models stishovite_relaxed.set_composition([1.0]) @@ -261,21 +231,6 @@ def misfit_cell(args, scalar_args): # ("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 - elif publication == "Wang_2012": - fP = fP_Wang - fP2 = fP2_Wang - elif publication == "Nishihara_2005": - fP = fP_Nishihara - fP2 = fP2_Nishihara - else: - fP = 1.0 - fP2 = 0.0 PTV = data["PTV"][phase][publication] abc_obs = data["abc"][phase][publication] @@ -287,7 +242,7 @@ def misfit_cell(args, scalar_args): axis=0, ) - P_actual = PTV[:, 0] * (fP + fP2 * PTV[:, 0]) + P_actual = PTV[:, 0] cell_parameters_model = stishovite_relaxed.evaluate( ["cell_parameters"], P_actual, PTV[:, 1] )[0] @@ -341,35 +296,14 @@ def misfit_elastic(args, scalar_args, cell_args): V0Q1overV0Q0, dDebye_0, P_tr_GPa, - fP_Zhang, - fP_Andrault, - fP_Wang, - fP_Nishihara, - fP2_Zhang, - fP2_Andrault, - fP2_Wang, - fP2_Nishihara, ) = scalar_args - (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = cell_args - (a11, a22, a33, a44, a55, a66, b11i, b33i, b44i, b66i, c44, c66, b112i, b332i) = ( - args - ) + (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_d, f_PsiI_22) = cell_args + (a11, a22, a33, a44, a55, a66, b11, b33, b44, b66, c44, c66, d44, d66) = args - b55i = b44i - b22i = b11i - b222i = b112i + b22 = b11 + b55 = b44 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) + d55 = d44 models = make_models( dVQ0, @@ -385,8 +319,8 @@ def misfit_elastic(args, scalar_args, cell_args): PsiI_33_a, PsiI_33_b, PsiI_33_c, - PsiI_33_b2, - PsiI_33_c2, + PsiI_33_d, + f_PsiI_22, a11, a22, a33, @@ -402,42 +336,60 @@ def misfit_elastic(args, scalar_args, cell_args): c44, c55, c66, - b112, - b222, - b332, + d44, + d55, + d66, ) - _, _, _, stishovite_relaxed = models + _, _, stishovite_unrelaxed, stishovite_relaxed = models stishovite_relaxed.set_composition([1.0]) misfits = [] + # Try to make b33 positive + # CN_model = stishovite_relaxed.evaluate( + # ["isentropic_stiffness_tensor"], [0., 10.e9, 20.e9], [298, 298, 298])[0] + # C33 = CN_model[:, 2, 2]/10.e9 + + # d2C33dP2 = (C33[2] - C33[1]) - (C33[1] - C33[0]) + # if d2C33dP2 < 0.: + # misfits.append(400.*d2C33dP2*d2C33dP2) + + # b33 = stishovite_relaxed.endmembers[0][0].anisotropic_params["b"][2][2] + # b11 = stishovite_relaxed.endmembers[0][0].anisotropic_params["b"][0][0] + + # misfits.append(400.*np.power(b33-b11, 2.)) + 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) + # Fit compliance data for Zhang + # Part I: Fit data under 25 GPa + mask = P_for_CN_new < 25.0e9 + P_actual = P_for_CN_new[mask] CN_model = stishovite_relaxed.evaluate( - ["isentropic_stiffness_tensor"], P_actual, T_for_CN + ["isentropic_stiffness_tensor"], P_actual, T_for_CN[mask] )[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] - misfits.extend(list(chis)) - if plot: + chis = (CN_GPa[mask] - CN_model / 1.0e9) / CN_err_GPa[mask] + chis[:, 2, 2] = chis[:, 2, 2] * 0.0 # don't fit C33 + chis = np.ravel(np.triu(chis)) + misfits.extend(list(chis)) + + if plot: + for i, j in [ + (0, 0), + (1, 1), + (2, 2), + (0, 1), + (0, 2), + (1, 2), + (3, 3), + (4, 4), + (5, 5), + ]: + axi = 2 if i < 2 and j < 2: axi = 0 @@ -460,6 +412,58 @@ def misfit_elastic(args, scalar_args, cell_args): ax[axi].legend() plt.show() + P_actual = P_for_CN_new[mask] + CN_model = stishovite_relaxed.evaluate( + ["isentropic_stiffness_tensor"], J2009_pressures, J2009_temperatures + )[0] + chis = (J2009_stiffness_matrices - CN_model) / 10.0e9 + chis = np.ravel(np.triu(chis)) + misfits.extend(list(chis)) + + if True: + # Part 2: Fit splitting data (>50 GPa) + chi_scale = 4.0 + split_mask = CN_GPa[:, 0, 0] != CN_GPa[:, 1, 1] + P_actual = P_for_CN_new[split_mask] + CN_model = stishovite_relaxed.evaluate( + ["isentropic_stiffness_tensor"], P_actual, T_for_CN[split_mask] + )[0] + sort = np.argsort(P_actual) + + # (11 and 22), (13 and 23), and (44 and 55) split + for i1, j1, i2, j2 in [(0, 0, 1, 1), (0, 2, 1, 2), (3, 3, 4, 4)]: + obs_split = CN_GPa[split_mask, i1, j1] - CN_GPa[split_mask, i2, j2] + unc_split = np.sqrt( + np.power(CN_err_GPa[split_mask, i1, j1], 2.0) + + np.power(CN_err_GPa[split_mask, i2, j2], 2.0) + ) + model_split = CN_model[:, i1, j1] / 1.0e9 - CN_model[:, i2, j2] / 1.0e9 + chis = (obs_split - model_split) / unc_split * chi_scale + misfits.extend(list(chis)) + + # Value of 44, 55 and 13, 23 + if True: + for i, j in [(3, 3), (4, 4), (0, 2), (1, 2)]: + obs = CN_GPa[split_mask, i, j] + unc = CN_err_GPa[split_mask, i, j] + model = CN_model[:, i, j] / 1.0e9 + chis = (obs - model) / unc + misfits.extend(list(chis)) + + # Chen 2024 data + K_mod, G_mod = stishovite_relaxed.evaluate( + ["isentropic_bulk_modulus_vrh", "isentropic_shear_modulus_vrh"], + C2024_pressures, + C2024_temperatures, + ) + misfits.extend(list((C2024_KS - K_mod) / C2024_KS_unc / 2.0)) + misfits.extend(list((C2024_G - G_mod) / C2024_G_unc / 2.0)) + # Finally, impose a weak penalty on large 0.5*(C11 + C22) - C12 + # stishovite_unrelaxed.set_composition([0.5, 0.5]) + # CN = stishovite_unrelaxed.evaluate(["isentropic_stiffness_tensor"], + # 100.e9, 298.15)[0] + # misfits.extend([(CN[0, 1] - 0.5*(CN[0, 0] + CN[1, 1])) / 100.e9]) + misfits = np.array(misfits) chisqr = np.sum(np.power(misfits, 2.0)) @@ -471,8 +475,8 @@ def misfit_elastic(args, scalar_args, cell_args): def misfit_scalar_and_cell(args): - scalar_args = args[:16] - cell_args = args[16:] + scalar_args = args[:8] + cell_args = args[8:] misfits = np.concatenate( (misfit_scalar(scalar_args), misfit_cell(cell_args, scalar_args)) ) @@ -489,7 +493,6 @@ def misfit_cell_and_elastic(args, scalar_args): cell_args = args[:7] elastic_args = args[7:] misfits_1 = misfit_cell(cell_args, scalar_args) - modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) misfits_2 = misfit_elastic(elastic_args, scalar_args, cell_args) misfits = np.concatenate((misfits_1, misfits_2)) @@ -503,12 +506,11 @@ def misfit_cell_and_elastic(args, scalar_args): def misfit_all(args): - scalar_args = args[:16] - cell_args = args[16:23] - elastic_args = args[23:] + scalar_args = args[:8] + cell_args = args[8:15] + elastic_args = args[15:] misfits_1 = misfit_scalar(scalar_args) misfits_2 = misfit_cell(cell_args, scalar_args) - modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) misfits_3 = misfit_elastic(elastic_args, scalar_args, cell_args) misfits = np.concatenate((misfits_1, misfits_2, misfits_3)) @@ -516,9 +518,9 @@ def misfit_all(args): if chisqr < min_misfit_combined[0]: min_misfit_combined[0] = chisqr - print(repr(args[:16])) - print(repr(args[16:23])) - print(repr(args[23:])) + print(repr(args[:8])) + print(repr(args[8:15])) + print(repr(args[15:])) print(chisqr) return misfits @@ -585,7 +587,6 @@ def lsq(args): 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, @@ -595,8 +596,35 @@ def lsq(args): ) if False: + + def misfit(all_args): + misfits = misfit_all(all_args) + return np.sum(np.power(misfits, 2.0)) + + sol = minimize( + misfit, + all_args, + bounds=scalar_bounds + cell_bounds + elastic_bounds, + # method="Nelder-Mead", + method="Powell", + ) + + sol = minimize( + misfit, + sol.x, + bounds=scalar_bounds + cell_bounds + elastic_bounds, + method="Nelder-Mead", + # method="Powell", + ) + + if False: + + def misfit(elastic_args, cell_args, scalar_args): + misfits = misfit_elastic(elastic_args, cell_args, scalar_args) + return np.sum(np.power(misfits, 2.0)) + sol = minimize( - misfit_elastic, + misfit, elastic_args, bounds=elastic_bounds, args=(scalar_args, cell_args), @@ -604,8 +632,13 @@ def lsq(args): ) if False: + + def misfit(cell_args, scalar_args): + misfits = misfit_cell(cell_args, scalar_args) + return np.sum(np.power(misfits, 2.0)) + sol = minimize( - misfit_cell, + misfit, cell_args, bounds=cell_bounds, args=(scalar_args), diff --git a/contrib/anisotropic_stishovite/stishovite_model.py b/contrib/anisotropic_stishovite/stishovite_model.py index 9a289208..a6e85e39 100644 --- a/contrib/anisotropic_stishovite/stishovite_model.py +++ b/contrib/anisotropic_stishovite/stishovite_model.py @@ -4,134 +4,12 @@ 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, - fP_Wang, - fP_Nishihara, - fP2_Zhang, - fP2_Andrault, - fP2_Wang, - fP2_Nishihara, - ) = scalar_args - (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = cell_args - (a11, a22, a33, a44, a55, a66, b11i, b33i, b44i, b66i, c44, c66, b112i, b332i) = ( - 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) @@ -233,8 +111,8 @@ def make_models( PsiI_33_a, PsiI_33_b, PsiI_33_c, - PsiI_33_b2, - PsiI_33_c2, + PsiI_33_d, + f_PsiI_22, a11, a22, a33, @@ -250,9 +128,9 @@ def make_models( c44, c55, c66, - b112, - b222, - b332, + d44, + d55, + d66, ): scalar_prms = make_scalar_model( @@ -267,19 +145,13 @@ def make_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) - ) + a = params["a"] + b = params["b"] + c = params["c"] + d = params["d"] + + dPsidf = a + b * np.tanh(c * f + d) + Psi = a * f + b / c * np.log(np.cosh(c * f + d) / np.cosh(d)) dPsidPth = np.zeros((6, 6)) return (Psi, dPsidf, dPsidPth) @@ -299,12 +171,17 @@ def psi_func(f, Pth, params): 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)), + "b": np.zeros((6, 6)), + "c": np.ones((6, 6)) * PsiI_33_c, + "d": np.ones((6, 6)) * PsiI_33_d, } + anisotropic_parameters["c"][3][3] = c44 + anisotropic_parameters["c"][4][4] = c55 + anisotropic_parameters["c"][5][5] = c66 + anisotropic_parameters["d"][3][3] = d44 + anisotropic_parameters["d"][4][4] = d55 + anisotropic_parameters["d"][5][5] = d66 """ beta_T/beta_RT = d(Psi*I)/df F_ij = exp(Psi I) @@ -340,55 +217,34 @@ def psi_func(f, Pth, params): # 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, + "b": PsiI_33_b, } # The following assumes that the compression of the a and - # b axes is the same + # b axes potentially differs 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, + "a": (1.0 - f_PsiI_22) * (1.0 - PsiI_33_a), + "b": -0.5 * PsiI_33_b, } 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, + "a": f_PsiI_22 * (1.0 - PsiI_33_a), + "b": -0.5 * PsiI_33_b, } 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["b"][0, 0] = b11 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["b"][1, 1] = b22 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["b"][2, 2] = b33 anisotropic_parameters["a"][3, 3] = a44 - anisotropic_parameters["b1"][3, 3] = b44 - anisotropic_parameters["c1"][3, 3] = c44 + anisotropic_parameters["b"][3, 3] = b44 anisotropic_parameters["a"][4, 4] = a55 - anisotropic_parameters["b1"][4, 4] = b55 - anisotropic_parameters["c1"][4, 4] = c55 + anisotropic_parameters["b"][4, 4] = b55 anisotropic_parameters["a"][5, 5] = a66 - anisotropic_parameters["b1"][5, 5] = b66 - anisotropic_parameters["c1"][5, 5] = c66 + anisotropic_parameters["b"][5, 5] = b66 # Fill the rest - for p in ["a", "b1", "b2"]: + for p in ["a", "b"]: anisotropic_parameters[p][0, 1] = 0.5 * ( (PsiI_11[p] + PsiI_22[p] - PsiI_33[p]) @@ -419,12 +275,6 @@ def psi_func(f, Pth, params): 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 @@ -438,7 +288,7 @@ def psi_func(f, Pth, params): 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 p in ["a", "b"]: for i, j in [(0, 1), (3, 4)]: l1 = copy(anti_anisotropic_parameters[p][j, :]) l2 = copy(anti_anisotropic_parameters[p][i, :]) @@ -449,10 +299,6 @@ def psi_func(f, Pth, params): 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, @@ -500,13 +346,11 @@ def fn(lnV, Pth, X, params): return stishovite_Q1, scalar_stv, stishovite_anisotropic, stishovite_relaxed -def get_models(): - - scalar_args = all_args[:16] - cell_args = all_args[16:23] - elastic_args = all_args[23:] +def get_models(all_args): - modify_Zhang_elasticity(scalar_args, cell_args, elastic_args) + scalar_args = all_args[:8] + cell_args = all_args[8:15] + elastic_args = all_args[15:] ( dVQ0, @@ -517,35 +361,16 @@ def get_models(): V0Q1overV0Q0, dDebye_0, P_tr_GPa, - fP_Zhang, - fP_Andrault, - fP_Wang, - fP_Nishihara, - fP2_Zhang, - fP2_Andrault, - fP2_Wang, - fP2_Nishihara, ) = scalar_args - (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_b2, PsiI_33_c2) = cell_args - (a11, a22, a33, a44, a55, a66, b11i, b33i, b44i, b66i, c44, c66, b112i, b332i) = ( + (a0Q1, b0Q1, PsiI_33_a, PsiI_33_b, PsiI_33_c, PsiI_33_d, f_PsiI_22) = cell_args + (a11, a22, a33, a44, a55, a66, b11, b33, b44, b66, c44, c66, d44, d66) = ( elastic_args ) - b55i = b44i - b22i = b11i - b222i = b112i + b22 = b11 + b55 = b44 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) + d55 = d44 models = make_models( dVQ0, @@ -561,8 +386,8 @@ def get_models(): PsiI_33_a, PsiI_33_b, PsiI_33_c, - PsiI_33_b2, - PsiI_33_c2, + PsiI_33_d, + f_PsiI_22, a11, a22, a33, @@ -578,16 +403,17 @@ def get_models(): c44, c55, c66, - b112, - b222, - b332, + d44, + d55, + d66, ) # 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) + # b axes potentially differs + PsiI_11_a = (1.0 - f_PsiI_22) * (1.0 - PsiI_33_a) PsiI_11_b = -0.5 * PsiI_33_b - PsiI_11_b2 = -0.5 * PsiI_33_b2 + PsiI_22_a = f_PsiI_22 * (1.0 - PsiI_33_a) + PsiI_22_b = -0.5 * PsiI_33_b sm = models[1] Q0, Q1 = sm.solution_model.interaction_endmembers @@ -605,13 +431,20 @@ def get_models(): ] table = [["", "$Q = 0$", "$Q = 1$"]] for i, p in enumerate(properties): - table.append([pps[i], f"{Q0.params[p]:.4e}", f"{Q1.params[p]:.4e}"]) + if p in ["K_0", "Kprime_0", "grueneisen_0", "q_0"]: + table.append([pps[i], f"{Q0.params[p]:.3e}", f"({Q1.params[p]:.3e})"]) + elif p == "Debye_0": + table.append( + [pps[i], f"({Q0.params[p]:.3e}; SLB22)", f"{Q1.params[p]:.3e}"] + ) + else: + table.append([pps[i], f"{Q0.params[p]:.3e}", f"{Q1.params[p]:.3e}"]) table[1][1] = "-" table.append(["$a_0$", "-", a0Q1]) table.append(["$b_0$", "-", b0Q1]) - print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".4e")) + print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".3e")) print("b_xs", sm.b_xs) @@ -621,60 +454,59 @@ def get_models(): table = [headers] table.extend( [ - ["$a$", PsiI_11_a, "as above", PsiI_33_a, a11, a22, a33, a44, a55, a66], [ - "$b_1$", + "$a$", + PsiI_11_a, + f"({PsiI_22_a:.3e})", + PsiI_33_a, + a11, + a22, + a33, + a44, + a55, + a66, + ], + [ + "$b$", PsiI_11_b, - "as above", - PsiI_33_b, + f"({PsiI_22_b:.3e})", + f"({PsiI_33_b:.3e})", b11, - "as above", + f"({b22:.3e})", b33, b44, - "as above", + f"({b55:.3e})", b66, ], [ - "$c_1$", + "$c$", PsiI_33_c, - "as above", - "as above", - "as above", - "as above", - "as above", + f"({PsiI_33_c:.3e})", + f"({PsiI_33_c:.3e})", + f"({PsiI_33_c:.3e})", + f"({PsiI_33_c:.3e})", + f"({PsiI_33_c:.3e})", c44, - "as above", + f"({c55:.3e})", 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", - "-", - "-", - "-", + "$d$", + PsiI_33_d, + f"({PsiI_33_d:.3e})", + f"({PsiI_33_d:.3e})", + f"({PsiI_33_d:.3e})", + f"({PsiI_33_d:.3e})", + f"({PsiI_33_d:.3e})", + d44, + f"({d55:.3e})", + d66, ], ] ) table = [ - [f"{item:.4e}" if type(item) is np.float64 else item for item in row] + [f"{item:.3e}" if type(item) is np.float64 else item for item in row] for row in table ] print( @@ -682,23 +514,8 @@ def get_models(): list(map(list, zip(*table))), headers="firstrow", tablefmt="latex_raw", - floatfmt=".4e", + floatfmt=".3e", ) ) - print("Corrections to pressure:") - print(f"${fP_Zhang:.4f}{fP2_Zhang*1e9:+.4f}P$ (Zhang)") - print(f"${fP_Andrault:.4f}{fP2_Andrault*1e9:+.4f}P$ (Andrault)") - print(f"${fP_Wang:.4f}{fP2_Wang*1e9:+.4f}P$ (Wang)") - print(f"${fP_Nishihara:.4f}{fP2_Nishihara*1e9:+.4f}P$ (Nishihara)") - return ( - models, - fP_Zhang, - fP_Andrault, - fP_Wang, - fP_Nishihara, - fP2_Zhang, - fP2_Andrault, - fP2_Wang, - fP2_Nishihara, - ) + return models diff --git a/contrib/anisotropic_stishovite/stishovite_model_plots.py b/contrib/anisotropic_stishovite/stishovite_model_plots.py index eddb1984..a4dc69dc 100644 --- a/contrib/anisotropic_stishovite/stishovite_model_plots.py +++ b/contrib/anisotropic_stishovite/stishovite_model_plots.py @@ -1,14 +1,39 @@ import numpy as np +import itertools +import scipy.interpolate as interpolate import matplotlib.pyplot as plt import matplotlib import matplotlib.cm as cm -from matplotlib.gridspec import GridSpec +from matplotlib.lines import Line2D +from matplotlib.legend_handler import HandlerTuple from copy import deepcopy from scipy.optimize import minimize -from stishovite_data import common_data +from scipy.integrate import cumulative_simpson +from stishovite_data import common_data, get_data, other_data from stishovite_model import get_models +from stishovite_parameters import all_args 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_data import ( + CN_GPa, + P_for_CN_new, + KNR_GPa, + CN_err_GPa, + SN_invGPa, + V_for_CN, + P_for_CN, + J2009_pressures, + J2009_temperatures, + J2009_dPsidfs, + J2009_stiffness_matrices, + C2024_pressures, + C2024_temperatures, + C2024_Vp, + C2024_Vp_unc, + C2024_Vs, + C2024_Vs_unc, + C2024_Vphi, + C2024_Vphi_unc, +) from stishovite_fit_eos import transition_pressure from burnman.tools.plot import plot_projected_elastic_properties from burnman.utils.math import is_positive_definite @@ -16,38 +41,435 @@ CN_GPa_orig = deepcopy(CN_GPa) -models = get_models() -_, stishovite_scalar, stishovite_unrelaxed, stishovite_relaxed = models[0] -( - fP_Zhang, - fP_Andrault, - fP_Wang, - fP_Nishihara, - fP2_Zhang, - fP2_Andrault, - fP2_Wang, - fP2_Nishihara, -) = models[1:] +models = get_models(all_args) +_, stishovite_scalar, stishovite_unrelaxed, stishovite_relaxed = models stishovite_relaxed.set_composition([1.0]) stishovite_unrelaxed.set_composition([0.5, 0.5]) -plot_isotropic_velocities = False +plot_isotropic_velocities = True plot_seismic_properties = False low_res_seismic_properties = False -plot_relaxed = True +plot_cell_properties = False +plot_relaxed = False +plot_relaxed_SN = False plot_Q_V_abc = False plot_G = False plot_lnabc = False +plot_Fischer_Andrault_splitting = False show_plots_one_by_one = False +if plot_Fischer_Andrault_splitting: + + data = common_data() + + PTV = np.concatenate( + ( + 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"]["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"]["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"]["Andrault_2003"], + data["abc_err"]["poststv"]["Andrault_2003"], + data["abc_err"]["stv"]["Fischer_2018"], + data["abc_err"]["poststv"]["Fischer_2018"], + ) + ) + + ln_abc = np.log(abc) + P = PTV[:, 0] + T = PTV[:, 1] + d_ba = ln_abc[:, 1] - ln_abc[:, 0] + + mask = np.abs(d_ba) < 1.0e-10 + + fig = plt.figure(figsize=(6, 4)) + ax = [fig.add_subplot(1, 1, 1)] + + scatter = plt.scatter(P / 1.0e9, T, c=d_ba, cmap="rainbow", s=50, edgecolor="k") + ax[0].scatter(P[mask] / 1.0e9, T[mask], c="white", s=50, edgecolor="k") + cbar = fig.colorbar(scatter, ax=ax[0], pad=0.1, label="$\\ln (b/a)$") + + args = deepcopy(all_args) + for ddebye_0 in [args[6] + 0.6, args[6] + 5.6, args[6] + 10.6]: + args[6] = ddebye_0 + models = get_models(args) + _, stishovite_scalar, stishovite_unrelaxed, stishovite_relaxed = models + debye_0 = stishovite_scalar.endmembers[0][0].params["Debye_0"] + temperatures = np.linspace(5.0, 4000.0, 41) + pressures = np.array( + [transition_pressure(stishovite_scalar, T) for T in temperatures] + ) + ax[0].plot( + pressures / 1.0e9, temperatures, label=f"$\\Theta_D$: {debye_0:.1f} K" + ) + + ax[0].legend() + ax[0].set_xlabel("Pressure (GPa)") + ax[0].set_ylabel("Temperature (K)") + ax[0].set_xlim(0.0, 125.0) + ax[0].set_ylim(0.0, 4000.0) + fig.set_tight_layout(True) + fig.savefig("figures/Fischer_Andrault_splitting.pdf") + plt.show() + + +if plot_cell_properties: + data = get_data() + models = get_models(all_args) + _, stishovite_scalar, stishovite_unrelaxed, stishovite_relaxed = models + + marker_cycle = itertools.cycle(Line2D.filled_markers) + colour_cycle = itertools.cycle(plt.rcParams["axes.prop_cycle"].by_key()["color"]) + + factor = 1.00 + KTR_GPa = factor * KNR_GPa + beta = np.sum(SN_invGPa[:, :3, :3], axis=1) + lnV_for_CN = np.log(V_for_CN) + + fig = plt.figure(figsize=(10, 8)) + ax = [fig.add_subplot(2, 2, i) for i in range(1, 5)] + + stishovite_relaxed.set_composition([1.0]) + J2009_V = stishovite_relaxed.evaluate(["V"], J2009_pressures, J2009_temperatures)[0] + + i = 0 + a_ln = ax[1].scatter( + -np.log(J2009_V / J2009_V[0]), + (J2009_dPsidfs[:, i, 0] + J2009_dPsidfs[:, i, 1] + J2009_dPsidfs[:, i, 2]) + - 1.0 / 3.0, + label="a, b (J09)", + marker="^", + ) + i = 2 + c_ln = ax[1].scatter( + -np.log(J2009_V / J2009_V[0]), + (J2009_dPsidfs[:, i, 0] + J2009_dPsidfs[:, i, 1] + J2009_dPsidfs[:, i, 2]) + - 1.0 / 3.0, + label="c (J09)", + marker="^", + ) + + beta_ab = (beta[:, 0] + beta[:, 1]) / 2.0 + ax[1].scatter( + lnV_for_CN[0] - lnV_for_CN, + beta[:, 0] * KNR_GPa - 1.0 / 3.0, + label="a, b (Z21)", + c=a_ln.get_facecolor(), + ) + ax[1].scatter( + lnV_for_CN[0] - lnV_for_CN, + beta[:, 1] * KNR_GPa - 1.0 / 3.0, + c=a_ln.get_facecolor(), + ) + ax[1].scatter( + lnV_for_CN[0] - lnV_for_CN, + beta[:, 2] * KNR_GPa - 1.0 / 3.0, + label="c (Z21)", + c=c_ln.get_facecolor(), + ) + + def dlncdlnV_linear(lnV, grad=1.2): + return beta[0, 2] * KNR_GPa[0] + grad * (lnV_for_CN[0] - lnV) + + stishovite_relaxed.set_composition([1.0]) + stishovite_relaxed.set_state(1.0e5, 298.15) + lnaxis0 = np.log(stishovite_relaxed.cell_parameters[:3]) + + pressures_model = np.linspace(1.0e5, 150.0e9, 301) + temperatures_model = 300.0 + 0.0 * pressures_model + stishovite_relaxed.set_composition([1.0]) + + prps = stishovite_relaxed.evaluate( + [ + "V", + "cell_parameters", + "isothermal_compressibility_tensor", + "isothermal_bulk_modulus_reuss", + ], + pressures_model, + temperatures_model, + ) + V_model = prps[0] + cell_parameters_model = prps[1] + isothermal_compressibility_model = prps[2] + bulk_moduli_model = prps[3] + + lnV_model = np.log(V_model) + for grad, linestyle in [(0.0, "--"), (1.2, ":")]: + label = f"$d^2\\Psi_3/df^2$ = {grad}" + ax[1].plot( + lnV_for_CN[0] - lnV_model, + dlncdlnV_linear(lnV_model, grad) - 1.0 / 3.0, + c=c_ln.get_facecolor(), + linestyle=linestyle, + label=label, + ) + + cs = np.exp( + lnaxis0[2] + - cumulative_simpson( + dlncdlnV_linear(lnV_model, grad), x=-lnV_model, initial=0 + ) + ) + ax[3].plot( + V_model * 1.0e6, + cs * 1.0e2, + c=c_ln.get_facecolor(), + linestyle=linestyle, + label=label, + ) + + ax[1].set_ylabel("$(\\beta_{{axis}}/\\beta_{{RT}}) - 1/3$") + ax[1].set_xlabel("$-\\ln (V/V_{{0}})$") + ax[1].set_xlim(-0.01, np.max(lnV_for_CN[0] - lnV_for_CN) + 0.1) + + marker = next(marker_cycle) + custom_legend = [] + labels = [] + for study in [ + "Andrault_2003", + # "Murakami_2003", + "Nishihara_2005", + "Wang_2012", + # "Grocholski_2013", + # "Fischer_2018", # No RT data + # "Sun_2019", + # "Zhang_2021", + "Zhang_2023", + ]: + marker = next(marker_cycle) + colour = next(colour_cycle) + + custom_legend.append( + ( + Line2D( + [0], + [0], + marker=marker, + color=colour, + markeredgecolor=colour, + markerfacecolor="none", + linestyle="", + ), + Line2D([0], [0], marker=marker, color=colour, linestyle=""), + ) + ) + labels.append(f"{study.split('_')[0]} et al. ({study.split('_')[1]})") + + for phase in ["stv", "poststv"]: + try: + V = data[study][phase]["V"] + a = data[study][phase]["a"] + b = data[study][phase]["b"] + c = data[study][phase]["c"] + + P = data[study][phase]["P"] + T = data[study][phase]["T"] + + mask = T < 310.0 + if phase == "poststv": + ax[0].scatter( + P[mask] / 1.0e9, V[mask] * 1.0e6, marker=marker, c=colour + ) + ax[2].scatter( + V[mask] * 1.0e6, a[mask] * 1.0e2, marker=marker, c=colour + ) + ax[2].scatter( + V[mask] * 1.0e6, b[mask] * 1.0e2, marker=marker, c=colour + ) + ax[2].scatter( + V[mask] * 1.0e6, + np.sqrt(a[mask] * b[mask]) * 1.0e2, + marker=marker, + c=colour, + s=3, + ) + ax[3].scatter( + V[mask] * 1.0e6, c[mask] * 1.0e2, marker=marker, c=colour + ) + else: + ax[0].scatter( + P[mask] / 1.0e9, + V[mask] * 1.0e6, + marker=marker, + edgecolors=colour, + facecolors="none", + ) + ax[2].scatter( + V[mask] * 1.0e6, + a[mask] * 1.0e2, + marker=marker, + edgecolors=colour, + facecolors="none", + ) + ax[2].scatter( + V[mask] * 1.0e6, + b[mask] * 1.0e2, + marker=marker, + edgecolors=colour, + facecolors="none", + ) + ax[3].scatter( + V[mask] * 1.0e6, + c[mask] * 1.0e2, + marker=marker, + edgecolors=colour, + facecolors="none", + ) + except KeyError: + pass + + ax[1].plot( + lnV_model[0] - lnV_model, + bulk_moduli_model * isothermal_compressibility_model[:, 0, 0] - 1.0 / 3.0, + linestyle="-", + c=a_ln.get_facecolor(), + ) + + ax[1].plot( + lnV_model[0] - lnV_model, + bulk_moduli_model * isothermal_compressibility_model[:, 1, 1] - 1.0 / 3.0, + linestyle="-", + c=a_ln.get_facecolor(), + ) + + ax[1].plot( + lnV_model[0] - lnV_model, + bulk_moduli_model * isothermal_compressibility_model[:, 2, 2] - 1.0 / 3.0, + linestyle="-", + c=c_ln.get_facecolor(), + label="model", + ) + + ax[2].plot( + V_model * 1.0e6, + cell_parameters_model[:, 0] * 1.0e2, + linestyle="-", + label="model", + c=a_ln.get_facecolor(), + ) + ax[2].plot( + V_model * 1.0e6, + cell_parameters_model[:, 1] * 1.0e2, + linestyle="-", + label="model", + c=a_ln.get_facecolor(), + ) + ax[3].plot( + V_model * 1.0e6, + cell_parameters_model[:, 2] * 1.0e2, + linestyle="-", + label="model", + c=c_ln.get_facecolor(), + ) + + ln = ax[0].plot(P_for_CN / 1.0e9, V_for_CN * 1.0e6, linestyle="--", c="purple") + custom_legend.append(ln[0]) + labels.append("Zhang et al. (2021), BLS") + + ln = ax[0].plot(P_for_CN_new / 1.0e9, V_for_CN * 1.0e6, c="purple") + custom_legend.append(ln[0]) + labels.append("Zhang et al. (2021), BLS cal.-free") + + ln = ax[0].plot(pressures_model / 1.0e9, V_model * 1.0e6, label="model", c="black") + custom_legend.append(ln[0]) + labels.append("model") + + ax[0].set_xlabel("P (GPa)") + ax[0].set_ylabel("V (cm$^3$/mol)") + ax[0].legend( + custom_legend, + labels, + handler_map={tuple: HandlerTuple(ndivide=None)}, + fontsize="small", + markerscale=0.8, + labelspacing=0.5, + ) + + for i in range(2, 4): + ax[i].set_xlim(14.1, 10.9) + ax[i].set_xlabel("V (cm$^3$/mol)") + ax[2].set_ylabel("a,b (cm/mol$^{{1/3}}$)") + ax[3].set_ylabel("c (cm/mol$^{{1/3}}$)") + ax[2].set_ylim(2.8 * 0.88, 2.8 * 1.01) + ax[3].set_ylim(1.786 * 0.88, 1.786 * 1.01) + + i = 1 + ax_twin = ax[i].twiny() + Ps = np.linspace(0.0, 120.0e9, 7) + Vs = stishovite_relaxed.evaluate(["V"], Ps, Ps * 0.0 + 298.15)[0] + ax_twin.set_xticks(-np.log(Vs / Vs[0])) + ax_twin.set_xticklabels(["%d" % P for P in Ps / 1.0e9]) + ax_twin.set_xlim(ax[i].get_xlim()) + ax_twin.set_xlabel("P (model, GPa)") + + for i in range(2, 4): + ax_twin = ax[i].twiny() + ax_twin.set_xticks(Vs * 1.0e6) + ax_twin.set_xticklabels(["%d" % P for P in Ps / 1.0e9]) + ax_twin.set_xlim(ax[i].get_xlim()) + ax_twin.set_xlabel("P (model, GPa)") + + ax[1].set_ylim(-1.0, 1.0) + ax[1].legend(fontsize=8) + fig.set_tight_layout(True) + fig.savefig("figures/stv_cell_properties.pdf") + plt.show() + 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) + fig = plt.figure() + ax = [fig.add_subplot(1, 1, 1)] + + for T, colour in [(300, "grey"), (1073, "black")]: + mask = np.abs(C2024_temperatures - T) < 1.0 + ax[0].errorbar( + seismic_model.depth(C2024_pressures[mask]) / 1.0e3, + C2024_Vp[mask] / 1.0e3, + C2024_Vp_unc[mask] / 1.0e3, + c=colour, + label=f"$V_{{VRH}}$ ({T:.0f} K; Chen et al., 2024)", + ) + ax[0].errorbar( + seismic_model.depth(C2024_pressures[mask]) / 1.0e3, + C2024_Vs[mask] / 1.0e3, + C2024_Vs_unc[mask] / 1.0e3, + c=colour, + ) + ax[0].errorbar( + seismic_model.depth(C2024_pressures[mask]) / 1.0e3, + C2024_Vphi[mask] / 1.0e3, + C2024_Vphi_unc[mask] / 1.0e3, + c=colour, + ) + + stishovite_relaxed.set_composition([1.0]) v = stishovite_relaxed.evaluate( [ "isentropic_shear_modulus_voigt", @@ -76,9 +498,6 @@ 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" ) @@ -99,9 +518,33 @@ ) 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() + T_labels = np.linspace(100.0, 3000.0, 30) + spl = interpolate.splrep(temperatures, depths, s=0, k=4) + depth_spline = interpolate.BSpline(*spl, extrapolate=False) + + ax_twin = ax[0].twiny() + ax_twin.set_xticks(depth_spline(T_labels) / 1.0e3) + xlabels = [f"{int(T)}" if np.abs(T % 200.0) < 1.0 else "" for T in T_labels] + ax_twin.set_xticklabels(xlabels) + ax_twin.set_xlim(ax[0].get_xlim()) + ax_twin.set_xlabel("Temperatures (K)", labelpad=10.0) + + P_labels = np.linspace(10, 140, 14) + spl = interpolate.splrep(pressures, depths, s=0, k=4) + depth_spline = interpolate.BSpline(*spl, extrapolate=False) + + ax_twin = ax[0].twiny() + ax_twin.tick_params(axis="x", direction="in", pad=-15) + ax_twin.set_xticks(depth_spline(P_labels * 1.0e9) / 1.0e3) + xlabels = [f"{int(P)}" if np.abs(P % 20.0) < 1.0 else "" for P in P_labels] + ax_twin.set_xticklabels(xlabels) + ax_twin.set_xlim(ax[0].get_xlim()) + ax_twin.set_xlabel("Pressures (GPa)", labelpad=-30.0) + + ax[0].set_ylim(0, 18) + ax[0].set_xlabel("Depth (km)") + ax[0].set_ylabel("Velocity (km/s)") + ax[0].legend(fontsize=8.0) fig.set_tight_layout(True) fig.savefig("figures/stv_isotropic_velocities.pdf") if show_plots_one_by_one: @@ -111,13 +554,28 @@ T = 2200.0 # for delta_P in [0.1e9]: # for delta_P in [-40.0e9, -1.0e9, -0.1e9, 0.1e9, 1.0e9, 40.0e9]: - # for delta_P in [-40.0e9, -30.e9, -20.e9, -10.e9, -1.0e9, -0.1e9, 0.1e9, 1.0e9, 10.e9, 20.e9, 30.e9, 40.0e9]: - for delta_P in [-10.0e9]: + for delta_P in [ + -40.0e9, + -30.0e9, + -20.0e9, + -10.0e9, + -1.0e9, + -0.1e9, + 0.1e9, + 1.0e9, + 10.0e9, + 20.0e9, + 30.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_unrelaxed.set_composition([0.5, 0.5]) + stishovite_relaxed.set_composition([1.0]) + stishovite_relaxed.set_state(P, T) stishovite_unrelaxed.set_state(P, T) @@ -164,6 +622,7 @@ plt.show() data = common_data() +data2 = other_data() PTV = np.concatenate( ( @@ -176,6 +635,10 @@ data["PTV"]["poststv"]["Andrault_2003"], data["PTV"]["stv"]["Fischer_2018"], data["PTV"]["poststv"]["Fischer_2018"], + data2["PTV"]["poststv"]["Murakami_2003"], + data2["PTV"]["poststv"]["Grocholski_2013"], + data2["PTV"]["stv"]["Zhang_2023"], + data2["PTV"]["poststv"]["Zhang_2023"], ) ) PTV_err = np.concatenate( @@ -189,6 +652,10 @@ data["PTV_err"]["poststv"]["Andrault_2003"], data["PTV_err"]["stv"]["Fischer_2018"], data["PTV_err"]["poststv"]["Fischer_2018"], + data2["PTV_err"]["poststv"]["Murakami_2003"], + data2["PTV_err"]["poststv"]["Grocholski_2013"], + data2["PTV_err"]["stv"]["Zhang_2023"], + data2["PTV_err"]["poststv"]["Zhang_2023"], ) ) abc = np.concatenate( @@ -202,6 +669,10 @@ data["abc"]["poststv"]["Andrault_2003"], data["abc"]["stv"]["Fischer_2018"], data["abc"]["poststv"]["Fischer_2018"], + data2["abc"]["poststv"]["Murakami_2003"], + data2["abc"]["poststv"]["Grocholski_2013"], + data2["abc"]["stv"]["Zhang_2023"], + data2["abc"]["poststv"]["Zhang_2023"], ) ) abc_err = np.concatenate( @@ -215,40 +686,69 @@ data["abc_err"]["poststv"]["Andrault_2003"], data["abc_err"]["stv"]["Fischer_2018"], data["abc_err"]["poststv"]["Fischer_2018"], + data2["abc_err"]["poststv"]["Murakami_2003"], + data2["abc_err"]["poststv"]["Grocholski_2013"], + data2["abc_err"]["stv"]["Zhang_2023"], + data2["abc_err"]["poststv"]["Zhang_2023"], + ) +) + +is_stv = np.concatenate( + ( + data["abc"]["stv"]["Ito_1974"] > 0.0, + data["abc"]["stv"]["Wang_2012"] > 0.0, + data["abc"]["stv"]["Nishihara_2005"] > 0.0, + data["abc"]["stv"]["Zhang_2021"] > 0.0, + data["abc"]["poststv"]["Zhang_2021"] < 0.0, + data["abc"]["stv"]["Andrault_2003"] > 0.0, + data["abc"]["poststv"]["Andrault_2003"] < 0.0, + data["abc"]["stv"]["Fischer_2018"] > 0.0, + data["abc"]["poststv"]["Fischer_2018"] < 0.0, + data2["abc"]["poststv"]["Murakami_2003"] < 0.0, + data2["abc"]["poststv"]["Grocholski_2013"] < 0.0, + data2["abc"]["stv"]["Zhang_2023"] > 0.0, + data2["abc"]["poststv"]["Zhang_2023"] < 0.0, ) ) +is_stv = is_stv[:, 0] id = np.concatenate( ( np.ones(len(data["abc_err"]["stv"]["Ito_1974"])) * 1, - np.ones(len(data["abc_err"]["stv"]["Wang_2012"])) * 2, + np.ones(len(data["abc_err"]["stv"]["Wang_2012"])) * 4, np.ones(len(data["abc_err"]["stv"]["Nishihara_2005"])) * 3, - np.ones(len(data["abc_err"]["stv"]["Zhang_2021"])) * 4, - np.ones(len(data["abc_err"]["poststv"]["Zhang_2021"])) * 4, - np.ones(len(data["abc_err"]["stv"]["Andrault_2003"])) * 5, - np.ones(len(data["abc_err"]["poststv"]["Andrault_2003"])) * 5, - np.ones(len(data["abc_err"]["stv"]["Fischer_2018"])) * 6, - np.ones(len(data["abc_err"]["poststv"]["Fischer_2018"])) * 6, + np.ones(len(data["abc_err"]["stv"]["Zhang_2021"])) * 6, + np.ones(len(data["abc_err"]["poststv"]["Zhang_2021"])) * 6, + np.ones(len(data["abc_err"]["stv"]["Andrault_2003"])) * 2, + np.ones(len(data["abc_err"]["poststv"]["Andrault_2003"])) * 2, + np.ones(len(data["abc_err"]["stv"]["Fischer_2018"])) * 5, + np.ones(len(data["abc_err"]["poststv"]["Fischer_2018"])) * 5, + np.ones(len(data2["abc_err"]["poststv"]["Murakami_2003"])) * 7, + np.ones(len(data2["abc_err"]["poststv"]["Grocholski_2013"])) * 8, + np.ones(len(data2["abc_err"]["stv"]["Zhang_2023"])) * 9, + np.ones(len(data2["abc_err"]["poststv"]["Zhang_2023"])) * 9, ) ) Ito_mask = id == 1 -Wang_mask = id == 2 +Andrault_mask = id == 2 Nishihara_mask = id == 3 -Zhang_mask = id == 4 -Andrault_mask = id == 5 -Fischer_mask = id == 6 - -PTV[Zhang_mask, 0] *= fP_Zhang + PTV[Zhang_mask, 0] * fP2_Zhang -PTV[Andrault_mask, 0] *= fP_Andrault + PTV[Andrault_mask, 0] * fP2_Andrault -PTV[Wang_mask, 0] *= fP_Wang + PTV[Wang_mask, 0] * fP2_Wang -PTV[Nishihara_mask, 0] *= fP_Nishihara + PTV[Nishihara_mask, 0] * fP2_Nishihara -P_for_CN_orig = deepcopy(P_for_CN) -P_for_CN *= fP_Zhang + P_for_CN * fP2_Zhang - -print( - f"Transition pressure at 3000 K: {transition_pressure(stishovite_scalar, 3000.0) / 1.0e9} GPa" -) +Wang_mask = id == 4 +Fischer_mask = id == 5 +Zhang_mask = id == 6 +Murakami_mask = id == 7 +Grocholski_mask = id == 8 +Zhang_mask2 = id == 9 + +id_labels = ["I74", "A03", "N05", "W12", "F18", "Z21"] +id_masks = [id == i for i in range(1, 6)] + +P_for_CN_orig = deepcopy(P_for_CN_new) + +for T in [298.15, 1000.0, 2000.0, 3000.0]: + print( + f"Transition pressure at {T} K: {transition_pressure(stishovite_scalar, T) / 1.0e9} GPa" + ) stishovite_relaxed.set_composition([1]) print(f"Consistent?: {check_anisotropic_eos_consistency(stishovite_relaxed)}") @@ -257,34 +757,43 @@ fig_lnabc = plt.figure(figsize=(10, 4)) ax_lnabc = [fig_lnabc.add_subplot(1, 2, i) for i in range(1, 3)] - colors = [ - "tab:olive", - "tab:cyan", - "tab:orange", - "tab:blue", - "tab:purple", - "tab:red", - ] labels = [ + "Zhang et al. (2023)", "Zhang et al. (2021)", "Fischer et al. (2018)", + "Grocholski et al. (2013)", "Wang et al. (2012)", "Nishihara et al. (2005)", + "Murakami et al. (2003)", "Andrault et al. (2003)", "Ito et al. (1974)", ] + marker_cycle = itertools.cycle(Line2D.filled_markers) + colour_cycle = itertools.cycle(plt.rcParams["axes.prop_cycle"].by_key()["color"]) + for i, mask in enumerate( - [Zhang_mask, Fischer_mask, Wang_mask, Nishihara_mask, Andrault_mask, Ito_mask] + [ + Zhang_mask2, + Zhang_mask, + Fischer_mask, + Grocholski_mask, + Wang_mask, + Nishihara_mask, + Murakami_mask, + Andrault_mask, + Ito_mask, + ] ): - c = colors[i] + c = next(colour_cycle) + marker = next(marker_cycle) 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) + s = 15 + ax_lnabc[0].scatter(lnV, np.log(abc[mask, 0]), color=c, marker=marker, s=s) ax_lnabc[0].errorbar( lnV, np.log(abc[mask, 0]), @@ -293,7 +802,7 @@ ls="none", color=c, ) - ax_lnabc[0].scatter(lnV, np.log(abc[mask, 1]), color=c) + ax_lnabc[0].scatter(lnV, np.log(abc[mask, 1]), color=c, marker=marker, s=s) ax_lnabc[0].errorbar( lnV, np.log(abc[mask, 1]), @@ -304,9 +813,11 @@ ) ax_lnabc[0].scatter( - lnV, 0.5 * np.log(abc[mask, 0] * abc[mask, 1]), s=5, color=c + lnV, 0.5 * np.log(abc[mask, 0] * abc[mask, 1]), s=2, color=c, marker=marker + ) + ax_lnabc[1].scatter( + lnV, np.log(abc[mask, 2]), color=c, label=labels[i], marker=marker, s=s ) - 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]), @@ -322,7 +833,7 @@ 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}}$)") handles, labels = ax_lnabc[1].get_legend_handles_labels() - ax_lnabc[1].legend(handles[::-1], labels[::-1]) + ax_lnabc[1].legend(handles[::-1], labels[::-1], fontsize=8) fig_lnabc.set_tight_layout(True) fig_lnabc.savefig("figures/stv_lnabc.pdf") if show_plots_one_by_one: @@ -380,9 +891,13 @@ if plot_Q_V_abc: + marker_cycle = itertools.cycle(Line2D.filled_markers) + colour_cycle = itertools.cycle(plt.rcParams["axes.prop_cycle"].by_key()["color"]) + 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)] + # 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)] + ax_Q = [fig_Q.add_subplot(1, 1, 1)] fig_V = plt.figure(figsize=(8, 4)) ax_V = [fig_V.add_subplot(1, 2, i) for i in range(1, 3)] @@ -399,12 +914,13 @@ temperatures = np.array([4000.0, 3000.0, 2000.0, 1000.0, 298.15]) stishovite_relaxed.set_composition([1.0]) + stishovite_relaxed.set_state(1.0e5, 300) 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) + pressures = np.linspace(Pmin, 150.0e9, 151) Ts = pressures * 0.0 + T (volumes, cell_params, molar_fractions) = stishovite_relaxed.evaluate( ["V", "cell_parameters", "molar_fractions"], pressures, Ts @@ -414,29 +930,91 @@ ["V", "cell_parameters"], pressures, Ts ) ax_Q[0].plot( - pressures / 1.0e9, molar_fractions[:, 0] - molar_fractions[:, 1], c=color + pressures / 1.0e9, + molar_fractions[:, 0] - molar_fractions[:, 1], + c=color, + label=f"{T} K", ) 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_V[0].plot(pressures / 1.0e9, volumes * 1.0e6, c=color, label=f"{T} K") + ax_V[1].plot( + pressures / 1.0e9, (volumes - volumes_Q0) * 1.0e6, c=color, label=f"{T} K" + ) + ax_abc[0].plot( + pressures / 1.0e9, cell_params[:, 0] * 1.0e2, c=color, label=f"{T} K" + ) 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) + ax_abc[1].plot( + pressures / 1.0e9, cell_params[:, 2] * 1.0e2, c=color, label=f"{T} K" + ) - 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 + temperatures = np.linspace(298.15, 4000.0, 41) + pressures = np.array( + [transition_pressure(stishovite_scalar, T) for T in temperatures] ) + volumes = stishovite_scalar.evaluate(["V"], pressures, temperatures)[0] + ax_V[0].plot(pressures / 1.0e9, volumes * 1.0e6, c="black", label="transition") + + for i in range(4): + marker = next(marker_cycle) + + custom_legend = [] + labels = [] + + for i in range(5): # don't plot Zhang et al (2021) + label = id_labels[i] + mask = id_masks[i] + marker = next(marker_cycle) + s = 15 + + # reversing the order makes for a better choice of label colours + scatter_V = ax_V[0].scatter( + PTV[mask, 0][::-1] / 1.0e9, + PTV[mask, 2][::-1] * 1.0e6, + c=PTV[mask, 1][::-1], + cmap=cmap, + norm=norm, + s=s, + marker=marker, + label=label, + ) + ax_abc[0].scatter( + PTV[mask, 0][::-1] / 1.0e9, + abc[mask, 0][::-1] * 1.0e2, + c=PTV[mask, 1][::-1], + cmap=cmap, + norm=norm, + s=s, + marker=marker, + label=label, + ) + ax_abc[0].scatter( + PTV[mask, 0][::-1] / 1.0e9, + abc[mask, 1][::-1] * 1.0e2, + c=PTV[mask, 1][::-1], + cmap=cmap, + norm=norm, + s=s, + marker=marker, + ) + scatter_c = ax_abc[1].scatter( + PTV[mask, 0][::-1] / 1.0e9, + abc[mask, 2][::-1] * 1.0e2, + c=PTV[mask, 1][::-1], + cmap=cmap, + norm=norm, + s=s, + marker=marker, + label=label, + ) + + ax_Q[0].legend(fontsize=8.0) + ax_abc[0].legend(fontsize=8.0) + ax_abc[1].legend(fontsize=8.0) + ax_V[0].legend(fontsize=8.0) + ax_V[1].legend(fontsize=8.0) ax_Q[0].set_ylim(0.0, 1.5) @@ -454,7 +1032,7 @@ def Q_misfit(args): return np.sum(np.power((angles * args[0] - Q_model) / angles_err, 2.0)) sol = minimize(Q_misfit, [0.0]) - + print(f"scaling parameter for Q to theta={1./sol.x[0]}") # instantiate a second axes that shares the same x-axis color = "tab:blue" ax_Q2 = ax_Q[0].twinx() @@ -464,12 +1042,12 @@ def Q_misfit(args): 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_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)") + fig_V.colorbar(scatter_V, ax=ax_V[1], label="Temperature (K)") + fig_abc.colorbar(scatter_c, ax=ax_abc[1], label="Temperature (K)") for i in range(2): ax_V[i].set_xlabel("Pressure (GPa)") @@ -482,6 +1060,12 @@ def Q_misfit(args): 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}}$)") + ax_V[0].set_xlim(0.0, 150.0) + ax_V[1].set_xlim(0.0, 150.0) + ax_abc[0].set_xlim(0.0, 150.0) + ax_abc[1].set_xlim(0.0, 150.0) + ax_Q[0].set_xlim(0.0, 150.0) + fig_Q.set_tight_layout(True) fig_V.set_tight_layout(True) fig_abc.set_tight_layout(True) @@ -532,7 +1116,7 @@ def Q_misfit(args): 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_{SR}$") + ax_relaxed[1].scatter(P_for_CN_new / 1.0e9, KNR_GPa, label="$K_{SR}$") ax_Q0[1].plot(pressures / 1.0e9, K_NR_Q0 / 1.0e9, label="$K_{SR} (Q0)$") ax_Q0[1].plot(pressures / 1.0e9, K_NR_Q1 / 1.0e9, label="$K_{SR} (Q1)$") @@ -549,14 +1133,26 @@ def Q_misfit(args): (2, 5, 5), ): c = ax_relaxed[axi].plot( - pressures / 1.0e9, C_N[:, i, j] / 1.0e9, label=f"$C_{{S{i+1}{j+1}}}$" + pressures / 1.0e9, + C_N[:, i, j] / 1.0e9, + label=f"$\\mathbb{{C}}_{{S{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 + J2009_pressures / 1.0e9, + J2009_stiffness_matrices[:, i, j] / 1.0e9, + label=f"$\\mathbb{{C}}_{{S{i+1}{j+1}}}$ (J09)", + color=color, + marker="^", + ) + ax_relaxed[axi].scatter( + P_for_CN_new / 1.0e9, + CN_GPa[:, i, j], + label=f"$\\mathbb{{C}}_{{S{i+1}{j+1}}}$ (Z21)", + color=color, ) # Plot original data @@ -580,7 +1176,7 @@ def Q_misfit(args): ) for i in range(3): - ax_relaxed[i].legend() + ax_relaxed[i].legend(fontsize=8) ax_Q0[i].legend() ax_relaxed[i].set_xlabel("Pressure (GPa)") @@ -597,5 +1193,116 @@ def Q_misfit(args): if show_plots_one_by_one: plt.show() +if plot_relaxed_SN: + 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, S_N, K_NR) = stishovite_relaxed.evaluate( + [ + "molar_fractions", + "isentropic_compliance_tensor", + "isentropic_bulk_modulus_reuss", + ], + pressures, + Ts, + ) + (S_N_u, K_NR_u) = stishovite_relaxed.unrelaxed.evaluate( + ["isentropic_compliance_tensor", "isentropic_bulk_modulus_reuss"], + pressures, + Ts, + molar_fractions, + ) + + stishovite_unrelaxed.set_composition([0.5, 0.5]) + (S_N_Q0, K_NR_Q0) = stishovite_unrelaxed.evaluate( + ["isentropic_compliance_tensor", "isentropic_bulk_modulus_reuss"], pressures, Ts + ) + stishovite_unrelaxed.set_composition([1.0, 0.0]) + (S_N_Q1, K_NR_Q1) = stishovite_unrelaxed.evaluate( + ["isentropic_compliance_tensor", "isentropic_bulk_modulus_reuss"], pressures, Ts + ) + + c = ax_relaxed[1].plot(pressures / 1.0e9, 1.0e9 / K_NR, label="$1/K_{SR}$") + ax_relaxed[1].plot( + pressures / 1.0e9, 1.0e9 / K_NR_u, linestyle=":", c=c[0].get_color() + ) + ax_relaxed[1].scatter(P_for_CN_new / 1.0e9, 1.0 / KNR_GPa, label="$1/K_{SR}$") + + ax_Q0[1].plot(pressures / 1.0e9, 1.0e9 / K_NR_Q0, label="$1/K_{SR} (Q0)$") + ax_Q0[1].plot(pressures / 1.0e9, 1.0e9 / K_NR_Q1, label="$1/K_{SR} (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, K_NR * S_N[:, i, j], label=f"$S_{{S{i+1}{j+1}}}$" + ) + color = c[0].get_color() + ax_relaxed[axi].plot( + pressures / 1.0e9, K_NR_Q0 * S_N_Q0[:, i, j], linestyle=":", c=color + ) + ax_relaxed[axi].scatter( + P_for_CN_new / 1.0e9, + KNR_GPa * SN_invGPa[:, i, j], + label=f"S_{{{i+1}{j+1}}} (Z21)", + color=color, + ) + ax_relaxed[axi].scatter( + J2009_pressures / 1.0e9, + J2009_dPsidfs, + label=f"S_{{{i+1}{j+1}}} (J09)", + color=color, + ) + + # Plot original data + ax_relaxed[axi].scatter( + P_for_CN_orig / 1.0e9, + KNR_GPa * SN_invGPa[:, i, j], + s=5, + alpha=0.5, + color=color, + ) + + ax_Q0[axi].plot( + pressures / 1.0e9, K_NR_Q0 * S_N_Q0[:, i, j], label=f"{i+1}{j+1} (Q0)" + ) + ax_Q0[axi].plot( + pressures / 1.0e9, K_NR_Q1 * S_N_Q1[:, i, j], label=f"{i+1}{j+1} (Q1)" + ) + + for i in range(3): + ax_relaxed[i].legend(fontsize=8) + ax_Q0[i].legend() + + ax_relaxed[i].set_xlabel("Pressure (GPa)") + ax_Q0[i].set_xlabel("Pressure (GPa)") + + ax_relaxed[i].set_ylabel("Compliance (1/GPa)") + ax_Q0[i].set_ylabel("Compliance (1/GPa)") + + fig_relaxed.set_tight_layout(True) + fig_Q0.set_tight_layout(True) + + fig_relaxed.savefig("figures/stv_relaxed_ST.pdf") + + if show_plots_one_by_one: + plt.show() + + if not show_plots_one_by_one: plt.show() diff --git a/contrib/anisotropic_stishovite/stishovite_parameters.py b/contrib/anisotropic_stishovite/stishovite_parameters.py index 9d93e02a..77feddf3 100644 --- a/contrib/anisotropic_stishovite/stishovite_parameters.py +++ b/contrib/anisotropic_stishovite/stishovite_parameters.py @@ -9,59 +9,18 @@ "$V_0 (Q=1) / V_0 (Q=0)$", "$\\Theta_0 (Q=1) - \\Theta_0 (Q=0)$", "$P_{{tr}} (T = T_{{ref}})$", - "$f_{{P}}$ (Zhang)", - "$f_{{P}}$ (Andrault)", - "$f_{{P}}$ (Wang)", - "$f_{{P}}$ (Nishihara)", - "$f_{{P2}}$ (Zhang)", - "$f_{{P2}}$ (Andrault)", - "$f_{{P2}}$ (Wang)", - "$f_{{P2}}$ (Nishihara)", ] -# 76 -scalar_args = [ - -1.15714715e-02, - -2.49554064e00, - -7.69118482e-04, - -2.08731542e-01, - -6.08850674e-01, - 9.93834483e-01, - 1.76227286e01, - 4.99610469e01, - 8.91475407e-01, - 9.89762819e-01, - 9.87845187e-01, - 9.95720410e-01, - 2.92917650e-13, - -9.08667999e-13, - -2.20268875e-13, - 1.55926921e-12, -] - - cell_param_names = [ "$a_0 (Q=1)$", "$b_0 (Q=1)$", "$\\Psi_{{3a}}$", "$\\Psi_{{3b}}$", "$\\Psi_{{3c}}$", - "$\\Psi_{{3b2}}$", - "$\\Psi_{{3c2}}$", + "$\\Psi_{{3d}}$", + "$f\\Psi_{{2}}$", ] -# 1117 -cell_args = [ - 2.72779892e-02, - 2.85698440e-02, - 2.10385858e-01, - -6.02813854e-01, - 1.10903361e00, - 8.86818639e-05, - -1.09619314e01, -] - - elastic_param_names = [ "$\\Psi_{{11a}}$", "$\\Psi_{{22a}}$", @@ -75,62 +34,99 @@ "$\\Psi_{{66b}}$", "$\\Psi_{{44c}}$", "$\\Psi_{{66c}}$", - "$\\Psi_{{11b2}}$", - "$\\Psi_{{33b2}}$", + "$\\Psi_{{44d}}$", + "$\\Psi_{{66d}}$", ] # let the b parameters be equal to # frel = -0.14 # b1_input = b1 * c1 * exp(c1 * frel) - 1) # 197 -elastic_args = [ - 3.24807928e-01, - 8.37663793e-01, - 4.93000185e-01, - 1.11959128e00, - 1.20914946e00, - 9.49237642e-01, - -4.29973046e-02, - 7.10161257e-04, - 3.85475084e-01, - 1.13312302e-01, - 9.74598492e-01, - 9.97390953e-01, - 5.85741559e-01, - 2.13037695e-01, -] + +scalar_args = np.array( + [ + 6.97812599e-03, + -3.18199894e00, + 2.02568387e-01, + 2.58439405e-01, + -4.82534573e-05, + 9.97227979e-01, + 7.24798635e00, + 5.05759030e01, + ] +) +cell_args = np.array( + [ + 0.02754879, + 0.02839605, + 0.19781129, + -0.13650963, + 7.93215475, + 0.03176523, + 0.46676305, + ] +) +elastic_args = np.array( + [ + 0.56655959, + 0.71555477, + 0.49784435, + 1.33005083, + 1.45154562, + 0.97764314, + 0.17631275, + -0.0516391, + -0.70121757, + -0.18867047, + 3.6433892, + 4.63773363, + 0.31665864, + 0.09046939, + ] +) + + +""" +# Removing the Fischer transition constraint +scalar_args = np.array([ 3.99584627e-03, -1.87172492e-01, 1.17470694e-01, -1.31671294e-01, + -1.99994808e+00, 9.97305028e-01, 1.69419908e+01, 4.97481580e+01]) +cell_args = np.array([ 0.02755226, 0.0284011 , 0.1912243 , -0.14349027, 7.76739196, + 0.02422077, 0.46590064]) +elastic_args = np.array([ 0.54560767, 0.66594729, 0.47720876, 1.33951091, 1.4542009 , + 0.98069674, 0.2092898 , -0.07160769, -0.72256214, -0.16372882, + 3.56540327, 5.38951401, 0.32085463, 0.12083754]) +""" + scalar_and_cell_args = np.concatenate((scalar_args, cell_args)) all_args = np.concatenate((scalar_args, cell_args, elastic_args)) all_param_names = scalar_param_names + cell_param_names + elastic_param_names scalar_bounds = ( - (-8.0e-3, -4.0e-3), # dVQ0 - (-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, + (-8.0e-3, 7.0e-3), # dVQ0 + (-10.0, 10.0), # dKQ0 + (0.0, 2.0), # dKpQ0 4.0292, + (-0.6, 0.3), # dgrQ0 1.55674, + (-2.0, 0.0), # dqQ0 2.2096, (0.99, 0.999), # V0Q1overV0Q0 - (0.0, 80.0), # dDebye_0 + (0.0, 100.0), # dDebye_0 (46, 53), # P_tr_GPa, - (0.9, 1.1), # s - (0.9, 1.1), # s - (0.9, 1.1), # s - (0.9, 1.1), # s - (-0.1, 0.1), # s - (-0.1, 0.1), - (-0.1, 0.1), # s - (-0.1, 0.1), -) # s +) 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 + (None, None), # PsiI_33_a + (None, None), # PsiI_33_b + (None, None), # PsiI_33_c + (None, None), # PsiI_33_d + (None, None), # f +) elastic_bounds = tuple((None, None) for i in range(14)) + + +print(f"Total number of arguments: {len(all_args)}") +print(f"Scalar arguments: {len(scalar_args)}") +print(f"Cell arguments: {len(cell_args)}") +print(f"Elastic arguments: {len(elastic_args)}")