pyCRS : Basic usage for PropPred and FastSigma

Basic usage

This basic example provides a minimal script input a molecule (paracetamol) as a smiles string, calculate a single property, and then estimate the default (COSMO-RS) sigma profile.

[show/hide code]
import pyCRS

mol = pyCRS.Input.read_smiles("CC(=O)Nc1ccc(O)cc1")  # paracetamol

print("available properties:", pyCRS.PropPred.available_properties)

pyCRS.PropPred.estimate(mol, "hfusion")
print("hfusion value:", mol.properties["hfusion"], pyCRS.PropPred.units["hfusion"])

pyCRS.FastSigma.estimate(mol, method="COSMO-RS", display=False)

sigma_profiles = mol.get_sigma_profile()
print("Total sigma profile:")
print(sigma_profiles["Total Profile"])
print("H-Bonding:")
print(sigma_profiles["H-Bonding Profile"])

The output produced is the following:

available properties: ['acentricfactor', 'autoignitiontemp', 'boilingpoint', 'critcompress', 'criticalpressure', 'criticaltemp', 'criticalvol', 'density', 'dielectricconstant', 'dipolemoment', 'entropygas', 'entropystd', 'flashpoint', 'gformstd', 'gidealgas', 'hcombust', 'hformstd', 'hfusion', 'hidealgas', 'hsublimation', 'liquidviscosity', 'lowflamlimper', 'meltingpoint', 'molarvol', 'parachor', 'radgyration', 'refractiveindex', 'solubilityparam', 'synacc', 'tpp', 'tpt', 'upflamlimper', 'vaporpressure', 'vdwarea', 'vdwvol']
hfusion value: 28.084352493286133 kJ/mol
Total sigma profile:
[0.0, 0.0, 0.0, 0.002353191375733, 0.050697326660157, 0.24713134765625, 0.5523872375488279, 0.840805053710938, 1.301651000976562, 1.316818237304688, 1.408416748046875, 1.173583984375, 1.027912139892578, 1.160888671875, 0.979301452636719, 0.8482666015625, 0.5888519287109371, 5.276319718325397, 11.539728505193898, 13.300330108724362, 12.906459205297642, 10.834775504652777, 8.539289410607651, 7.653812772468701, 7.964459592741903, 6.8641550757498555, 9.304299127480718, 7.66864582217209, 11.452185796196918, 12.639082653211547, 13.748599090727417, 11.745518829550344, 3.4796992518586065, 3.1053283341103315, 2.864667892456055, 3.6309814453125, 4.3218994140625, 3.9415283203125, 2.786376953125, 0.968017578125, 0.48760986328125, 0.319122314453125, 0.047515869140625, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
H-Bonding:
[0.0, 0.0, 0.0, 0.002353191375733, 0.050697326660157, 0.246719360351562, 0.552310943603516, 0.840805053710938, 1.301651000976562, 1.313278198242188, 1.3948974609375, 1.173583984375, 1.027912139892578, 1.15216064453125, 0.979301452636719, 0.8482666015625, 0.5888519287109371, 0.065513610839843, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.277877807617188, 2.451416015625, 3.6309814453125, 4.3218994140625, 3.9354248046875, 2.78125, 0.968017578125, 0.48760986328125, 0.319122314453125, 0.047515869140625, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Temperature-dependent properties

This example calculates the vapor pressure and produces a plot of vapor pressure against temperature.

[show/hide code]
import pyCRS
import matplotlib.pyplot as plt

mol = pyCRS.Input.read_smiles("CCCCCCO")

prop_name = "vaporpressure"
pyCRS.PropPred.estimate(mol, temperatures=[290, 295, 300, 305, 310, 315])
print("Results:", mol.properties_tdep[prop_name])

x, y = mol.get_tdep_values(prop_name)
unit = pyCRS.PropPred.units[prop_name]
plt.plot(x, y, "-o")
plt.ylabel(f"vapor pressure ({unit})")
plt.xlabel("Temperature (K)")
# plt.savefig('./pyCRS_PropPred_Tdep.png')
plt.show()

The output shows the format of the results: (temperature, vapor pressure) pairs

Results: [(290.0, 0.0006316340979498092), (295.0, 0.0009549864170162676), (300.0, 0.0014201952225525484), (305.0, 0.002079200184963613), (310.0, 0.0029991100667307634), (315.0, 0.004265490872465904)]

Finally, the plot produced is the following:

/scm-uploads/doc.trunk/COSMO-RS/_images/pyCRS_PropPred_Tdep.png

Fig. 4 The estimated vapor pressure versus temperature for 1-Hexanol

Estimating multiple properties

Estimating multiple properties is as simple as supplying a list of property names to the PropPred interface. All properties are estimated by default if no property argument is supplied. In this example, we first estimate a few properties, and then estimate all properties.

[show/hide code]
import pyCRS


def print_props(mol):
    for prop, value in mol.properties.items():
        unit = pyCRS.PropPred.units[prop]
        print(f"{prop:<20s}: {value:.3f} {unit}")

    for prop, value in mol.properties_tdep.items():
        print(f"{prop:<20s}:")
        unit = pyCRS.PropPred.units[prop]
        propunit = f"{prop} ({unit})"
        print("T (K)".rjust(30) + f"{propunit:>30s}")

        for t, v in value:
            print(f"{t:>30.3f}{v:>30.8g}")


mol = pyCRS.Input.read_smiles("CCCCCCO")

props = ["meltingpoint", "boilingpoint", "density", "flashpoint", "vaporpressure"]
pyCRS.PropPred.estimate(mol, props, temperatures=[298.15, 308.15, 318.15, 328.15])
print("Results (temp-independent) :", mol.properties)
print("Results (temp-dependent)   :", mol.properties_tdep)

# we can also estimate all properties by supplying the property name 'all' or simply omitting this argument
pyCRS.PropPred.estimate(mol, temperatures=[298.15, 308.15, 318.15, 328.15])
print_props(mol)

The output produced is the following:

Results (temp-independent) : {'boilingpoint': 435.7771752780941, 'density': 0.7918196941677842, 'flashpoint': 342.2705857793571, 'meltingpoint': 231.1412353515625, 'molarvol': 0.1289491355419159}
Results (temp-dependent)   : {'vaporpressure': [(298.1499938964844, 0.00122854727137622), (308.1499938964844, 0.0026233569814824815), (318.1499938964844, 0.005288582928457778), (328.1499938964844, 0.010122673317257832)]}
boilingpoint        : 435.777 K
criticalpressure    : 34.349 bar
criticaltemp        : 878.101 K
criticalvol         : 0.404 L/mol
density             : 0.792 kg/L (298.15 K)
dielectricconstant  : 10.951
entropygas          : 439.885 J/(mol K)
flashpoint          : 342.271 K
gidealgas           : -131.869 kJ/mol
hcombust            : -3678.121 kJ/mol
hformstd            : -384.388 kJ/mol
hfusion             : 18.505 kJ/mol
hidealgas           : -316.821 kJ/mol
hsublimation        : 80.980 kJ/mol
meltingpoint        : 231.141 K
molarvol            : 0.129 L/mol
parachor            : 289.059
solubilityparam     : 10.129 √(cal/cm^3)
synacc              : 6.747
tpt                 : 230.404 K
vdwarea             : 171.059 Ų
vdwvol              : 120.519 ų
liquidviscosity     :
                         T (K)        liquidviscosity (Pa-s)
                       298.150                  0.0044653385
                       308.150                   0.003363708
                       318.150                  0.0025843814
                       328.150                  0.0020210327
vaporpressure       :
                         T (K)           vaporpressure (bar)
                       298.150                  0.0012285473
                       308.150                   0.002623357
                       318.150                  0.0052885829
                       328.150                   0.010122673

Calculating sigma profiles with all models

This example demonstrates how to calculate COSMO-RS sigma profiles with both available models (SG1 and FS1). We’ll use 4-Methylphenol for this example. The variable ref_sp in the script is the \(\sigma\)-profile calculated in the AMS COSMO-RS.

[show/hide code]
import pyCRS
import matplotlib.pyplot as plt

chdens = [-0.025 + 0.001 * x for x in range(51)]
ref_sp = [
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.17050209,
    0.73517288,
    0.98383886,
    0.86397332,
    1.15571158,
    0.67051286,
    0.79380670,
    0.71422825,
    0.68854784,
    0.72126747,
    0.72414167,
    0.83293801,
    2.33617486,
    5.42197092,
    8.58230745,
    9.84294559,
    8.81993395,
    8.80932110,
    6.79578176,
    6.57764296,
    6.07802490,
    9.19651976,
    11.04721719,
    12.80255519,
    11.17672544,
    12.28788009,
    10.72121222,
    2.77829981,
    1.11136819,
    1.58235813,
    1.13444125,
    1.81316980,
    2.44468560,
    2.02558727,
    0.01243933,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
    0.00000000,
]

mol = pyCRS.Input.read_smiles("Cc1ccc(O)cc1")

pyCRS.FastSigma.estimate(mol, method="COSMO-RS", model="FS1")
sp_fs1 = mol.get_sigma_profile()["Total Profile"]

pyCRS.FastSigma.estimate(mol, method="COSMO-RS", model="SG1")
sp_sg1 = mol.get_sigma_profile()["Total Profile"]

plt.plot(chdens, ref_sp, "--", label="Reference $\sigma$-profile")
plt.plot(chdens, sp_fs1, label="FS1 $\sigma$-profile")
plt.plot(chdens, sp_sg1, label="SG1 $\sigma$-profile")
plt.ylabel("Area ($\AA^2$)")
plt.xlabel("$\sigma$")
plt.grid()
plt.legend()
# plt.savefig('./pyCRS_PropPred_SigmaProfile.png')
plt.show()

Finally, the plot produced shows the various \(\sigma\)-profiles produced.

/scm-uploads/doc.trunk/COSMO-RS/_images/pyCRS_PropPred_SigmaProfile.png

Fig. 5 The sigma profile of 4-Methylphenol