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