Source code for qcdb.parker
from __future__ import absolute_import
from __future__ import print_function
from .vecutil import *
from .physconst import *
from .cov_radii import *
BOND_FACTOR = 1.2 # fudge factor for bond length threshold
_expected_bonds = {
'H': 1,
'C': 4,
'N': 3,
'O': 2,
'F': 1,
'P': 3,
'S': 2,
}
[docs]def xyz2mol(self):
"""Returns a string of Molecule formatted for mol2.
Written by Trent M. Parker 9 Jun 2014
"""
factor = 1.0 if self.PYunits == 'Angstrom' else psi_bohr2angstroms
bonds = self.bond_profile()
N = 0
for i in range(self.natom()):
if self.Z(i):
N += 1
# header
text = '%s\n' % (self.tagline)
text += ' Generated by xyz2mol\n\n'
text += '%3i%3i 0 0 0 0 0 0 0 0999 V2000\n' % (N, len(bonds))
# coordinates
for i in range(self.natom()):
[x, y, z] = self.atoms[i].compute()
if self.Z(i):
text += ' %9.4f %9.4f %9.4f %-2s 0 0 0 0 0\n' % \
(x * factor, y * factor, z * factor, self.symbol(i))
# bonds
for p in range(len(bonds)):
text += '%3i%3i%3i' % (bonds[p][0] + 1, bonds[p][1] + 1, bonds[p][2])
text += ' 0 0 0\n'
text += 'M END\n'
return text
[docs]def missing_bonds(bonds, bond_tree, at_types):
"""Determine number of bonds missing for each atom"""
n_missing = []
for i in range(len(at_types)):
n_bonds_i = 0
for p in range(len(bonds)):
at1 = bonds[p][0]
at2 = bonds[p][1]
if (at1 == i or at2 == i):
bond_order = bonds[p][2]
n_bonds_i += bond_order
n_expect_i = _expected_bonds[at_types[i]]
n_missing.append(n_expect_i - n_bonds_i)
return n_missing
[docs]def missing_neighbors(bond_tree, n_missing):
"""Determine number of neighboring atoms missing bonds for each atom"""
missing_neighbors = []
for i in range(len(bond_tree)):
N_neighbors = len(bond_tree[i])
missing = 0
for a in range(N_neighbors):
j = bond_tree[i][a]
if n_missing[j] > 0:
missing += 1
missing_neighbors.append(missing)
return missing_neighbors
[docs]def bond_profile(self):
"""Obtain bonding topology of molecule"""
# determine bond topology from covalent radii
bonds = []
for i in range(self.natom()):
for j in range(i + 1, self.natom()):
dist = norm(sub(self.xyz(j), self.xyz(i))) * psi_bohr2angstroms
# TOOD check bohr/ang progress
bonded_dist = BOND_FACTOR * (psi_cov_radii[self.symbol(i)] + psi_cov_radii[self.symbol(j)])
if bonded_dist > dist:
bonds.append([i, j, 1])
# determine bond order from number of bonds
N_atoms = self.natom()
N_bonds = len(bonds)
at_types = [self.symbol(i) for i in range(self.natom())]
bond_tree = [[] for i in range(N_atoms)]
for i in range(N_bonds):
at1 = bonds[i][0]
at2 = bonds[i][1]
bond_tree[at1].append(at2)
bond_tree[at2].append(at1)
# determine bond order for all bonds from bond tree and element types
n_missing = missing_bonds(bonds, bond_tree, at_types)
n_neighbors_missing = missing_neighbors(bond_tree, n_missing)
# add double / triple bonds if only one neighbor missing bonds
N_left = math.floor(sum(n_missing) / 2)
N_left_previous = N_left + 1
N_iter = 0
while N_left > 0:
N_iter += 1
if N_left == N_left_previous:
neighbor_min += 1
else:
neighbor_min = 1
N_left_previous = N_left
# add a multiple bond to a deficient atom with the fewest number of deficient neighbors
BREAK_LOOP = False
for i in range(N_atoms):
if n_missing[i] > 0 and n_neighbors_missing[i] == neighbor_min:
N_neighbors = len(bond_tree[i])
for a in range(N_neighbors):
j = bond_tree[i][a]
if n_missing[j] > 0:
for p in range(N_bonds):
at1 = bonds[p][0]
at2 = bonds[p][1]
if (at1 == i and at2 == j) or (at1 == j and at2 == i):
bonds[p][2] += 1
n_missing[i] += -1
n_missing[j] += -1
n_neighbors_missing[i] += -1
n_neighbors_missing[j] += -1
N_left = math.floor(sum(n_missing) / 2)
BREAK_LOOP = True
if BREAK_LOOP:
break
if BREAK_LOOP:
break
# recalculate incomplete bond topology
n_missing = missing_bonds(bonds, bond_tree, at_types)
n_neighbors_missing = missing_neighbors(bond_tree, n_missing)
# break cycle if takes more than given number of iterations
max_iter = 100
if N_iter > max_iter:
print("""Error: multiple bond determination not complete""")
print(""" %i bonds unaccounted for""" % (N_left))
break
# bond order is number of bonds between each bonded atom pair
bond_order = []
for p in range(N_bonds):
bond_order.append(bonds[p][2])
for p in range(len(bond_order)):
bonds[p][2] = bond_order[p]
return bonds