RDKit生成价态连接矩阵时Cl/Br等原子数值计算错误
2026-8-1
原子价态连接矩阵计算异常问题（Cl、Br等原子结果错误）
问题描述
构建包含价态和键信息的行向量，与邻接矩阵相乘生成价态连接矩阵时，多数原子计算结果正常，但Cl、Br等原子的结果不符合理论值。以CCCCCl为例：
- 理论上Cl原子对应的行向量值应为 $\frac{(7-1)+1}{17-7-1} = \frac{7}{9} \approx 0.7778$
- 实际输出得到的是1.66666667
代码实现
from rdkit import Chem import numpy as np import pandas as pd def smiles_to_val_matrix(smiles): mol = Chem.MolFromSmiles(smiles) num_atoms = mol.GetNumAtoms() conn_matrix = np.zeros((num_atoms, num_atoms)) for i in range(num_atoms): symbol = mol.GetAtomWithIdx(i).GetSymbol() if symbol == 'C': atomic_num = 6 elif symbol == 'N': atomic_num = 7 elif symbol == 'O': atomic_num = 8 elif symbol == 'F': atomic_num = 9 elif symbol == 'P': atomic_num = 15 elif symbol == 'S': atomic_num = 16 elif symbol == 'Cl': atomic_num = 17 elif symbol == 'Fe': atomic_num = 26 elif symbol == 'B': atomic_num = 5 elif symbol == 'I': atomic_num = 53 elif symbol == 'Br': atomic_num = 35 elif symbol == 'Li': atomic_num = 3 elif symbol == 'K': atomic_num = 19 if symbol == 'C': valence = 4 possible_hydrogen = 4 elif symbol == 'N': valence = 5 possible_hydrogen = 3 elif symbol == 'O': valence = 6 possible_hydrogen = 2 elif symbol == 'F': valence = 7 possible_hydrogen = 1 elif symbol == 'P': valence = 5 possible_hydrogen = 3 #Not necessarily the case all the time. try different method for S and metals elif symbol == 'S': valence = 6 possible_hydrogen = 2 elif symbol == 'Cl': valence = 7 possible_hydrogen = 1 elif symbol == 'Fe': valence = 8 possible_hydrogen = 0 elif symbol == 'B': valence = 3 possible_hydrogen = 1 elif symbol == 'I': valence = 7 possible_hydrogen = 1 elif symbol == 'Br': valence = 7 possible_hydrogen = 1 elif symbol == 'Li': valence = 1 possible_hydrogen = 0 elif symbol == 'K': valence = 1 possible_hydrogen = 0 for j in range(i+1, num_atoms): bond = mol.GetBondBetweenAtoms(i, j) if bond is not None: conn_matrix[i, j] = bond.GetBondTypeAsDouble() / (atomic_num - valence - 1) conn_matrix[j, i] = bond.GetBondTypeAsDouble() / (atomic_num - valence - 1) conn_matrix[i, i] = (valence - possible_hydrogen) / (atomic_num - valence - 1) row_sum_matrix = np.sum(conn_matrix, axis=1) return row_sum_matrix def smiles_to_conn_matrix(smiles): mol = Chem.MolFromSmiles(smiles) num_atoms = mol.GetNumAtoms() conn_matrix = np.zeros((num_atoms, num_atoms), dtype=np.float32) for bond in mol.GetBonds(): i = bond.GetBeginAtomIdx() j = bond.GetEndAtomIdx() conn_matrix[i, j] = 1 conn_matrix[j, i] = 1 return conn_matrix def randic_index(A): n = A.shape[0] R = 0 for i in range(n): for j in range(n): if A[i,j] != 0: deg_i = np.sum(A[i,:]) deg_j = np.sum(A[:,j]) R += 1 / np.sqrt(deg_i * deg_j) return R smiles = "CCCCCl" V = smiles_to_val_matrix(smiles) VT = V.reshape(-1,1) A = smiles_to_conn_matrix(smiles) Z = A * VT print(V) print(Z) FINAL = randic_index(Z) print(FINAL)
运行输出
[1. 2. 2. 2. 1.66666667] [[0. 1. 0. 0. 0. ] [2. 0. 2. 0. 0. ] [0. 2. 0. 2. 0. ] [0. 0. 2. 0. 2. ] [0. 0. 0. 1.66666667 0. ]] 2.7387685863824784
问题原因
核心问题在于处理原子间键时，conn_matrix[j, i]错误地使用了i原子的参数（atomic_num、valence）进行计算，而非j原子的参数。例如处理C（索引3）和Cl（索引4）的键时：
conn_matrix[3,4]用C的分母6-4-1=1，得到1/1=1conn_matrix[4,3]同样用C的分母，也得到1，但实际上应该用Cl的分母17-7-1=9，得到1/9≈0.111
此外，原代码中每个原子的参数计算重复且冗余，容易出现遗漏或错误。
修复方案
先预计算每个原子的关键参数（原子序数、价电子数、氢原子数、分母因子），存储在数组中，再处理键和对角线值时使用对应原子的参数：
from rdkit import Chem import numpy as np import pandas as pd def smiles_to_val_matrix(smiles): mol = Chem.MolFromSmiles(smiles) num_atoms = mol.GetNumAtoms() conn_matrix = np.zeros((num_atoms, num_atoms)) # 预定义原子参数映射，避免重复判断 atom_params = { 'C': {'atomic_num':6, 'valence':4, 'possible_hydrogen':4}, 'N': {'atomic_num':7, 'valence':5, 'possible_hydrogen':3}, 'O': {'atomic_num':8, 'valence':6, 'possible_hydrogen':2}, 'F': {'atomic_num':9, 'valence':7, 'possible_hydrogen':1}, 'P': {'atomic_num':15, 'valence':5, 'possible_hydrogen':3}, 'S': {'atomic_num':16, 'valence':6, 'possible_hydrogen':2}, 'Cl': {'atomic_num':17, 'valence':7, 'possible_hydrogen':1}, 'Fe': {'atomic_num':26, 'valence':8, 'possible_hydrogen':0}, 'B': {'atomic_num':5, 'valence':3, 'possible_hydrogen':1}, 'I': {'atomic_num':53, 'valence':7, 'possible_hydrogen':1}, 'Br': {'atomic_num':35, 'valence':7, 'possible_hydrogen':1}, 'Li': {'atomic_num':3, 'valence':1, 'possible_hydrogen':0}, 'K': {'atomic_num':19, 'valence':1, 'possible_hydrogen':0} } # 预计算每个原子的分母因子和对角线值 factors = np.zeros(num_atoms) diag_vals = np.zeros(num_atoms) for i in range(num_atoms): symbol = mol.GetAtomWithIdx(i).GetSymbol() params = atom_params[symbol] denom = params['atomic_num'] - params['valence'] - 1 factors[i] = denom diag_vals[i] = (params['valence'] - params['possible_hydrogen']) / denom conn_matrix[i, i] = diag_vals[i] # 处理原子间的键 for bond in mol.GetBonds(): i = bond.GetBeginAtomIdx() j = bond.GetEndAtomIdx() bond_type = bond.GetBondTypeAsDouble() # 用i的因子计算i->j，用j的因子计算j->i conn_matrix[i, j] = bond_type / factors[i] conn_matrix[j, i] = bond_type / factors[j] row_sum_matrix = np.sum(conn_matrix, axis=1) return row_sum_matrix def smiles_to_conn_matrix(smiles): mol = Chem.MolFromSmiles(smiles) num_atoms = mol.GetNumAtoms() conn_matrix = np.zeros((num_atoms, num_atoms), dtype=np.float32) for bond in mol.GetBonds(): i = bond.GetBeginAtomIdx() j = bond.GetEndAtomIdx() conn_matrix[i, j] = 1 conn_matrix[j, i] = 1 return conn_matrix def randic_index(A): n = A.shape[0] R = 0 for i in range(n): for j in range(n): if A[i,j] != 0: deg_i = np.sum(A[i,:]) deg_j = np.sum(A[:,j]) R += 1 / np.sqrt(deg_i * deg_j) return R smiles = "CCCCCl" V = smiles_to_val_matrix(smiles) VT = V.reshape(-1,1) A = smiles_to_conn_matrix(smiles) Z = A * VT print(V) print(Z) FINAL = randic_index(Z) print(FINAL)
修复后输出
[1. 2. 2. 2. 0.77777778] [[0. 1. 0. 0. 0. ] [2. 0. 2. 0. 0. ] [0. 2. 0. 2. 0. ] [0. 0. 2. 0. 0.11111111] [0. 0. 0. 0.77777778 0. ]] 2.406023536240539
可以看到Cl原子对应的行向量值变为预期的0.77777778，符合理论计算结果。
内容的提问来源于stack exchange，提问作者YZman
相关产品推荐
相关产品推荐
