Python rdkit.Chem.AddHs() 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 CalculateBasakSIC1(mol):
    """
    #################################################################
    Obtain the structural information content with order 1

    proposed by Basak.

    ---->SIC1
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC1(mol)
    if nAtoms <= 1:
        BasakSIC = 0.0
    else:
        BasakSIC = IC / numpy.log2(nAtoms)

    return BasakSIC 

def prepare_smiles_and_mol(self, mol):
        """Prepare `smiles` and `mol` used in following preprocessing.

        This method is called before `get_input_features` is called, by parser
        class.
        This method may be overriden to support custom `smile`/`mol` extraction

        Args:
            mol (mol): mol instance

        Returns (tuple): (`smiles`, `mol`)
        """
        # Note that smiles expression is not unique.
        # we obtain canonical smiles which is unique in `mol`
        canonical_smiles = Chem.MolToSmiles(mol, isomericSmiles=False,
                                            canonical=True)
        mol = Chem.MolFromSmiles(canonical_smiles)
        if self.add_Hs:
            mol = Chem.AddHs(mol)
        if self.kekulize:
            Chem.Kekulize(mol)
        return canonical_smiles, mol 

def CalculateHydrogenNumber(mol):

    """
    #################################################################
    Calculation of Number of Hydrogen in a molecule

    ---->nhyd

    Usage:

        result=CalculateHydrogenNumber(mol)

        Input: mol is a molecule object.

        Output: result is a numeric value.
    #################################################################
    """
    i = 0
    Hmol = Chem.AddHs(mol)
    for atom in Hmol.GetAtoms():
        if atom.GetAtomicNum() == 1:
            i = i + 1

    return i 

def CalculateAllAtomNumber(mol):
    """
    #################################################################
    Calculation of all atom counts in a molecule

    ---->nta

    Usage:

        result=CalculateAllAtomNumber(mol)

        Input: mol is a molecule object.

        Output: result is a numeric value.
    #################################################################
    """

    return Chem.AddHs(mol).GetNumAtoms() 

def CalculateBasakCIC6(mol):
    """
    #################################################################
    Obtain the complementary information content with order 6 proposed

    by Basak.

    ---->CIC6
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC6(mol)
    if nAtoms <= 1:
        BasakCIC = 0.0
    else:
        BasakCIC = numpy.log2(nAtoms) - IC

    return BasakCIC 

def CalculateBasakSIC0(mol):
    """
    #################################################################
    Obtain the structural information content with order 0

    proposed by Basak

    ---->SIC0
    #################################################################
    """

    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC0(mol)
    if nAtoms <= 1:
        BasakSIC = 0.0
    else:
        BasakSIC = IC / numpy.log2(nAtoms)

    return BasakSIC 

def CalculateBasakCIC0(mol):

    """
    #################################################################
    Obtain the complementary information content with order 0

    proposed by Basak

    ---->CIC0
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC0(mol)
    if nAtoms <= 1:
        BasakCIC = 0.0
    else:
        BasakCIC = numpy.log2(nAtoms) - IC

    return BasakCIC 

def CalculateBasakCIC5(mol):
    """
    #################################################################
    Obtain the complementary information content with order 5 proposed

    by Basak.

    ---->CIC5
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC5(mol)
    if nAtoms <= 1:
        BasakCIC = 0.0
    else:
        BasakCIC = numpy.log2(nAtoms) - IC

    return BasakCIC 

def featurize(self, x):
        # check if type(x) = list
        if isinstance(x, pd.Series):
            x = x.tolist()
        if not isinstance(x, list):
            x = [x]
        # check input format, assume SMILES if not RDKit-MOL
        if not isinstance(x[0], Chem.rdchem.Mol):
            x_mol = []
            for z in x:
                x_mol.append(Chem.MolFromSmiles(z))
                if x_mol[-1] is None:
                    raise ValueError('can not convert Mol from SMILES %s' % z)
        else:
            x_mol = x
        
        # convert to counting dictionary
        mol = [Chem.AddHs(z) for z in x_mol]
        d_list = [dict(Counter([atom.GetSymbol() for atom in z.GetAtoms()])) for z in mol]

        self.output = self._cal.transform(d_list)
        
        return self.output 

def CalculateBasakSIC3(mol):
    """
    #################################################################
    Obtain the structural information content with order 3 proposed

    by Basak.

    ---->SIC3
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC3(mol)
    if nAtoms <= 1:
        BasakSIC = 0.0
    else:
        BasakSIC = IC / numpy.log2(nAtoms)

    return BasakSIC 

def CalculateBasakSIC4(mol):
    """
    #################################################################
    Obtain the structural information content with order 4 proposed

    by Basak.

    ---->SIC4
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC4(mol)
    if nAtoms <= 1:
        BasakSIC = 0.0
    else:
        BasakSIC = IC / numpy.log2(nAtoms)

    return BasakSIC 

def CalculateBasakSIC5(mol):
    """
    #################################################################
    Obtain the structural information content with order 5 proposed

    by Basak.

    ---->SIC5
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC5(mol)
    if nAtoms <= 1:
        BasakSIC = 0.0
    else:
        BasakSIC = IC / numpy.log2(nAtoms)

    return BasakSIC 

def CalculateBasakSIC6(mol):
    """
    #################################################################
    Obtain the structural information content with order 6 proposed

    by Basak.

    ---->SIC6
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC6(mol)
    if nAtoms <= 1:
        BasakSIC = 0.0
    else:
        BasakSIC = IC / numpy.log2(nAtoms)

    return BasakSIC 

def CalculateBasakCIC1(mol):
    """
    #################################################################
    Obtain the complementary information content with order 1 proposed

    by Basak.

    ---->CIC1
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC1(mol)
    if nAtoms <= 1:
        BasakCIC = 0.0
    else:
        BasakCIC = numpy.log2(nAtoms) - IC

    return BasakCIC 

def CalculateBasakCIC3(mol):
    """
    #################################################################
    Obtain the complementary information content with order 3 proposed

    by Basak.

    ---->CIC3
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC3(mol)
    if nAtoms <= 1:
        BasakCIC = 0.0
    else:
        BasakCIC = numpy.log2(nAtoms) - IC

    return BasakCIC 

def CalculateBasakCIC4(mol):
    """
    #################################################################
    Obtain the complementary information content with order 4 proposed

    by Basak.

    ---->CIC4
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    nAtoms = Hmol.GetNumAtoms()
    IC = CalculateBasakIC4(mol)
    if nAtoms <= 1:
        BasakCIC = 0.0
    else:
        BasakCIC = numpy.log2(nAtoms) - IC

    return BasakCIC 

def _extra_docs(self):
        # method: to_smiles
        self._to_smiles_core_names = ("rdkit.Chem.MolToSmarts",)
        self._to_smiles_core_docs = ("http://rdkit.org/docs/source/rdkit.Chem.inchi.html?highlight=inchi#rdkit.Chem.inchi.MolToInchi",)
        # method: to_smarts
        self._to_smarts_core_names = ("rdkit.Chem.MolToSmarts",)
        self._to_smarts_core_docs = ("http://rdkit.org/docs/source/rdkit.Chem.inchi.html?highlight=inchi#rdkit.Chem.inchi.MolToInchi",)
        # method: to_inchi
        self._to_inchi_core_names = ("rdkit.Chem.MolToInchi",)
        self._to_inchi_core_docs = ("http://rdkit.org/docs/source/rdkit.Chem.inchi.html?highlight=inchi#rdkit.Chem.inchi.MolToInchi",)
        # method: hydrogens
        self._hydrogens_core_names = ("rdkit.Chem.AddHs","rdkit.Chem.RemoveHs")
        self._hydrogens_core_docs = ("http://rdkit.org/docs/source/rdkit.Chem.rdmolops.html?highlight=addhs#rdkit.Chem.rdmolops.AddHs",
                                    "http://rdkit.org/docs/source/rdkit.Chem.rdmolops.html?highlight=addhs#rdkit.Chem.rdmolops.RemoveHs")
        #
        self._to_xyz_core_names = ("rdkit.Chem.AllChem.MMFFOptimizeMolecule","rdkit.Chem.AllChem.UFFOptimizeMolecule")
        self._to_xyz_core_docs =(
            "http://rdkit.org/docs/source/rdkit.Chem.rdForceFieldHelpers.html?highlight=mmff#rdkit.Chem.rdForceFieldHelpers.MMFFOptimizeMolecule",
            "http://rdkit.org/docs/source/rdkit.Chem.rdForceFieldHelpers.html?highlight=mmff#rdkit.Chem.rdForceFieldHelpers.UFFOptimizeMolecule"
        ) 

def _CalculateElementSumSquareCharge(mol, AtomicNum=6):
    """
    #################################################################
    **Internal used only**

    Ths sum of square Charges on all atoms with atomicnumber equal to n
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    GMCharge.ComputeGasteigerCharges(Hmol, iter_step)
    res = []
    for atom in Hmol.GetAtoms():
        if atom.GetAtomicNum() == AtomicNum:
            res.append(float(atom.GetProp("_GasteigerCharge")))
    if res == []:
        return 0
    else:
        return round(sum(numpy.square(res)), 3) 

def CalculateAllMaxNCharge(mol):
    """
    #################################################################
    Most negative charge on all atoms

    -->Qmin

    Usage:

        result=CalculateAllMaxNCharge(mol)

        Input: mol is a molecule object.

        Output: result is a numeric value.
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    GMCharge.ComputeGasteigerCharges(Hmol, iter_step)
    res = []
    for atom in Hmol.GetAtoms():
        res.append(float(atom.GetProp("_GasteigerCharge")))
    if res == []:
        return 0
    else:
        return round(min(res), 3) 

def CalculateAllMaxPCharge(mol):
    """
    #################################################################
    Most positive charge on ALL atoms

    -->Qmax

    Usage:

        result=CalculateAllMaxPCharge(mol)

        Input: mol is a molecule object.

        Output: result is a numeric value.
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    GMCharge.ComputeGasteigerCharges(Hmol, iter_step)
    res = []
    for atom in Hmol.GetAtoms():
        res.append(float(atom.GetProp("_GasteigerCharge")))
    if res == []:
        return 0
    else:
        return round(max(res), 3) 

def _CalculateElementMaxNCharge(mol, AtomicNum=6):
    """
    #################################################################
    **Internal used only**

    Most negative charge on atom with atomic number equal to n
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    GMCharge.ComputeGasteigerCharges(Hmol, iter_step)
    res = []
    for atom in Hmol.GetAtoms():
        if atom.GetAtomicNum() == AtomicNum:
            res.append(float(atom.GetProp("_GasteigerCharge")))
    if res == []:
        return 0
    else:
        return round(min(res), 3) 

def _CalculateElementMaxPCharge(mol, AtomicNum=6):
    """
    #################################################################
    **Internal used only**

    Most positive charge on atom with atomic number equal to n
    #################################################################
    """
    Hmol = Chem.AddHs(mol)
    GMCharge.ComputeGasteigerCharges(Hmol, iter_step)
    res = []
    for atom in Hmol.GetAtoms():
        if atom.GetAtomicNum() == AtomicNum:
            res.append(float(atom.GetProp("_GasteigerCharge")))

    if res == []:
        return 0
    else:
        return round(max(res), 3) 

def __init__(self, mol, atom_count=None, MW=None, Tb=None):
        if type(mol) == Chem.rdchem.Mol:
            self.rdkitmol = mol
        else:
            self.rdkitmol = Chem.MolFromSmiles(mol)
        if atom_count is None:
            self.rdkitmol_Hs = Chem.AddHs(self.rdkitmol)
            self.atom_count = len(self.rdkitmol_Hs.GetAtoms())
        else:
            self.atom_count = atom_count
        if MW is None:
            self.MW = rdMolDescriptors.CalcExactMolWt(self.rdkitmol_Hs)
        else:
            self.MW = MW
            
        self.counts, self.success, self.status = smarts_fragment(J_BIGGS_JOBACK_SMARTS_id_dict, rdkitmol=self.rdkitmol)
            
        if Tb is not None:
            self.Tb_estimated = self.Tb(self.counts)
        else:
            self.Tb_estimated = Tb 

def setUp(self):
    """
        Set up tests.
        """
    smiles = '[H]C([H])([H])[H]'
    from rdkit import Chem
    mol = Chem.MolFromSmiles(smiles)
    mol = Chem.AddHs(mol)
    engine = conformers.ConformerGenerator(max_conformers=1)
    self.mol = engine.generate_conformers(mol)
    assert self.mol.GetNumConformers() > 0 

def water_octanol_partition_coefficient_scores(mols, norm=False):
        scores = [MolecularMetrics._avoid_sanitization_error(lambda: Crippen.MolLogP(Chem.AddHs(mol, 1))) if mol is not None else None
                  for mol in mols]
        scores = np.array(list(map(lambda x: -3 if x is None else x, scores)))
        scores = np.clip(MolecularMetrics.remap(scores, -2.12178879609, 6.0429063424), 0.0, 1.0) if norm else scores

        return scores 

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 hydrogens(self, action='add', **kwargs):
        """
        This function adds/removes hydrogens to/from a prebuilt molecule object.

        Parameters
        ----------
        action: str
            Either 'add' or 'remove', to add hydrogns or remove them from the rdkit molecule.

        kwargs:
            The arguments that can be passed to the rdkit functions:
            - `Chem.AddHs`: documentation at http://rdkit.org/docs/source/rdkit.Chem.rdmolops.html?highlight=addhs#rdkit.Chem.rdmolops.AddHs
            - `Chem.RemoveHs`: documentation at http://rdkit.org/docs/source/rdkit.Chem.rdmolops.html?highlight=addhs#rdkit.Chem.rdmolops.RemoveHs

        Notes
        -----
            - The rdkit or pybel molecule object must be created in advance.
            - Only rdkit or pybel molecule object will be modified in place.
            - If you remove hydrogens from molecules, the atomic 3D coordinates might not be accurate for the conversion to xyz representation.

        """

        # molecule must exist
        _ = self._check_original_molecule()

        if action == 'add':
            if self.pybel_molecule:
                self.pybel_molecule.addh()
            # Note: just if not elif
            if self.rdkit_molecule:
                self.rdkit_molecule = Chem.AddHs(self.rdkit_molecule, **kwargs)
        elif action == 'remove':
            if self.pybel_molecule:
                self.pybel_molecule.removeh()
            if self.rdkit_molecule:
                self.rdkit_molecule = Chem.RemoveHs(self.rdkit_molecule, **kwargs)
        else:
            raise ValueError("The parameter 'action' must be either of 'add' or 'remove'.") 

def embed_molecule(self, mol):
        """Generate conformers, possibly with pruning.

        Parameters
        ----------
        mol : RDKit Mol
            Molecule.
        """
        logging.debug("Adding hydrogens for %s" % mol.GetProp("_Name"))
        mol = Chem.AddHs(mol)  # add hydrogens
        logging.debug("Hydrogens added to %s" % mol.GetProp("_Name"))
        logging.debug("Sanitizing mol for %s" % mol.GetProp("_Name"))
        Chem.SanitizeMol(mol)
        logging.debug("Mol sanitized for %s" % mol.GetProp("_Name"))
        if self.max_conformers == -1 or type(self.max_conformers) is not int:
            self.max_conformers = self.get_num_conformers(mol)
        n_confs = self.max_conformers * self.pool_multiplier
        if self.first_conformers == -1:
            self.first_conformers = self.max_conformers
        logging.debug(
            "Embedding %d conformers for %s" % (n_confs, mol.GetProp("_Name"))
        )
        AllChem.EmbedMultipleConfs(
            mol,
            numConfs=n_confs,
            maxAttempts=10 * n_confs,
            pruneRmsThresh=-1.0,
            randomSeed=self.seed,
            ignoreSmoothingFailures=True,
        )
        logging.debug("Conformers embedded for %s" % mol.GetProp("_Name"))
        return mol 

def find_largest_number_of_atoms_atomic_number_and_ms_peak(
    mol_list, add_hs_to_molecule=False):
  """Finds the greatest number of molecules and the largest peak in the spectra.

  Used to determine the largest number of atoms in a molecule, and the
  largest mass/charge ratio found in the spectra.

  Args:
    mol_list : list of rdkit.Mol objects
    add_hs_to_molecule : whether or not to add hydrogens to the molecule.

  Returns:
    max_atoms : Largest number of atoms found in a molecule in mol_list
    max_atom_type: Largest atomic number in all the molecules from the mol_list
    max_peak_loc : Greatest mass/charge peak found in samples in sdf_file

  Raises:
    ValueError: If first rdkit.Mol in mol_list does not have
        MASS SPECTRAL PEAKS as one of the properties.
        This should already be present as one of the tags in the original
        sdf file the rdkit.Mols were parsed from.
  """
  if not mol_list[0].HasProp(ms_constants.SDF_TAG_MASS_SPEC_PEAKS):
    raise ValueError('first molecule in list does not contain SDF tag'
                     '\'{}\''.format(ms_constants.SDF_TAG_MASS_SPEC_PEAKS))

  # Add hydrogens to atoms:
  if add_hs_to_molecule:
    mol_list = [Chem.rdmolops.AddHs(mol) for mol in mol_list]

  max_atoms = max(mol.GetNumAtoms() for mol in mol_list)
  max_atom_type = max(
      [at.GetAtomicNum() for mol in mol_list for at in mol.GetAtoms()])

  max_peak_loc = max(
      feature_utils.get_largest_mass_spec_peak_loc(mol) for mol in mol_list)

  return max_atoms, max_atom_type, max_peak_loc