diff --git a/src/pyEQL/Benchmarking_plot b/src/pyEQL/Benchmarking_plot new file mode 100644 index 00000000..a6013347 --- /dev/null +++ b/src/pyEQL/Benchmarking_plot @@ -0,0 +1,74 @@ +import numpy as np +import pandas as pd +from importlib import resources +import matplotlib.pyplot as plt +import pyEQL +import pyEQL.database as db + +def _load_db_csv(name: str): + '''load dataset from pyEQL database''' + return pd.read_csv(resources.files(db).joinpath(name)) +ac = _load_db_csv("CRC_activity_coefficient.csv") +cond = _load_db_csv("CRC_conductivity.csv") + +def ionic_strength(df): + '''calculate ionic strength of the selected data''' + z={"Na[+1]":1,"K[+1]":1,"Ca[+2]":2,"Mg[+2]":2,"Cl[-1]":-1,"SO4[-2]":-2} + I=np.zeros(len(df)) + for k,v in z.items(): + I+=df[k].values*v*v + return 0.5*I + +def salt_label(ions): + '''generate salt label as plot title''' + a=ions[0].split("[")[0] + b=ions[1].split("[")[0] + if a in ["Ca","Mg"] and b=="Cl":return a+"Cl2" + if a in ["Na","K"] and b=="SO4":return a+"2SO4" + return a+b + +def pyEQL_sim(ions, conc, prop): + '''pyEQL model result calculation''' + s=pyEQL.Solution([[i,f"{c} mol/kg"] for i,c in zip(ions,conc)]) + if prop == "activity_coefficient": + nu = [abs(int(i.split("[")[1].split("]")[0])) for i in ions] + gammas = [s.get_activity_coefficient(i).magnitude for i in ions] + prod = 1.0 + for g, n in zip(gammas, nu): + prod *= g ** n + return prod ** (1 / sum(nu)) + elif prop == "conductivity": + return s.conductivity.to("mS/cm").magnitude + else: + raise ValueError(f"Unknown property: {prop}") + +def filter_df(df,ions): + '''filter selected experimental data''' + m=(df[ions[0]]>0)&(df[ions[1]]>0) + return df[m] + +# main benchmark +def benchmark(df,ions,prop,n=200): + '''plot experimental data against pyEQL calculation''' + df = filter_df(df,ions) + Is = ionic_strength(df) + y = df["mean_activity_coefficient"] if prop=="activity_coefficient" else df["conductivity"] + + I_sim = np.linspace(Is.min(),Is.max()+0.1,n) + base_conc = df[ions].mean().values + y_sim = [] + for I_s in I_sim: + scale = I_s / np.mean(Is) if np.mean(Is) != 0 else 1.0 + conc = base_conc * scale + y_sim.append(pyEQL_sim(ions, conc, prop)) + y_sim = np.array(y_sim) + + ylabel_map={"activity_coefficient":"Mean Activity Coefficient","conductivity":"Conductivity (mS/cm)"} + plt.plot(I_sim,y_sim,color="r",lw=2,label="pyEQL") + plt.scatter(Is,y,s=18,color="k",marker="x",label="Experiment") + plt.xlabel("Ionic Strength (m)") + plt.ylabel(ylabel_map[prop]) + plt.title(salt_label(ions)) + plt.legend() + plt.tight_layout() + plt.show() diff --git a/src/pyEQL/database/CRC_activity_coefficient.csv b/src/pyEQL/database/CRC_activity_coefficient.csv new file mode 100644 index 00000000..5d884239 --- /dev/null +++ b/src/pyEQL/database/CRC_activity_coefficient.csv @@ -0,0 +1,74 @@ +Na[+1],K[+1],Ca[+2],Mg[+2],Cl[-1],SO4[-2],mean_activity_coefficient +0.0,0.0,0.1,0.0,0.2,0.0,0.518 +0.0,0.0,0.2,0.0,0.4,0.0,0.472 +0.0,0.0,0.3,0.0,0.6,0.0,0.455 +0.0,0.0,0.4,0.0,0.8,0.0,0.448 +0.0,0.0,0.5,0.0,1.0,0.0,0.448 +0.0,0.0,0.6,0.0,1.2,0.0,0.453 +0.0,0.0,0.7,0.0,1.4,0.0,0.46 +0.0,0.0,0.8,0.0,1.6,0.0,0.47 +0.0,0.0,0.9,0.0,1.8,0.0,0.484 +0.0,0.0,1.0,0.0,2.0,0.0,0.5 +0.0,0.1,0.0,0.0,0.1,0.0,0.77 +0.0,0.2,0.0,0.0,0.2,0.0,0.718 +0.0,0.3,0.0,0.0,0.3,0.0,0.688 +0.0,0.4,0.0,0.0,0.4,0.0,0.666 +0.0,0.5,0.0,0.0,0.5,0.0,0.649 +0.0,0.6,0.0,0.0,0.6,0.0,0.637 +0.0,0.7,0.0,0.0,0.7,0.0,0.626 +0.0,0.8,0.0,0.0,0.8,0.0,0.618 +0.0,0.9,0.0,0.0,0.9,0.0,0.61 +0.0,1.0,0.0,0.0,1.0,0.0,0.604 +0.0,0.2,0.0,0.0,0.0,0.1,0.441 +0.0,0.4,0.0,0.0,0.0,0.2,0.36 +0.0,0.6,0.0,0.0,0.0,0.3,0.316 +0.0,0.8,0.0,0.0,0.0,0.4,0.286 +0.0,1.0,0.0,0.0,0.0,0.5,0.264 +0.0,1.2,0.0,0.0,0.0,0.6,0.246 +0.0,1.4,0.0,0.0,0.0,0.7,0.232 +0.0,0.0,0.0,0.1,0.2,0.0,0.529 +0.0,0.0,0.0,0.2,0.4,0.0,0.489 +0.0,0.0,0.0,0.3,0.6,0.0,0.477 +0.0,0.0,0.0,0.4,0.8,0.0,0.475 +0.0,0.0,0.0,0.5,1.0,0.0,0.481 +0.0,0.0,0.0,0.6,1.2,0.0,0.491 +0.0,0.0,0.0,0.7,1.4,0.0,0.506 +0.0,0.0,0.0,0.8,1.6,0.0,0.522 +0.0,0.0,0.0,0.9,1.8,0.0,0.544 +0.0,0.0,0.0,1.0,2.0,0.0,0.57 +0.0,0.0,0.0,0.1,0.0,0.1,0.15 +0.0,0.0,0.0,0.2,0.0,0.2,0.107 +0.0,0.0,0.0,0.3,0.0,0.3,0.0874 +0.0,0.0,0.0,0.4,0.0,0.4,0.0756 +0.0,0.0,0.0,0.5,0.0,0.5,0.0675 +0.0,0.0,0.0,0.6,0.0,0.6,0.0616 +0.0,0.0,0.0,0.7,0.0,0.7,0.0571 +0.0,0.0,0.0,0.8,0.0,0.8,0.0536 +0.0,0.0,0.0,0.9,0.0,0.9,0.0508 +0.0,0.0,0.0,1.0,0.0,1.0,0.0485 +0.1,0.0,0.0,0.0,0.1,0.0,0.778 +0.2,0.0,0.0,0.0,0.2,0.0,0.735 +0.3,0.0,0.0,0.0,0.3,0.0,0.71 +0.4,0.0,0.0,0.0,0.4,0.0,0.693 +0.5,0.0,0.0,0.0,0.5,0.0,0.681 +0.6,0.0,0.0,0.0,0.6,0.0,0.673 +0.7,0.0,0.0,0.0,0.7,0.0,0.667 +0.8,0.0,0.0,0.0,0.8,0.0,0.662 +0.9,0.0,0.0,0.0,0.9,0.0,0.659 +1.0,0.0,0.0,0.0,1.0,0.0,0.657 +0.2,0.0,0.0,0.0,0.0,0.1,0.445 +0.4,0.0,0.0,0.0,0.0,0.2,0.365 +0.6,0.0,0.0,0.0,0.0,0.3,0.32 +0.8,0.0,0.0,0.0,0.0,0.4,0.289 +1.0,0.0,0.0,0.0,0.0,0.5,0.266 +1.2,0.0,0.0,0.0,0.0,0.6,0.248 +1.4,0.0,0.0,0.0,0.0,0.7,0.233 +1.6,0.0,0.0,0.0,0.0,0.8,0.221 +1.8,0.0,0.0,0.0,0.0,0.9,0.21 +2.0,0.0,0.0,0.0,0.0,1.0,0.201 +0.0,0.0,2.0,0.0,4.0,0.0,0.784 +0.0,2.0,0.0,0.0,2.0,0.0,0.573 +0.0,5.0,0.0,0.0,5.0,0.0,0.593 +2.0,0.0,0.0,0.0,2.0,0.0,0.668 +5.0,0.0,0.0,0.0,5.0,0.0,0.874 +4.0,0.0,0.0,0.0,0.0,2.0,0.155 diff --git a/src/pyEQL/database/CRC_conductivity.csv b/src/pyEQL/database/CRC_conductivity.csv new file mode 100644 index 00000000..5bbac325 --- /dev/null +++ b/src/pyEQL/database/CRC_conductivity.csv @@ -0,0 +1,60 @@ +Na[+1],K[+1],Ca[+2],Mg[+2],Cl[-1],SO4[-2],conductivity +0.0005,0,0,0,0.0005,0,0.06222 +0.001,0,0,0,0.001,0,0.12368 +0.005,0,0,0,0.005,0,0.60295 +0.01,0,0,0,0.01,0,1.1845 +0.02,0,0,0,0.02,0,2.314 +0.05,0,0,0,0.05,0,5.5505 +0.1,0,0,0,0.1,0,10.669 +0,0.0005,0,0,0.0005,0,0.07387 +0,0.001,0,0,0.001,0,0.14688 +0,0.005,0,0,0.005,0,0.7169 +0,0.01,0,0,0.01,0,1.412 +0,0.02,0,0,0.02,0,2.7654 +0,0.05,0,0,0.05,0,6.665 +0,0.1,0,0,0.1,0,12.89 +0,0,0.00025,0,0.0005,0,0.06593 +0,0,0.0005,0,0.001,0,0.1303 +0,0,0.0025,0,0.005,0,0.62095 +0,0,0.005,0,0.01,0,1.203 +0,0,0.01,0,0.02,0,2.3118 +0,0,0.025,0,0.05,0,5.421 +0,0,0.05,0,0.1,0,10.241 +0,0,0,0.00025,0.0005,0,0.062775 +0,0,0,0.0005,0.001,0,0.12415 +0,0,0,0.0025,0.005,0,0.59125 +0,0,0,0.005,0.01,0,1.1449 +0,0,0,0.01,0.02,0,2.1998 +0,0,0,0.025,0.05,0,5.1515 +0,0,0,0.05,0.1,0,9.705 +0.0005,0,0,0,0,0.00025,0.06284 +0.001,0,0,0,0,0.0005,0.12409 +0.005,0,0,0,0,0.0025,0.58545 +0.01,0,0,0,0,0.005,1.1238 +0.02,0,0,0,0,0.01,2.1346 +0.05,0,0,0,0,0.025,4.885 +0.1,0,0,0,0,0.05,8.994 +0,0,0.045045,0,0.09009009,0,8.94305451 +0,0,0.09009009,0,0.18018018,0,17.3340686 +0,0,0.18018018,0,0.36036036,0,32.4599756 +0,0,0,0.05263,0.105263,0,9.49509491 +0,0,0,0.105263,0.210526,0,18.3277413 +0,0,0,0.210526,0.421052,0,34.4473211 +0,0,0,0.04166666,0,0.04166666,4.52673129 +0,0,0,0.08333333,0,0.08333333,8.3910141 +0,0,0,0.166666667,0,0.166666667,14.6842747 +0,0,0,0.416666667,0,0.4166667,30.251814 +0,0.067114,0,0,0.067114,0,9.05346259 +0,0.134228,0,0,0.134228,0,17.3340686 +0,0.268456,0,0,0.268456,0,32.5703837 +0,0.67114,0,0,0.67114,0,79.3834098 +0,0.05747126,0,0,0,0.02873563,6.40366866 +0,0.11494253,0,0,0,0.05747126,12.365705 +0,0.22988506,0,0,0,0.11494253,23.1856969 +0,0.57471264,0,0,0,0.28735632,52.9958786 +0.085470085,0,0,0,0.085470085,0,9.05346259 +0.17094017,0,0,0,0.17094017,0,17.6652929 +0.34188034,0,0,0,0.34188034,0,33.3432403 +0.070422535,0,0,0,0,0.035211268,6.51407674 +0.14084507,0,0,0,0,0.070422535,12.365705 +0.28169014,0,0,0,0,0.14084507,21.8607999 \ No newline at end of file