{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmarking models for oxygen fugacity related variables\n",
    "\n",
    "This notebook benchmarks the models for oxygen fugacity related variables in VolFe where possible.\n",
    "\n",
    "## Python set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import VolFe as vf\n",
    "import math\n",
    "import pandas as pd\n",
    "import Thermobar as tb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for FMQ\n",
    "\n",
    "option = FMQbuffer, function = FMQ\n",
    "\n",
    "### 'ONeill87' O'Neill (1987) AmMin 72(1-2):67-75\n",
    "\n",
    "Appendix A. Supplementary material - Supplementary data 2. Tab = Table S6 S redox calculator (sulfate,KSOg2=ONeill22.xlsx).\n",
    "\n",
    "Cell R12: 0.13 [DFMQ]\n",
    "\n",
    "Cell S12: -8.29 [logfO2]\n",
    "\n",
    "Matches to 2 decimal places. Note spreadsheet uses +273 to convert to K, rather than 273.15 used in VolFe so T in spreadsheet = 1200.15 'C\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-8.294378372874451"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PT = {\"P\":1000.}\n",
    "PT[\"T\"]=1200.\n",
    "D = 0.13 # DFMQ\n",
    "\n",
    "my_models = [[\"FMQbuffer\", \"ONeill87\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "math.log10(vf.Dbuffer2fO2(PT,D,'FMQ',my_models))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Frost91' Frost (1991) in \"Oxide Minerals: Petrologic and Magnetic Significance\"\n",
    "\n",
    "Comparison with results from Thermobar (Wieser et al., 2022) - values match perfectly.\n",
    "\n",
    "Define input conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "PT = {\"P\":1000.}\n",
    "PT[\"T\"]=1200.\n",
    "D = 0.13 # DFMQ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run using VolFe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-8.096212368054847"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"FMQbuffer\", \"Frost91\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "math.log10(vf.Dbuffer2fO2(PT,D,'FMQ',my_models))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run using Thermobar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-8.096212368054847"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tb.calculate_logfo2_from_buffer_pos(T_K=PT[\"T\"]+273.15,P_kbar=PT['P']/1000.,fo2_offset=D,buffer='QFM')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for NNO\n",
    "\n",
    "option = NNObuffer, function = NNO\n",
    "\n",
    "### 'Frost91' Frost (1991) in \"Oxide Minerals: Petrologic and Magnetic Significance\"\n",
    "\n",
    "Comparison with results from Thermobar (Wieser et al., 2022) - values match perfectly.\n",
    "\n",
    "Define input conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "PT = {\"P\":1000.}\n",
    "PT[\"T\"]=1200.\n",
    "D = 0.13 # NNO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run using VolFe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-7.401725893493535"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [[\"NNObuffer\", \"Frost91\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "math.log10(vf.Dbuffer2fO2(PT,D,'NNO',my_models))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run using Thermobar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-7.401725893493535"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tb.calculate_logfo2_from_buffer_pos(T_K=PT[\"T\"]+273.15,P_kbar=PT['P']/1000.,fo2_offset=D,buffer='NNO')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for fO2-Fe3+/FeT\n",
    "\n",
    "option = fO2, function = fO2\n",
    "\n",
    "The following models do not have material available in the original papers for benchmarking:\n",
    "- 'Kress91A' Eq. (A-5, A-6) in Kress and Carmichael (1991) CMP 108:82-92 doi:10.1007/BF00307328 (although the value is similar to 'Kress91' as expected)\n",
    "- 'Borisov18' Eq. (4) from Borisov et al. (2018) CMP 173:98 doi:10.1007/s00410-018-1524-8   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Kress91' Eq. (7) in Kress and Carmichael (1991) CMP 108:82-92\n",
    "\n",
    "Comparison with results from Thermobar (Wieser et al., 2022) - values match to two signifcant figures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_analysis = {\n",
    "        \"Sample\": \"Hawaiian basalt\",\n",
    "        \"SiO2\": 51.29,  # wt%\n",
    "        \"TiO2\": 2.50,  # wt%\n",
    "        \"Al2O3\": 13.70,  # wt%\n",
    "        \"FeOT\": 11.04,  # wt%\n",
    "        \"MnO\": 0.02,  # wt%\n",
    "        \"MgO\": 6.70,  # wt%\n",
    "        \"CaO\": 11.03,  # wt%\n",
    "        \"Na2O\": 2.27,  # wt%\n",
    "        \"K2O\": 0.43,  # wt%\n",
    "        \"P2O5\": 0.,  # wt%\n",
    "        \"H2O\": 3.,  # wt%\n",
    "        \"CO2ppm\": 1000.,  # ppm\n",
    "        \"STppm\": 0.,  # ppm\n",
    "        \"Xppm\": 0.0,  # ppm\n",
    "        \"Fe3FeT\": 0.1,\n",
    "    }\n",
    "\n",
    "my_analysis = pd.DataFrame(my_analysis, index=[0])\n",
    "\n",
    "PT = {\"P\":1000.}\n",
    "PT[\"T\"]=1200.\n",
    "\n",
    "melt_wf=vf.melt_comp(0.,my_analysis)\n",
    "melt_wf['CO2'] = my_analysis.loc[0.,\"CO2ppm\"]/1000000.\n",
    "melt_wf[\"H2OT\"] = my_analysis.loc[0,\"H2O\"]/100.\n",
    "melt_wf[\"H2O\"] = my_analysis.loc[0,\"H2O\"]/100.\n",
    "melt_wf['ST'] = my_analysis.loc[0.,\"STppm\"]/1000000.\n",
    "melt_wf['CT'] = (melt_wf['CO2']/vf.species.loc['CO2','M'])*vf.species.loc['C','M']\n",
    "melt_wf['HT'] = (melt_wf['H2OT']/vf.species.loc['H2O','M'])*(2.*vf.species.loc['H','M'])\n",
    "melt_wf['XT'] = my_analysis.loc[0.,\"Xppm\"]/1000000.\n",
    "melt_wf[\"Fe3FeT\"] = my_analysis.loc[0.,\"Fe3FeT\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.2695394895766"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [['fO2','Kress91']]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "math.log10(vf.f_O2(PT,melt_wf,models=my_models))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overwriting Fe3Fet_Liq to that specified in the function input\n",
      "(1,)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\ehughes\\AppData\\Local\\Temp\\ipykernel_25252\\3582113947.py:3: FutureWarning: Calling float on a single element Series is deprecated and will raise a TypeError in the future. Use float(ser.iloc[0]) instead\n",
      "  math.log10(result['fo2_calc'])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-9.274068975767225"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "liq_comps=vf.melt_pysulfsat(melt_wf)\n",
    "result = tb.convert_fe_partition_to_fo2(liq_comps=liq_comps, T_K=pd.Series(PT['T'])+273.15, P_kbar=PT['P']/1000., model=\"Kress1991\", Fe3Fet_Liq=melt_wf[\"Fe3FeT\"],renorm=False)\n",
    "math.log10(result['fo2_calc'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the value for 'Kress91A' is similar, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-9.206761031816585"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [['fO2','Kress91A']]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "math.log10(vf.f_O2(PT,melt_wf,models=my_models))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'ONeill18' Eq. (9a) in O'Neill et al. (2018) EPSL 504:152-162\n",
    "\n",
    "Appendix A. Supplementary material - Supplementary data 2. Tab = Table S6 S redox calculator (sulfate,KSOg2=ONeill22.xlsx).\n",
    "\n",
    "Cell S12: -8.29 [logfO2]\n",
    "\n",
    "Matches to 2 decimal places. Note spreadsheet uses +273 to convert to K, rather than 273.15 used in VolFe so T in spreadsheet = 1200.15 'C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-8.289923385626128"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_models = [['fO2','ONeill18']]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "\n",
    "math.log10(vf.f_O2(PT,melt_wf,models=my_models))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "volfe-dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}