OpenBabel构造小分子二聚体

from openbabel import pybel, openbabel

mol1 = pybel.readstring('smi','CC(C)C')
mol2 = pybel.readstring('smi','COC')

def have_h(atom):
    '''
    check if the provided atom has at least one hydrogen
    '''
    for neighbour_atom in openbabel.OBAtomAtomIter(atom.OBAtom):
        if neighbour_atom.GetAtomicNum() == 1:
            return True
    return False

def copy_mol(dst,src):
    '''
    copy atoms and related bond from `src` to `dst`
    '''
    dst_ = dst.clone
    shift = len(dst.atoms)
    for idx, atom in enumerate(src.atoms):
        new = dst_.OBMol.NewAtom()
        new.Duplicate(atom.OBAtom)
        new.SetIdx(idx + 1 + shift)
        new.SetId(idx + 1 + shift)
    for bond in openbabel.OBMolBondIter(src.OBMol):
        dst_.OBMol.AddBond(
            bond.GetBeginAtomIdx() + shift,
            bond.GetEndAtomIdx() + shift,
            bond.GetBondOrder()
        )
    return dst_

def lower_h(atom):
    '''
    make the implicit hydrogen of provided atom decrease by 1
    '''
    atom.OBAtom.SetImplicitHCount(
        atom.OBAtom.GetImplicitHCount() - 1
    )

def join(mol1,mol2):
    '''
    create a dimer of mol1 by using mol2 as a linker
    each atom having at least one hydrogen will be selected
    as an anchor for further linkage.
    the first and the last atom of mol2 will be used.
    '''
    mol = copy_mol(copy_mol(mol1,mol1),mol2)
    n_mole = len(mol1.atoms)
    n_link = len(mol2.atoms)
    mol1_ = mol1.clone
    mol1_.addh()
    ret = set()
    for idx, atom in enumerate([i for i in mol1_.atoms if i.atomicnum > 1]):
        if have_h(atom):
            mol_ = mol.clone
            mol_.OBMol.AddBond(idx+1,2*n_mole+1,1)
            mol_.OBMol.AddBond(idx+n_mole+1,2*n_mole+n_link,1)
            lower_h(mol_.atoms[idx])
            lower_h(mol_.atoms[2*n_mole])
            lower_h(mol_.atoms[idx+n_mole])
            lower_h(mol_.atoms[2*n_mole+n_link-1])
            ret.add(mol_.write('can')[:-3])
    return ret