"""
Ligand atom typing
==================
This module contains the :func:`assign_sybyl_type` function to analyze
all :class:`propka.atom.Atom` in terms of SYBYL atom types (see
:data:`ALL_SYBYL_TYPES`).
"""
from propka.calculations import squared_distance
from propka.vector_algebra import Vector
#: SYBYL atom types
ALL_SYBYL_TYPES = [
'C.3', # carbon sp3
'H', # hydrogen
'C.2', # carbon sp2
'H.spc', # hydrogen in Single Point Charge (SPC) water model
'C.1', # carbon sp
'H.t3p', # hydrogen in Transferable intermolecular Potential (TIP3P) water model
'C.ar', # carbon aromatic
'LP', # lone pair
'C.cat', # carbocation (C+) used only in a guadinium group
'Du', # dummy atom
'N.3', # nitrogen sp3
'Du.C', # dummy carbon
'N.2', # nitrogen sp2
'Any', # any atom
'N.1', # nitrogen sp
'Hal', # halogen
'N.ar', # nitrogen aromatic
'Het', # heteroatom = N, O, S, P
'N.am', # nitrogen amide
'Hev', # heavy atom (non hydrogen)
'N.pl3', # nitrogen trigonal planar
'Li', # lithium
'N.4', # nitrogen sp3 positively charged
'Na', # sodium
'O.3', # oxygen sp3
'Mg', # magnesium
'O.2', # oxygen sp2
'Al', # aluminum
'O.co2', # oxygen in carboxylate and phosphate groups
'Si', # silicon
'O.spc', # oxygen in Single Point Charge (SPC) water model
'K', # potassium
'O.t3p', # oxygen in Transferable Intermolecular Potential (TIP3P) water model
'Ca', # calcium
'S.3', # sulfur sp3
'Cr.th', # chromium (tetrahedral)
'S.2', # sulfur sp2
'Cr.oh', # chromium (octahedral)
'S.O', # sulfoxide sulfur
'Mn', # manganese
'S.O2', # sulfone sulfur
'Fe', # iron
'P.3', # phosphorous sp3
'Co.oh', # cobalt (octahedral)
'F', # fluorine
'Cu', # copper
'Cl', # chlorine
'Zn', # zinc
'Br', # bromine
'Se', # selenium
'I', # iodine
'Mo', # molybdenum
'Sn'] # tin
#: PROPKA input types
PROPKA_INPUT_TYPES = ['1P', '1N', '2P', '2N', 'C3', 'H', 'C2', 'Hsp', 'C1',
'Ht3', 'Car', 'LP', 'Cca', 'Du', 'N3', 'DuC', 'N2',
'Any', 'N1', 'Hal', 'Nar', 'Het', 'Nam', 'Hev', 'Npl',
'Li', 'N4', 'Na', 'O3', 'Mg', 'O2', 'Al', 'Oco', 'Si',
'Osp', 'K', 'Ot3', 'Ca', 'S3', 'Crt', 'S2', 'Cro', 'SO',
'Mn', 'SO2', 'Fe', 'P3', 'Coo', 'F', 'Cu', 'Cl', 'Zn',
'Br', 'Se', 'I', 'Mo', 'Sn']
MAX_C_DOUBLE_BOND = 1.3
MAX_C_TRIPLE_BOND = 1.2
MAX_C_DOUBLE_BOND_SQUARED = MAX_C_DOUBLE_BOND*MAX_C_DOUBLE_BOND
MAX_C_TRIPLE_BOND_SQUARED = MAX_C_TRIPLE_BOND*MAX_C_TRIPLE_BOND
PLANARITY_MARGIN = 0.20
[docs]
def assign_sybyl_type(atom):
"""Assign Sybyl type to atom.
Args:
atom: atom to assign
"""
# check if we already have assigned a name to this atom
if atom.sybyl_assigned:
return
# find some properties of the atom
ring_atoms = is_ring_member(atom)
aromatic = is_aromatic_ring(ring_atoms)
planar = is_planar(atom)
bonded_elements = {}
for i, bonded_atom in enumerate(atom.bonded_atoms):
bonded_elements[i] = bonded_atom.element
# Aromatic carbon/nitrogen
if aromatic:
for ring_atom in ring_atoms:
if ring_atom.element in ['C', 'N']:
set_type(ring_atom, ring_atom.element+'.ar')
return
# check for amide
if atom.element in ['O', 'N', 'C']:
o_atom = None
n_atom = None
c_atom = None
# oxygen, nitrogen
if atom.element in ['O', 'N']:
for bonded_elem in atom.get_bonded_elements('C'):
for bonded_atom in bonded_elem.bonded_atoms:
if (bonded_atom.element == 'N' and atom.element == 'O'):
o_atom = atom
c_atom = bonded_elem
n_atom = bonded_atom
elif (bonded_atom.element == 'O' and atom.element == 'N'):
n_atom = atom
c_atom = bonded_elem
o_atom = bonded_atom
# carbon
if atom.element == 'C':
nitrogens = atom.get_bonded_elements('N')
oxygens = atom.get_bonded_elements('O')
if len(nitrogens) == 1 and len(oxygens) == 1:
c_atom = atom
n_atom = nitrogens[0]
o_atom = oxygens[0]
if c_atom and n_atom and o_atom:
# make sure that the Nitrogen is not aromatic and that it has two
# heavy atom bonds
if (not is_aromatic_ring(is_ring_member(n_atom))
and len(n_atom.get_bonded_heavy_atoms()) == 2):
set_type(n_atom, 'N.am')
set_type(c_atom, 'C.2')
set_type(o_atom, 'O.2')
return
if atom.element == 'C':
# check for carboxyl
if (len(atom.bonded_atoms) == 3
and list(bonded_elements.values()).count('O') == 2):
index1 = list(bonded_elements.values()).index('O')
index2 = list(bonded_elements.values()).index('O', index1+1)
if (len(atom.bonded_atoms[index1].bonded_atoms) == 1
and len(atom.bonded_atoms[index2].bonded_atoms) == 1):
set_type(atom.bonded_atoms[index1], 'O.co2-')
set_type(atom.bonded_atoms[index2], 'O.co2')
set_type(atom, 'C.2')
return
# sp carbon
if len(atom.bonded_atoms) <= 2:
for bonded_atom in atom.bonded_atoms:
if (squared_distance(atom, bonded_atom)
< MAX_C_TRIPLE_BOND_SQUARED):
set_type(atom, 'C.1')
set_type(bonded_atom, bonded_atom.element + '.1')
if atom.sybyl_assigned:
return
# sp2 carbon
if planar:
set_type(atom, 'C.2')
# check for N.pl3
for bonded_atom in atom.bonded_atoms:
if bonded_atom.element == 'N':
if (len(bonded_atom.bonded_atoms) < 3
or is_planar(bonded_atom)):
set_type(bonded_atom, 'N.pl3')
return
# sp3 carbon
set_type(atom, 'C.3')
return
# Nitrogen
if atom.element == 'N':
# check for planar N
if len(atom.bonded_atoms) == 1:
if is_planar(atom.bonded_atoms[0]):
set_type(atom, 'N.pl3')
return
if planar:
set_type(atom, 'N.pl3')
return
set_type(atom, 'N.3')
return
# Oxygen
if atom.element == 'O':
set_type(atom, 'O.3')
if len(atom.bonded_atoms) == 1:
# check for carboxyl
if atom.bonded_atoms[0].element == 'C':
the_carbon = atom.bonded_atoms[0]
if (len(the_carbon.bonded_atoms) == 3
and the_carbon.count_bonded_elements('O') == 2):
[oxy1, oxy2] = the_carbon.get_bonded_elements('O')
if (len(oxy1.bonded_atoms) == 1
and len(oxy2.bonded_atoms) == 1):
set_type(oxy1, 'O.co2-')
set_type(oxy2, 'O.co2')
set_type(the_carbon, 'C.2')
return
# check for X=O
if (squared_distance(atom, atom.bonded_atoms[0])
< MAX_C_DOUBLE_BOND_SQUARED):
set_type(atom, 'O.2')
if atom.bonded_atoms[0].element == 'C':
set_type(atom.bonded_atoms[0], 'C.2')
return
# Sulphur
if atom.element == 'S':
# check for SO2
if list(bonded_elements.values()).count('O') == 2:
index1 = list(bonded_elements.values()).index('O')
index2 = list(bonded_elements.values()).index('O', index1+1)
set_type(atom.bonded_atoms[index1], 'O.2')
set_type(atom.bonded_atoms[index2], 'O.2')
set_type(atom, 'S.o2')
return
# check for SO4
if list(bonded_elements.values()).count('O') == 4:
no_o2 = 0
for i in range(len(atom.bonded_atoms)):
if len(atom.bonded_atoms[i].bonded_atoms) == 1 and no_o2 < 2:
set_type(atom.bonded_atoms[i], 'O.2')
no_o2 += 1
else:
set_type(atom.bonded_atoms[i], 'O.3')
set_type(atom, 'S.3')
return
# Phosphorus
if atom.element == 'P':
set_type(atom, 'P.3')
# check for phosphate group
bonded_oxygens = atom.get_bonded_elements('O')
for o_atom in bonded_oxygens:
set_type(o_atom, 'O.3')
return
element = atom.element.capitalize()
set_type(atom, element)
[docs]
def is_ring_member(atom):
"""Determine if atom is a member of a ring.
Args:
atom: atom to test
Returns:
list of atoms
"""
return identify_ring(atom, atom, 0, [])
[docs]
def identify_ring(this_atom, original_atom, number, past_atoms):
"""Identify the atoms in a ring
Args:
this_atom: atom to test
original_atom: some other atom
number: number of atoms
past_atoms: atoms that have already been found
Returns:
list of atoms
"""
number += 1
past_atoms = past_atoms + [this_atom]
return_atoms = []
if number > 10:
return return_atoms
for atom in this_atom.get_bonded_heavy_atoms():
if atom == original_atom and number > 2:
return past_atoms
if atom not in past_atoms:
these_return_atoms = identify_ring(atom, original_atom, number,
past_atoms)
if len(these_return_atoms) > 0:
if (len(return_atoms) > len(these_return_atoms)
or len(return_atoms) == 0):
return_atoms = these_return_atoms
return return_atoms
[docs]
def set_type(atom, type_):
"""Set atom type..
Args:
atom: atom to set
type_: type value to set
"""
atom.sybyl_type = type_
atom.sybyl_assigned = True
[docs]
def is_planar(atom):
"""Finds out if atom forms a plane together with its bonded atoms.
Args:
atom: atom to test
Returns:
Boolean
"""
atoms = [atom] + atom.bonded_atoms
return are_atoms_planar(atoms)
[docs]
def are_atoms_planar(atoms):
"""Test whether a group of atoms are planar.
Args:
atoms: list of atoms
Returns:
Boolean
"""
if len(atoms) == 0:
return False
if len(atoms) < 4:
return False
vec1 = Vector(atom1=atoms[0], atom2=atoms[1])
vec2 = Vector(atom1=atoms[0], atom2=atoms[2])
norm = vec1.cross(vec2).rescale(1.0)
margin = PLANARITY_MARGIN
for atom in atoms[3:]:
vec = Vector(atom1=atoms[0], atom2=atom).rescale(1.0)
if abs(vec.dot(norm)) > margin:
return False
return True
[docs]
def is_aromatic_ring(atoms):
"""Determine whether group of atoms form aromatic ring.
Args:
atoms: list of atoms to test
Returns:
Boolean
"""
if len(atoms) < 5:
return False
for i in range(len(atoms)):
if not are_atoms_planar(atoms[i:]+atoms[:i]):
return False
return True