Python rdkit.Chem.GetFormalCharge() Examples

def write_xtb_input_file(fragment, fragment_name):
    number_of_atoms = fragment.GetNumAtoms()
    charge = Chem.GetFormalCharge(fragment)
    symbols = [a.GetSymbol() for a in fragment.GetAtoms()] 
    for i,conf in enumerate(fragment.GetConformers()):
        file_name = fragment_name+"+"+str(i)+".xyz"
        with open(file_name, "w") as file:
            file.write(str(number_of_atoms)+"\n")
            file.write("title\n")
            for atom,symbol in enumerate(symbols):
                p = conf.GetAtomPosition(atom)
                line = " ".join((symbol,str(p.x),str(p.y),str(p.z),"\n"))
                file.write(line)
            if charge !=0:
                file.write("$set\n")
                file.write("chrg "+str(charge)+"\n")
                file.write("$end") 

def charge(self):
        r'''Charge of a chemical, computed with RDKit from a chemical's SMILES.
        If RDKit is not available, holds None.

        Examples
        --------
        >>> Chemical('sodium ion').charge
        1
        '''
        try:
            if not self.rdkitmol:
                return charge_from_formula(self.formula)
            else:
                return Chem.GetFormalCharge(self.rdkitmol)
        except:
            return charge_from_formula(self.formula) 

def run(self, mol):
        charge = Chem.GetFormalCharge(mol)
        if not charge == 0:
            chargestring = '+%s' % charge if charge > 0 else '%s' % charge
            self.log.info('Not an overall neutral system (%s)', chargestring) 

def mol2xyz(mol, multiplicity=1):
    charge = Chem.GetFormalCharge(mol)
    xyz_string = "\n{} {}\n".format(charge, multiplicity)
    for atom in mol.GetAtoms():
        pos = mol.GetConformer().GetAtomPosition(atom.GetIdx())
        xyz_string += "{} {} {} {}\n".format(atom.GetSymbol(), pos.x, pos.y, pos.z)
    return xyz_string 

def get_mol(smiles):
    mol = Chem.MolFromSmiles(smiles)
    Chem.Kekulize(mol, clearAromaticFlags=True)
    charge = Chem.GetFormalCharge(mol)
    mol = Chem.AddHs(mol)
    return mol 

def test_smiles_from_adjacent_matrix(smiles):

    charged_fragments = True
    quick = True

    # Cut apart the smiles
    mol = get_mol(smiles)
    atoms = get_atoms(mol)
    charge = Chem.GetFormalCharge(mol)
    adjacent_matrix = Chem.GetAdjacencyMatrix(mol)

    #
    mol = Chem.RemoveHs(mol)
    canonical_smiles = Chem.MolToSmiles(mol)

    # Define new molecule template from atoms
    new_mol = x2m.get_proto_mol(atoms)

    # reconstruct the molecule from adjacent matrix, atoms and total charge
    new_mol = x2m.AC2mol(new_mol, adjacent_matrix, atoms, charge, charged_fragments, quick)
    new_mol = Chem.RemoveHs(new_mol)
    new_mol_smiles = Chem.MolToSmiles(new_mol)

    assert new_mol_smiles == canonical_smiles

    return 

def test_smiles_from_coord_vdw(smiles):

    # The answer
    mol = Chem.MolFromSmiles(smiles)
    charge = Chem.GetFormalCharge(mol)
    canonical_smiles = Chem.MolToSmiles(mol, isomericSmiles=False)

    # generate forcefield coordinates
    atoms, coordinates = generate_structure_from_smiles(smiles)

    # Generate molobj from atoms, charge and coordinates
    mol = x2m.xyz2mol(atoms, coordinates, charge=charge)

    # For this test, remove chira. clean and canonical
    Chem.Kekulize(mol)
    mol = Chem.RemoveHs(mol)
    Chem.RemoveStereochemistry(mol)
    smiles = Chem.MolToSmiles(mol, isomericSmiles=False)

    # Please look away. A small hack that removes the explicit hydrogens
    mol = Chem.MolFromSmiles(smiles)
    smiles = Chem.MolToSmiles(mol)

    assert smiles == canonical_smiles

    return 

def test_smiles_from_coord_huckel(smiles):

    # The answer
    mol = Chem.MolFromSmiles(smiles)
    charge = Chem.GetFormalCharge(mol)
    canonical_smiles = Chem.MolToSmiles(mol, isomericSmiles=False)

    # generate forcefield coordinates
    atoms, coordinates = generate_structure_from_smiles(smiles)

    # Generate molobj from atoms, charge and coordinates
    mol = x2m.xyz2mol(atoms, coordinates, charge=charge, use_huckel=True)

    # For this test, remove chira. clean and canonical
    Chem.Kekulize(mol)
    mol = Chem.RemoveHs(mol)
    Chem.RemoveStereochemistry(mol)
    smiles = Chem.MolToSmiles(mol, isomericSmiles=False)

    # Please look away. A small hack that removes the explicit hydrogens
    mol = Chem.MolFromSmiles(smiles)
    smiles = Chem.MolToSmiles(mol)

    assert smiles == canonical_smiles

    return 

def reionize(self, mol):
        """If molecule with multiple acid groups is partially ionized, ensure strongest acids ionize first.

        The algorithm works as follows:

        - Use SMARTS to find the strongest protonated acid and the weakest ionized acid.
        - If the ionized acid is weaker than the protonated acid, swap proton and repeat.

        :param mol: The molecule to reionize.
        :type mol: :rdkit:`Mol <Chem.rdchem.Mol-class.html>`
        :return: The reionized molecule.
        :rtype: :rdkit:`Mol <Chem.rdchem.Mol-class.html>`
        """
        log.debug("Running Reionizer")
        while True:
            ppos, poccur = self._strongest_protonated(mol)
            ipos, ioccur = self._weakest_ionized(mol)
            if ioccur and poccur and ppos < ipos:
                log.info(
                    "Moved proton from %s to %s",
                    self.acid_base_pairs[ppos].name,
                    self.acid_base_pairs[ipos].name,
                )
                patom = mol.GetAtomWithIdx(poccur[-1])
                patom.SetFormalCharge(patom.GetFormalCharge() - 1)
                patom.SetNumExplicitHs(max(0, patom.GetNumExplicitHs() - 1))
                iatom = mol.GetAtomWithIdx(ioccur[-1])
                iatom.SetFormalCharge(iatom.GetFormalCharge() + 1)
                iatom.SetNumExplicitHs(iatom.GetTotalNumHs() + 1)
            else:
                Chem.SanitizeMol(mol)
                return mol 

def run(self, mol):
        charge = Chem.GetFormalCharge(mol)
        if not charge == 0:
            chargestring = "+%s" % charge if charge > 0 else "%s" % charge
            self.log.info("Not an overall neutral system (%s)", chargestring) 

def uncharge(self, mol):
        """Neutralize molecule by adding/removing hydrogens.

        :param mol: The molecule to uncharge.
        :type mol: rdkit.Chem.rdchem.Mol
        :return: The uncharged molecule.
        :rtype: rdkit.Chem.rdchem.Mol
        """
        log.debug('Running Uncharger')
        mol = copy.deepcopy(mol)

        # Neutralize positive charges
        pos_remainder = 0
        neg_count = 0
        for atom in mol.GetAtoms():
            # Remove hydrogen from positive atoms and reduce formal change until neutral or no more hydrogens
            while atom.GetFormalCharge() > 0 and atom.GetNumExplicitHs() > 0:
                atom.SetNumExplicitHs(atom.GetNumExplicitHs() - 1)
                atom.SetFormalCharge(atom.GetFormalCharge() - 1)
                log.info('Removed positive charge')
            chg = atom.GetFormalCharge()
            if chg > 0:
                # Record number of non-neutralizable positive charges
                pos_remainder += chg
            elif chg < 0:
                # Record total number of negative charges
                neg_count += -chg

        # Choose negative charges to leave in order to balance non-neutralizable positive charges
        neg_skip = self._get_neg_skip(mol, pos_remainder)

        # Neutralize remaining negative charges
        for atom in mol.GetAtoms():
            log.info(atom.GetIdx())
            if atom.GetIdx() in neg_skip:
                continue
            # Make sure to stop when neg_count <= pos_remainder, as it is possible that neg_skip is not large enough
            while atom.GetFormalCharge() < 0 and neg_count > pos_remainder:
                atom.SetNumExplicitHs(atom.GetNumExplicitHs() + 1)
                atom.SetFormalCharge(atom.GetFormalCharge() + 1)
                neg_count -= 1
                log.info('Removed negative charge')
        return mol 

def disconnect(self, mol):
        """Break covalent bonds between metals and organic atoms under certain conditions.

        The algorithm works as follows:

        - Disconnect N, O, F from any metal.
        - Disconnect other non-metals from transition metals + Al (but not Hg, Ga, Ge, In, Sn, As, Tl, Pb, Bi, Po).
        - For every bond broken, adjust the charges of the begin and end atoms accordingly.

        :param mol: The input molecule.
        :type mol: :rdkit:`Mol <Chem.rdchem.Mol-class.html>`
        :return: The molecule with metals disconnected.
        :rtype: :rdkit:`Mol <Chem.rdchem.Mol-class.html>`
        """
        log.debug("Running MetalDisconnector")
        # TODO: In future, maybe insert zero-order complex/ionic/dative bonds instead of disconnecting?
        # Remove bonds that match SMARTS
        for smarts in [self._metal_nof, self._metal_non]:
            pairs = mol.GetSubstructMatches(smarts)
            edmol = Chem.EditableMol(mol)
            orders = []
            for i, j in pairs:
                # TODO: Could get the valence contributions of the bond instead of GetBondTypeAsDouble?
                orders.append(int(mol.GetBondBetweenAtoms(i, j).GetBondTypeAsDouble()))
                edmol.RemoveBond(i, j)
            # Adjust neighbouring charges accordingly
            mol = edmol.GetMol()
            for n, (i, j) in enumerate(pairs):
                chg = orders[n]
                atom1 = mol.GetAtomWithIdx(i)
                atom1.SetFormalCharge(atom1.GetFormalCharge() + chg)
                atom2 = mol.GetAtomWithIdx(j)
                atom2.SetFormalCharge(atom2.GetFormalCharge() - chg)
                log.info(
                    "Removed covalent bond between %s and %s",
                    atom1.GetSymbol(),
                    atom2.GetSymbol(),
                )
        # Ionize a free neutral metal with a protonated carboxylic acid
        # TODO: Extend this to other acids?
        # TODO: Move to charge module?
        fms = mol.GetSubstructMatches(self._free_metal)
        carbs = mol.GetSubstructMatches(self._carboxylic)
        if len(fms) == len(carbs) > 0:
            for fm in fms:
                atom = mol.GetAtomWithIdx(fm[0])
                atom.SetFormalCharge(atom.GetFormalCharge() + 1)
            for carb in carbs:
                atom = mol.GetAtomWithIdx(carb[2])
                atom.SetFormalCharge(atom.GetFormalCharge() - 1)
        log.debug(Chem.MolToSmiles(mol))
        Chem.SanitizeMol(mol)
        return mol


# ==============================================================================
# from .fragment import PREFER_ORGANIC, LargestFragmentChooser, FragmentRemover
# ============================================================================== 

def uncharge(self, mol):
        """Neutralize molecule by adding/removing hydrogens. Attempts to preserve zwitterions.

        :param mol: The molecule to uncharge.
        :type mol: :rdkit:`Mol <Chem.rdchem.Mol-class.html>`
        :return: The uncharged molecule.
        :rtype: :rdkit:`Mol <Chem.rdchem.Mol-class.html>`
        """
        log.debug("Running Uncharger")
        mol = copy.deepcopy(mol)
        # Get atom ids for matches
        p = [x[0] for x in mol.GetSubstructMatches(self._pos_h)]
        q = [x[0] for x in mol.GetSubstructMatches(self._pos_quat)]
        n = [x[0] for x in mol.GetSubstructMatches(self._neg)]
        a = [x[0] for x in mol.GetSubstructMatches(self._neg_acid)]
        # Neutralize negative charges
        if q:
            # Surplus negative charges more than non-neutralizable positive charges
            neg_surplus = len(n) - len(q)
            if a and neg_surplus > 0:
                # zwitterion with more negative charges than quaternary positive centres
                while neg_surplus > 0 and a:
                    # Add hydrogen to first negative acid atom, increase formal charge
                    # Until quaternary positive == negative total or no more negative acid
                    atom = mol.GetAtomWithIdx(a.pop(0))
                    atom.SetNumExplicitHs(atom.GetNumExplicitHs() + 1)
                    atom.SetFormalCharge(atom.GetFormalCharge() + 1)
                    neg_surplus -= 1
                    log.info("Removed negative charge")
        else:
            #
            for atom in [mol.GetAtomWithIdx(x) for x in n]:
                while atom.GetFormalCharge() < 0:
                    atom.SetNumExplicitHs(atom.GetNumExplicitHs() + 1)
                    atom.SetFormalCharge(atom.GetFormalCharge() + 1)
                    log.info("Removed negative charge")
        # Neutralize positive charges
        for atom in [mol.GetAtomWithIdx(x) for x in p]:
            # Remove hydrogen and reduce formal change until neutral or no more hydrogens
            while atom.GetFormalCharge() > 0 and atom.GetNumExplicitHs() > 0:
                atom.SetNumExplicitHs(atom.GetNumExplicitHs() - 1)
                atom.SetFormalCharge(atom.GetFormalCharge() - 1)
                log.info("Removed positive charge")
        return mol