Skip to content

[ BondOrder_System] add sanitize_guanidine_Catom, and formatting code #231

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
wanghan-iapcm merged 6 commits into from
Dec 16, 2021
Merged

[ BondOrder_System] add sanitize_guanidine_Catom, and formatting code #231

Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
1d1084d
Merge pull request #216 from amcadmus/master
wanghan-iapcm Nov 14, 2021
6569780
Merge branch 'master' into devel
Nov 17, 2021
f242cbe
Merge pull request #221 from amcadmus/master
wanghan-iapcm Nov 17, 2021
c377f85
Merge remote-tracking branch 'upstream/devel' into devel
Nov 23, 2021
306172d
Merge pull request #228 from amcadmus/master
wanghan-iapcm Nov 23, 2021
9685e2e
add sanitize_guanidine_Catom, and formatting code
TablewareBox Dec 5, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
108 changes: 89 additions & 19 deletions dpdata/rdkit/sanitize.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,24 +3,29 @@
from rdkit.Chem.rdchem import Atom, Bond, Mol, BondType
import os
import time

# openbabel
try:
from openbabel import openbabel

USE_OBABEL = True
except ModuleNotFoundError as e:
USE_OBABEL = False


def get_explicit_valence(atom, verbose=False):
exp_val_calculated_from_bonds = int(sum([bond.GetBondTypeAsDouble() for bond in atom.GetBonds()]))
try:
exp_val = atom.GetExplicitValence()
if exp_val != exp_val_calculated_from_bonds:
if verbose:
print(f"Explicit valence given by GetExplicitValence() and sum of bond order are inconsistent on {atom.GetSymbol()}{atom.GetIdx()+1}, using sum of bond order.")
print(
f"Explicit valence given by GetExplicitValence() and sum of bond order are inconsistent on {atom.GetSymbol()}{atom.GetIdx() + 1}, using sum of bond order.")
return exp_val_calculated_from_bonds
except:
return exp_val_calculated_from_bonds


def regularize_formal_charges(mol, sanitize=True, verbose=False):
"""
Regularize formal charges of atoms
Expand All @@ -37,6 +42,7 @@ def regularize_formal_charges(mol, sanitize=True, verbose=False):
else:
return mol


def assign_formal_charge_for_atom(atom, verbose=False):
"""
assigen formal charge according to 8-electron rule for element B,C,N,O,S,P,As
Expand All @@ -48,12 +54,13 @@ def assign_formal_charge_for_atom(atom, verbose=False):
elif atom.GetSymbol() == "C":
atom.SetFormalCharge(valence - 4)
if valence == 3:
print(f"Detect a valence of 3 on #C{atom.GetIdx()+1}, the formal charge of this atom will be assigned to -1")
print(
f"Detect a valence of 3 on #C{atom.GetIdx() + 1}, the formal charge of this atom will be assigned to -1")
elif valence > 4:
raise ValueError(f"#C{atom.GetIdx()+1} has a valence larger than 4")
raise ValueError(f"#C{atom.GetIdx() + 1} has a valence larger than 4")
elif atom.GetSymbol() == "N":
if valence > 4:
raise ValueError(f"#N{atom.GetIdx()+1} has a valence larger than 4")
raise ValueError(f"#N{atom.GetIdx() + 1} has a valence larger than 4")
else:
atom.SetFormalCharge(valence - 3)
elif atom.GetSymbol() == "O":
Expand All @@ -64,34 +71,35 @@ def assign_formal_charge_for_atom(atom, verbose=False):
elif valence == 3:
atom.SetFormalCharge(1)
elif valence > 6:
raise ValueError(f"#S{atom.GetIdx()+1} has a valence larger than 6")
raise ValueError(f"#S{atom.GetIdx() + 1} has a valence larger than 6")
else:
atom.SetFormalCharge(0)
elif atom.GetSymbol() == "P" or atom.GetSymbol() == "As":
if valence == 5:
atom.SetFormalCharge(0)
elif valence > 5:
raise ValueError(f"#{atom.GetSymbol()}{atom.GetIdx()+1} has a valence larger than 5")
raise ValueError(f"#{atom.GetSymbol()}{atom.GetIdx() + 1} has a valence larger than 5")
else:
atom.SetFormalCharge(valence - 3)


# print bond and atom information (for debugger)
def print_bonds(mol):
for bond in mol.GetBonds():
begin_atom = bond.GetBeginAtom()
end_atom = bond.GetEndAtom()
print(f'{begin_atom.GetSymbol()}{begin_atom.GetIdx()+1} {end_atom.GetSymbol()}{end_atom.GetIdx()+1} {bond.GetBondType()}')
print(
f'{begin_atom.GetSymbol()}{begin_atom.GetIdx() + 1} {end_atom.GetSymbol()}{end_atom.GetIdx() + 1} {bond.GetBondType()}')


def print_atoms(mol):
for atom in mol.GetAtoms():
print(f'{atom.GetSymbol()}{atom.GetIdx()+1} {atom.GetFormalCharge()} {get_explicit_valence(atom)}')
print(f'{atom.GetSymbol()}{atom.GetIdx() + 1} {atom.GetFormalCharge()} {get_explicit_valence(atom)}')


def is_terminal_oxygen(O_atom):
if len(O_atom.GetNeighbors()) == 1:
return True
else:
return False
return len(O_atom.GetNeighbors()) == 1


def get_terminal_oxygens(atom):
terminal_oxygens = []
Expand All @@ -101,6 +109,21 @@ def get_terminal_oxygens(atom):
terminal_oxygens.append(nei)
return terminal_oxygens


def is_terminal_NR2(N_atom):
return len(N_atom.GetNeighbors()) == 3


def get_terminal_NR2s(atom):
terminal_NR2s = []
for nei in atom.GetNeighbors():
if nei.GetSymbol() == "N":
if is_terminal_NR2(nei):
terminal_NR2s.append(nei)
terminal_NR2s.sort(key=lambda N_atom: len([atom for atom in N_atom.GetNeighbors() if atom.GetSymbol() == 'H']))
return terminal_NR2s


def sanitize_phosphate_Patom(P_atom, verbose=True):
if P_atom.GetSymbol() == "P":
terminal_oxygens = get_terminal_oxygens(P_atom)
Expand All @@ -116,11 +139,13 @@ def sanitize_phosphate_Patom(P_atom, verbose=True):
bond.SetBondType(Chem.rdchem.BondType.SINGLE)
terminal_oxygens[ii].SetFormalCharge(-1)


def sanitize_phosphate(mol):
for atom in mol.GetAtoms():
sanitize_phosphate_Patom(atom)
return mol


def sanitize_sulfate_Satom(S_atom, verbose=True):
if S_atom.GetSymbol() == "S":
terminal_oxygens = get_terminal_oxygens(S_atom)
Expand All @@ -136,11 +161,13 @@ def sanitize_sulfate_Satom(S_atom, verbose=True):
bond = mol.GetBondBetweenAtoms(S_atom.GetIdx(), terminal_oxygens[ii].GetIdx())
bond.SetBondType(Chem.rdchem.BondType.DOUBLE)


def sanitize_sulfate(mol):
for atom in mol.GetAtoms():
sanitize_sulfate_Satom(atom)
return mol


def sanitize_carboxyl_Catom(C_atom, verbose=True):
if C_atom.GetSymbol() == "C":
terminal_oxygens = get_terminal_oxygens(C_atom)
Expand All @@ -157,11 +184,40 @@ def sanitize_carboxyl_Catom(C_atom, verbose=True):
bond2.SetBondType(Chem.rdchem.BondType.DOUBLE)
terminal_oxygens[1].SetFormalCharge(0)


def sanitize_carboxyl(mol):
for atom in mol.GetAtoms():
sanitize_carboxyl_Catom(atom)
return mol


def sanitize_guanidine_Catom(C_atom, verbose=True):
if C_atom.GetSymbol() == "C":
terminal_NR2s = get_terminal_NR2s(C_atom)
mol = C_atom.GetOwningMol()
if len(terminal_NR2s) == 3:
if verbose:
print("Guanidyl group detected, sanitizing it...")
# set two C-N and one C=N+
bond1 = mol.GetBondBetweenAtoms(C_atom.GetIdx(), terminal_NR2s[0].GetIdx())
bond1.SetBondType(Chem.rdchem.BondType.SINGLE)
terminal_NR2s[0].SetFormalCharge(-1)

bond2 = mol.GetBondBetweenAtoms(C_atom.GetIdx(), terminal_NR2s[1].GetIdx())
bond2.SetBondType(Chem.rdchem.BondType.SINGLE)
terminal_NR2s[1].SetFormalCharge(0)

bond3 = mol.GetBondBetweenAtoms(C_atom.GetIdx(), terminal_NR2s[2].GetIdx())
bond3.SetBondType(Chem.rdchem.BondType.DOUBLE)
terminal_NR2s[2].SetFormalCharge(1)


def sanitize_guanidine(mol):
for atom in mol.GetAtoms():
sanitize_guanidine_Catom(atom)
return mol


def sanitize_nitro_Natom(N_atom, verbose=True):
if N_atom.GetSymbol() == "N":
terminal_oxygens = get_terminal_oxygens(N_atom)
Expand All @@ -178,17 +234,20 @@ def sanitize_nitro_Natom(N_atom, verbose=True):
bond2.SetBondType(Chem.rdchem.BondType.DOUBLE)
terminal_oxygens[1].SetFormalCharge(0)


def sanitize_nitro(mol):
for atom in mol.GetAtoms():
sanitize_nitro_Natom(atom)
return mol


def is_terminal_nitrogen(N_atom):
if N_atom.GetSymbol() == 'N' and len(N_atom.GetNeighbors()) == 1:
return True
else:
return False


def sanitize_nitrine_Natom(atom, verbose=True):
if atom.GetSymbol() == "N" and len(atom.GetNeighbors()) == 2:
mol = atom.GetOwningMol()
Expand All @@ -213,7 +272,8 @@ def sanitize_nitrine_Natom(atom, verbose=True):

bond = mol.GetBondBetweenAtoms(atom.GetIdx(), N_non_terminal.GetIdx())
bond.SetBondType(Chem.rdchem.BondType.DOUBLE)
atom.SetFormalCharge(1)
atom.SetFormalCharge(1)


def contain_hetero_aromatic(mol):
flag = False
Expand All @@ -223,6 +283,7 @@ def contain_hetero_aromatic(mol):
break
return flag


# for carbon with explicit valence > 4
def regularize_carbon_bond_order(atom, verbose=True):
if atom.GetSymbol() == "C" and get_explicit_valence(atom) > 4:
Expand All @@ -240,6 +301,7 @@ def regularize_carbon_bond_order(atom, verbose=True):
if bond.GetIdx() != double_bond_idx:
bond.SetBondType(Chem.rdchem.BondType.SINGLE)


# for nitrogen with explicit valence > 4
def regularize_nitrogen_bond_order(atom, verbose=True):
mol = atom.GetOwningMol()
Expand All @@ -255,6 +317,7 @@ def regularize_nitrogen_bond_order(atom, verbose=True):
def sanitize_mol(mol, verbose=False):
for atom in mol.GetAtoms():
sanitize_carboxyl_Catom(atom, verbose)
sanitize_guanidine_Catom(atom, verbose)
sanitize_phosphate_Patom(atom, verbose)
sanitize_sulfate_Satom(atom, verbose)
sanitize_nitro_Natom(atom, verbose)
Expand All @@ -272,6 +335,7 @@ def mol_edit_log(mol, i, j):
edited = mol.GetProp("edit")
mol.SetProp("edit", edited + ",%d_%d" % (i, j))


def kekulize_aromatic_heterocycles(mol_in, assign_formal_charge=True, sanitize=True):
mol = Chem.RWMol(mol_in)
rings = Chem.rdmolops.GetSymmSSSR(mol)
Expand Down Expand Up @@ -345,7 +409,7 @@ def hetero_priority(idx, mol):
elif bAllAr and not bAllC:
HAr.append(ring)

if len(HAr) == 0:
if len(HAr) == 0:
# no hetrerocycles
return mol_in
else:
Expand All @@ -364,7 +428,7 @@ def hetero_priority(idx, mol):
fuseCAr[i] = j
break
if i > 1:
if (fuseCAr[i] == fuseCAr[i-1]) & (fuseCAr[i] >= 0):
if (fuseCAr[i] == fuseCAr[i - 1]) & (fuseCAr[i] >= 0):
fuseDouble.append(i)
atom = mol.GetAtomWithIdx(ring[i])
if atom.GetSymbol() != 'C':
Expand All @@ -376,7 +440,7 @@ def hetero_priority(idx, mol):
hasDouble.append(i)
bond = mol.GetBondBetweenAtoms(ring[i], ring[(i + 1) % lring])

if (fuseCAr[0] == fuseCAr[lring-1]) & (fuseCAr[0] >= 0):
if (fuseCAr[0] == fuseCAr[lring - 1]) & (fuseCAr[0] >= 0):
fuseDouble.append(0)

if (len(hetero) > 0) | (len(hasDouble) > 0):
Expand Down Expand Up @@ -446,6 +510,7 @@ def hetero_priority(idx, mol):
except Exception as e:
raise RuntimeError(f"Manual kekulization for aromatic heterocycles failed, below are errors:\n\t {e}")


def convert_by_obabel(mol, cache_dir=os.path.join(os.getcwd(), '.cache'), obabel_path="obabel"):
if not os.path.exists(cache_dir):
os.mkdir(cache_dir)
Expand All @@ -464,6 +529,7 @@ def convert_by_obabel(mol, cache_dir=os.path.join(os.getcwd(), '.cache'), obabel
mol_obabel = Chem.MolFromMolFile(mol_file_out, removeHs=False, sanitize=False)
return mol_obabel


def super_sanitize_mol(mol, name=None, verbose=True):
if name is None:
if mol.HasProp("_Name"):
Expand All @@ -484,11 +550,13 @@ def super_sanitize_mol(mol, name=None, verbose=True):
except Exception as e:
try:
if verbose:
print("Hermite procedure failed, maybe due to unsupported representation of hetero aromatic rings, re-try with obabel")
print(
"Hermite procedure failed, maybe due to unsupported representation of hetero aromatic rings, re-try with obabel")
print("=====Stage 2: re-try with obabel=====")
mol = convert_by_obabel(mol)
mol = sanitize_mol(mol, verbose)
mol = kekulize_aromatic_heterocycles(mol, assign_formal_charge=False, sanitize=False) # aromatic heterocycles
mol = kekulize_aromatic_heterocycles(mol, assign_formal_charge=False,
sanitize=False) # aromatic heterocycles
mol = regularize_formal_charges(mol, sanitize=False)
mol_copy = deepcopy(mol)
Chem.SanitizeMol(mol_copy)
Expand All @@ -501,6 +569,7 @@ def super_sanitize_mol(mol, name=None, verbose=True):
print(name, "Failed!")
return None


class Sanitizer(object):
def __init__(self, level='medium', raise_errors=True, verbose=False):
'''
Expand All @@ -526,7 +595,7 @@ def _check_level(self, level):
else:
if level == 'high' and not USE_OBABEL:
raise ModuleNotFoundError("obabel not installed, high level sanitizer cannot work")

def _handle_exception(self, error_info):
if self.raise_errors:
raise SanitizeError(error_info)
Expand Down Expand Up @@ -561,6 +630,7 @@ def sanitize(self, mol):
self._handle_exception(error_info)
return mol


class SanitizeError(Exception):
def __init__(self, content="Sanitization Failed."):
self.content = content
Expand Down