Eutectic systems¶

import os
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


# the ethanol water system
# experimental numbers for this eutectic are 0.86 mole fraction ethanol at a temperature of -124.3 degrees C
files = ["Ethanol.coskf","Water.coskf"]
tm    = [158.5,273.15] #K
hfus  = [1.2,1.43] #kcal/mol
initial_t_range = [100,300] #K -- the temperature range over which the eutectic search is done
steps           = 20 # number of steps to take within the temperature range

# another eutectic system
# files = ["L-Menthol.coskf","Camphor.coskf"]
# tm    = [316.2,451.5] #K
# hfus  = [2.84,1.63] #kcal/mol
# initial_t_range = [100,460] #K
# steps           = 20

# if we know the eutectic temperature is bounded to within a range of <= estimate_precision, we simply use a linear interpolation between two x,T pairs
estimate_precision = 1.0 #K

class Eutectic:

    def __init__(self,files,tm,hfus,steps=10):
        self.files = files
        self.tm    = tm
        self.hfus  = hfus
        self.steps = steps

        if not (len(self.files)==len(self.tm)==len(self.hfus)):
            print("Error.  Inputs must be the same length.")


    def calc_xt_curves(self,t_range):

        # initialize settings object
        settings = Settings()
        settings.input.property._h = 'SOLUBILITY'
        # optionally, change to the COSMOSAC2013 method
        settings.input.method = 'COSMOSAC2013'

        # make compounds
        compounds = [Settings() for i in range(len(files))]
        for i,file in enumerate(self.files):
            compounds[i]._h = os.path.join( database_path, file )
            compounds[i].meltingpoint = self.tm[i]
            compounds[i].hfusion      = self.hfus[i]

        comp1_fracs = []
        for i,frac1 in enumerate([0.0,1.0]):
            compounds[0].frac1 = frac1
            compounds[1].frac1 = 1.0-frac1

            settings.input.temperature = " ".join([str(t) for t in t_range] + [str(self.steps)])

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

            if i==0:
                comp1_fracs.append(res["molar fraction"][0])
            else:
                comp1_fracs.append(1.0-res["molar fraction"][1])

        return comp1_fracs


    def calc_eutectic(self,t_range,history=[[],[],[]]):

        comp1_fracs = self.calc_xt_curves(t_range)

        # the temperatures used in the calculation
        temps = [t_range[0]+(t_range[1]-t_range[0])/self.steps * i for i in range(self.steps+1)]

        history[0].extend(comp1_fracs[0])
        history[1].extend(comp1_fracs[1])
        history[2].extend(temps)
        # the difference between compound 1's mole fraction in the two calculations
        # when these mole fractions are the same, we've found the eutectic
        diffs = comp1_fracs[0] - comp1_fracs[1]
        #find where the sign changes (intersection of SLE lines)
        for i in range(self.steps):
            if diffs[i]*diffs[i+1] < 0:

                if temps[i+1]-temps[i] < estimate_precision:
                    # use linear combination of t's
                    tot = abs(diffs[i])+abs(diffs[i+1])
                    w1  = tot-abs(diffs[i]) # same as abs(diffs[i+1])
                    w2  = tot-abs(diffs[i+1])
                    return (((w1*comp1_fracs[0][i]+w2*comp1_fracs[1][i+1])/tot, (w1*temps[i]+w2*temps[i+1])/tot)), history
                else:
                    return self.calc_eutectic( [temps[i],temps[i+1]], history )

        return None

eutectic_calc = Eutectic(files,tm,hfus,steps=steps)
eutectic, history = eutectic_calc.calc_eutectic(initial_t_range)

if not eutectic:
    print("No eutectic point found in the temperature range")
else:
    print("Found eutectic point:")
    print("x_1".rjust(10),"T (K)".rjust(10), "T (C)".rjust(10))
    x,t = eutectic
    print(str(x.round(5)).rjust(10),str(t.round(5)).rjust(10), str((-273.15+t).round(5)).rjust(10))

    # plot the solubility curves and eutectic point
    h_s1 = sorted( list(zip(history[0],history[2])), key=lambda x:x[0])
    h_s2 = sorted( list(zip(history[1],history[2])), key=lambda x:x[0])

    h_s1 = [x for x in h_s1 if x[0]>=eutectic[0]]
    h_s2 = [x for x in h_s2 if x[0]<=eutectic[0]]

    # adjust the melting point back to the correct value for high or low solubility
    for i in range(len(h_s1)):
        if h_s1[i][0] > 0.9999:
            h_s1[i] = (h_s1[i][0],tm[0])
    for i in range(len(h_s2)):
        if h_s2[i][0] < 0.0001:
            h_s2[i] = (h_s2[i][0],tm[1])


    plt.plot([x[0] for x in h_s1],[x[1] for x in h_s1],label="x_1 (solvent compound 1)")
    plt.plot([x[0] for x in h_s2],[x[1] for x in h_s2],label="x_1 (solvent compound 2)")
    plt.plot(eutectic[0],eutectic[1],'o',label="Eutectic point")
    if eutectic[0]<0.5:
        plt.annotate( "  "+str(tuple([xt.round(3) for xt in eutectic])),eutectic,va='center',ha='left')
    else:
        plt.annotate( str(tuple([xt.round(3) for xt in eutectic]))+"  ",eutectic,va='center',ha='right')

    plt.xlabel("Mole fraction compound 1")
    plt.ylabel("Melting point of mixture (K)")
    plt.legend(loc='upper right')
    plt.grid()
    plt.show()

finish()