1# -*- coding: utf-8 -*- 2# MolMod is a collection of molecular modelling tools for python. 3# Copyright (C) 2007 - 2019 Toon Verstraelen <Toon.Verstraelen@UGent.be>, Center 4# for Molecular Modeling (CMM), Ghent University, Ghent, Belgium; all rights 5# reserved unless otherwise stated. 6# 7# This file is part of MolMod. 8# 9# MolMod is free software; you can redistribute it and/or 10# modify it under the terms of the GNU General Public License 11# as published by the Free Software Foundation; either version 3 12# of the License, or (at your option) any later version. 13# 14# MolMod is distributed in the hope that it will be useful, 15# but WITHOUT ANY WARRANTY; without even the implied warranty of 16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17# GNU General Public License for more details. 18# 19# You should have received a copy of the GNU General Public License 20# along with this program; if not, see <http://www.gnu.org/licenses/> 21# 22# -- 23"""Bond length database and bond type definitions 24 25 An object ``bonds`` of the class ``BondData`` is created upon importing this 26 module. It loads information about average bond lengths from a csv file 27 when it is initialized. Missing data points in the csv file are estimated 28 by adding Van der Waals radii of the two atoms of the given bond. 29 30 This module also defines a few constants for different bond types: 31 32 ============= ===== 33 Name Value 34 ============= ===== 35 BOND_SINGLE 0 36 BOND_DOUBLE 1 37 BOND_TRIPLE 2 38 BOND_HYBRID 3 39 BOND_HYDROGEN 4 40 ============= ===== 41""" 42 43 44from __future__ import division 45 46from builtins import range 47import pkg_resources 48 49from molmod.periodic import periodic 50import molmod.units as units 51 52 53__all__ = [ 54 "BOND_SINGLE", "BOND_DOUBLE", "BOND_TRIPLE", "BOND_HYBRID", 55 "BOND_HYDROGEN", "bond_types", "BondData", "bonds" 56] 57 58 59BOND_SINGLE = 1 60BOND_DOUBLE = 2 61BOND_TRIPLE = 3 62BOND_HYBRID = 4 63BOND_HYDROGEN = 5 64 65bond_types = [ 66 BOND_SINGLE, BOND_DOUBLE, BOND_TRIPLE, BOND_HYBRID, BOND_HYDROGEN 67] 68 69 70class BondData(object): 71 """Database with bond lengths""" 72 bond_tolerance = 1.2 73 74 def __init__(self): 75 """ 76 This object is created when importing this module. There is no need 77 to do it a second time manually. 78 """ 79 self.lengths = dict([bond_type, {}] for bond_type in bond_types) 80 self._load_bond_data() 81 self._approximate_unkown_bond_lengths() 82 self.max_length = max( 83 max(lengths.values()) 84 for lengths 85 in self.lengths.values() 86 if len(lengths) > 0 87 ) 88 89 def _load_bond_data(self): 90 """Load the bond data from the given file 91 92 It's assumed that the uncommented lines in the data file have the 93 following format: 94 symbol1 symbol2 number1 number2 bond_length_single_a bond_length_double_a bond_length_triple_a bond_length_single_b bond_length_double_b bond_length_triple_b ..." 95 where a, b, ... stand for different sources. 96 """ 97 98 def read_units(unit_names): 99 """convert unit_names into conversion factors""" 100 tmp = { 101 "A": units.angstrom, 102 "pm": units.picometer, 103 "nm": units.nanometer, 104 } 105 return [tmp[unit_name] for unit_name in unit_names] 106 107 def read_length(BOND_TYPE, words, col): 108 """Read the bondlengths from a single line in the data file""" 109 nlow = int(words[2]) 110 nhigh = int(words[3]) 111 for i, conversion in zip(range((len(words) - 4) // 3), conversions): 112 word = words[col + 3 + i*3] 113 if word != 'NA': 114 self.lengths[BOND_TYPE][frozenset([nlow, nhigh])] = float(word)*conversion 115 return 116 117 with pkg_resources.resource_stream(__name__, 'data/bonds.csv') as f: 118 for line in f: 119 words = line.decode('utf-8').split() 120 if (len(words) > 0) and (words[0][0] != "#"): 121 if words[0] == "unit": 122 conversions = read_units(words[1:]) 123 else: 124 read_length(BOND_SINGLE, words, 1) 125 read_length(BOND_DOUBLE, words, 2) 126 read_length(BOND_TRIPLE, words, 3) 127 128 def _approximate_unkown_bond_lengths(self): 129 """Completes the bond length database with approximations based on VDW radii""" 130 dataset = self.lengths[BOND_SINGLE] 131 for n1 in periodic.iter_numbers(): 132 for n2 in periodic.iter_numbers(): 133 if n1 <= n2: 134 pair = frozenset([n1, n2]) 135 atom1 = periodic[n1] 136 atom2 = periodic[n2] 137 #if (pair not in dataset) and hasattr(atom1, "covalent_radius") and hasattr(atom2, "covalent_radius"): 138 if (pair not in dataset) and (atom1.covalent_radius is not None) and (atom2.covalent_radius is not None): 139 dataset[pair] = (atom1.covalent_radius + atom2.covalent_radius) 140 #print "%3i %3i %s %30s %30s" % (n1, n2, dataset.get(pair), atom1, atom2) 141 142 def bonded(self, n1, n2, distance): 143 """Return the estimated bond type 144 145 Arguments: 146 | ``n1`` -- the atom number of the first atom in the bond 147 | ``n2`` -- the atom number of the second atom the bond 148 | ``distance`` -- the distance between the two atoms 149 150 This method checks whether for the given pair of atom numbers, the 151 given distance corresponds to a certain bond length. The best 152 matching bond type will be returned. If the distance is a factor 153 ``self.bond_tolerance`` larger than a tabulated distance, the 154 algorithm will not relate them. 155 """ 156 if distance > self.max_length * self.bond_tolerance: 157 return None 158 159 deviation = 0.0 160 pair = frozenset([n1, n2]) 161 result = None 162 for bond_type in bond_types: 163 bond_length = self.lengths[bond_type].get(pair) 164 if (bond_length is not None) and \ 165 (distance < bond_length * self.bond_tolerance): 166 if result is None: 167 result = bond_type 168 deviation = abs(bond_length - distance) 169 else: 170 new_deviation = abs(bond_length - distance) 171 if deviation > new_deviation: 172 result = bond_type 173 deviation = new_deviation 174 return result 175 176 def get_length(self, n1, n2, bond_type=BOND_SINGLE): 177 """Return the length of a bond between n1 and n2 of type bond_type 178 179 Arguments: 180 | ``n1`` -- the atom number of the first atom in the bond 181 | ``n2`` -- the atom number of the second atom the bond 182 183 Optional argument: 184 | ``bond_type`` -- the type of bond [default=BOND_SINGLE] 185 186 This is a safe method for querying a bond_length. If no answer can be 187 found, this get_length returns None. 188 """ 189 dataset = self.lengths.get(bond_type) 190 if dataset == None: 191 return None 192 return dataset.get(frozenset([n1, n2])) 193 194 195bonds = BondData() 196