molvs/charge.py
def reionize(self, mol):
        """Enforce charges on certain atoms, then perform competitive reionization.

        First, charge corrections are applied to ensure, for example, that free metals are correctly ionized. Then, if
        a molecule with multiple acid groups is partially ionized, ensure the 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.Chem.rdchem.Mol
        :return: The reionized molecule.
        :rtype: rdkit.Chem.rdchem.Mol
        """
        log.debug('Running Reionizer')

        start_charge = Chem.GetFormalCharge(mol)

        # Apply forced charge corrections
        for cc in self.charge_corrections:
            for match in mol.GetSubstructMatches(cc.smarts):
                atom = mol.GetAtomWithIdx(match[0])
                log.info('Applying charge correction %s (%s %+d)', cc.name, atom.GetSymbol(), cc.charge)
                atom.SetFormalCharge(cc.charge)

        current_charge = Chem.GetFormalCharge(mol)
        charge_diff = Chem.GetFormalCharge(mol) - start_charge
        # If molecule is now neutral, assume everything is now fixed
        # But otherwise, if charge has become more positive, look for additional protonated acid groups to ionize
        if not current_charge == 0:
            while charge_diff > 0:
                ppos, poccur = self._strongest_protonated(mol)
                if ppos is None:
                    break
                log.info('Ionizing %s to balance previous charge corrections', self.acid_base_pairs[ppos].name)
                patom = mol.GetAtomWithIdx(poccur[-1])
                patom.SetFormalCharge(patom.GetFormalCharge() - 1)
                if patom.GetNumExplicitHs() > 0:
                    patom.SetNumExplicitHs(patom.GetNumExplicitHs() - 1)
                # else:
                patom.UpdatePropertyCache()
                charge_diff -= 1

        already_moved = set()
        while True:
            ppos, poccur = self._strongest_protonated(mol)
            ipos, ioccur = self._weakest_ionized(mol)
            if ioccur and poccur and ppos < ipos:
                if poccur[-1] == ioccur[-1]:
                    # Bad! H wouldn't be moved, resulting in infinite loop.
                    log.warning('Aborted reionization due to unexpected situation')
                    break

                key = tuple(sorted([poccur[-1], ioccur[-1]]))
                if key in already_moved:
                    log.warning('Aborting reionization to avoid infinite loop due to it being ambiguous where to put a Hydrogen')
                    break
                already_moved.add(key)

                log.info('Moved proton from %s to %s', self.acid_base_pairs[ppos].name, self.acid_base_pairs[ipos].name)

                # Remove hydrogen from strongest protonated
                patom = mol.GetAtomWithIdx(poccur[-1])
                patom.SetFormalCharge(patom.GetFormalCharge() - 1)
                # If no implicit Hs to autoremove, and at least 1 explicit H to remove, reduce explicit count by 1
                if patom.GetNumImplicitHs() == 0 and patom.GetNumExplicitHs() > 0:
                    patom.SetNumExplicitHs(patom.GetNumExplicitHs() - 1)
                    # TODO: Remove any chiral label on patom?
                patom.UpdatePropertyCache()

                # Add hydrogen to weakest ionized
                iatom = mol.GetAtomWithIdx(ioccur[-1])
                iatom.SetFormalCharge(iatom.GetFormalCharge() + 1)
                # Increase explicit H count if no implicit, or aromatic N or P, or non default valence state
                if (iatom.GetNoImplicit() or
                        ((patom.GetAtomicNum() == 7 or patom.GetAtomicNum() == 15) and patom.GetIsAromatic()) or
                        iatom.GetTotalValence() not in list(Chem.GetPeriodicTable().GetValenceList(iatom.GetAtomicNum()))):
                    iatom.SetNumExplicitHs(iatom.GetNumExplicitHs() + 1)
                iatom.UpdatePropertyCache()
            else:
                break

        # TODO: Canonical ionization position if multiple equivalent positions?

        Chem.SanitizeMol(mol)
        return mol
def reionize(self, mol):
        """Enforce charges on certain atoms, then perform competitive reionization.

        First, charge corrections are applied to ensure, for example, that free metals are correctly ionized. Then, if
        a molecule with multiple acid groups is partially ionized, ensure the 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.Chem.rdchem.Mol
        :return: The reionized molecule.
        :rtype: rdkit.Chem.rdchem.Mol
        """
        log.debug('Running Reionizer')

        start_charge = Chem.GetFormalCharge(mol)

        # Apply forced charge corrections
        for cc in self.charge_corrections:
            for match in mol.GetSubstructMatches(cc.smarts):
                atom = mol.GetAtomWithIdx(match[0])
                log.info('Applying charge correction %s (%s %+d)', cc.name, atom.GetSymbol(), cc.charge)
                atom.SetFormalCharge(cc.charge)

        current_charge = Chem.GetFormalCharge(mol)
        charge_diff = Chem.GetFormalCharge(mol) - start_charge
        # If molecule is now neutral, assume everything is now fixed
        # But otherwise, if charge has become more positive, look for additional protonated acid groups to ionize
        if not current_charge == 0:
            while charge_diff > 0:
                ppos, poccur = self._strongest_protonated(mol)
                if ppos is None:
                    break
                log.info('Ionizing %s to balance previous charge corrections', self.acid_base_pairs[ppos].name)
                patom = mol.GetAtomWithIdx(poccur[-1])
                patom.SetFormalCharge(patom.GetFormalCharge() - 1)
                if patom.GetNumExplicitHs() > 0:
                    patom.SetNumExplicitHs(patom.GetNumExplicitHs() - 1)
                # else:
                patom.UpdatePropertyCache()
                charge_diff -= 1

        already_moved = set()
        while True:
            ppos, poccur = self._strongest_protonated(mol)
            ipos, ioccur = self._weakest_ionized(mol)
            if ioccur and poccur and ppos < ipos:
                if poccur[-1] == ioccur[-1]:
                    # Bad! H wouldn't be moved, resulting in infinite loop.
                    log.warning('Aborted reionization due to unexpected situation')
                    break

                key = tuple(sorted([poccur[-1], ioccur[-1]]))
                if key in already_moved:
                    log.warning('Aborting reionization to avoid infinite loop due to it being ambiguous where to put a Hydrogen')
                    break
                already_moved.add(key)

                log.info('Moved proton from %s to %s', self.acid_base_pairs[ppos].name, self.acid_base_pairs[ipos].name)

                # Remove hydrogen from strongest protonated
                patom = mol.GetAtomWithIdx(poccur[-1])
                patom.SetFormalCharge(patom.GetFormalCharge() - 1)
                # If no implicit Hs to autoremove, and at least 1 explicit H to remove, reduce explicit count by 1
                if patom.GetNumImplicitHs() == 0 and patom.GetNumExplicitHs() > 0:
                    patom.SetNumExplicitHs(patom.GetNumExplicitHs() - 1)
                    # TODO: Remove any chiral label on patom?
                patom.UpdatePropertyCache()

                # Add hydrogen to weakest ionized
                iatom = mol.GetAtomWithIdx(ioccur[-1])
                iatom.SetFormalCharge(iatom.GetFormalCharge() + 1)
                # Increase explicit H count if no implicit, or aromatic N or P, or non default valence state
                if (iatom.GetNoImplicit() or
                        ((patom.GetAtomicNum() == 7 or patom.GetAtomicNum() == 15) and patom.GetIsAromatic()) or
                        iatom.GetTotalValence() not in list(Chem.GetPeriodicTable().GetValenceList(iatom.GetAtomicNum()))):
                    iatom.SetNumExplicitHs(iatom.GetNumExplicitHs() + 1)
                iatom.UpdatePropertyCache()
            else:
                break

        # TODO: Canonical ionization position if multiple equivalent positions?

        Chem.SanitizeMol(mol)
        return mol
Enforce charges on certain atoms, then perform competitive reionization.

        First, charge corrections are applied to ensure, for example, that free metals are correctly ionized. Then, if
        a molecule with multiple acid groups is partially ionized, ensure the 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.Chem.rdchem.Mol
        :return: The reionized molecule.
        :rtype: rdkit.Chem.rdchem.Mol