Binodal and Spinodal Curves¶
Binodal and spinodal curves are useful for understanding phase stability. The binodal curve (or coexistence curve) defines the temperatures and compositions at which phase separation is thermodynamically favorable. The spinodal curve is located within the binodal curve and indicates the limit of local phase stability. Compositions between the spinodal and binodal curves – while not thermodynamically stable – are robust against small fluctuations (i.e., the free energy surface is locally convex for points in this region).
Python code (Binary mixture)¶
[show/hide code]
import os
import numpy as np
import matplotlib.pyplot as plt
from scm.plams import *
################## Note: Be sure to add the path to your own AMSCRS directory here ##################
database_path = os.getcwd()
if not os.path.exists(database_path):
raise OSError(f'The provided path does not exist. Exiting.')
init()
#suppress plams output
config.log.stdout = 0
files = ["Nitrobenzene.coskf","Hexane.coskf"]
nring = [6,0]
# problem paramters
temp_range = [203.15,243.15] #K -- the temperature range over which the curves are calculated
steps = 40# number of steps to take within the temperature range
def binmix_at_T(files,temp):
# initialize settings object
settings = Settings()
settings.input.property._h = 'BINMIXCOEF'
settings.input.property.nfrac = 100
# make compounds
compounds = [Settings() for i in range(len(files))]
for i,(file,nr) in enumerate(zip(files,nring)):
compounds[i]._h = os.path.join( database_path, file )
compounds[i].nring = nr
settings.input.temperature = temp
# optionally, change to the COSMOSAC2013 method
settings.input.method = 'COSMOSAC2013'
# we'll also tighten the convergence threshold for better numerical accuracy
settings.input.Technical.sacconv = 1e-10
# add the compounds to the settings object
settings.input.compound = compounds
# create a job that can be run by COSMO-RS
my_job = CRSJob(settings=settings)
# run the job
out = my_job.run()
# convert all the results into a python dict
res = out.get_results()
return res
def calc_binodal_at_T(res):
# this is the miscibility gap, so we can use the result calculated by the program
if res["showmiscgap"]:
return res["xlle"][:2]
else:
return None
def calc_spinodal_at_T(res):
# here, we'll look for points with d^2(G_mix)/dx^2 = 0
# we'll calculate a numerical second derivative for every point
spinodal = []
gmix = res["Gibbs energy of mixing"]
frac1 = res["molar fraction"][0]
second_deriv = np.zeros(len(gmix))
# initial values for endpoints (assuming convexity close to pure compounds)
second_deriv[0] = 0.0001
second_deriv[-1] = 0.0001
for i in range(1,len(gmix)-1):
delta1 = frac1[i]-frac1[i-1]
delta2 = frac1[i+1]-frac1[i]
d1 = (gmix[i]-gmix[i-1])/delta1
d2 = (gmix[i+1]-gmix[i])/delta2
second_deriv[i] = 2*(d2-d1)/(delta1+delta2)
for i in range(len(second_deriv)-1):
if second_deriv[i]*second_deriv[i+1] < 0:
dist1 = abs(second_deriv[i])
dist2 = abs(second_deriv[i+1])
tot = dist1+dist2
zero = (dist2*frac1[i]+dist1*frac1[i+1])/tot
spinodal.append(zero)
return spinodal if spinodal else None
temps = [temp_range[0] + (temp_range[1]-temp_range[0])/steps * i for i in range(steps+1)]
bin_left_points = []
bin_right_points = []
spin_left_points = []
spin_right_points = []
print("Temperature".ljust(15),"Binodal points".ljust(25),"Spinodal points")
for temp in temps:
res = binmix_at_T(files,temp)
binodal = calc_binodal_at_T(res)
spinodal = calc_spinodal_at_T(res)
bin_str = '(' + ",".join([ '{0:<10.5g}'.format(x) for x in binodal]) + ')' if binodal is not None else "--"
spin_str = '(' + ",".join([ '{0:<10.5g}'.format(x) for x in spinodal]) + ')' if spinodal is not None else "--"
print( '{0:.5g}'.format(temp).ljust(15), bin_str.ljust(25) ,spin_str)
if binodal is not None:
bin_left_points.append((binodal[0],temp))
bin_right_points.append((binodal[1],temp))
if spinodal is not None:
spin_left_points.append((spinodal[0],temp))
spin_right_points.append((spinodal[-1],temp))
bin_points = bin_left_points + list(reversed(bin_right_points))
spin_points = spin_left_points + list(reversed(spin_right_points))
plt.plot([x[0] for x in bin_points],[y[1] for y in bin_points],label="Binodal curve")
plt.plot([x[0] for x in spin_points],[y[1] for y in spin_points], label="Spinodal curve")
plt.xlabel("Mole fraction compound 1")
plt.ylabel("Temperature (K)")
plt.legend(loc='upper right')
plt.grid()
plt.show()
finish()
This code produces the following output:
Python code (Ternary mixture)¶
Note
This example uses the python package ternary, but this is only required for plotting. This package can be installed in amspython using pip as follows: amspython -m pip install python-ternary. Users may choose to remove the plotting features of the code and not install ternary.
[show/hide code]
import os
import numpy as np
import matplotlib.pyplot as plt
from scm.plams import *
try:
import ternary
except ImportError:
print ("Cannot find ternary package.")
print ("Try to install with:")
print ("amspython -m pip install python-ternary")
exit()
################## Note: Be sure to add the path to your own AMSCRS directory here ##################
database_path = os.getcwd()
if not os.path.exists(database_path):
raise OSError(f'The provided path does not exist. Exiting.')
init()
#suppress plams output
config.log.stdout = 0
files = ["Water.coskf","Chloroform.coskf","Acetic_acid.coskf"]
nring = [0,0,0]
nfrac = 50 # the nfrac parameter for ternary mixtures
temp_range = [273.15,273.15] #K -- the temperature range over which the curves are calculated
steps = 0 # number of steps to take within the temperature range
def ternmix_at_T(files,temp,nfrac=20):
# initialize settings object
settings = Settings()
settings.input.property._h = 'TERNARYMIX'
settings.input.property.nfrac = nfrac
# make compounds
compounds = [Settings() for i in range(len(files))]
for i,(file,nr) in enumerate(zip(files,nring)):
compounds[i]._h = os.path.join( database_path, file )
compounds[i].nring = nr
settings.input.temperature = temp
# optionally, change to the COSMOSAC2013 method
settings.input.method = 'COSMOSAC2013'
# we'll also tighten the convergence threshold for better numerical accuracy
settings.input.Technical.sacconv = 1e-10
settings.input.Technical.rsconv = 1e-10
# add the compounds to the settings object
settings.input.compound = compounds
# create a job that can be run by COSMO-RS
my_job = CRSJob(settings=settings)
# run the job
out = my_job.run()
# convert all the results into a python dict
res = out.get_results()
return res
def calc_binodal_at_T(res):
# this is the miscibility gap, so we can use the result calculated by the program
points_l = []
points_r = []
if res["nxll"]>0:
for i in range(len(res["xll"])//6):
points_l.append(res["xll"][6*i:6*i+3])
points_r.append(res["xll"][6*i+3:6*i+6])
return points_l + list(reversed(points_r))
else:
return None
def calc_spinodal_at_T(res):
# in this function, we search for mole fraction values where the determinant of the hessian of G_mix with respect to mole fractions is 0
# using the first two mole fractions as degrees of freedom
spinodal = []
gmix_res = res["Gibbs energy of mixing"]
fracs = res["molar fraction"][:2]
step_size = abs(fracs[0][1]-fracs[0][0])
tot_steps = round(1.0/step_size)+1
gmix = np.empty((tot_steps,tot_steps))
gmix.fill(None)
det_mat = np.empty((tot_steps,tot_steps))
det_mat.fill(None)
for x1,x2,gmix_val in zip(fracs[0],fracs[1],gmix_res):
idx1 = int(round(x1/step_size,0))
idx2 = int(round(x2/step_size,0))
gmix[idx1,idx2] = gmix_val
# calculate second derivatives and determinants
for i in range(tot_steps):
for j in range(tot_steps-i):
# these finite difference expressions work because the ternary mixture always has a constant step size
if 0<i<tot_steps-1 and 0<j<tot_steps-1:
d_xx = (gmix[i+1,j]-2*gmix[i,j]+gmix[i-1,j])/(step_size**2)
d_yy = (gmix[i,j+1]-2*gmix[i,j]+gmix[i,j-1])/(step_size**2)
d_xy = (gmix[i+1,j+1]-gmix[i+1,j-1]-gmix[i-1,j+1]+gmix[i-1,j-1])/(4*step_size**2)
det_mat[i,j] = d_xx*d_yy-d_xy**2
for i in range(tot_steps):
for j in range(tot_steps-i):
# compare right and below
for (i1,j1),(i2,j2) in [[(i,j),(i+1,j)],[(i,j),(i,j+1)]]:
if i2 < tot_steps and j2 < tot_steps and not np.isnan(det_mat[i1,j1]) and not np.isnan(det_mat[i2,j2]) and det_mat[i1,j1]*det_mat[i2,j2] < 0:
x1 = step_size*np.array([i1,j1,tot_steps-1-i1-j1])
x2 = step_size*np.array([i2,j2,tot_steps-1-i2-j2])
spin = (abs(det_mat[i2,j2])*x1+abs(det_mat[i1,j1])*x2)/(abs(det_mat[i1,j1])+abs(det_mat[i2,j2]))
spinodal.append(spin)
return spinodal
## Make ternary figure
figure, tax = ternary.figure(scale=1.0)
tax.boundary(linewidth=2.0)
tax.gridlines(color="black", multiple=0.05)
# Set Axis labels and Title
fontsize = 10
tax.bottom_axis_label("$x_1$", fontsize=fontsize,offset=0.2)
tax.right_axis_label( "$x_2$", fontsize=fontsize,offset=0.2)
tax.left_axis_label( "$x_3$", fontsize=fontsize,offset=0.2)
tax.ticks(axis='lbr',multiple=0.1, linewidth=1, tick_formats="%.1f",offset=0.02,fontsize=fontsize)
# Remove default Matplotlib Axes
tax.clear_matplotlib_ticks()
temps = []
if steps > 0:
temps = [temp_range[0] + (temp_range[1]-temp_range[0])/steps * i for i in range(steps+1)]
else:
temps = [temp_range[0]] if isinstance(temps,list) else [temp_range]
for temp in temps:
res = ternmix_at_T(files,temp,nfrac)
binodal = calc_binodal_at_T(res)
spinodal = calc_spinodal_at_T(res)
# the sorting index might need to be changed here if the curves are in a different position
binodal.sort( key = lambda x: x[1])
spinodal.sort(key = lambda x: x[1])
tax.plot(binodal, label="Binodal T="+str(temp)+" K")
tax.plot(spinodal,label="Spinodal T="+str(temp)+" K")
tax.get_axes().axis('off')
ternary.plt.legend(bbox_to_anchor=(1.1,1.0), loc='right')
ternary.plt.show()
finish()
This code produces the following output:
