空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

上机作业1: 不对称催化中的非共价作用（学生版）

机器学习及其在化学中的应用

CuiChang2022

更新于 2024-10-14

推荐镜像 :mlgpu:0929

推荐机型 :c2_m4_cpu

242

数据集

nci_birman(v1)

背景概述

对映选择性(enantioselectivity)是有机化学重要的研究主题之一. 在有机反应中, 手性催化剂的对映选择性往往与反应过渡态中一些重要的非共价相互作用(non-covalent interactions, NCIs)有关. 以Birman等报道的酰基转移反应为例:

alt image.png

在Birman等假设的机理中, 手性催化剂的对映选择性由 $π$ - $π$ 堆积或 $π$ -阳离子堆积作用介导: 特定手性的过渡态因位阻作用破坏了芳环的平行取向, 不利于堆积作用的形成.

alt image.png

本次上机作业中, 我们将以线性回归模型定量地探索所用手性催化剂的性质对对映选择性的影响, 试图探索:

哪些性质参数会影响对映选择性?
如何优化性质参数可以提升对映选择性 (例如减少还是增加某参数)?

为此, 我们将考察如下图所示的模型体系:

alt image.png

其中, 底物的芳香环(A, a, b)按照特定的几何约束条件与催化剂(D, c, d)排列成复合物, 而芳环与底物的其余分子骨架(包含醇羟基和烷基)在图中箭头所示位置相连.

参考文献: J. Am. Chem. Soc. 2017, 139, 20, 6803–6806

代码

文本

数据探索与清洗

代码

文本

数据的初步探索

在数据集nci_birman.csv中, 包含如下三组(用作特征的)参数:

$π$ - $π$ 堆积作用参数: d_pi_d, d_pi_D, e_pi_d, e_pi_D. 根据特定的电子结构计算方法(B97-D/def2TZVP水平), 对复合物进行几何优化, 给出几何结构信息与互作能(复合物能量减去各组分能量和). 其中:
- 每组复合物中, 底物芳环由于具体朝向的不同, 存在2种可能构象(比如下图的例子), 分别以后缀_d和_D作区分;
- 在计算中, 对底物与催化剂分子的两个芳环施加约束条件, 使二者完全“正对彼此”(环平面平行, 且两环中心的连线与环平面垂直). 此时, 我们找到使电子能量(前缀e_, 单位kcal/mol)最低的、也就是使复合物最稳定的环间距(前缀d_, 单位Å)
烷基几何参数: L_Alk, B1_Alk, B5_Alk, 分别代表长度、最小宽度、最大宽度.
芳环几何参数: L_Ar, B1_Ar, B5_Ar, 分别代表长度、最小宽度、最大宽度. 后两组参数称为Sterimol参数, 其具体含义可从下图中做简单的理解:

alt image.png

我们的线性回归模型将以上述参数作为特征(在化学信息学中, 这些特征也常称为描述符, descriptor), 试图预测特定催化剂针对特定底物的对映选择性. 在数据集中, 选择性以对映比 $er$ 描述, 这是一个实验可测的量. 通过热力学公式, 我们可以将其与两种对映体的活化自由能之差 $ΔΔ G^{‡}$ 之间建立定量联系:

$ΔΔ G^{‡} = - RT ln (er) .$

代码

文本

[ ]

import pandas as pd

nci_birman = pd.read_csv("/bohr/mlchem-2024fall-qmtv/v1/nci_birman.csv", sep="\t")

nci_birman.iloc[:5, :]

代码

文本

数据的清洗与加工

对于堆积相互作用参数, 我们不妨对其进行降维: 取两个可能构象_d和_D的参数加权和: $D π_{w} E π_{w} := c_{d} D π_{d} + c_{D} D π_{D}, := c_{d} E π_{d} + c_{D} E π_{D},$ 其中, 权重 $c_{X} (X = d, D)$ 取为(按电子能量计算的)Boltzmann因子: $c_{X} = \frac{exp ( - \frac{E π _{X}}{RT} )}{exp ( - \frac{E π _{d}}{RT} ) + exp ( - \frac{E π _{D}}{RT} ) .}$

代码

文本

任务1: 对数据nci_birman进行进一步加工处理.

编写函数prepare_data(), 完成数据的初步清洗.
- 输入:
  - 直接读取自csv的原始数据集in_df;
  - 用于进行加权和与 $ΔΔ G^{‡}$ 计算的温度temperature, 以 $K$ 为单位.
- 返回: 处理好的特征和标签. 其中:
  - 特征为pd.DataFrame格式, 包括:
    - 降维后的(加权)参数 $D π_{w}, E π_{w}$ , 以及二者的交叉项(乘积) $D π_{w} E π_{w}$ , 分别命名为d_pi_w, e_pi_w与de_pi_w;
    - 所有原有的Sterimol几何参数, 名称不变.
  - 标签为pd.Series格式, 为自由能差 $ΔΔ G^{‡} = - RT ln (er)$ , 单位取为kcal/mol, 命名为delta_delta_G. 这里给出单位换算因子: $1 kcal/mol = 4.184 kJ/mol$ .

代码

文本

提示

注意 $er (%)$ 一栏的单位, 例如表格中若该栏的值为 $1$ , 则实际代表 $er = 1%$ , 或说 $ΔΔ G^{‡} = - RT ln 0.01$ .

代码

文本

[ ]

from typing import Tuple

import numpy as np

R = 8.314 # ideal gas constant

KCAL_TO_KJ = 4.184 # unit conversion

def prepare_data(in_df: pd.DataFrame, temperature: float=298.0) -> Tuple[pd.DataFrame, pd.Series]:

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

代码

文本

完成该函数后, 请运行下面两个代码块先做初步检查.

代码

文本

[ ]

### SANITY CHECK CELL 1.1 ###

### EXPECTED OUTPUT ###

'''

d_pi_w e_pi_w de_pi_w L_Alk B1_Alk B5_Alk L_Ar B1_Ar B5_Ar

0 3.54 -10.24 -36.2496 4.36 2.92 3.35 6.38 1.77 3.15

1 3.54 -10.24 -36.2496 4.35 2.09 3.34 6.38 1.77 3.15

2 3.54 -10.24 -36.2496 4.38 1.73 3.33 6.38 1.77 3.15

3 3.54 -10.24 -36.2496 3.08 1.70 2.20 6.38 1.77 3.15

4 3.49 -11.88 -41.4612 3.08 1.70 2.20 6.38 1.88 4.52

'''

X, y_true = prepare_data(nci_birman)

X.iloc[:5, :]

代码

文本

[ ]

### SANITY CHECK CELL 1.2 ###

### EXPECTED OUTPUT ###

'''

0 2.409283

1 2.009314

2 1.939568

3 1.770392

4 1.770392

Name: delta_delta_G, dtype: float64

'''

y_true[:5]

代码

文本

完成初步检查后, 可以运行下面的代码块对该函数进行测试. 测试通过情况将关系到该任务的得分. 请勿修改该代码块中的任何内容.

代码

文本

[ ]

### TEST CELL 1 ###

### CAUTION: DO NOT MODIFY THIS CELL. ###

test_in_df = pd.concat([

pd.DataFrame({

"d_pi_d": [4.89, 4.91, 6.48, 7.00],

"d_pi_D": [7.00, 7.21, 7.22, 7.20],

"e_pi_d": [-3.48, -2.98, -1.80, -1.45],

"e_pi_D": [-1.48, -2.00, -1.96, -2.26]

}),

nci_birman.iloc[:4, 4:]

], axis=1)

test_labels = np.array([2.409283, 2.009314, 1.939568, 1.770392])

test_weighted_features = np.array([

[4.959643, -3.413988, -16.932160],

[5.279007, -2.822771, -14.901427],

[6.899685, -1.890743, -13.045529],

[7.159408, -2.095601, -15.003264]

])

def test_prepare_data_names(in_df):

X_, y_true_ = prepare_data(in_df)

assert set(X.columns) == {"d_pi_w", "e_pi_w", "de_pi_w", "L_Alk", "B1_Alk", "B5_Alk", "L_Ar", "B1_Ar", "B5_Ar"}, "wrong column names."

assert y_true_.name == "delta_delta_G", "wrong label name."

def test_prepare_data_labels(in_df, labels):

_, y_true_ = prepare_data(in_df)

np.testing.assert_allclose(

y_true_, labels, rtol=1.0e-3,

err_msg="wrong label values."

)

def test_prepare_data_weighted_features(in_df, weighted_features):

X_, _ = prepare_data(in_df)

np.testing.assert_allclose(

X_.loc[:, ["d_pi_w", "e_pi_w", "de_pi_w"]], weighted_features, rtol=1.0e-3,

err_msg="wrong weighted feature values."

)

try:

print("names:", end= " ")

test_prepare_data_names(test_in_df)

print("PASSED")

except Exception as err:

print(err)

try:

print("labels:", end=" ")

test_prepare_data_labels(test_in_df, test_labels)

print("PASSED")

except Exception as err:

print(err)

try:

print("weighted features:", end=" ")

test_prepare_data_weighted_features(test_in_df, test_weighted_features)

print("PASSED")

except Exception as err:

print(err)

代码

文本

拆分与归一化

原数据没有直接分出训练集和测试集, 我们可以用sklearn.model_selection.train_test_split对数据集进行手拆分, 训练 : 测试 = 23 : 5.

代码

文本

[ ]

from sklearn.model_selection import train_test_split

test_size = 5

random_state = 0

X_train, X_test, y_true_train, y_true_test = train_test_split(X, y_true, test_size=test_size, random_state=random_state)

代码

文本

随后, 进行数据的归一化(重标度)处理, 以使各个特征维度量纲一致.

任务2: 在各个数据集上分别对特征作均值-方差归一化:

$x_{n}^{'} := \frac{x _{n} - Mean _{n = 1}^{N} ( x )}{Std _{n = 1}^{N} ( x )}, n = 1, \dots, N,$

返回(归一化后的)训练集与测试集. 注意: 测试集上归一化时, 正确的做法是使用训练集上算出的均值和方差, 而非在测试集上重新计算.

任务2.1: 编写函数normalize(), 完成数据特征的归一化.
- 输入:
  - 训练集X_train;
  - 其它数据集other_datasets, 以列表打包起来.
- 返回: 经过归一化处理的X_train_normalized以及其它数据集normalized_datasets(以列表打包起来). 为了和StandardScaler保持一致, 请确保返回值均为np.array格式, 例如对某DataFrame类型的变量df使用格式转换np.array(df).
任务2.2: 在以下代码块的注释区回答问题: 为什么均值-方差归一化需要在训练集和测试集上分开进行? 如果先归一化再拆分, 会导致什么后果?

代码

文本

提示

归一化(在sklearn中叫做标准化, standardization)既可以像上机实习1演示的那样手动计算, 也可以通过函数sklearn.preprocessing.StandardScaler实现.
- 注意: 若选择pandas库手动计算...
  - 请在计算std时设置参数ddof=0, 以和StandardScaler方法保持一致.
  - 与上机实习稍有不同的是, 参与运算的是每一列的均值或标准差, 因此需要对mean()或std()等函数添加参数axis.
警惕机器学习模型在数据集上“作弊”.

代码

文本

[ ]

from sklearn.preprocessing import StandardScaler

from typing import List, Union

def normalize(

X_train: Union[pd.DataFrame, np.array],

other_datasets: Union[List[pd.DataFrame], List[np.array]]

) -> Tuple[np.array, List[np.array]]:

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

### 任务2.2答题区(另起一行时请记得加注释符号#) ###

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

代码

文本

完成该函数后, 请运行下面2个代码块先做初步检查.

代码

文本

[ ]

### SANITY CHECK CELL 2.1 ###

### EXPECTED OUTPUT ###

'''

array([ 2.41352831e-15, -2.02736378e-16, -5.79246795e-16, 9.26794873e-16,

-1.13918536e-15, -9.36448986e-16, 8.83351363e-16, -3.53340545e-15,

-1.68464276e-15]

'''

X_train, other = normalize(X_train, [X_test])

X_test = other[0]

X_train.mean(axis=0)

代码

文本

[ ]

### SANITY CHECK CELL 2.2 ###

### EXPECTED OUTPUT ###

'''

array([1., 1., 1., 1., 1., 1., 1., 1., 1.])

'''

X_train.std(axis=0)

代码

文本

完成初步检查后, 可以运行下面的代码块对该函数进行测试. 测试通过情况将关系到该任务的得分. 请勿修改该代码块中的任何内容.

代码

文本

[ ]

### TEST CELL 2 ###

### CAUTION: DO NOT MODIFY THIS CELL. ###

test_X_train = np.random.randn(100, 5)

test_other_datasets = np.random.randn(10, 50, 5)

def test_normalize_train(X_train):

X_train, _ = normalize(X_train, [])

np.testing.assert_allclose(

X_train.mean(axis=0), np.zeros(5), atol=1.0e-5,

err_msg="invalid X_train normalization: non-zero mean."

)

np.testing.assert_allclose(

X_train.std(axis=0), np.ones(5), rtol=1.0e-3,

err_msg="invalid X_train normalization: non-unit std."

)

def test_normalize_consistency(X_train, other_datasets):

X_train_mean, X_train_std = X_train.mean(axis=0), X_train.std(axis=0)

X_train, normalized_datasets = normalize(X_train, [*other_datasets])

for dataset, normalized_dataset in zip(other_datasets, normalized_datasets):

normalized_mean = (dataset.mean(axis=0) - X_train_mean) / X_train_std

normalized_std = dataset.std(axis=0) / X_train_std

np.testing.assert_allclose(

normalized_dataset.mean(axis=0), normalized_mean, atol=1.0e-5,

err_msg="inconsistent normalization: didn't use the training mean correctly."

)

np.testing.assert_allclose(

normalized_dataset.std(axis=0), normalized_std, rtol=1.0e-3,

err_msg="inconsistent normalization: didn't use the training std correctly."

)

try:

print("training set normalized:", end=" ")

test_normalize_train(test_X_train)

print("PASSED")

except Exception as err:

print(err)

try:

print("normalization consistency:", end=" ")

test_normalize_consistency(test_X_train, test_other_datasets)

print("PASSED")

except Exception as err:

print(err)

代码

文本

模型训练与评估

我们在9个特征与23个训练样本上训练一个线性回归模型, 并评估其在5个测试样本上的预测表现.

任务3: 搭建模型, 完成训练与评估.

任务3.1: 编写函数train_model(), 实现模型训练.
- 输入: 训练集X_train, y_true_train.
- 返回: 一个训练好的线性回归模型sklearn.linear_model.LinearRegression.
任务3.2: 编写函数evaluate_model(), 实现预测图、 $R^{2}$ 与RMSE值的报告(作图用的函数plot_prediction()已经给出).
- 输入:
  - 模型model与数据集X, y_true;
  - 评估模式mode, 可在plot和metrics中二选一.
- 返回: 如果评估模式为plot, 则作预测图、不返回任何内容; 如果评估模式为metrics, 则返回RMSE与 $R^{2}$ 值.
任务3.3: 调用你编写的函数evaluate_model(), 完成下述任务:
- 将训练集上的预测图保存为文件pred.png, 以备提交;
- 将测试集上的RMSE值与 $R^{2}$ 值打印出来.

代码

文本

[ ]

from sklearn.linear_model import LinearRegression

def train_model(X_train, y_true_train) -> object:

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

代码

文本

完成任务3.1后, 请运行下面的代码块做初步检查. 可以看到, 训练集上的RMSE值应当与0很接近.

代码

文本

[ ]

### SANITY CHECK CELL 3.1 ###

### EXPECTED OUTPUT ###

'''

0.08926981854646464

'''

from sklearn.metrics import root_mean_squared_error as RMSE

lr = train_model(X_train, y_true_train)

RMSE(lr.predict(X_train), y_true_train)

代码

文本

检查完成后, 你可以继续完成任务3.2: 结合给出的函数plot_prediction()(已经给出, 取自上机实习1), 编写函数evaluate_model().

代码

文本

[ ]

from matplotlib import pyplot as plt

from sklearn.metrics import r2_score

def plot_prediction(y_true: np.array, y_pred: np.array):

r2 = r2_score(y_true, y_pred)

# text annotation setup

plt.title(r"True values vs predicted values ($R^2$ = " + f"{r2:.4f}" + ")")

plt.xlabel("True values")

plt.ylabel("Predicted values")

# plot the scatter and line

plt.scatter(y_true, y_pred, c="red", marker="o")

plt.plot(y_true, y_true, "b--")

# show the plot!

plt.show()

代码

文本

[ ]

def evaluate_model(model, X, y_true, mode: str):

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

代码

文本

完成该函数后, 你可以继续完成任务3.3: 运行下面的代码块. 我们先在训练集上作出预测图, 再报告测试集上的RMSE与 $R^{2}$ .

代码

文本

[ ]

evaluate_model(lr, X_train, y_true_train, mode="plot")

RMSE, r2 = evaluate_model(lr, X_test, y_true_test, mode="metrics")

print(f"RMSE: {RMSE}")

print(f"r2: {r2}")

代码

文本

结果分析

代码

文本

任务4: 运行下面的代码块, 根据参数正负, 考察各个因素对对映选择性的影响, 在注释的答题区中写下你的分析.

$π$ - $π$ 堆积作用如何影响对映选择性? (从几何结构与能量两方面考虑)
烷基的几何参数如何影响对映选择性? 和芳基几何参数相比, 总体看, 谁的影响更大?

代码

文本

提示

9列特征分别为:
- d_pi_w, e_pi_w, de_pi_w, 表征 $π$ - $π$ 堆积结构的距离或能量参数;
- L_Alk, B1_Alk, B5_Alk, 表征烷基的几何参数;
- L_Ar, B1_Ar, B5_Ar, 表征芳基的几何参数.
$ΔΔ G^{‡}$ 越大, 表明两种对映异构体的动力学活性相差越大, 于是对映选择性越好.

代码

文本

[ ]

lr.coef_

### 任务4答题区(另起一行时请记得加注释符号#) ###

### BEGIN YOUR SOLUTION ###

### END YOUR SOLUTION ###

代码

文本

机器学习及其在化学中的应用

已赞1

本文被以下合集收录

2023机器学习及其在化学中的应用

昌珺涵

更新于 2024-10-13

23 篇175 人关注

机器学习及其在化学中的应用

微信用户pYt0

更新于 2024-09-20

10 篇3 人关注