Python rdkit.Chem.Conformer() Examples

def xyz2AC_vdW(atoms, xyz):

    # Get mol template
    mol = get_proto_mol(atoms)

    # Set coordinates
    conf = Chem.Conformer(mol.GetNumAtoms())
    for i in range(mol.GetNumAtoms()):
        conf.SetAtomPosition(i, (xyz[i][0], xyz[i][1], xyz[i][2]))
    mol.AddConformer(conf)

    AC = get_AC(mol)

    return AC, mol 

def initialize_conformer(self, conf):
        """Retrieve atom coordinates and instantiate shells generator.

        Parameters
        ----------
        conf : RDKit Conformer
            Conformer to fingerprint
        """
        self.conf = conf
        self.atom_coords = coords_from_atoms(self.atoms, self.conf)
        self.shells_gen = ShellsGenerator(
            self.conf,
            self.atoms,
            radius_multiplier=self.radius_multiplier,
            atom_coords=self.atom_coords,
            include_disconnected=self.include_disconnected,
            bound_atoms_dict=self.bound_atoms_dict,
        ) 

def coords_from_atoms(atoms, conf):
    """Build `dict` matching atom id to coordinates.

    Parameters
    ----------
    atoms : list of int
        Atom ids
    conf : RDKit Conformer
        Conformer from which to fetch coordinates

    Returns
    -------
    dict : Dict matching atom id to 1-D array of coordinates.
    """
    coordinates = [
        np.array(conf.GetAtomPosition(int(x)), dtype=np.float64) for x in atoms
    ]
    return dict(zip(atoms, coordinates)) 

def get_lowest_confs(label: str,
                     confs: dict or list,
                     n: int = 10,
                     e: float = 5.0,
                     energy: str = 'FF energy',
                     ) -> list:
    """
    Get the most stable conformer

    Args:
        label (str): The species' label.
        confs (dict, list): Entries are either conformer dictionaries or a length two list of xyz coordinates and energy
        n (int, optional): Number of lowest conformers to return.
        e (float, optional): The energy threshold above the lowest energy conformer in kJ/mol
                             below which all conformers will be returned.
        energy (str, optional): The energy attribute to search by. Currently only 'FF energy' is supported.

    Raises:
        ConformerError: If n < 1, e < 0, both n and e are ``None``, or if no conformers are given.

    Returns:
        list: Conformer dictionaries.
    """
    if e is not None:
        if e < 0:
            raise ConformerError(f'e cannot be negative, got: {e}')
    elif n is not None:
        if n < 1:
            raise ConformerError(f'n cannot be lower than 1, got: {n}')
    else:
        raise ConformerError(f'Either n or e must be specified')
    if not confs or confs is None:
        raise ConformerError(f'get_lowest_confs() got no conformers for {label}')

    if isinstance(confs[0], list):
        conformer_list = list()
        for entry in confs:
            if entry[1] is not None:
                conformer_list.append({'xyz': entry[0], energy: entry[1]})
    elif isinstance(confs[0], dict):
        conformer_list = [conformer for conformer in confs if energy in conformer and conformer[energy] is not None]
    else:
        raise ConformerError(f'confs could either be a list of dictionaries or a list of lists. '
                             f'Got a list of {type(confs[0])}s for {label}')

    conformer_list.sort(key=lambda conformer: conformer[energy], reverse=False)
    if e is not None:
        min_e = min([conf[energy] for conf in conformer_list])
    lowest_confs = [conformer_list[0]]
    for index in range(len(conformer_list)):
        if (e is not None and conformer_list[index][energy] > min_e + e) or (n is not None and len(lowest_confs) >= n):
            break
        if index > 0 and not any([converter.compare_confs(lowest_conf['xyz'], conformer_list[index]['xyz'])
                                  for lowest_conf in lowest_confs]):
            lowest_confs.append(conformer_list[index])
    return lowest_confs 

def embed_rdkit(label, mol, num_confs=None, xyz=None):
    """
    Generate unoptimized conformers in RDKit. If ``xyz`` is not given, random conformers will be generated.

    Args:
        label (str): The species' label.
        mol (RMG Molecule or RDKit RDMol): The molecule object with connectivity and bond order information.
        num_confs (int, optional): The number of random 3D conformations to generate.
        xyz (dict, optional): The 3D coordinates.

    Returns:
        RDMol: An RDKIt molecule with embedded conformers.
    """
    if num_confs is None and xyz is None:
        raise ConformerError(f'Either num_confs or xyz must be set when calling embed_rdkit() for {label}')
    if isinstance(mol, RDMol):
        rd_mol = mol
    elif isinstance(mol, Molecule):
        rd_mol = converter.to_rdkit_mol(mol=mol, remove_h=False)
    else:
        raise ConformerError(f'Argument mol can be either an RMG Molecule or an RDKit RDMol object. '
                             f'Got {type(mol)} for {label}')
    if num_confs is not None:
        Chem.AllChem.EmbedMultipleConfs(rd_mol, numConfs=num_confs, randomSeed=1, enforceChirality=True)
        # Chem.AllChem.EmbedMultipleConfs(rd_mol, numConfs=num_confs, randomSeed=15, enforceChirality=False)
    elif xyz is not None:
        rd_conf = Chem.Conformer(rd_mol.GetNumAtoms())
        for i in range(rd_mol.GetNumAtoms()):
            rd_conf.SetAtomPosition(i, xyz['coords'][i])
        rd_mol.AddConformer(rd_conf)
    return rd_mol 

def _process_molecule_rdkit(jmol):
        from rdkit import Chem

        # Handle errors
        if abs(jmol.molecular_charge) > 1.0e-6:
            raise InputError("RDKit does not currently support charged molecules.")

        if not jmol.connectivity:  # Check for empty list
            raise InputError("RDKit requires molecules to have a connectivity graph.")

        # Build out the base molecule
        base_mol = Chem.Mol()
        rw_mol = Chem.RWMol(base_mol)
        for sym in jmol.symbols:
            rw_mol.AddAtom(Chem.Atom(sym.title()))

        # Add in connectivity
        bond_types = {1: Chem.BondType.SINGLE, 2: Chem.BondType.DOUBLE, 3: Chem.BondType.TRIPLE}
        for atom1, atom2, bo in jmol.connectivity:
            rw_mol.AddBond(atom1, atom2, bond_types[bo])

        mol = rw_mol.GetMol()

        # Write out the conformer
        natom = len(jmol.symbols)
        conf = Chem.Conformer(natom)
        bohr2ang = ureg.conversion_factor("bohr", "angstrom")
        for line in range(natom):
            conf.SetAtomPosition(
                line,
                (
                    bohr2ang * jmol.geometry[line, 0],
                    bohr2ang * jmol.geometry[line, 1],
                    bohr2ang * jmol.geometry[line, 2],
                ),
            )

        mol.AddConformer(conf)
        Chem.rdmolops.SanitizeMol(mol)

        return mol 

def xyz2AC_huckel(atomicNumList,xyz,charge):
    """

    args
        atomicNumList - atom type list
        xyz - coordinates
        charge - molecule charge

    returns
        ac - atom connectivity
        mol - rdkit molecule

    """
    mol = get_proto_mol(atomicNumList)

    conf = Chem.Conformer(mol.GetNumAtoms())
    for i in range(mol.GetNumAtoms()):
        conf.SetAtomPosition(i,(xyz[i][0],xyz[i][1],xyz[i][2]))
    mol.AddConformer(conf)

    num_atoms = len(atomicNumList)
    AC = np.zeros((num_atoms,num_atoms)).astype(int)

    mol_huckel = Chem.Mol(mol)
    mol_huckel.GetAtomWithIdx(0).SetFormalCharge(charge) #mol charge arbitrarily added to 1st atom    

    passed,result = rdEHTTools.RunMol(mol_huckel)
    opop = result.GetReducedOverlapPopulationMatrix()
    tri = np.zeros((num_atoms, num_atoms))
    tri[np.tril(np.ones((num_atoms, num_atoms), dtype=bool))] = opop #lower triangular to square matrix
    for i in range(num_atoms):
        for j in range(i+1,num_atoms):
            pair_pop = abs(tri[j,i])   
            if pair_pop >= 0.15: #arbitry cutoff for bond. May need adjustment
                AC[i,j] = 1
                AC[j,i] = 1

    return AC, mol 

def AtomListToSubMol(mol, amap, includeConformer=False):
    """
    Parameters
    ----------
        mol: rdkit.Chem.rdchem.Mol
            Molecule
        amap: array-like
            List of atom indices (zero-based)
        includeConformer: bool (default=True)
            Toogle to include atoms coordinates in submolecule.

    Returns
    -------
        submol: rdkit.Chem.rdchem.RWMol
            Submol determined by specified atom list
    """
    if not isinstance(amap, list):
        amap = list(amap)
    submol = Chem.RWMol(Chem.Mol())
    for aix in amap:
        submol.AddAtom(mol.GetAtomWithIdx(aix))
    for i, j in combinations(amap, 2):
        bond = mol.GetBondBetweenAtoms(i, j)
        if bond:
            submol.AddBond(amap.index(i),
                           amap.index(j),
                           bond.GetBondType())
    if includeConformer:
        for conf in mol.GetConformers():
            new_conf = Chem.Conformer(len(amap))
            for i in range(len(amap)):
                new_conf.SetAtomPosition(i, conf.GetAtomPosition(amap[i]))
                new_conf.SetId(conf.GetId())
                new_conf.Set3D(conf.Is3D())
            submol.AddConformer(new_conf)
    return submol 

def xyz2AC(atomicNumList,xyz):
    import numpy as np
    mol = get_proto_mol(atomicNumList)

    conf = Chem.Conformer(mol.GetNumAtoms())
    for i in range(mol.GetNumAtoms()):
        conf.SetAtomPosition(i,(xyz[i][0],xyz[i][1],xyz[i][2]))
    mol.AddConformer(conf)

    dMat = Chem.Get3DDistanceMatrix(mol)
    pt = Chem.GetPeriodicTable()

    num_atoms = len(atomicNumList)
    AC = np.zeros((num_atoms,num_atoms)).astype(int)

    for i in range(num_atoms):
        a_i = mol.GetAtomWithIdx(i)
        Rcov_i = pt.GetRcovalent(a_i.GetAtomicNum())*1.30
        for j in range(i+1,num_atoms):
            a_j = mol.GetAtomWithIdx(j)
            Rcov_j = pt.GetRcovalent(a_j.GetAtomicNum())*1.30
            if dMat[i,j] <= Rcov_i + Rcov_j:
                AC[i,j] = 1
                AC[j,i] = 1

    return AC,mol 

def run(self, conf=None, mol=None, return_substruct=False):
        """Generate fingerprint from provided conformer or mol and conf id.

        Parameters
        ----------
        conf : RDKit Conformer or int, optional
            Input conformer or conformer in `mol`.
        mol : RDKit Mol, optional
            Input molecule object, with at least one conformer. If `conf` not
            specified, first conformer is used.
        return_substruct : bool, optional
            Return dict mapping substructure to fingerprint indices. Keys are
            indices, values are list of substructures, represented as a tuple
            of atom indices where the first index is the central atom and the
            remaining indices (within the sphere) are sorted.
        """
        if mol is None:  # mol not provided; get from conf
            try:
                mol = conf.GetOwningMol()
            except AttributeError:  # conf is int ID; use existing mol
                mol = self.mol
        else:
            if not isinstance(conf, Chem.Conformer):
                try:
                    conf = mol.GetConformer(conf)
                except TypeError:  # conf isn't ID either. Fall back to first
                    conf = mol.GetConformer(0)

        if mol is not self.mol:
            self.reset_mol()
            self.initialize_mol(mol)
        elif conf is not self.conf:
            self.reset_conf()

        self.initialize_conformer(conf)

        for i in iter(self):
            pass 

def determine_chirality(conformers, label, mol, force=False):
    """
    Determines the Cahn–Ingold–Prelog (CIP) chirality (R or S) of atoms in the conformer,
    as well as the CIP chirality of double bonds (E or Z).

    Args:
        conformers (list): Entries are conformer dictionaries.
        label (str): The species' label.
        mol (RMG Molecule or RDKit RDMol): The molecule object with connectivity and bond order information.
        force (bool, optional): Whether to override data, ``True`` to override, default is ``False``.

    Returns:
        list: Conformer dictionaries with updated with 'chirality'. ``conformer['chirality']`` is a dictionary.
              Keys are either a 1-length tuple of atom indices (for chiral atom centers) or a 2-length tuple of atom
              indices (for chiral double bonds), values are either 'R' or 'S' for chiral atom centers
              (or 'NR' or 'NS' for chiral nitrogen centers), or 'E' or 'Z' for chiral double bonds.
              All atom indices are 0-indexed.
    """
    chiral_nitrogen_centers = identify_chiral_nitrogen_centers(mol)
    new_mol, elements_to_insert = replace_n_with_c_in_mol(mol, chiral_nitrogen_centers)
    for conformer in conformers:
        if 'chirality' not in conformer:
            # keys are either 1-length atom indices (for chiral atom centers)
            # or 2-length atom indices (for chiral double bonds)
            # values are either 'R', 'S', 'NR', 'NS', 'E', or 'Z'
            conformer['chirality'] = dict()
        elif conformer['chirality'] != dict() and not force:
            # don't override data
            continue
        new_xyz = replace_n_with_c_in_xyz(label, mol, conformer['xyz'], chiral_nitrogen_centers, elements_to_insert)
        rd_mol = embed_rdkit(label, new_mol, xyz=new_xyz)
        Chem.rdmolops.AssignStereochemistryFrom3D(rd_mol, 0)
        for i, rd_atom in enumerate(rd_mol.GetAtoms()):
            rd_atom_props_dict = rd_atom.GetPropsAsDict()
            if '_CIPCode' in list(rd_atom_props_dict.keys()):
                if mol.atoms[i].is_nitrogen():
                    # this is a nitrogen site in the original molecule, mark accordingly
                    conformer['chirality'][(i,)] = 'N' + rd_atom_props_dict['_CIPCode']
                else:
                    conformer['chirality'][(i,)] = rd_atom_props_dict['_CIPCode']
        for rd_bond in rd_mol.GetBonds():
            stereo = str(rd_bond.GetStereo())
            if stereo in ['STEREOE', 'STEREOZ']:
                # possible values are 'STEREOANY', 'STEREOCIS', 'STEREOE', 'STEREONONE', 'STEREOTRANS', and 'STEREOZ'
                rd_atoms = [rd_bond.GetBeginAtomIdx(), rd_bond.GetEndAtomIdx()]  # indices of atoms bonded by this bond
                conformer['chirality'][tuple(rd_atom for rd_atom in rd_atoms)] = stereo[-1]
    return conformers