rdkitで特定の原子に構造を付加するスクリプト

from rdkit import Chem
from rdkit.Chem import AllChem

# PDBファイルから分子を読み込む
mol = Chem.MolFromPDBFile('path_to_your_pdb_file.pdb')

# 原子を識別する。ここでは、例として原子番号5をターゲットとします。
target_atom = mol.GetAtomWithIdx(4)  # 0-based index

# ターゲット原子に結合している水素の一つを削除
for neighbor in target_atom.GetNeighbors():
    if neighbor.GetSymbol() == 'H':
        mol.RemoveBond(target_atom.GetIdx(), neighbor.GetIdx())
        mol.RemoveAtom(neighbor.GetIdx())
        break

# COOHを持つフラグメントを作成
cooh = Chem.MolFromSmiles("COOH")
cooh = Chem.AddHs(cooh)
AllChem.EmbedMolecule(cooh)

# 結合
combo = Chem.CombineMols(mol, cooh)

# 新しい結合を作成
combo.AddBond(target_atom.GetIdx(), mol.GetNumAtoms(), Chem.BondType.SINGLE)

# 最適化
AllChem.MMFFOptimizeMolecule(combo)

# 新しいPDBファイルとして保存
Chem.MolToPDBFile(combo, 'modified_structure.pdb')

from rdkit import Chem
from rdkit.Chem import AllChem

# 2つのPDBファイルから分子を読み込む
mol1 = Chem.MolFromPDBFile('path_to_your_first_pdb_file.pdb')
mol2 = Chem.MolFromPDBFile('path_to_your_second_pdb_file.pdb')

# 結合させたい原子を識別する
# 例: mol1の原子番号5、mol2の原子番号7 (0-based index)
atom1 = mol1.GetAtomWithIdx(4)
atom2 = mol2.GetAtomWithIdx(6)

# 分子を結合
combo = Chem.CombineMols(mol1, mol2)

# 新しい結合を作成
combo.AddBond(atom1.GetIdx(), atom2.GetIdx() + mol1.GetNumAtoms(), Chem.BondType.SINGLE)

# 最適化
AllChem.MMFFOptimizeMolecule(combo)

# 新しいPDBファイルとして保存
Chem.MolToPDBFile(combo, 'combined_structure.pdb')

from rdkit import Chem
from rdkit.Chem import AllChem

# 2つのPDBファイルから分子を読み込む
mol1 = Chem.MolFromPDBFile('path_to_your_first_pdb_file.pdb')
mol2 = Chem.MolFromPDBFile('path_to_your_second_pdb_file.pdb')

# 結合させたい原子を識別する
atom1 = mol1.GetAtomWithIdx(4)
atom2 = mol2.GetAtomWithIdx(6)

# 分子を結合
combo = Chem.CombineMols(mol1, mol2)

# 新しい結合を作成
combo.AddBond(atom1.GetIdx(), atom2.GetIdx() + mol1.GetNumAtoms(), Chem.BondType.SINGLE)

# 結合原子の近傍の原子を特定（例：距離2の原子まで考慮）
affected_atoms = set([atom1.GetIdx(), atom2.GetIdx() + mol1.GetNumAtoms()])
for bond in atom1.GetBonds():
    affected_atoms.add(bond.GetOtherAtom(atom1).GetIdx())
for bond in atom2.GetBonds():
    affected_atoms.add(bond.GetOtherAtom(atom2).GetIdx() + mol1.GetNumAtoms())

# ターゲットとなる原子群のみを最適化
ff = AllChem.MMFFGetMoleculeForceField(combo)
for i in range(combo.GetNumAtoms()):
    if i not in affected_atoms:
        ff.AddFixedPoint(i)
ff.Minimize()

# 最適化後の構造を保存
Chem.MolToPDBFile(combo, 'locally_optimized_structure.pdb')