Skip to content

Commit

Permalink
Update bond definitions (duartegroup#151)
Browse files Browse the repository at this point in the history 
* Migrate to bonds to using covalent radii

* Update maximum valence

* Update changelog

* Add tests
  • Loading branch information
@t-young31
t-young31 committed Oct 1, 2022
1 parent 811c34f commit ce7e8c1
Show file tree
Hide file tree
Showing 16 changed files with 259 additions and 204 deletions.
25 changes: 12 additions & 13 deletions autode/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,23 +1,22 @@
from autode.wrappers.keywords import (KeywordsSet, OptKeywords,
HessianKeywords, SinglePointKeywords,
Keywords, GradientKeywords)
from autode.reactions.reaction import Reaction
from autode.reactions.multistep import MultiStepReaction
from autode.atoms import Atom
from autode.species.molecule import (Reactant, Product, Molecule, Species)
from autode.species.complex import NCIComplex
from autode.config import Config
from autode.calculation import Calculation
from autode import methods
from autode import geom
from autode import pes
from autode import utils
from autode import neb
from autode import mol_graphs
from autode import hessians

from autode.reactions.reaction import Reaction
from autode.reactions.multistep import MultiStepReaction
from autode.atoms import Atom
from autode.species.molecule import (Reactant, Product, Molecule, Species)
from autode.species.complex import NCIComplex
from autode.config import Config
from autode.calculation import Calculation
from autode.wrappers.keywords import (KeywordsSet, OptKeywords,
HessianKeywords, SinglePointKeywords,
Keywords, GradientKeywords)
"""
So, you want to bump the version.. make sure the following steps are followed
Bumping the version number requires following the release proceedure:
- Update changelog (doc/changelog.rst)
Expand All @@ -33,7 +32,7 @@
- Merge when tests pass
"""

__version__ = '1.2.3'
__version__ = '1.3.0'


__all__ = [
Expand Down
201 changes: 153 additions & 48 deletions autode/atoms.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -296,17 +296,17 @@ def maximal_valance(self) -> int:
Returns:
(int): Maximal valance
"""
max_valances = {'H': 1, 'B': 4, 'C': 4, 'N': 4, 'O': 3, 'F': 1,
'Si': 4, 'P': 6, 'S': 6, 'Cl': 4, 'Br': 4, 'I': 6}

if self.label in max_valances:
return max_valances[self.label]

else:
logger.warning(f'Could not find a valid valance for {self}. '
f'Guessing at 6')
if self.is_metal:
return 6

if self.label in _max_valances:
return _max_valances[self.label]

logger.warning(f'Could not find a valid valance for {self}. '
f'Guessing at 6')
return 6

@property
def vdw_radius(self) -> Distance:
"""
Expand All @@ -332,6 +332,24 @@ def vdw_radius(self) -> Distance:

return Distance(radius, 'Å')

@property
def covalent_radius(self) -> Distance:
"""
Covalent radius for this atom. Example:
.. code-block:: Python
>>> import autode as ade
>>> ade.Atom('H').covalent_radius
Distance(0.31 Å)
-----------------------------------------------------------------------
Returns:
(autode.values.Distance): Van der Waals radius
"""
radius = Distance(_covalent_radii_pm[self.atomic_number-1], units='pm')
return radius.to('Å')

def is_pi(self, valency: int) -> bool:
"""
Determine if this atom is a 'π-atom' i.e. is unsaturated. Only
Expand Down Expand Up @@ -518,6 +536,16 @@ def mass(self) -> Mass:
"""Dummy atoms do not have any weight/mass"""
return Mass(0.0)

@property
def vdw_radius(self) -> Distance:
"""Dummy atoms have no radius"""
return Distance(0., units='Å')

@property
def covalent_radius(self) -> Distance:
"""Dummy atoms have no radius"""
return Distance(0., units='Å')


class Atoms(list):

Expand Down Expand Up @@ -626,6 +654,78 @@ def contain_metals(self) -> bool:
"""
return any(atom.label in metals for atom in self)

def idxs_are_present(self, *args: int) -> bool:
"""Are all these indexes present in this set of atoms"""
return set(args).issubset(set(range(len(self))))

def eqm_bond_distance(self,
i: int,
j: int) -> Distance:
"""
Equilibrium distance between two atoms. If known then use the
experimental dimer distance, otherwise estimate if from the
covalent radii of the two atoms. Example
Example:
.. code-block:: Python
>>> import autode as ade
>>> mol = ade.Molecule(atoms=[ade.Atom('H'), ade.Atom('H')])
>>> mol.distance(0, 1)
Distance(0.0 Å)
>>> mol.eqm_bond_distance(0, 1)
Distance(0.741 Å)
-----------------------------------------------------------------------
Returns:
(autode.values.Distance): Equlirbium distance
"""
if not self.idxs_are_present(i, j):
raise ValueError(f'Cannot calculate the equilibrium distance '
f'between {i}-{j}. At least one atom not present')

if i == j:
return Distance(0.0, units='Å')

symbols = f"{self[i].atomic_symbol}{self[j].atomic_symbol}"

if symbols in _bond_lengths:
return Distance(_bond_lengths[symbols], units='Å')

# TODO: Something more accurate here
return self[i].covalent_radius + self[j].covalent_radius

def distance(self,
i: int,
j: int) -> Distance:
"""
Distance between two atoms (Å), indexed from 0.
.. code-block:: Python
>>> import autode as ade
>>> mol = ade.Molecule(atoms=[ade.Atom('H'), ade.Atom('H', x=1.0)])
>>> mol.distance(0, 1)
Distance(1.0 Å)
-----------------------------------------------------------------------
Arguments:
i (int): Atom index of the first atom
j (int): Atom index of the second atom
Returns:
(autode.values.Distance): Distance
Raises:
(ValueError):
"""
if not self.idxs_are_present(i, j):
raise ValueError(f'Cannot calculate the distance between {i}-{j}. '
f'At least one atom not present')

return Distance(np.linalg.norm(self[i].coord - self[j].coord))

def vector(self,
i: int,
j: int) -> np.ndarray:
Expand Down Expand Up @@ -843,43 +943,11 @@ def weight(self) -> Mass:

return sum(atom.mass for atom in self.atoms)

def _idxs_are_present(self, *args):
"""
Are a set of indexes present in the collection of atoms?
-----------------------------------------------------------------------
Arguments:
args (int):
Returns:
(bool):
"""
return set(args).issubset(set(range(self.n_atoms)))

def distance(self,
i: int,
j: int) -> Distance:
"""
Distance between two atoms (Å), indexed from 0.
-----------------------------------------------------------------------
Arguments:
i (int): Atom index of the first atom
j (int): Atom index of the second atom
Returns:
(autode.values.Distance): Distance
Raises:
(ValueError):
"""
if not self._idxs_are_present(i, j):
raise ValueError(f'Cannot calculate the distance between {i}-{j}. '
f'At least one atom not present')

value = np.linalg.norm(self.atoms[i].coord - self.atoms[j].coord)
def distance(self, i: int, j: int) -> Distance:
return self.atoms.distance(i, j)

return Distance(value)
def eqm_bond_distance(self, i: int, j: int) -> Distance:
return self.atoms.eqm_bond_distance(i, j)

def angle(self,
i: int,
Expand Down Expand Up @@ -916,7 +984,7 @@ def angle(self,
Raises:
(ValueError): If any of the atom indexes are not present
"""
if not self._idxs_are_present(i, j, k):
if not self.atoms.idxs_are_present(i, j, k):
raise ValueError(f'Cannot calculate the angle between {i}-{j}-{k}.'
f' At least one atom not present')

Expand Down Expand Up @@ -973,7 +1041,7 @@ def dihedral(self,
(ValueError): If any of the atom indexes are not present in the
molecule
"""
if not self._idxs_are_present(w, x, y, z):
if not self.atoms.idxs_are_present(w, x, y, z):
raise ValueError(f'Cannot calculate the dihedral angle involving '
f'atoms {z}-{w}-{x}-{y}. At least one atom not '
f'present')
Expand Down Expand Up @@ -1147,7 +1215,6 @@ def transition_metals(cls, row: int):
'Rh': [0, 1, 2, 3, 4, 5, 6]
}


# Atomic weights in amu from:
# IUPAC-CIAWW's Atomic weights of the elements: Review 2000
atomic_weights = {"H": 1.00794, "He": 4.002602, "Li": 6.941, "Be": 9.012182,
Expand Down Expand Up @@ -1185,7 +1252,6 @@ def transition_metals(cls, row: int):
'Ts': 294.0, 'Og': 294.0
}


# van der Walls radii from https://books.google.no/books?id=bNDMBQAAQBAJ
vdw_radii = {'H': 1.1, 'He': 1.4, 'Li': 1.82, 'Be': 1.53, 'B': 1.92, 'C': 1.7, 'N': 1.55, 'O': 1.52, 'F': 1.47, 'Ne': 1.54, 'Na': 2.27, 'Mg': 1.73, 'Al': 1.84,
'Si': 2.1, 'P': 1.8, 'S': 1.8, 'Cl': 1.75, 'Ar': 1.88, 'K': 2.75, 'Ca': 2.31, 'Sc': 2.15, 'Ti': 2.11, 'V': 2.07, 'Cr': 2.06, 'Mn': 2.05, 'Fe': 2.04,
Expand Down Expand Up @@ -1218,3 +1284,42 @@ def transition_metals(cls, row: int):
'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm',
'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn', 'Nh', 'Fl',
'Mc', 'Lv']

# Covalent radii in picometers from https://en.wikipedia.org/wiki/Covalent_radius
_covalent_radii_pm = [
31., 28.,
128., 96., 84., 76., 71., 66., 57., 58.,
166., 141., 121., 111., 107., 105., 102., 106.,
102., 203., 176., 170., 160., 153., 139., 161., 152., 150., 124., 132., 122., 122., 120., 119., 120., 116.,
220., 195., 190., 175., 164., 154., 147., 146., 142., 139., 145., 144., 142., 139., 139., 138., 139., 140.,
244., 215.,
207., 204., 203., 201., 199., 198., 198., 196., 194., 192., 192., 189., 190., 187.,
175., 187., 170., 162., 151., 144., 141., 136., 136., 132., 145., 146., 148., 140., 150., 150.
]

# Experimental bond lengths from https://cccbdb.nist.gov/diatomicexpbondx.asp
_bond_lengths = {
"HH": 0.741,
"FF": 1.412,
"ClCl": 1.988,
"II": 2.665
}


_max_valances = {'H': 1,
'He': 0,
'B': 4,
'C': 4,
'N': 4,
'O': 3,
'F': 1,
'Si': 4,
'P': 6,
'S': 6,
'Cl': 4,
'Br': 4,
'I': 6,
'Xe': 6,
'Al': 4,
}

Loading

0 comments on commit ce7e8c1

Please sign in to comment.