【Pythonコード付】RDkitで化合物の構造最適化

import rdkit
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.rdForceFieldHelpers import MMFFGetMoleculeForceField
import py3Dmol
import matplotlib.pyplot as plt

def optimize_molecule(smiles, force_field_type='UFF', max_steps=1000, force_tolerance=1e-3):
    """
    分子の構造最適化を行い、エネルギーの変化をプロットする関数。

    Args:
        mol (Mol): RDKitの分子オブジェクト。
        force_field_type (str): 使用する力場の種類（'UFF', 'MMFF'）。デフォルトは 'UFF'。
        max_steps (int): 最適化の最大ステップ数。デフォルトは1000。
        force_tolerance (float): 収束条件としての勾配の閾値。デフォルトは1e-3。
    """
    # 分子を読み込む
    mol = Chem.AddHs(Chem.MolFromSmiles(smiles))
    
    # 3D座標を生成する
    AllChem.EmbedMolecule(mol, randomSeed=42)
    
    if force_field_type == 'UFF':
        
        # 力場の作成（デフォルトはUFFを使用）
        force_field = AllChem.UFFGetMoleculeForceField(mol)
        
    elif force_field_type == 'MMFF':
        
        # MMFF用のプロパティを取得
        props = AllChem.MMFFGetMoleculeProperties(mol)

        # MMFF力場を作成
        force_field = MMFFGetMoleculeForceField(mol, props)

    # エネルギーの変化を記録するリスト
    energies = []

    # 構造最適化を実行し、エネルギーを記録
    for step in range(max_steps):
        # 最適化の1ステップ実行
        converged = force_field.Minimize(maxIts=1, forceTol=force_tolerance)

        # 現在のエネルギーを記録
        energy = force_field.CalcEnergy()
        energies.append(energy)
        print(f"Step {step}: Energy = {energy:.6f} kcal/mol")

        # 勾配（力）の確認
        max_force = max(force_field.CalcGrad())
        print(f"Maximum Force: {max_force:.6f} kcal/mol/Å")

        # 勾配が閾値以下であれば終了
        if max_force < force_tolerance:
            print("Optimization converged.")
            break

    # 最適化後のエネルギー
    final_energy = force_field.CalcEnergy()
    print(f"Final Energy: {final_energy:.6f} kcal/mol")

    # エネルギーの推移をプロット
    plt.figure(figsize=(8, 6))
    plt.plot(energies, marker='o', linestyle='-', color='b')
    plt.title(f'Energy Convergence During {force_field_type} Optimization')
    plt.xlabel('Optimization Step')
    plt.ylabel('Energy (kcal/mol)')
    plt.grid(True)
    plt.show()

    # 分子をスティックスタイルで表示
    mb = Chem.MolToMolBlock(mol)
    viewer = py3Dmol.view(width=400, height=400)
    viewer.addModel(mb, 'mol')
    viewer.setStyle({'stick': {}})
    viewer.zoomTo()
    viewer.show()

#計算する化合物のSMILESを定義
smiles = 'C1CCCCC1C'    

# 最適化の実行（MMFF力場使用）
optimize_molecule(smiles=smiles, force_field_type='MMFF')

# 分子を読み込む
mol = Chem.AddHs(Chem.MolFromSmiles(smiles))
    
# 3D座標を生成する
AllChem.EmbedMolecule(mol, randomSeed=42)

if force_field_type == 'UFF':
        
    # 力場の作成（デフォルトはUFFを使用）
    force_field = AllChem.UFFGetMoleculeForceField(mol)
        
elif force_field_type == 'MMFF':
        
    # MMFF用のプロパティを取得
    props = AllChem.MMFFGetMoleculeProperties(mol)

    # MMFF力場を作成
    force_field = MMFFGetMoleculeForceField(mol, props)

# エネルギーの変化を記録するリスト
    energies = []

    # 構造最適化を実行し、エネルギーを記録
    for step in range(max_steps):
        # 最適化の1ステップ実行
        converged = force_field.Minimize(maxIts=1, forceTol=force_tolerance)

        # 現在のエネルギーを記録
        energy = force_field.CalcEnergy()
        energies.append(energy)
        print(f"Step {step}: Energy = {energy:.6f} kcal/mol")

        # 勾配（力）の確認
        max_force = max(force_field.CalcGrad())
        print(f"Maximum Force: {max_force:.6f} kcal/mol/Å")

        # 勾配が閾値以下であれば終了
        if max_force < force_tolerance:
            print("Optimization converged.")
            break

    # 最適化後のエネルギー
    final_energy = force_field.CalcEnergy()
    print(f"Final Energy: {final_energy:.6f} kcal/mol")

    # エネルギーの推移をプロット
    plt.figure(figsize=(8, 6))
    plt.plot(energies, marker='o', linestyle='-', color='b')
    plt.title(f'Energy Convergence During {force_field_type} Optimization')
    plt.xlabel('Optimization Step')
    plt.ylabel('Energy (kcal/mol)')
    plt.grid(True)
    plt.show()

    # 分子をスティックスタイルで表示
    mb = Chem.MolToMolBlock(mol)
    viewer = py3Dmol.view(width=400, height=400)
    viewer.addModel(mb, 'mol')
    viewer.setStyle({'stick': {}})
    viewer.zoomTo()
    viewer.show()