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

[show/hide code]

import os
import numpy as np
from scm.plams import Settings, init, finish, CRSJob, config

########  Note: Ensure to configure the database path to either the installed ADFCRS-2018 directory or your own specified directory  ########

database_path = os.path.join(os.environ["SCM_PKG_ADFCRSDIR"], "ADFCRS-2018")
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()

The output produced is the following

Moment         Water.coskf     Hexane.coskf     Ethanol.coskf     Acetone.coskf
MOM_0               43.011           160.38            90.019            103.28
MOM_1           0.00026666         0.002145        0.00089555         0.0011418
MOM_2            0.0062556         0.001061         0.0046302          0.004566
MOM_3          -3.8253e-07       1.1557e-07        1.5947e-05         2.883e-05
MOM_hb_acc        0.078179                0           0.06562          0.068401
MOM_hb_don        0.078272                0          0.034516                 0

References

1

A. Klamt, COSMO-RS From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design, Elsevier. Amsterdam (2005), ISBN 0-444-51994-7.