RDKitで拘束を入れて構造最適化をする
RDKitの構造最適化
前回の記事で、RDKit上で複数分子を取り扱い、移動や回転を行ってみました。ただし、ここでの分子操作はあくまでも幾何的に行ったものであり化学的に正しいとは限りません。そのため、こういった化学構造を初期構造として量子化学計算を行うとうまく収束しないことが多いです。うまく構造を収束させるために低レベル基底や汎関数での事前の最適化や分子力場による構造最適化を行います。
RDKitにも分子力場による構造最適化の機能は備わっています。しかし、よくある公式リファレンスのやり方だと構造の拘束がないので、二量体のように特定の原子間の位置関係を制御した状態で最適化することはできません。
・・・ですが、実はRDKitにも構造最適化時に特定の原子の座標や結合長、角度の拘束を与えるメソッドがあります。今回は前回紹介した二量体の構造をいい感じに最適化してみます。
二量体の準備
今回は酢酸の二量体を題材に扱います。酢酸は日常生活でもよく利用・接種する物質の一つですが、
酢酸のようなカルボン酸は液相・気相で二量体を形成することが知られています。
二量体を作成するにはまず二つの酢酸を準備する必要があります。これは下記のようなコードで実現できます。下記では一つ目の酢酸のカルボニルのOと二つ目の酢酸のカルボキシル基のHをx軸に沿って対向するように回転し、両分子の重心が6Å離れるようにしています。この操作は後の構造最適化をうまく進めるために行っていますが、構造最適化の条件によってはこういった座標操作は行わなくても良い構造ができる場合もあります。なお、分子の回転操作と平行移動についてはこちらを参照ください。
def rotate_molecule(position_df, target_atom_id, vector):
def rotation_matrix(frm, to):
if np.array_equal(frm, to):
return np.identity(3)
if np.array_equal(frm, -to):
return -np.identity(3)
s = frm + to
return 2.0 * np.outer(s, s) / np.dot(s, s) - np.identity(3)
position_df[['x', 'y', 'z']] -= position_df[['x', 'y', 'z']].mean()
atom_vector = position_df.iloc[target_atom_id].values
print(atom_vector)
atom_vector = atom_vector/np.linalg.norm(atom_vector)
rot_matrix = rotation_matrix(atom_vector[0], vector)
for i, value in enumerate(position_df[['x', 'y', 'z']].values):
position_df.at[i, 'x'] = np.dot(rot_matrix, value)[0]
position_df.at[i, 'y'] = np.dot(rot_matrix, value)[1]
position_df.at[i, 'z'] = np.dot(rot_matrix, value)[2]
return position_df
from rdkit.Geometry import Point3D
import pandas as pd
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
import py3Dmol
smiles = "CC(=O)O"
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
AllChem.EmbedMolecule(mol)
mols = Chem.CombineMols(mol, mol)
conf = mols.GetConformer()
for c in mols.GetConformers():
position = pd.DataFrame(c.GetPositions(), columns=['x', 'y', 'z'])
mol1_pos = position.iloc[:len(mol.GetAtoms())].copy().reset_index(drop=True)
mol2_pos = position.iloc[len(mol.GetAtoms()):].copy().reset_index(drop=True)
mol1_pos = rotate_molecule(mol1_pos, 7, np.array([1, 0, 0]))
mol2_pos = rotate_molecule(mol2_pos, 2, np.array([-1, 0, 0]))
rot_position = pd.concat([mol1_pos-3, mol2_pos+3]).reset_index(drop=True)
for c in mols.GetConformers():
for i, p in enumerate(c.GetPositions()):
rot_p = rot_position.iloc[i].values
conf.SetAtomPosition(i, Point3D(rot_p[0], rot_p[1], rot_p[2]))
v = py3Dmol.view(width=400, height=400)
mb = Chem.MolToMolBlock(mols)
v.addModel(mb, 'sdf')
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()
拘束を入れた構造最適化
カルボン酸の二量体はカルボキシル基の相互作用により６員環のような構造を形成することが知られてます。二つの分子から６員環の構造を形成するには、以下のように分子間の位置関係に構造に制約を与える必要があります。
- カルボキシル基のプロトンとカルボニルのOの距離
- カルボキシル基のプロトンとプロトンと結合しているO、及びカルボニルのOとの角度
- カルボキシル基のCとプロトンと結合しているO、及びカルボニルのOとの角度
RDKitにはUFFAddDistanceConstraint()のように構造最適化を行うメソッドがあり、
これを組み合わせることで任意の拘束を与えることが可能です。これらのメソッドは基本的に拘束対象の原子と拘束の範囲、拘束エネルギーを指定します。例えばUFFAddDistanceConstraint()の場合、
UFFAddDistanceConstraint(原子i, 原子j, (相対距離か絶対距離か), 最小距離, 最大距離, 拘束エネルギー)
のように指定します。そのため、どの原子間に拘束を与えるべきかは事前に原子のインデックスを確認するなどしておく必要があります。今回の酢酸の場合、各原子のインデックスは次のようになっています。これらを踏まえて、拘束を定義します。これを間違えるとおかしな構造になるので注意が必要です。
これらの拘束を与えた状態で構造最適化を行うコードは下記のとおりです。距離や角度、拘束エネルギーについては特に吟味していないので適当な値です。
mols.UpdatePropertyCache()
Chem.GetSymmSSSR(mols)
mols.GetRingInfo().NumRings()
ff = AllChem.UFFGetMoleculeForceField(mols)
# カルボキシル基のプロトンとカルボニルのOの距離の拘束
ff.UFFAddDistanceConstraint(7, 10, False, 1.8, 1.8, 50.0)
ff.UFFAddDistanceConstraint(2, 15, False, 1.8, 1.8, 50.0)
# カルボキシル基のプロトンとプロトンと結合しているO、及びカルボニルのOとの角度
ff.UFFAddAngleConstraint(3, 7, 10, False, 180, 180, 300.0)
ff.UFFAddAngleConstraint(2, 15, 11, False, 180, 180, 300.0)
# カルボキシル基のCとプロトンと結合しているO、及びカルボニルのOとの角度
ff.UFFAddAngleConstraint(1, 2, 15, False, 120, 180, 300.0)
ff.Initialize()
ff.Minimize(maxIts=1000)
v = py3Dmol.view(width=400, height=400)
mb = Chem.MolToMolBlock(mols)
v.addModel(mb, 'sdf')
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()
それっぽい構造ができてますね。
Discussion