Benchmarking visualization

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()

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

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