Python rdkit.Chem.AllChem.EmbedMolecule() Examples

def generate_structure_from_smiles(smiles):

    # Generate a 3D structure from smiles

    mol = Chem.MolFromSmiles(smiles)
    mol = Chem.AddHs(mol)

    status = AllChem.EmbedMolecule(mol)
    status = AllChem.UFFOptimizeMolecule(mol)

    conformer = mol.GetConformer()
    coordinates = conformer.GetPositions()
    coordinates = np.array(coordinates)

    atoms = get_atoms(mol)

    return atoms, coordinates 

def make3D(self, forcefield="mmff94", steps=50):
        """Generate 3D coordinates.

        Optional parameters:
           forcefield -- default is "uff". See the forcefields variable
                         for a list of available forcefields.
           steps -- default is 50

        Once coordinates are generated, a quick
        local optimization is carried out with 50 steps and the
        UFF forcefield. Call localopt() if you want
        to improve the coordinates further.
        """
        forcefield = forcefield.lower()
        success = AllChem.EmbedMolecule(self.Mol,
                                        useExpTorsionAnglePrefs=True,
                                        useBasicKnowledge=True,
                                        enforceChirality=True,
                                        )
        if success == -1:
            raise Exception("Embedding failed!")

        self.localopt(forcefield, steps)
        self._clear_cache() 

def get_3d_mols(smiles_strs):
    """
    Convert SMILES strings to Mol objects with explicit hydrogens and 3D coordinates

    Args:
        smiles_strs (iterable of str): List of SMILES strings to convert

    Returns:
        tuple (mols, is_valid):
            mols (ndarray of Mol): Mol objects for valid SMILES strings only
            is_valid (ndarray of bool): True for each input SMILES string that was valid according to RDKit
            
    """
    log.debug('Converting SMILES to RDKit Mols')
    nsmiles = len(smiles_strs)
    mols = [None]*nsmiles
    for i, smi in enumerate(smiles_strs):
        try:
            mols[i] = Chem.MolFromSmiles(smi)
        except TypeError:
            pass

    log.debug('Adding hydrogens to mols')
    mols = [Chem.AddHs(m) if m is not None else None for m in mols]
    log.debug('Computing 3D coordinates')
    for i, m in enumerate(mols):
        if m is not None:
            try:
                AllChem.EmbedMolecule(m)
            except RuntimeError:
                # This sometimes fails in the RDKit code. Give up on this molecule.
                mols[i] = None
    is_valid = np.array([(m is not None) for m in mols], dtype=bool)
    mols = np.array(mols)[is_valid]
    return mols, is_valid 

def make_input_files(mol, charge = 0.0):

    new_mol = AllChem.AddHs(mol)
    new_mol = replace_invalid_atoms(new_mol)
    
    symbols = list(set([a.GetSymbol() for a in new_mol.GetAtoms()]))
    num_atoms = len(new_mol.GetAtoms())

    # Do we need to generate a geometry? Only if this is the first time
    if charge == 0.0:
        AllChem.EmbedMolecule(new_mol)
        conf = new_mol.GetConformer()
        with open(os.path.join(dftb_root, 'geom.gen'), 'w') as fid:
            fid.write('{} C\n'.format(num_atoms))
            fid.write('  {}\n'.format(' '.join(symbols)))
            for (i, a) in enumerate(new_mol.GetAtoms()):
                fid.write('%i\t%i\t%f\t%f\t%f\n' % (i + 1, symbols.index(a.GetSymbol()) + 1, 
                    conf.GetAtomPosition(a.GetIdx()).x, conf.GetAtomPosition(a.GetIdx()).y, conf.GetAtomPosition(a.GetIdx()).z))

    maxangularmomentum = ''
    for symbol in symbols:
        try:
            maxangularmomentum += '    {} = "{}"\n'.format(symbol, L_max[symbol])
        except KeyError:
            raise ValueError('Cannot run DFTB+ on this molecule because symbol {} does not have a defined MaxAngularMomentum! Could it be substituted for one of the elements for which we have parameters? Check {}'.format(symbol, __file__))

    with open(dftb_in, 'w') as fid:
        if charge == 0.0:
            fid.write(template_input_file % (charge, maxangularmomentum))
        else:
            # Use optimized neutral geometry
            fid.write(template_input_file.replace('geom.gen', 'geom.out.gen').replace('MaxSteps = 200', 'MaxSteps = 0') % (charge, maxangularmomentum))

    return new_mol 

def rdkit_optimize(self, addHs=True):
        if addHs:
            self.mol = Chem.AddHs(self.mol)
        AllChem.EmbedMolecule(self.mol, useExpTorsionAnglePrefs=True,useBasicKnowledge=True)
        AllChem.UFFOptimizeMolecule(self.mol) 

def load_data_from_smiles(x_smiles, labels, add_dummy_node=True, one_hot_formal_charge=False):
    """Load and featurize data from lists of SMILES strings and labels.

    Args:
        x_smiles (list[str]): A list of SMILES strings.
        labels (list[float]): A list of the corresponding labels.
        add_dummy_node (bool): If True, a dummy node will be added to the molecular graph. Defaults to True.
        one_hot_formal_charge (bool): If True, formal charges on atoms are one-hot encoded. Defaults to False.

    Returns:
        A tuple (X, y) in which X is a list of graph descriptors (node features, adjacency matrices, distance matrices),
        and y is a list of the corresponding labels.
    """
    x_all, y_all = [], []

    for smiles, label in zip(x_smiles, labels):
        try:
            mol = MolFromSmiles(smiles)
            try:
                mol = Chem.AddHs(mol)
                AllChem.EmbedMolecule(mol, maxAttempts=5000)
                AllChem.UFFOptimizeMolecule(mol)
                mol = Chem.RemoveHs(mol)
            except:
                AllChem.Compute2DCoords(mol)

            afm, adj, dist = featurize_mol(mol, add_dummy_node, one_hot_formal_charge)
            x_all.append([afm, adj, dist])
            y_all.append([label])
        except ValueError as e:
            logging.warning('the SMILES ({}) can not be converted to a graph.\nREASON: {}'.format(smiles, e))

    return x_all, y_all 

def optimize_conformer(idx, m, prop, algo="MMFF"):
    print("Calculating {}: {} ...".format(idx, Chem.MolToSmiles(m)))

    mol = Chem.AddHs(m)

    if algo == "ETKDG":
        # Landrum et al. DOI: 10.1021/acs.jcim.5b00654
        k = AllChem.EmbedMolecule(mol, AllChem.ETKDG())

        if k != 0:
            return None, None

    elif algo == "UFF":
        # Universal Force Field
        AllChem.EmbedMultipleConfs(mol, 50, pruneRmsThresh=0.5)
        try:
            arr = AllChem.UFFOptimizeMoleculeConfs(mol, maxIters=2000)
        except ValueError:
            return None, None

        if not arr:
            return None, None

        else:
            arr = AllChem.UFFOptimizeMoleculeConfs(mol, maxIters=20000)
            idx = np.argmin(arr, axis=0)[1]
            conf = mol.GetConformers()[idx]
            mol.RemoveAllConformers()
            mol.AddConformer(conf)

    elif algo == "MMFF":
        # Merck Molecular Force Field
        AllChem.EmbedMultipleConfs(mol, 50, pruneRmsThresh=0.5)
        try:
            arr = AllChem.MMFFOptimizeMoleculeConfs(mol, maxIters=2000)
        except ValueError:
            return None, None

        if not arr:
            return None, None

        else:
            arr = AllChem.MMFFOptimizeMoleculeConfs(mol, maxIters=20000)
            idx = np.argmin(arr, axis=0)[1]
            conf = mol.GetConformers()[idx]
            mol.RemoveAllConformers()
            mol.AddConformer(conf)

    mol = Chem.RemoveHs(mol)

    return mol, prop 

def construct_distance_matrix(mol, out_size=-1, contain_Hs=False):
    """Construct distance matrix

    Args:
        mol (Chem.Mol):
        out_size (int):
        contain_Hs (bool):

    Returns (numpy.ndarray): 2 dimensional array which represents distance
        between atoms

    """
    if mol is None:
        raise MolFeatureExtractionError('mol is None')

    N = mol.GetNumAtoms()
    if out_size < 0:
        size = N
    elif out_size >= N:
        size = out_size
    else:
        raise MolFeatureExtractionError('out_size {} is smaller than number '
                                        'of atoms in mol {}'
                                        .format(out_size, N))

    if contain_Hs:
        mol2 = mol
    else:
        mol2 = AllChem.AddHs(mol)

    conf_id = AllChem.EmbedMolecule(mol2)
    if not contain_Hs:
        mol2 = AllChem.RemoveHs(mol2)

    try:
        dist_matrix = rdmolops.Get3DDistanceMatrix(mol2, confId=conf_id)
    except ValueError as e:
        logger = getLogger(__name__)
        logger.info('construct_distance_matrix failed, type: {}, {}'
                    .format(type(e).__name__, e.args))
        logger.debug(traceback.format_exc())
        raise MolFeatureExtractionError

    if size > N:
        dists = numpy.zeros((size, size), dtype=numpy.float32)
        a0, a1 = dist_matrix.shape
        dists[:a0, :a1] = dist_matrix
    else:
        dists = dist_matrix
    return dists.astype(numpy.float32) 

def _to_xyz_rdkit(self, optimizer, **kwargs):
        """
        The internal function creates and stores the xyz coordinates for a pre-built molecule object.
        """
        # add hydrogens >> commented out and left for the users to take care of it using hydrogens method.
        # self.hydrogens('add')

        # embeding and optimization
        if optimizer:
            AllChem.EmbedMolecule(self.rdkit_molecule)
            if optimizer ==  'MMFF':
                if AllChem.MMFFHasAllMoleculeParams(self.rdkit_molecule):
                    AllChem.MMFFOptimizeMolecule(self.rdkit_molecule, **kwargs)
                else:
                    msg = "The MMFF parameters are not available for all of the molecule's atoms."
                    raise ValueError(msg)
            elif optimizer == 'UFF':
                if AllChem.UFFHasAllMoleculeParams(self.rdkit_molecule):
                    AllChem.UFFOptimizeMolecule(self.rdkit_molecule, **kwargs)
                else:
                    msg = "The UFF parameters are not available for all of the molecule's atoms."
                    raise ValueError(msg)
            else:
                msg = "The '%s' is not a legit value for the optimizer parameter."%str(optimizer)
                raise ValueError(msg)

        # get conformer
        try:
            conf = self.rdkit_molecule.GetConformer()
        except ValueError:
            msg = "The conformation has not been built yet (maybe due to the 2D representation of the creator).\n" \
                  "You should set the optimizer value if you wish to embed and optimize the 3D geometry."
            raise ValueError(msg)
        geometry = conf.GetPositions()
        atoms_list = self.rdkit_molecule.GetAtoms()
        atomic_nums = np.array([i.GetAtomicNum() for i in atoms_list])
        atomic_symbols = np.array([i.GetSymbol() for i in atoms_list])
        self._xyz = XYZ(geometry, atomic_nums.reshape(-1,1), atomic_symbols.reshape(-1,1))

        if optimizer=='UFF':
            self._UFF_args = update_default_kwargs(self._default_UFF_args, kwargs,
                                                 self._to_xyz_core_names[1], self._to_xyz_core_docs[1])
        elif optimizer=='MMFF':
            self._MMFF_args = update_default_kwargs(self._default_MMFF_args, kwargs,
                                                 self._to_xyz_core_names[0], self._to_xyz_core_docs[0])