# -*- coding: utf-8 -*- """ This code is part of the publication "Re-evaluating fitting ambiguities and physicochemical assumptions of the Debye-Hückel theory" by Benjamin Janotta*, Maximilian Schalenbach*, Hermann Tempel, Rüdiger-A. Eichel *Corresponding Authors The code includes measurement data, models, and adapted codes of other publications. We refer the reader/user to the respective publications. The sources are given in the text at the appropriate lines as well as the main text of this publication. Script used in Spyder (Python 3.9) """ """ Parametrization for Figure 6B: r_+ r_- epsilon_0 Eq. 18 for electrostatic? BiMSA: 0.95 1.81 78.14 no EDH-H+MAL (C=0.2) 3.84 3.60 78.14 no EDH-H (C=0.1) 2.35 2.35 78.14 no DHFULL+HS 3.2 2.8 78.14 no DHFULL,eps_r+Born+HS 3.4 2.8 78.14 yes DHLL+HS+Born (epsilon_w=95) 3.2 2.151 95 no """ import numpy as np # Version 1.26.4 import scipy as sp # Version 1.11.2 import matplotlib import matplotlib.pyplot as plt from BiMSA_adapted import BiMSA ## import constants # vacuum permittivity, electron charge, Avogadro number, ideal gas constant from scipy.constants import epsilon_0,e,N_A,R F = sp.constants.physical_constants['Faraday constant'][0] #Faraday constant kB = sp.constants.k #Boltzmann constant def get_nacl_activity_data(): """ Returns experimental NaCl activity data (concentration and activity coefficients). """ molar_mass_nacl = 58.44e-3 pure_water_density = 0.99704 molality = np.array([0, 1e-3, 1e-2, 5e-2, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6]) activity_coefficients_molal = np.array([ 10, 9.65, 9.01, 8.23, 7.78, 7.35, 7.1, 6.93, 6.81, 6.73, 6.67, 6.62, 6.59, 6.57, 6.54, 6.55, 6.57, 6.62, 6.68, 6.88, 7.14, 7.46, 7.83, 8.26, 8.74, 9.28, 9.86 ]) * 1e-1 molality_density = np.array([0.0986, 0.4971, 0.9865, 1.993, 2.9296, 3.4865, 3.9873, 5.4952, 5.8023]) density_values = np.array([1.0011, 1.0169, 1.0357, 1.072, 1.1034, 1.1212, 1.1366, 1.1803, 1.1888]) interpolated_density = np.interp(molality, molality_density, density_values) activity_coeff_molar = activity_coefficients_molal * (1 + molality * molar_mass_nacl) * pure_water_density / interpolated_density concentration_molar = molality * interpolated_density / (1 + molality * molar_mass_nacl) return concentration_molar, activity_coeff_molar def calculate_dielectric_constant(pure_solvent_epsilon, salt_concentration): """ Calculates the dielectric constant of an electrolyte solution based on salt concentration. equation from: R. Buchner, G.T. Hefter, P.M. May, Dielectric Relaxation of Aqueous NaCl Solutions, J. Phys. Chem. A 103 (1999) 1–9. https://doi.org/10.1021/jp982977k. parameterization (15.2;3.64) for 25°C. epsilon_r0 = 78.14 """ return pure_solvent_epsilon-15.2*(salt_concentration/1000)+3.64*(salt_concentration/1000)**(3/2) def compute_free_ions(activity_func, salt_concentration, ion_radii): """ Computes the fraction of free ions using the mass action law and activity coefficients. Data for calculation from molecular dynamics simulation from: T. Driesner, T.M. Seward, I.G. Tironi, Molecular dynamics simulation study of ionic hydration and ion association in dilute and 1 molal aqueous sodium chloride solutions from ambient to supercritical conditions, Geochimica et Cosmochimica Acta 62 (1998) 3095–3107. https://doi.org/10.1016/S0016-7037(98)00207-5. """ min_radius, max_radius = ion_radii total_concentration = 0.985e3 # mol/m³ c_na_cl = 0.633 * total_concentration c_na_cl_aq = 0.163 * total_concentration c_na2_cl = 1/3 * 0.185 * total_concentration c_na_cl2 = 1/3 * 0.185 * total_concentration c_na2_cl2_aq = 0.019 * total_concentration c_total_ions = c_na_cl + c_na2_cl + c_na_cl2 log_gamma = activity_func(min_radius, max_radius, c_total_ions) gamma = np.exp(log_gamma[0]) activity_na_cl = c_na_cl * gamma activity_na2_cl = c_na2_cl * gamma k2 = c_na_cl_aq / (activity_na_cl ** 2) k3 = activity_na2_cl / (activity_na_cl ** 3) k4 = (c_na2_cl2_aq / 2) / (activity_na_cl ** 4) def optimize_concentration_and_gamma(variables, args): """ Optimization function for ionic concentrations and activity coefficients. """ c_na, gamma = variables k2, k3, k4, c_na_total = args c_ions_total = (2 * k3 * gamma**2 * c_na**3 + c_na) log_gamma_new = activity_func(min_radius, max_radius, c_ions_total) gamma_new = np.exp(log_gamma_new[0]) return [c_na + k2 * gamma_new**2 * c_na**2 + k3 * 3 * gamma_new**2 * c_na**3 + k4 * 2 * gamma_new**4 * c_na**4 - c_na_total, gamma - gamma_new] c_na_guess = salt_concentration[0] gamma_guess = 1 c_na_values = np.zeros(len(salt_concentration)) gamma_values = np.zeros(len(salt_concentration)) for idx, c_na_total in enumerate(salt_concentration): solution = sp.optimize.fsolve(optimize_concentration_and_gamma, (c_na_guess, gamma_guess), args=[k2, k3, k4, c_na_total]) c_na_guess, gamma_guess = solution gamma_values[idx] = gamma_guess c_na_values[idx] = c_na_guess alpha_single_ions = c_na_values / salt_concentration alpha_all_ions = (c_na_values + 2 * k3 * gamma_values**2 * c_na_values**3) / salt_concentration return alpha_single_ions, alpha_all_ions, gamma_values def mean_spherical_approximation(radii,permittivity,salt_concentration,ion_valences,T): ''' calulation of the electrostatic (El) and Hard Sphere (HS) contributions for the gamma: activity coefficient; An: anion; Cat:Cation Mean Spherical Approximation (MSA) radii: function for acitivity coefficients of ions permittivity salt_concentration: list of concentrations # should be monotonically increasing. ion_valences: valence of [cation,anion] T: temperature implementation based on equations in: W. Verweij, J.-P. Simonin, Implementing the Mean Spherical Approximation Model in the Speciation Code CHEAQS Next at High Salt Concentrations, Journal of Solution Chemistry 49 (2020) 1319–1327. https://doi.org/10.1007/s10953-020-01008-9. ''' beta = 1/(kB*T) ln_y_An = np.zeros(len(salt_concentration)) ln_y_Cat = np.zeros(len(salt_concentration)) ln_y_An_HS = np.zeros(len(salt_concentration)) ln_y_Cat_HS = np.zeros(len(salt_concentration)) ln_y_An_El = np.zeros(len(salt_concentration)) ln_y_Cat_El = np.zeros(len(salt_concentration)) # epsilon_c_dep = False sigmaAnion = radii[1]#3.62e-10 # m # Cl- radii[1]#radii[1]# sigmaCation = radii[0]#2.887e-10 # m Na+ radii[0] z_Anion = ion_valences[1] z_Cation = ion_valences[0] def MSA_Electrostatic_Contribution(inputs): #optimizer for electrostatic contribution Gamma,ln_y_AnEl,ln_y_CatEl,eta,Omega= inputs F = np.zeros(5) F[0] = 4*Gamma**2 - beta * e**2/ (permittivity)* (rhoAnion*((z_Anion-eta*sigmaAnion**2)/(1+Gamma*sigmaAnion))**2 +rhoCation*((z_Cation-eta*sigmaCation**2)/(1+Gamma*sigmaCation))**2) F[1] = eta - ( np.pi/(Omega*2*x)*(rhoAnion*z_Anion*sigmaAnion/(1+Gamma*sigmaAnion)+rhoCation*z_Cation*sigmaCation/(1+Gamma*sigmaCation))) F[2] = Omega - ( 1 + np.pi/(2*x) * (rhoAnion*sigmaAnion**3/(1+Gamma*sigmaAnion)+rhoCation*sigmaCation**3/(1+Gamma*sigmaCation))) F[3] = ln_y_AnEl + beta*e**2/(permittivity*4*np.pi) * (Gamma*z_Anion**2/(1+Gamma*sigmaAnion) + eta*sigmaAnion*((2*z_Anion-eta*sigmaAnion**2)/(1+Gamma*sigmaAnion)+eta/3*sigmaAnion**2)) F[4] = ln_y_CatEl + beta*e**2/(permittivity*4*np.pi) * (Gamma*z_Cation**2/(1+Gamma*sigmaCation) + eta*sigmaCation*((2*z_Cation-eta*sigmaCation**2)/(1+Gamma*sigmaCation)+eta/3*sigmaCation**2)) return F # intitial guesses for optimization of MSA electrostatic contribution kappa = np.sqrt(beta *e**2/(permittivity)*(z_Anion**2*salt_concentration[0]+z_Cation**2*salt_concentration[0])) ln_y_CatEl0 = 0 ln_y_AnEl0 = 0 Gamma0 = kappa/2 eta0 = -4.414229*1e8 Omega0 = 1 initGuess = Gamma0,ln_y_AnEl0,ln_y_CatEl0,eta0,Omega0 # go through optimization for each concentration in salt_concentration for idx, ci in enumerate(salt_concentration): rhoAnion = ci*N_A # ions/m3 rhoCation = ci*N_A #### Hard Sphere Contribution ### X0 = np.pi/6 * (rhoAnion*sigmaAnion**0+rhoCation*sigmaCation**0) X1 = np.pi/6 * (rhoAnion*sigmaAnion**1+rhoCation*sigmaCation**1) X2 = np.pi/6 * (rhoAnion*sigmaAnion**2+rhoCation*sigmaCation**2) X3 = np.pi/6 * (rhoAnion*sigmaAnion**3+rhoCation*sigmaCation**3) x = 1-X3 F1 = 3*X2/x F2 = 3*X1/x+3*X2**2/(X3*x**2)+3*X2**2/X3**2*np.log(x) F3 = (X0-X2**3/X3**2)/x+(3*X1*X2-X2**3/X3**2)/x**2+2*X2**3/(X3*x**3)-2*X2**3/X3**3*np.log(x) ln_gamma_HS_An = -np.log(x) + sigmaAnion*F1+sigmaAnion**2*F2+sigmaAnion**3*F3 ln_yHSCation = -np.log(x) + sigmaCation*F1+sigmaCation**2*F2+sigmaCation**3*F3 ### Electrical contribution ### sol=sp.optimize.leastsq(MSA_Electrostatic_Contribution,initGuess) initGuess = sol[0] #print(sol[0]) Gamma,ln_y_AnEl,ln_y_CatEl,eta,Omega = initGuess ln_y_An_HS[idx] = ln_gamma_HS_An ln_y_Cat_HS[idx] = ln_yHSCation ln_y_An_El[idx] = ln_y_AnEl ln_y_Cat_El[idx] = ln_y_CatEl ln_y_An[idx] = ln_y_AnEl ln_y_Cat[idx] = ln_y_CatEl ln_yEl = (ln_y_An+ln_y_Cat)/2 ln_yHS = (ln_y_An_HS+ln_y_Cat_HS)/2 #return ln(y) for electrostatic contr. and HS contr. return ln_yEl,ln_yHS def DHLL(radii,permittivity,salt_concentration,ion_valences,T): ''' Debye Hückel limiting law ;implementation based on: P. Debye, E. Hückel, Zur Theorie der Elektrolyte., Physikalische Zeitschrift 1923 (1923) 185–206. ''' sum_zi2_ci = salt_concentration * np.sum(ion_valences**2) kappa = np.sqrt(F**2/(permittivity*R*T)*sum_zi2_ci) RT_ln_y = -F**2 /(4*np.pi*permittivity ) *(kappa*ion_valences[0]**2)/(2*N_A) ln_y = RT_ln_y/(R*T) return ln_y def EDH(radii,permittivity,salt_concentration,ion_valences,T): ''' extended Debye Hückel ;implementation based on: Erich Hückel, Zur Theorie konzentrierter wässriger Lösungen starker Elektrolyte., Physikalische Zeitschrift 1925 (1925) 93–147 ''' sum_zi2_ci = salt_concentration * np.sum(ion_valences**2) r_a = radii[1] r_c = radii[0] r_ca = (r_c+r_a)/2 kappa = np.sqrt(F**2/(permittivity*R*T)*sum_zi2_ci) RT_ln_y = -F**2 /(4*np.pi*permittivity ) *(ion_valences[0]**2*kappa)/(2*N_A*(1+kappa*r_ca)) ln_y = RT_ln_y/(R*T) return ln_y def DHFULL(radii,permittivity,salt_concentration,ion_valences,T): ''' full Debye Hückel ;implementation based on: P. Debye, E. Hückel, Zur Theorie der Elektrolyte., Physikalische Zeitschrift 1923 (1923) 185–206. ''' sum_zi2_ci = salt_concentration * np.sum(ion_valences**2) r_a = radii[1] r_c = radii[0] kappa = np.sqrt(F**2/(permittivity*R*T)*sum_zi2_ci) sigma_a = 3/(kappa*r_a)**3 *(1+kappa*r_a-1/(1+kappa*r_a)-2*np.log(1+kappa*r_a)) sigma_c = 3/(kappa*r_c)**3 *(1+kappa*r_c-1/(1+kappa*r_c)-2*np.log(1+kappa*r_c)) sumnz2sig = salt_concentration*(ion_valences[0]**2*sigma_c+ion_valences[1]**2*sigma_a) sumnz2 = salt_concentration*np.sum(ion_valences**2) Xi_c = 3/(kappa*r_c)**3*(3/2+np.log(1+kappa*r_c)-2*(1+kappa*r_c)+(1+kappa*r_c)**2/2) Xi_a = 3/(kappa*r_a)**3*(3/2+np.log(1+kappa*r_a)-2*(1+kappa*r_a)+(1+kappa*r_a)**2/2) RT_ln_y_c = - F**2 /(4*np.pi*permittivity ) *(ion_valences[0]**2*kappa)/(6*N_A)*(2*Xi_c+sumnz2sig/sumnz2) RT_ln_y_a = - F**2 /(4*np.pi*permittivity ) *(ion_valences[0]**2*kappa)/(6*N_A)*(2*Xi_a+sumnz2sig/sumnz2) ln_y = (RT_ln_y_a+RT_ln_y_c)/(2*R*T) return ln_y def BornTerm(radii,epsilon_r0,salt_concentration,ion_valences,T): ''' Born-0 ;implementation based on: G.M. Silva, X. Liang, G.M. Kontogeorgis, How to account for the concentration dependency of relative permittivity in the Debye–Hückel and Born equations, Fluid Phase Equilibria 566 (2023) 113671. https://doi.org/10.1016/j.fluid.2022.113671 ''' epsilon_r = calculate_dielectric_constant(epsilon_r0,salt_concentration) r_c,r_a = radii RT_ln_y_a = F**2*ion_valences[1]**2/(8*np.pi*epsilon_0*r_a*N_A)*(1/epsilon_r-1/epsilon_r0) RT_ln_y_c = F**2*ion_valences[0]**2/(8*np.pi*epsilon_0*r_c*N_A)*(1/epsilon_r-1/epsilon_r0) return (RT_ln_y_a+RT_ln_y_c)/(2*R*T) def DHFULLBORN(radii,ion_valences,epsilon_r0,salt_concentration,T=298.15): # usage of concentration dependence of permittivity for DH theory + Born epsilon_r = calculate_dielectric_constant(epsilon_r0,salt_concentration) permittivity = np.array([epsilon_r*epsilon_0]) ln_yDH = DHFULL(radii,permittivity,salt_concentration,ion_valences,T) ln_yBorn = BornTerm(radii,epsilon_r0,salt_concentration,ion_valences,T) ln_y = ln_yDH+ln_yBorn return ln_y ############################################################################# ## functions for Mass action law. def DHLL_MAL(r_c,r_a,salt_concentration,T=298.15): # simplified function for DHLL epsilon_r = 78.14 ion_valences = np.array([1,-1]) radii = r_c,r_a permittivity = np.array([epsilon_r*epsilon_0]) gamma = DHLL(radii,permittivity,salt_concentration,ion_valences,T) return gamma def EDH_MAL(r_c,r_a,salt_concentration,T=298.15): epsilon_r = 78.14 ion_valences = np.array([1,-1]) radii = r_c,r_a permittivity = np.array([epsilon_r*epsilon_0]) gamma = EDH(radii,permittivity,salt_concentration,ion_valences,T) return gamma def Huckel(r_c,r_a,c,C=0.2): #C=0.1#10.5e-2 gammaEl = EDH_MAL(r_c,r_a,c)# return gammaEl+2*C*c/1000 ############################################################################## # Conductivity def DHO3(T,epsilon,ion_valences,viscosity,salt_concentration,conductivities,diameters): ''' implementation of DHO3 from S. Naseri Boroujeni, B. Maribo-Mogensen, X. Liang, G.M. Kontogeorgis, New Electrical Conductivity Model for Electrolyte Solutions Based on the Debye-Hückel-Onsager Theory, J. Phys. Chem. B 127 (2023) 9954–9975. https://doi.org/10.1021/acs.jpcb.3c03381. ''' cond_c,cond_a = conductivities diam_c,diam_a = diameters #lambdai0 =21 kappa = np.sqrt(F**2/(epsilon*R*T)*2*salt_concentration) Tetta = e**2/(epsilon*4*np.pi*kB*T) q = 0.5 diam = (diam_c+diam_a)/2 zeta_c = F**2*abs(ion_valences[0])/(6*np.pi*viscosity*N_A*cond_c) zeta_a = F**2*abs(ion_valences[1])/(6*np.pi*viscosity*N_A*cond_a) dnubynu_c = -zeta_c*(kappa)/(1+kappa*diam) dnubynu_a = -zeta_a*(kappa)/(1+kappa*diam) dkbyk = -Tetta*q/(1+np.sqrt(q)) * (kappa)/(3*(1+kappa*diam)*(1+kappa*diam*np.sqrt(q)+kappa**2*diam**2/6)) lambda_c = (1+dnubynu_c)*(1+dkbyk) lambda_a = (1+dnubynu_a)*(1+dkbyk) Lambda = cond_c*lambda_c+lambda_a*cond_a return Lambda def model1(T,epsilon,ion_valences,viscosity,salt_concentration,conductivities,diameters): ''' implementation of model1 based on S. Naseri Boroujeni, B. Maribo-Mogensen, X. Liang, G.M. Kontogeorgis, New Electrical Conductivity Model for Electrolyte Solutions Based on the Debye-Hückel-Onsager Theory, J. Phys. Chem. B 127 (2023) 9954–9975. https://doi.org/10.1021/acs.jpcb.3c03381. ''' cond_c,cond_a = conductivities diam_c,diam_a = diameters diam = (diam_c+diam_a)/2 zeta_c = F**2*abs(ion_valences[0])/(6*np.pi*viscosity*N_A*cond_c) zeta_a = F**2*abs(ion_valences[1])/(6*np.pi*viscosity*N_A*cond_a) #lambdai0 =21 kappa = np.sqrt(F**2/(epsilon*R*T)*2*salt_concentration) kappaq = np.sqrt((salt_concentration*F*e)/(epsilon*kB*T)) Tetta = e**2/(epsilon*4*np.pi*kB*T) dnubynu_c = -zeta_c*(kappa)/(1+kappa*diam_c) dnubynu_a = -zeta_a*(kappa)/(1+kappa*diam_a) OMEGA = kappaq**2 /3 * (np.sinh(kappaq*diam)/(kappaq*diam)-kappaq*diam**2/Tetta*( np.cosh(kappaq*diam)/(kappaq*diam) -(np.sinh(kappaq*diam)/(kappaq**2*diam**2)) )) dkbyk_c = -OMEGA*Tetta*np.exp(kappa*(diam_c-diam))*np.exp(-kappaq*diam)/((kappa+kappaq)*(1+kappa*diam_c)) dkbyk_a = -OMEGA*Tetta*np.exp(kappa*(diam_a-diam))*np.exp(-kappaq*diam)/((kappa+kappaq)*(1+kappa*diam_a)) lambda2_c = (1+dnubynu_c)*(1+dkbyk_c) lambda2_a = (1+dnubynu_a)*(1+dkbyk_a) Lambda2 = cond_c*lambda2_c+lambda2_a*cond_a return Lambda2 def DHOLL(lam0,ion_valences,salt_concentration,epsi,q,viscosity,T=298.15): #implemented for 25 °C IonicStr = salt_concentration/2000*(ion_valences**2).sum() cond_DHO = lam0 - ( (2.801e6*q*lam0)/((epsi*T)**(3/2)*(1+np.sqrt(q))) + (41.25*(abs(ion_valences).sum()))/(viscosity*10*np.sqrt(epsi*T)))*np.sqrt(IonicStr) return cond_DHO def DHO_stokes(lam0,salt_concentration,epsi,q,viscosity,diameters,T): #lam0 = conductivity[1][0] ion_valences = np.array([1,-1]) kappa = np.sqrt(F**2/(epsilon*R*T)*2*salt_concentration) diam = (diameters[0]+diameters[1])/2 IonicStr = salt_concentration/2000*(ion_valences**2).sum() cond_DHO_stokes = lam0 - ( (2.801e6*q*lam0)/((epsi*T)**(3/2)*(1+np.sqrt(q))) + (41.25*(abs(ion_valences).sum()))/(viscosity*10*np.sqrt(epsi*T)))*np.sqrt(IonicStr)/(1+kappa*diam) return cond_DHO_stokes ########################################################################## ########################################################################## # Diagrams for NaCl for exemplary visualization. ########################################################################## ########################################################################## if __name__ == '__main__': pauling_diameter = { 'Li+': 1.20 , 'Na+': 1.90 , 'K+': 2.66 , 'F-': 2.72 , 'Cl-': 3.62 , 'Br−': 3.90 } ########################################################################## # DATA for NaCl measurement_concentrations, measurement_y = get_nacl_activity_data() T =298.15 # K c_list = np.logspace(-2, 3.73, 700 ) #mol/m^3 epsilon_r_0 = 78.14 epsilon_r_0_max = 95 epsilon = epsilon_r_0*epsilon_0 epsilon_max = epsilon_r_0_max*epsilon_0 ions = ['Na+','Cl-'] ion_valences = np.array([1,-1]) lam0 = 126.5 #conductivity NaCl q = 0.5 # for 1-1 electrolytes diameter_water = 2.75e-10 r_cmin = pauling_diameter[ions[0]]*1e-10/2 + .0*diameter_water r_cmax = pauling_diameter[ions[0]]*1e-10/2 + 1.05*diameter_water r_amin = pauling_diameter[ions[1]]*1e-10/2 + .0*diameter_water r_amax = pauling_diameter[ions[1]]*1e-10/2 + .65*diameter_water ########################################################################## # Data Driesner: cges_D = 0.985 #mol/l c_Driesner = np.array([cges_D,cges_D]) alpha_Driesner =np.array([ 0.633+2/3*0.185, 0.633]) alpha_Bjerrum = np.ones_like(c_list) c_Alejandre = np.array([0.494,2.375]) #J. Chem. Phys. 130, 174505 (2009), p8 alpha_Alejandre = np.array([0.99,0.55]) ########################################################################## # Evaluate models MSA_Max, HSMax = mean_spherical_approximation([r_cmax,r_amax], epsilon,c_list,ion_valences,T) MSA_Min, HSMin = mean_spherical_approximation([r_cmin,r_amin], epsilon,c_list,ion_valences,T) DHLLMax = DHLL([r_cmin,r_amin], epsilon,c_list,ion_valences,T) DHLLmin = DHLL([r_cmin,r_amin], epsilon,c_list,ion_valences,T) EDHMax = EDH([r_cmax,r_amax], epsilon,c_list,ion_valences,T) EDHmin = EDH([r_cmin,r_amin], epsilon,c_list,ion_valences,T) DHFULLMax = DHFULL([r_cmax,r_amax],epsilon,c_list,ion_valences,T) DHFULLmin = DHFULL([r_cmin,r_amin],epsilon,c_list,ion_valences,T) BornMax = BornTerm([r_cmax,r_amax],epsilon_r_0,c_list,ion_valences,T) Bornmin = BornTerm([r_cmin,r_amin],epsilon_r_0,c_list,ion_valences,T) AllMax= DHFULLBORN([r_cmax,r_amax],ion_valences,epsilon_r_0,c_list) Allmin= DHFULLBORN([r_cmin,r_amin],ion_valences,epsilon_r_0,c_list) huckel = Huckel(2.35e-10,2.35e-10,c_list,C=0.1) gammaBiMSA,alphaBiMSA= BiMSA([0.95e-10,1.81e-10],c_list) rc = 3.2e-10 ra = 2.8e-10 FULL= DHFULL([rc,ra],epsilon,c_list,ion_valences,T) ct,HS= mean_spherical_approximation([rc,ra],epsilon,c_list,ion_valences,T) rc = 3.2e-10 ra = 2.151e-10 FULLe= DHFULL([rc,ra],epsilon_max,c_list,ion_valences,T) ct,HSe= mean_spherical_approximation([rc,ra],epsilon_max,c_list,ion_valences,T) rc = 3.4e-10 ra = 2.151e-10 aasa,HS2 = mean_spherical_approximation([rc,ra],epsilon,c_list,ion_valences,T) BornL = DHFULLBORN([rc,ra],ion_valences,epsilon_r_0,c_list) #eps variation rc=pauling_diameter[ions[0]]*1e-10/2 ra=pauling_diameter[ions[1]]*1e-10/2 DHLLMaxe = DHLL([rc,ra], epsilon_max,c_list,ion_valences,T) DHLLmine = DHLL([rc,ra], epsilon,c_list,ion_valences,T) EDHMaxe = EDH([rc,ra], epsilon_max,c_list,ion_valences,T) EDHmine = EDH([rc,ra], epsilon,c_list,ion_valences,T) DHFULLMaxe = DHFULL([rc,ra],epsilon_max,c_list,ion_valences,T) DHFULLmine = DHFULL([rc,ra],epsilon,c_list,ion_valences,T) BornMaxe = BornTerm([rc,ra],epsilon_r_0_max,c_list,ion_valences,T) Bornmine = BornTerm([rc,ra],epsilon_r_0,c_list,ion_valences,T) rel_viscosity = np.array([ 1.000, 0.08728783 ,0.00844365, 0.00124189, 0])#polinomial A0, A1, A2, A3, A4 = rel_viscosity ####################### # conductivity # data from (and sources in): Thomas W. Chapman and John Newman, A compilation of selected thermodynamic and transport properties of binary electrolytes in aqueous solution, 1968. conductivity = [np.array([5.9441e-05, 1.1280e-04, 1.8980e-04, 4.2677e-04, 5.0050e-04, 6.7890e-04, 6.9190e-04, 7.1870e-04, 7.3840e-04, 9.2240e-04, 1.3440e-03, 1.3868e-03, 1.4466e-03, 1.5306e-03, 1.6271e-03, 2.0348e-03, 2.1253e-03, 2.2275e-03, 2.3388e-03, 2.6538e-03, 2.8988e-03, 2.9806e-03, 3.1862e-03, 3.7361e-03, 3.8776e-03, 1.0000e-03, 2.0000e-03, 5.0000e-03, 1.0000e-02, 2.0000e-02, 5.0000e-02, 1.0000e-01, 1.1000e-01, 1.2000e-01, 1.9000e-01, 2.6000e-01, 3.6000e-01, 4.7000e-01, 6.0000e-01, 6.6000e-01, 7.4000e-01, 8.6000e-01, 1.2300e+00, 5.354,4.5199,3.5037,3.1431,2.7439,2.3793,1.9279,1.5193]) , np.array([125.8 , 125.6 , 125.1 , 124.6 , 124.4 , 124.2 , 124.2 , 124. , 124. , 123.7 , 123.3 , 123.1 , 123.2 , 123.2 , 122.9 , 122.5 , 122. , 122.4 , 122.4 , 122. , 121.88, 121.9 , 121.7 , 121.2 , 121.2 , 124.72, 123.67, 121.58, 119.43, 116.64, 111.88, 107.67, 106. , 106. , 103. , 100. , 97. , 95. , 93. , 91. , 90. , 88. , 83. , 46.83,53.13,61.3,64.36,67.8,71.14,75.41,79.65,]) ] viscosity= 0.890e-3 #Pa*s viscosity H20 cond_Na = 50.08e-4 # Sm^2/mol cond_Cl = 76.31e-4# Sm^2/mol conductivities = cond_Na,cond_Cl c_lis=c_list/1000 viscosity_frac = A0/(A0 + A1 * c_lis + A2 * c_lis**2 + A3 * c_lis**3 + A4 * c_lis**4) viscosity_c = viscosity*viscosity_frac diametersmin = 2*np.array([r_amin,r_cmin]) diametersmax = 2*np.array([r_amax,r_cmax]) condDHO= DHOLL(lam0,ion_valences,c_list,epsilon_r_0,q,viscosity,T) condDHOe= DHOLL(lam0,ion_valences,c_list,epsilon_r_0_max,q,viscosity,T) # viscosity c dependent Lambda_DHO3_min = DHO3(T,epsilon,ion_valences,viscosity_c,c_list,conductivities,diametersmin) Lambda_DHO3_mine = DHO3(T,epsilon_max,ion_valences,viscosity_c,c_list,conductivities,diametersmin) Lambda_DHO3_max = DHO3(T,epsilon,ion_valences,viscosity_c,c_list,conductivities,diametersmax) # viscosity c const Lambda_DHO3_min_e0 = DHO3(T,epsilon,ion_valences,viscosity,c_list,conductivities,diametersmin) Lambda_DHO3_min_e0e = DHO3(T,epsilon_max,ion_valences,viscosity,c_list,conductivities,diametersmin) Lambda_DHO3_max_e0 = DHO3(T,epsilon,ion_valences,viscosity,c_list,conductivities,diametersmax) # viscosity c dependent condDHO_Stokes_min = DHO_stokes(lam0,c_list,epsilon_r_0,q,viscosity_c,diametersmin,T) condDHO_Stokes_mine = DHO_stokes(lam0,c_list,epsilon_r_0_max,q,viscosity_c,diametersmin,T) condDHO_Stokes_max = DHO_stokes(lam0,c_list,epsilon_r_0,q,viscosity_c,diametersmax,T) ######################################################################## ######################################################################## # plot properties matplotlib.rc('xtick', labelsize=18) matplotlib.rc('ytick', labelsize=18) font1 = {'family' : 'Arial', 'weight' : 'normal', 'size' : 18} matplotlib.rc('font', **font1) matplotlib.rc('lines', linewidth=3) matplotlib.rcParams["figure.figsize"] = (7.5,5) matplotlib.rcParams["figure.dpi"] = 600 from matplotlib.ticker import ScalarFormatter class ScalarFormatterClass(ScalarFormatter): def _set_format(self): self.format = "%1.1f" plotlimit = 6.14*1.202/(1+0.05844*6.14) #solubility limit ############################# ############################# Plot 3A: electrostatic contr. alpha = .7 matplotlib.rcParams['hatch.linewidth'] = 2 # previous pdf hatch linewidth # Fig2 fig = plt.figure() ax1 = fig.add_subplot(111) ax1.scatter(np.sqrt(measurement_concentrations), np.log(measurement_y), label = 'measurement',color='purple',zorder=2) ax1.fill_between(np.sqrt(c_list/1000), DHLLMax, DHLLMax, alpha=alpha,color='none', edgecolor='y',hatch='--', linewidth=4,label = 'DHLL',zorder=1) ax1.fill_between(np.sqrt(c_list/1000), EDHmin,EDHMax, alpha=alpha, color='none',edgecolor='b', hatch='||', linewidth=4,label = 'EDH',zorder=1) ax1.fill_between(np.sqrt(c_list/1000), DHFULLmin,DHFULLMax, alpha=alpha, color='none',edgecolor='r', hatch='//', linewidth=4,label = 'DHFULL',zorder=1) ax1.fill_between(np.sqrt(c_list/1000), MSA_Min,MSA_Max, alpha=alpha, color='none',edgecolor='g', hatch = r'\\', linewidth=4,label = 'MSA',zorder=1) plt.fill_between(np.sqrt(c_list/1000), Allmin[0]-Bornmin, AllMax[0]-BornMax, alpha=alpha, color= 'none',edgecolor = 'indianred', linewidth=4,label = r'DHFULL,$\epsilon_r (c)$',hatch = '--',zorder=1) ax1.set_xlim(0, max(np.sqrt(c_list/1000))) ax1.set_ylim(-1.49,1.5) ax1.grid('both') ax1.set_xticklabels([]) new_tick_locations = np.array([0,.5, 1, 1.5,2.0]) def tick_function(X): V = X**2 return ["%.2f" % z for z in V] ax2 = ax1.twiny() ax2.set_xlim(ax1.get_xlim()) ax2.set_xticks(new_tick_locations) ax2.set_xticklabels(tick_function(new_tick_locations)) ax2.set_xlabel(r'$\it{c}$ ($\mathrm{M}$)') yScalarFormatter = ScalarFormatterClass(useMathText=True) ax1.yaxis.set_major_formatter(yScalarFormatter) ax1.set_ylabel(r'ln$\it{(y_\pm^{\mathrm{el}})}$') ax1.legend(ncol=2,loc='upper left')#,edgecolor= '0.8',framealpha=0.2) plt.savefig("figures/myImagePDF.pdf", format="pdf", bbox_inches="tight") plt.show() ############################# Plot 3B: HS contr. alpha = .3 fig = plt.figure() ax1 = fig.add_subplot(111) ax1.yaxis.set_major_formatter(yScalarFormatter) ax1.scatter(np.sqrt(measurement_concentrations), np.log(measurement_y),label = 'measurement', color='purple') ax1.fill_between(np.sqrt(c_list/1000), HSMin, HSMax, alpha=alpha, color='dimgray', linewidth=4,label = 'HS (MSA)') ax1.set_xlim(0, max(np.sqrt(c_list/1000))) ax1.set_ylim(-1,2) ax1.legend(loc='upper left') ax1.grid() ax1.set_ylabel(r'$\mathrm{ln(}y_\pm^{\mathrm{HS}}\mathrm{)}$') ax1.set_xlabel(r'$\it{\sqrt{c}}$ $\it{(\sqrt{\mathrm{M}})}$') plt.show() ################## Plot 3C: Born contr. fig = plt.figure() ax1 = fig.add_subplot(111) ax1.scatter(np.sqrt(measurement_concentrations), np.log(measurement_y),label = 'measurement', color='purple') plt.fill_between(np.sqrt(c_list/1000), Bornmin,BornMax, alpha=alpha, color='darkgrey', linewidth=4,label = 'Born') ax1.set_xlim(0, max(np.sqrt(c_list/1000))) ax1.set_ylim(-.999,2) ax1.grid('both') ax1.set_xticklabels([])#axes.get_xaxis().set_visible(False) ax2 = ax1.twiny() ax2.set_xlim(ax1.get_xlim()) ax2.set_xticks(new_tick_locations) ax2.set_xticklabels(tick_function(new_tick_locations)) ax2.set_xlabel(r'$\it{c}$ ($\mathrm{M}$)') yScalarFormatter = ScalarFormatterClass(useMathText=True) ax1.yaxis.set_major_formatter(yScalarFormatter) ax1.set_ylabel(r'$\mathrm{ln}(y_\pm^{\mathrm{sol}})$') ax1.legend(loc='upper left')#ncol=2,loc='upper left')#,edgecolor= '0.8',framealpha=0.2) plt.show() ################## Plot 3D: DHFULL+HS+Born contr. plt.scatter(np.sqrt(measurement_concentrations), np.log(measurement_y), label = 'measurement',color='purple',zorder=2) plt.fill_between(np.sqrt(c_list/1000), Allmin[0]+HSMin, AllMax[0]+HSMax, alpha=alpha, color='indianred', linewidth=4,label = r'DHFULL,$\epsilon_r(c)$+HS+Born',hatch = '--',zorder=1) plt.xlim(0, np.sqrt(plotlimit)) plt.ylim(-1.,2) plt.legend(loc='upper left') plt.grid() plt.ylabel(r'$\mathrm{ln}(y_\pm)$') plt.xlabel(r'$\it{\sqrt{c}}$ $(\sqrt{\mathrm{M}})$') plt.show() ################## CONDUCTIVITY fig = plt.figure() ax1 = fig.add_subplot(111) ax2 = ax1.twiny() ax1.scatter(np.sqrt(conductivity[0]),conductivity[1],color = 'purple', label= 'measurement data',zorder=3) ax1.plot(np.sqrt(c_list/1000),condDHO, label = r'DHO-LL, $\eta_\mathrm{0}$',zorder=2,color='y') ax1.fill_between(np.sqrt(c_list/1000), condDHO_Stokes_min, condDHO_Stokes_max, alpha=.7, color='none',edgecolor='b', hatch='/',linewidth=4,label = r"DHO, $\eta_\mathrm{exp}$",zorder=1) ax1.fill_between(np.sqrt(c_list/1000), Lambda_DHO3_min_e0*1e4, Lambda_DHO3_max_e0*1e4, alpha=.7, color='none',edgecolor='chocolate', hatch='|',linewidth=4,label = r"DHO3, $\eta_\mathrm{0}$",zorder=1) ax1.fill_between(np.sqrt(c_list/1000), Lambda_DHO3_min*1e4, Lambda_DHO3_max*1e4, alpha=.7, color='none',edgecolor='maroon', hatch='\\',linewidth=4,label = r"DHO3, $\eta_\mathrm{exp}$",zorder=1) ax2.set_xlim(ax1.get_xlim()) ax2.set_xticks(new_tick_locations) ax2.set_xticklabels(tick_function(new_tick_locations)) ax1.set_xticklabels([]) ax1.set_xlim(0,np.sqrt(plotlimit)) ax1.set_ylim(0,) ax1.grid() ax1.set_ylabel(r'$\Lambda$ $\mathrm{(Scm^2/mol)}$') ax2.set_xlabel(r'$\it{c}$ ($\mathrm{M}$)') #plt.legend(bbox_to_anchor=(1.62, 1)) plt.show() ################## Conductivity EPS plt.scatter(np.sqrt(conductivity[0]),conductivity[1],color = 'purple', label= 'measurement data',zorder=3) plt.fill_between(np.sqrt(c_list/1000), condDHO, condDHOe, alpha=.7, color='none',edgecolor='y', hatch='--',linewidth=4,label = r"DHOLL, $\eta_0$",zorder=1) plt.fill_between(np.sqrt(c_list/1000), condDHO_Stokes_min, condDHO_Stokes_mine, alpha=.7, color='none',edgecolor='b', hatch='/',linewidth=4,label = r"DHO, $\eta_\mathrm{exp}$",zorder=1) plt.fill_between(np.sqrt(c_list/1000), Lambda_DHO3_min_e0*1e4, Lambda_DHO3_min_e0e*1e4, alpha=.7, color='none',edgecolor='chocolate', hatch='|',linewidth=4,label = r"DHO3, $\eta_\mathrm{0}$",zorder=1) plt.fill_between(np.sqrt(c_list/1000), Lambda_DHO3_min*1e4, Lambda_DHO3_mine*1e4, alpha=.7, color='none',edgecolor='maroon', hatch='\\',linewidth=4,label = r"DHO3, $\eta_\mathrm{exp}$",zorder=1) plt.xlim(0,np.sqrt(plotlimit)) plt.ylim(0,) plt.grid() plt.ylabel(r'$\Lambda$ $\mathrm{(Scm^2/mol)}$') plt.xlabel(r'$\it{\sqrt{c}}$ ($\mathrm{\sqrt{M}}$)') plt.legend(bbox_to_anchor=(1.62, 1)) plt.show() ###### ASSOCIATION: Fraction of free ions funcs = [DHLL_MAL,Huckel,Huckel] alpha_singles = np.zeros([len(c_list),len(funcs)]) alpha_allions = np.zeros([len(c_list),len(funcs)]) gamma_lists = np.zeros([len(c_list),len(funcs)]) rs= [(r_cmin,r_amin),(r_cmin,r_amin),(r_cmax,r_amax)] for i,func in enumerate(funcs): rss=rs[i] alpha_single_ions, alpha_all_ions,gamma_list = compute_free_ions(func,c_list,rss) alpha_singles[:,i] = alpha_single_ions alpha_allions[:,i] = alpha_all_ions gamma_lists[:,i] = gamma_list plt.show() fig = plt.figure() ax1 = fig.add_subplot(111) ax2 = ax1.twiny() ax1.set_ylabel(r'$\alpha$') ax1.set_xlim(0,np.sqrt(plotlimit)) ax1.set_ylim(0.001,1.05) #ax1.set_xlabel(r'$\it{\sqrt{c}}$ ($\it{\sqrt{M}}$)') ax2.set_xlabel(r'$\it{c}$ ($\mathrm{M}$)') #c_Alejandre = np.array([0.494,2.375]) #J. Chem. Phys. 130, 174505 (2009), p8 #alpha_Alejandre = np.array([0.99,0.55]) markers= ["x","o"] ax1.plot(np.sqrt(c_list/1000),alpha_Bjerrum,color = 'teal', label = r'Bjerrum',zorder=1) ax1.plot(np.sqrt(c_list/1000),alphaBiMSA,color='g', label = r'BiMSA',zorder=1) ax1.plot(np.sqrt(c_list/1000),alpha_singles[:,0],'--',color='y',zorder=1) ax1.plot(np.sqrt(c_list/1000),alpha_allions[:,0],color='y', label = r'MAL(DHLL)',zorder=1) ax1.fill_between(np.sqrt(c_list/1000), alpha_allions[:,1], alpha_allions[:,2], alpha=alpha*1.5,color='none', edgecolor='lightblue',hatch='\\', linewidth=4,label = 'MAL (EDH-H)',zorder=2) ax1.fill_between(np.sqrt(c_list/1000), alpha_singles[:,1], alpha_singles[:,2], alpha=alpha*1.5,color='none', edgecolor='lightblue',hatch='\\', linewidth=4,zorder=2, linestyle='--') ax1.scatter(np.sqrt(c_Driesner),alpha_Driesner,color = 'deeppink', label= 'MD Driesner et al.',zorder=2) ax1.scatter(np.sqrt(c_Alejandre[1]),alpha_Alejandre[1],color = 'purple',marker = markers[1],zorder=2, label= 'MD Alejandre et al.', s = 40) ax1.scatter(np.sqrt(c_Alejandre[0]),alpha_Alejandre[0],color = 'purple',marker = markers[0],zorder=2) ax2.set_xlim(ax1.get_xlim()) ax2.set_xticks(new_tick_locations) ax2.set_xticklabels(tick_function(new_tick_locations)) ax1.grid() ax1.legend(bbox_to_anchor=(1.02,0.94)) ax1.set_xticklabels([])#axes.get_xaxis().set_visible(False) plt.show() plt.plot(np.sqrt(c_list/1000), gammaBiMSA, linewidth=4,label = 'BiMSA',color = 'green',zorder=1) plt.plot(np.sqrt(c_list/1000), np.log(gamma_lists[:,1]),linewidth=4,label = 'EDH-H+MAL', color = 'lightblue',zorder=2) plt.plot(np.sqrt(c_list/1000), huckel,linewidth=4,label = r'EDH-H', color = 'blue',zorder=1) plt.plot(np.sqrt(c_list/1000), FULL+HS,linewidth=4,label = r'DHFULL+HS',color = 'red',zorder=1) plt.plot(np.sqrt(c_list/1000), BornL[0]+HS2,linewidth=4,label = r'DHFULL,$\epsilon_r (c)$+Born+HS',color = 'firebrick',zorder=2) plt.plot(np.sqrt(c_list/1000), DHLLMaxe+HSe+BornMaxe,linewidth=4,label = r'DHLL+HS+Born, $\epsilon_w=95$',color = 'k',zorder=1) plt.scatter(np.sqrt(measurement_concentrations), np.log(measurement_y), label = 'measurement',color='purple',zorder=3) plt.xlim(0, np.sqrt(plotlimit)) plt.ylim(-.8,.499) plt.legend(ncols=1,bbox_to_anchor=(1.02,0.94)) plt.grid() plt.ylabel(r'ln($\it{y_\pm}$)') plt.xlabel(r'$\it{\sqrt{c}}$ $\mathrm{(\sqrt{M})}$') plt.show() ###EPS Variation alpha=0.6 fig = plt.figure() ax1 = fig.add_subplot(111) ax1.scatter(np.sqrt(measurement_concentrations), np.log(measurement_y), label = 'measurement',color='purple',zorder=3) ax1.fill_between(np.sqrt(c_list/1000), DHLLmine, DHLLMaxe, alpha=alpha,color='none', edgecolor='y',hatch='\\', linewidth=4,label = 'DHLL',zorder=2) ax1.fill_between(np.sqrt(c_list/1000), EDHmine,EDHMaxe, alpha=alpha, color='none',edgecolor='b', hatch='||', linewidth=4,label = 'EDH',zorder=1) ax1.fill_between(np.sqrt(c_list/1000), DHFULLmine,DHFULLMaxe, alpha=alpha, color='none',edgecolor='r', hatch='//', linewidth=4,label = 'DHFULL',zorder=1) ax1.set_xlim(0, 1)#max(np.sqrt(c_list/1000))) ax1.set_ylim(-.79,.4) ax1.grid('both') #ax1.set_xticklabels([]) new_tick_locations = np.array([0,.2,.4,.6,.8, 1]) def tick_function(X): V = X**2 return ["%.2f" % z for z in V] ax2 = ax1.twiny() ax2.set_xlim(ax1.get_xlim()) ax2.set_xticks(new_tick_locations) ax2.set_xticklabels(tick_function(new_tick_locations)) ax1.set_xlabel(r'$\it{\sqrt{c}}$ $\it{(\sqrt{\mathrm{M}})}$') ax2.set_xlabel(r'$\it{c}$ ($\mathrm{M}$)') yScalarFormatter = ScalarFormatterClass(useMathText=True) ax1.yaxis.set_major_formatter(yScalarFormatter) ax1.set_ylabel(r'ln$\it{(y_\pm^{\mathrm{el}})}$') ax1.legend(ncol=2,loc='upper left')#,edgecolor= '0.8',framealpha=0.2) plt.show() fig = plt.figure() ax1 = fig.add_subplot(111) ax1.scatter(np.sqrt(measurement_concentrations), np.log(measurement_y),label = 'measurement', color='purple') plt.fill_between(np.sqrt(c_list/1000), Bornmine,BornMaxe, alpha=alpha, color='darkgrey', linewidth=4,label = 'Born') ax1.set_xlim(0, max(np.sqrt(c_list/1000))) ax1.set_ylim(-.999,2) ax1.grid('both') new_tick_locations = np.array([0,.5, 1, 1.5,2.0]) #ax1.set_xticklabels([])#axes.get_xaxis().set_visible(False) ax2 = ax1.twiny() ax2.set_xlim(ax1.get_xlim()) ax2.set_xticks(new_tick_locations) ax2.set_xticklabels(tick_function(new_tick_locations)) ax1.set_xlabel(r'$\it{\sqrt{c}}$ $\it{(\sqrt{\mathrm{M}})}$') ax2.set_xlabel(r'$\it{c}$ ($\mathrm{M}$)') yScalarFormatter = ScalarFormatterClass(useMathText=True) ax1.yaxis.set_major_formatter(yScalarFormatter) ax1.set_ylabel(r'$\mathrm{ln}(y_\pm^{\mathrm{sol}})$') ax1.legend(loc='upper left')#ncol=2,loc='upper left')#,edgecolor= '0.8',framealpha=0.2) plt.show()