ADF and COSMO-RS workflow¶
Note: This example requires AMS2023 or later.
This example uses ADF to generate the .coskf file for benzene. You can also modify it to instead use the Benzene.coskf from the ADFCRS-2018 database. Note that you first need to install the ADFCRS-2018 database.
The script then plots the sigma profile. This is not a necessary step but is done for illustration purposes only.
Then a solubility calculation is performed for benzene in water between 0 and 10 degrees C. The melting point and enthalpy of fusion can either be estimated using the property prediction tool, or the experimental numbers can be given (recommended).
Example usage: AMS2023+ required (ams_crs.py
)
This script is built on the PLAMS. For more detailed useage of COSMO-RS with PLAMS, refer to the documentation.
[show/hide code]
#!/usr/bin/env amspython
import os
import matplotlib.pyplot as plt
from scm.plams import Settings, Units, from_smiles, CRSJob, init
from scm.plams.recipes.adfcosmorscompound import ADFCOSMORSCompoundJob
def solubility():
# database can also be replaced with the output of "$AMSBIN/amspackages loc adfcrs" /ADFCRS-2018
database = CRSJob.database()
if not os.path.exists(database):
raise OSError(f"The provided path does not exist. Exiting.")
solute_smiles = "c1ccccc1"
solute_coskf = generate_coskf(solute_smiles, "adf_benzene") # generate files with ADF
# solute_coskf = os.path.abspath('plams_workdir/adf_benzene/adf_benzene.coskf') # to not rerun the ADF calculation
# solute_coskf = os.path.join(database, 'Benzene.coskf') # to load from database
# You can also estimate the solute properties with the Property Prediction tool. See the Property Prediction example
solute_properties = {"meltingpoint": 278.7, "hfusion": 9.91} # experimental values for benzene, hfusion in kJ/mol
solvent_coskf = os.path.join(database, "Water.coskf")
solvent_density = 1.0
s = Settings()
s.input.property._h = "solubility"
s.input.property.DensitySolvent = solvent_density
s.input.temperature = "273.15 283.15 10"
s.input.pressure = "1.01325 1.01325 10"
s.input.compound = [Settings(), Settings()]
s.input.compound[0]._h = solvent_coskf
s.input.compound[0].frac1 = 1.0
s.input.compound[1]._h = solute_coskf
s.input.compound[1].nring = 6 # number of ring atoms benzene
s.input.compound[1].meltingpoint = solute_properties["meltingpoint"]
s.input.compound[1].hfusion = solute_properties["hfusion"] * Units.convert(
1.0, "kJ/mol", "kcal/mol"
) # convert from kJ/mol to kcal/mol
job = CRSJob(name="benzene_in_water", settings=s)
job.run()
plot_results(job)
def generate_coskf(smiles, jobname=None):
molecule = from_smiles(smiles, nconfs=100, forcefield="uff")[0]
job = ADFCOSMORSCompoundJob(name=jobname, molecule=molecule)
job.run()
plot_sigma_profile(job.results)
return job.results.coskfpath()
def plot_results(job):
res = job.results.get_results("SOLUBILITY")
solubility_g_L = res["solubility g_per_L_solvent"][1]
temperatures = res["temperature"]
for temperature, sol_g_l in zip(temperatures, solubility_g_L):
print(f"{temperature:.2f} {sol_g_l:.4f}")
plt.plot(temperatures, solubility_g_L)
plt.xlabel("Temperature (K)")
plt.ylabel("Solubility (g/L solvent)")
plt.show()
def get_sigma_profile(coskf_file):
s = Settings()
s.input.property._h = "PURESIGMAPROFILE"
s.input.compound._h = coskf_file
job = CRSJob(name="sigma_profile", settings=s)
res = job.run()
return res.get_sigma_profile()
def plot_sigma_profile(results):
coskf_path = results.coskfpath()
sigma = get_sigma_profile(coskf_path)
xlabel = "σ (e/A**2)"
for profile in sigma:
if profile == xlabel:
continue
plt.plot(sigma[xlabel], sigma[profile], label=profile.split(".")[0])
plt.xlabel("σ (e/Å**2)")
plt.ylabel("p(σ)")
plt.legend()
plt.show()
def main():
# this line is not required in AMS2025+
init()
solubility()
if __name__ == "__main__":
main()
pyCRS example usage: AMS2024+ required (ams_pyCRS.py
)
This script utilizes pyCRS, a convenient python wrapper built on PLAMS and SQL database for various thermodynamic calculations. For more detatils, refer to the pyCRS Overview.
[show/hide code]
#!/usr/bin/env amspython
"""
- Solubility Calculation using pyCRS workflow with AMS2024+ and AMS2025+ Features
- Utilizes COSKFDatabase and CRSSystem from pyCRS (AMS2024+ Required)
- Use of densf2hbc in ADFCOSMORSCompoundJob for COSMO-SAC DHB-MESP method (AMS2025+ Required)
"""
import os
import matplotlib.pyplot as plt
from pyCRS.Database import COSKFDatabase
from pyCRS.CRSManager import CRSSystem
from scm.plams import Settings, from_smiles, CRSJob, CRSResults
from scm.plams.recipes.adfcosmorscompound import ADFCOSMORSCompoundJob
def solubility_pyCRS():
db = COSKFDatabase("my_coskf_db.db")
if db.get_compounds_id("water")[0] is None:
db.add_compound("Water.coskf")
if db.get_compounds_id("benzene")[0] is None:
solute_smiles = "c1ccccc1"
# Define molecular information for Benzene to be stored in the "Compound Data" section of the *.coskf* file
mol_info = {}
mol_info["IUPAC"] = "Benzene"
mol_info["Other Name"] = None
mol_info["CAS"] = "71-43-2"
mol_info["SMILES"] = solute_smiles
# Generate COSKF file with densf2hbc for DHB-MESP in COSMO-SAC
coskf_file = generate_coskf(solute_smiles, jobname="adf_benzene_densf", mol_info=mol_info, densf2hbc=True)
db.add_compound(coskf_file)
db.add_physical_property("benzene", "meltingpoint", 278.7)
db.add_physical_property("benzene", "hfusion", 9.91, unit="kJ/mol")
plot_sigma_profile_pyCRS(name="benzene", db=db)
crs = CRSSystem()
mixture = {}
mixture["water"] = 1.0
mixture["benzene"] = 0.0
temp = "273.15 283.15 10"
solvent_density = 1.0
additional_sett = Settings()
additional_sett.input.property.DensitySolvent = solvent_density
crs.add_Mixture(
mixture=mixture,
temperature=temp,
database="my_coskf_db.db",
problem_type="solubility",
additional_sett=additional_sett,
)
crs.runCRSJob()
plot_results(crs.outputs[0])
def generate_coskf(smiles, jobname=None, mol_info=None, densf2hbc=False):
molecule = from_smiles(smiles, nconfs=100, forcefield="uff")[0]
job = ADFCOSMORSCompoundJob(name=jobname, molecule=molecule, mol_info=mol_info, densf2hbc=densf2hbc)
job.run()
return job.results.coskfpath()
def plot_results(job):
if isinstance(job, CRSJob):
res = job.results.get_results("SOLUBILITY")
elif isinstance(job, CRSResults):
res = job.get_results("SOLUBILITY")
solubility_g_L = res["solubility g_per_L_solvent"][1]
temperatures = res["temperature"]
for temperature, sol_g_l in zip(temperatures, solubility_g_L):
print(f"{temperature:.2f} {sol_g_l:.4f}")
plt.plot(temperatures, solubility_g_L)
plt.xlabel("Temperature (K)")
plt.ylabel("Solubility (g/L solvent)")
plt.show()
def plot_sigma_profile_pyCRS(name, db):
crs = CRSSystem()
mixture = {name: 1.0}
crs.add_Mixture(mixture=mixture, problem_type="PURESIGMAPROFILE", database=db)
crs.runCRSJob()
sigma = crs.outputs[0].get_sigma_profile()
xlabel = "σ (e/A**2)"
for profile in sigma:
if profile == xlabel:
continue
plt.plot(sigma[xlabel], sigma[profile], label=profile.split(".")[0])
plt.xlabel("σ (e/Å**2)")
plt.ylabel("p(σ)")
plt.legend()
plt.show()
def main():
solubility_pyCRS()
if __name__ == "__main__":
main()
Results

273.15 1.6797
274.15 1.7258
275.15 1.7731
276.15 1.8215
277.15 1.8711
278.15 1.9220
279.15 1.9603
280.15 1.9823
281.15 2.0046
282.15 2.0273
283.15 2.0503
