Sigma Moments¶
Sigma moments are useful chemical descriptors derived from the sigma profile. They are analogous to moments of a statistical distribution and can be thought of as a way to reduce the high-dimensional information present in a sigma profile to a smaller number of descriptors that characterize that sigma profile. Sigma moments are known to be valuable descriptors in QSPR and are thought to represent the solvent space well [1].
The following script will calculate the first several sigma moments as well as a H-bond acceptor and H-bond donor moment for a few common molecules.
Python code¶
import os
import numpy as np
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
class SigmaMoments:
def __init__(self,filenames,hb_cutoff=0.00854):
self.filenames = filenames
self.hb_cutoff = hb_cutoff
def calculate_moments(self) -> dict:
self.moments = {}
self.calc_profiles_and_chdens()
self.calc_standard_moments()
self.calc_hb_moments()
return self.moments
def calc_profiles_and_chdens(self):
# initialize settings object
settings = Settings()
settings.input.property._h = 'PURESIGMAPROFILE'
# set the cutoff value for h-bonding
settings.parameters.sigmahbond = self.hb_cutoff
compounds = [Settings() for i in range(len(self.filenames))]
for i,filename in enumerate(filenames):
compounds[i]._h = os.path.join(database_path, filename)
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()
# retain profiles and charge density values
self.tot_profiles = res["profil"]
self.hb_profiles = res["hbprofil"]
self.chdens = res["chdval"]
def calc_standard_moments(self,max_power=3):
for i in range(max_power+1):
tmp_moms = []
for prof in self.tot_profiles:
tmp_moms.append( np.sum(prof*np.power(self.chdens,i)) )
self.moments["MOM_"+str(i)] = tmp_moms
def calc_hb_moments(self):
self.moments["MOM_hb_acc"] = []
self.moments["MOM_hb_don"] = []
zeros = np.zeros(len(self.chdens))
for prof in self.hb_profiles:
self.moments["MOM_hb_acc"].append(np.sum( prof * np.maximum(zeros,self.chdens-self.hb_cutoff) ))
self.moments["MOM_hb_don"].append(np.sum( prof * np.maximum(zeros,-self.chdens-self.hb_cutoff) ))
# the files we want to use to calculate sigma moments
filenames = ["Water.coskf", "Hexane.coskf","Ethanol.coskf","Acetone.coskf"]
sm = SigmaMoments(filenames)
moms = sm.calculate_moments()
max_mom_len = max([len(m) for m in moms])
print()
print( (" "*5).join(["Moment".ljust(max_mom_len)]+filenames))
lens = [len(fn) for fn in filenames]
for mom_name in moms:
print( (" "*5).join([mom_name.ljust(max_mom_len)]+[('{0:.5g}'.format(m)).rjust(l) for m,l in zip(moms[mom_name],lens)]))
finish()