Sigma Moments 的计算/土壤吸附

Sigma Moments 是从 sigma profile 中导出的非常有用的化学描述符，类似于统计分布的“矩”（Moment），并可以认为是将 sigma profile 中存在的高维信息，减少为表征该 sigma profile 的少数描述符。Sigma Moments 已知是 QSPR 中有价值的描述符，并被认为很好地代表了溶剂空间。例如，在以下文献中，发现 Sigma Moments 与土壤吸附系数具有良好的相关性：

• Prediction of soil sorption coefficients with a conductor-like screening model for real solvents, Environmental Toxicology and Chemistry, Volume21, Issue12, 2002, 2562-2566, DOI: 10.1002/etc.5620211206

以下脚本将计算一些常见分子的前几个 Sigma Moments 以及氢键受体和氢键供体矩。

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()

1. 该脚本目前计算的是Water、Hexane、Ethanol、Acetone四个分子，用户可以在“filenames = ”这一行自行修改为自己需要计算的分子，分子的.coskf文件与该脚本放置在同一个文件夹内（注意路径中不要包含中文、希腊字母、空格、括号等）
2. 将上述内容，创建为Python脚本（参考：Python脚本的专业编辑器）
3. 在如上命令行环境中，进入脚本所在文件夹（如果 AMSJobs 进入该文件夹，则 Help → Command-line输入 sh 回车即处于该文件夹下），输入命令并回车即可：

amspython main.py

其中包括 Sigma Moments 的四个分量，以及氢键受体和氢键供体矩。

物质 X 的土壤吸附采用 K_OC 定义：

那么物质 X 的 K_OC 计算，可以通过如下公式：

log K_OC = 0.0168(±8) M₀ -0.017(±2) M₂ - 0.040(±4) M₃ + 0.19(±5) M_acc - 0.27(±5) M_don + 0.37(±14)

将上面计算得到的，物质 X 的 Sigma Moments 数据带入上式，即得到其土壤吸附 log K_OC，求 10 的 log K_OC次方，即得到 K_OC。

其他固-液分配系数的计算，原则上可以沿用上述公式，但需要对上面公式中的常数，根据实验重新拟合。