Python scripting for COSMO-RS with PLAMS¶

General Information¶

Python and the PLAMS library can be used for scripting with COSMO-RS. Due to the speed of COSMO-RS calculations, these jobs can be run interactively from the python interpreter. Larger numbers of jobs or high-throughput calculations can also easily be automated with python scripting. All results are returned as a python object, meaning the properties calculated with COSMO-RS can immediately be post-processed or used directly in other python functions.

Executing the code from the command line¶

It is recommended to use the version of python that is shipped with AMS. This version ensures that all the necessary libraries (e.g., PLAMS) are properly imported and are mutually compatible. The best way to do this is to run the amspython program. That can be executed from the command line as follows:

$AMSBIN/amspython <your_program.py>

where <your_program.py> should of course be replaced by the name of your program.

Specifying a problem type¶

To run COSMO-RS, the user must first provide a problem type for the calculation. This can be done by first creating a Settings object and then specifying the .input.property._h attribute. For example, to set up an activity coefficient calculation, we do the following:

from scm.plams import *

settings = Settings()
settings.input.property._h = 'ACTIVITYCOEF'

For other problem types, the .input.property._h attribute must be set to other values. The other options for this value are summarized below:

Inputting Compounds¶

In PLAMS, each compound is also input as a Settings object. Additional information about the compounds required for the calculation (e.g., mole fraction) can be specified as an attribute of the compound’s Settings object. An example for a calculation with two compounds is given below.

# set the number of compounds
num_compounds = 2

compounds = [Settings() for i in range(num_compounds)]
compounds[0]._h = "Water.coskf"
compounds[1]._h = "1-Hexanol.coskf"

Specifying mole fractions, temperatures, and pressures¶

Mole fractions are attributes of the compound Settings object. There are two types of mole fractions used in COSMO-RS. frac1 is for standard specification of mole fractions in most problem types. frac2 is used when the problem type requires two distinct liquid phases (COMPOSITIONLINE or LOGP). Additionally, the temperature can be specified using the input.temperature attribute of the Settings object. An example of this is shown below:

#set compound mole fractions
compounds[0].frac1 = 0.3
compounds[1].frac1 = 0.7

#set temperature (range)
#to specify a range, use 3 numbers: (1) the lowest temperature,
#(2) the highest temperature, and (3) the steps taken between these temperatures
settings.input.temperature = "298.15"

To specify a temperature range, set the input.temperature object equal to a python str which contains the lower temperature, upper temperature, and number of steps taken between the temperatures. These values should simply be separated by spaces. For example, to specify that a calculation should go over the temperature range 298.15K to 398.15K with 10 temperature steps, do the following:

settings.input.temperature = "298.15 398.15 10"

Pressure works in much the same way. To input the system pressure (in bar), do the following:

settings.input.pressure = "1.5"

Running jobs¶

To run a job with COSMO-RS, first assign the input.compound attribute to the list of compound Settings objects used previously. Then, simply create the job using CRSJob(settings=<your previously defined Settings object>). Once a job is created, you can run it with the .run() function. An example of this is given below:

# specify the compounds as the compounds to be used in the calculation
settings.input.compound = compounds
# create a job that can be run by COSMO-RS
my_job = CRSJob(settings=settings)
# run the job
init()
out = my_job.run()
finish()

Reading the results of a job¶

Once a job has finished running, we can access the results directly in python. First, we can check to see which properties are available. We can do this using the get_prop_names() function on the output. For example, adding the line:

# check for the available properties
print( "Available properties:", out.get_prop_names() )

gives us the available properties as a python set for our calculation type (“ACTIVITYCOEF” in this case). The result of the print statement is the following:

Available properties: {'henrycnodim', 'property', 'deltag', 'henryc', 'nitems', 'gamma',
'ncomp', 'filename', 'temperature', 'frac1', 'G solute', 'mu gas', 'molmass', 'E gas',
'mu', 'usepolyunits', 'mu pure', 'method'}

We can also convert all of the calculation results to a python dict using the get_results() function. For example, to collect all of the results and then print the activity coefficient values (“gamma”), we write the following code:

# convert all the results into a python dict
res = out.get_results()
print( "Activity coef values:\n", res["gamma"] )

This results in the following program output:

Activity coef values:
 [[ 3.71486   ]
 [ 1.04484607]]

Here the two activity coefficient values are returned as elements in a numpy.ndarray. Properties with multiple values are always stored as a numpy array.

Putting all the previous code together, we have the following working example for calculating activity coefficients for 2 components:

from scm.plams import *
from os import path

##################  Be sure to add the path to your own ADFCRS directory here  ##################
#################################################################################################
database_path = "<Path to ADFCRS directory containing .coskf files>"
#################################################################################################

if not path.exists(database_path):
    raise OSError(f'The provided path does not exist. Exiting.')

# initialize settings object
settings = Settings()
settings.input.property._h = 'ACTIVITYCOEF'

# set the number of compounds
num_compounds = 2

compounds = [Settings() for i in range(num_compounds)]
compounds[0]._h = path.join( database_path, "Water.coskf" )
compounds[1]._h = path.join( database_path, "1-Hexanol.coskf" )

#set compound mole fractions
compounds[0].frac1 = 0.3
compounds[1].frac1 = 0.7

#set temperature (range)
#to specify a range, use 3 numbers: (1) the lowest temperature,
#(2) the highest temperature, and (3) the steps taken between these temperatures
settings.input.temperature = "298.15"

# specify the compounds as the compounds to be used in the calculation
settings.input.compound = compounds
# create a job that can be run by COSMO-RS
my_job = CRSJob(settings=settings)
# run the job
init()
out = my_job.run()
finish()

# check for the available properties
print( "Available properties:", out.get_prop_names() )

# convert all the results into a python dict
res = out.get_results()
print( "Activity coef values:\n", res["gamma"] )

Plotting results¶

2D graphs can also be generated to visualize the results with the plot function. The plot function takes as a first argument any (or multiple) of the following:

• a numpy.ndarray object. This can be passed to the function as a dictionary value after calling the get_results() function.
• the name of a property. This property is read from the results and plotted. For a list of available properties, use the get_prop_names() function.

Additionally, the plot function takes the following keyword arguments:

• x_axis. This can be the name of a property or a numpy.ndarray object. This represents the independent variable in the plot. This value must be one dimensional, meaning it cannot be indexed over both compounds and temperatures.
• x_label. This can be used to label the x axis in the plot.
• y_label. This can be used to label the y axis in the plot.
• plot_fig. This is set to True/False to indicate whether a plotted figure should be displayed. The default is True.

The results of plot are returned as a matplotlib.pyplot object and can be further modified.

To demonstrate the use of plot, we do an example in which we calculate the solubility of methane gas in 1-Octanol and Ethanol across the temperature range from 298.15K to 398.15K. We also include the vapor pressure of methane using the VPM1 model. The code is shown below:

from scm.plams import *
from os import path

##################  Be sure to add the path to your own ADFCRS directory here  ##################
#################################################################################################
database_path = "<Path to ADFCRS directory containing .coskf files>"
#################################################################################################

if not path.exists(database_path):
    raise OSError(f'The provided path does not exist. Exiting.')

# initialize settings object
settings = Settings()
settings.input.property._h = 'PURESOLUBILITY'

# this indicates we're calculating gas solubility
settings.input.property.isobar = ''

# set the number of compounds
num_compounds = 3

compounds = [Settings() for i in range(num_compounds)]
compounds[0]._h = path.join( database_path, "Methane.coskf" )
compounds[1]._h = path.join( database_path, "1-Octanol.coskf" )
compounds[2]._h = path.join( database_path, "Ethanol.coskf" )

#set compound mole fractions
#for pure solubility the solvent gets a mole fraction of 1
#and the solute does not have the frac1 attribute
compounds[1].frac1 = 1
compounds[2].frac1 = 1

# specify the vapor pressure equation for methane
compounds[0].vp_equation = "VPM1"
compounds[0].vp_params = "-1039.67755001 -0.183945615995 0.00061368649128 10.1113503603315 0.0"

#set temperature (range)
#to specify a range, use 3 numbers: (1) the lowest temperature,
#(2) the highest temperature, and (3) the steps taken between these temperatures
settings.input.temperature = "298.15 398.15 10"

#1 atm = 1.01325 bar
settings.input.pressure = "1.01325"

# specify the compounds as the compounds to be used in the calculation
settings.input.compound = compounds
# create a job that can be run by COSMO-RS
my_job = CRSJob(settings=settings)
# run the job
init()
out = my_job.run()
finish()

# convert all the results into a python dict
res = out.get_results()

#plot the solubilities in g/L solution
#the [1:] indicates that we're not plotting the values for methane (these are automatically set to 0)
out.plot( res["solubility g_per_L_solution"][1:], x_axis = "temperature", x_label="Temperature", y_label = "solubility g/L solution")

This code generates the following plot:

Fig. 3 The output of the plot function for a gas solubility calculation.

Examples¶

from scm.plams import *
from os import path

##################  Be sure to add the path to your own ADFCRS directory here  ##################
#################################################################################################
database_path = "<Path to ADFCRS directory containing .coskf files>"
#################################################################################################

if not path.exists(database_path):
    raise OSError(f'The provided path does not exist. Exiting.')

# initialize settings object
settings = Settings()
settings.input.property._h = 'LOGP'

# set the number of compounds
num_compounds = 3

compounds = [Settings() for i in range(num_compounds)]
compounds[0]._h = path.join( database_path, "1-Octanol.coskf" )
compounds[1]._h = path.join( database_path, "Water.coskf" )
compounds[2]._h = path.join( database_path, "Ibuprofen.coskf" )

#phase1 (octanol phase)
compounds[0].frac1 = 0.725
compounds[1].frac1 = 0.275

#phase2 (water phase)
compounds[0].frac2 = 0
compounds[1].frac2 = 1

#set temperature (range)
#to specify a range, use 3 numbers: (1) the lowest temperature,
#(2) the highest temperature, and (3) the steps taken between these temperatures
settings.input.temperature = "298.15"

# specify the compounds as the compounds to be used in the calculation
settings.input.compound = compounds
# create a job that can be run by COSMO-RS
my_job = CRSJob(settings=settings)
# run the job
init()
out = my_job.run()
finish()

# convert all the results into a python dict
res = out.get_results()

# print the logP of Ibuprofen
print ("logP of Ibuprofen:", res["logp"][2])

logP of Ibuprofen: [ 4.67381309]

Binary mixture¶

In this example, we calculate a binary mixture of water and 2-Hexanone and plot the vapor pressures as a function of composition. We also show how to change the method and calculate the binary mixture with the COSMO-SAC2013-Xiong model.

from scm.plams import *
from os import path

##################  Be sure to add the path to your own ADFCRS directory here  ##################
#################################################################################################
database_path = "<Path to ADFCRS directory containing .coskf files>"
#################################################################################################

if not path.exists(database_path):
    raise OSError(f'The provided path does not exist. Exiting.')

# initialize settings object
settings = Settings()
settings.input.property._h = 'BINMIXCOEF'

# let's also change to the COSMOSAC2013 method
settings.input.method = 'COSMOSAC2013'

# set the number of compounds
num_compounds = 2

compounds = [Settings() for i in range(num_compounds)]
compounds[0]._h = path.join( database_path, "Water.coskf" )
compounds[1]._h = path.join( database_path, "2-Hexanone.coskf" )

# use the vapor pressures from the VPM1 model
compounds[0].vp_equation = "VPM1"
compounds[0].vp_params = "-6093.40215895 -3.09584608667 0.000498622924643 34.47450247140318 0.0"
compounds[1].vp_equation = "VPM1"
compounds[1].vp_params = "-6474.348470271438 -6.057589837807771 0.003390587477679571 51.07134238467479 0.0"

#set temperature (range)
#to specify a range, use 3 numbers: (1) the lowest temperature,
#(2) the highest temperature, and (3) the steps taken between these temperatures
settings.input.temperature = "298.15"

# specify the compounds as the compounds to be used in the calculation
settings.input.compound = compounds
# create a job that can be run by COSMO-RS
my_job = CRSJob(settings=settings)
# run the job
init()
out = my_job.run()
finish()

# convert all the results into a python dict
res = out.get_results()

#plot all the pressures as a function of mole fraction of water
out.plot( "vapor pressure", "pressure", x_axis = res["molar fraction"][0], x_label="mole fraction water", y_label = "Pressure (bar)")

The code generates the following plot:

Fig. 4 A plot showing the total and partial vapor pressures for the water/2-Hexanone system.