Merged equiv_atoms

.travis.yml

@@ -1,7 +1,6 @@
 branches:
   only:
     - master
-    - equiv_atoms
     - joseph
 dist: xenial
 language: python

README.md

@@ -1,4 +1,4 @@
-[![DOI](https://zenodo.org/badge/196085570.svg)](https://zenodo.org/badge/latestdoi/196085570) [![Build Status](https://travis-ci.org/duartegroup/autodE.svg?branch=joseph)](https://travis-ci.org/duartegroup/autodE) [![codecov](https://codecov.io/gh/duartegroup/autodE/branch/joseph/graph/badge.svg)](https://codecov.io/gh/duartegroup/autodE/branch/joseph) [![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](https://opensource.org/licenses/MIT)
+[![DOI](https://zenodo.org/badge/196085570.svg)](https://zenodo.org/badge/latestdoi/196085570) [![Build Status](https://travis-ci.org/duartegroup/autodE.svg?branch=equiv_atoms)](https://travis-ci.org/duartegroup/autodE) [![codecov](https://codecov.io/gh/duartegroup/autodE/branch/equiv_atoms/graph/badge.svg)](https://codecov.io/gh/duartegroup/autodE/branch/joseph) [![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](https://opensource.org/licenses/MIT)
 
 ![alt text](autode/common/llogo.png)
 ***
@@ -9,12 +9,6 @@ profiles from just SMILES strings of reactant(s) and product(s).
 
 ***
 
-## Installation
-
-If the requirements are already satisfied to install **autodE** as a module
-```
-python setup.py install
-```
 
 ### Dependencies
 * [Python](https://www.python.org/) > v. 3.5
@@ -29,12 +23,19 @@ python setup.py install
 The Python dependencies are listed in requirements.txt best satisfied using a conda install (Miniconda or Anaconda) i.e.
 ```
 conda config --append channels conda-forge
-conda config --append channels conda-forge --append channels omnia
-conda activate autode_env
+conda config --append channels omnia
+conda install --file requirements.txt
 ```
 
 ***
 
+## Installation
+
+Once the requirements are satisfied to install **autodE** as a module
+```
+python setup.py install
+```
+
 ## Usage
 
 Broadly, **autodE** is invoked by first setting appropriate parameters in config.py, or specifying in a python script. 

autode/atoms.py

@@ -1,4 +1,50 @@
 from autode.log import logger
+import numpy as np
+
+
+class Atom:
+
+    def __repr__(self):
+        return f'[{self.label}, {self.coord[0]:.4f}, {self.coord[1]:.4f}, {self.coord[2]:.4f}]'
+
+    def translate(self, vec):
+        """Translate this atom by a vector (np.ndarray, length 3)"""
+        self.coord += vec
+        return None
+
+    def rotate(self, axis, theta, origin=None):
+        """Rotate this atom by theta radians (float) in an axis (np.ndarray, length 3)"""
+
+        # If specified shift so that the origin is at (0, 0, 0), apply the rotation, and shift back
+        if origin is not None:
+            self.translate(vec=-origin)
+
+        # Normalise the axis
+        axis = np.asarray(axis)
+        axis = axis / np.linalg.norm(axis)
+
+        # Compute the 3D rotation matrix using https://en.wikipedia.org/wiki/Euler–Rodrigues_formula
+        a = np.cos(theta / 2.0)
+        b, c, d = -axis * np.sin(theta / 2.0)
+        aa, bb, cc, dd = a * a, b * b, c * c, d * d
+        bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
+        rot_matrix = np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
+                               [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
+                               [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
+
+        # Apply the rotation
+        self.coord = np.matmul(rot_matrix, self.coord)
+
+        if origin is not None:
+            self.translate(vec=origin)
+
+        return None
+
+    def __init__(self, atomic_symbol, x, y, z):
+
+        self.label = atomic_symbol
+        self.coord = np.array([x, y, z])
+
 
 # A set of reasonable valances for anionic/neutral/cationic atoms
 valid_valances = {'H': [0, 1],
@@ -39,6 +85,9 @@
              'Au': 2.14, 'Hg': 2.23, 'Tl': 1.96, 'Pb': 2.02, 'Bi': 2.07, 'Po': 1.97, 'At': 2.02, 'Rn': 2.2, 'Fr': 3.48, 'Ra': 2.83, 'Ac': 2.47, 'Th': 2.45, 'Pa': 2.43,
              'U': 2.41, 'Np': 2.39, 'Pu': 2.43, 'Am': 2.44, 'Cm': 2.45, 'Bk': 2.44, 'Cf': 2.45, 'Es': 2.45, 'Fm': 2.45, 'Md': 2.46, 'No': 2.46, 'Lr': 2.46}
 
+pi_valencies = {'B': [1, 2], 'N': [1, 2], 'O': [1], 'C': [1, 2, 3], 'P': [1, 2, 3, 4], 'S': [1, 3, 4, 5],
+                'Si': [1, 2, 3]}
+
 
 def get_maximal_valance(atom_label):
     """Get the maximum valance of an atom
@@ -86,3 +135,24 @@ def get_vdw_radius(atom_label):
     else:
         logger.error(f'Couldn\'t find the VdV radii for {atom_label}. Guessing at 2.3')
         return 2.3
+
+
+def is_pi_atom(atom_label, valency):
+    """
+    Determine if an atom is a 'π-atom' i.e. is unsaturated and is a first or second row element
+
+    Arguments;
+        atom_label (str):
+        valency (int):
+
+    Returns:
+        (bool)
+    """
+
+    if atom_label not in pi_valencies.keys():
+        return False
+
+    if valency in pi_valencies[atom_label]:
+        return True
+
+    return False

autode/bond_lengths.py

@@ -1,115 +1,42 @@
 import numpy as np
 from autode.log import logger
-from autode.geom import xyz2coord
 
 
-def get_xyz_bond_list(xyzs, relative_tolerance=0.2):
-    """Determine the 'bonds' between atoms defined in a xyzs list.
-
-    Arguments:
-        xyzs (list(list)): e.g. [['C', 0.0, 0.0, 0.0], ...]
-
-    Keyword Arguments:
-        relative_tolerance (float): relative tolerance to consider a bond e.g. max bond length = 1.1 x avg. bond length (default: {0.2})
-
-    Returns:
-        list(tuple): list of bonds given as tuples of atom ids
-    """
-    logger.info(f'Getting bond list from xyzs. Maximum bond is {1 + relative_tolerance}x average')
-
-    bond_list = []
-
-    for i in range(len(xyzs)):
-        i_coords = xyz2coord(xyzs[i])
-
-        for j in range(len(xyzs)):
-            if i > j:
-                j_coords = xyz2coord(xyzs[j])
-                # Calculate the distance between the two points in Å
-                dist = np.linalg.norm(j_coords - i_coords)
-
-                # Get the possible keys e.g. CH and HC that might be in the avg_bond_lengths dictionary
-                atom_i_label, atom_j_label = xyzs[i][0], xyzs[j][0]
-
-                i_j_bond_length = get_avg_bond_length(atom_i_label, atom_j_label)
-
-                if dist < i_j_bond_length * (1.0 + relative_tolerance):
-                    bond_list.append((i, j))
-
-    if len(bond_list) == 0:
-        logger.warning('Bond list is empty')
-
-    return bond_list
-
-
-def get_bond_list_from_rdkit_bonds(rdkit_bonds_obj):
-    """For an RDKit bonds object get the standard xyz_bond_list
-
-    Arguments:
-        rdkit_bonds_obj (rdkit bonds object): rdkit bonds object
-    
-    Returns:
-        list(tuple): list of bonds given as tuples of atom ids
-    """
-    logger.info('Converting RDKit bonds to bond list')
-
-    bond_list = []
-
-    for bond in rdkit_bonds_obj:
-        atom_i = bond.GetBeginAtomIdx()
-        atom_j = bond.GetEndAtomIdx()
-        if atom_i > atom_j:
-            bond_list.append((atom_i, atom_j))
-        else:
-            bond_list.append((atom_j, atom_i))
-
-    if len(bond_list) == 0:
-        logger.warning('Bond list is empty')
-
-    return bond_list
-
-
-def get_ideal_bond_length_matrix(xyzs, bonds):
+def get_ideal_bond_length_matrix(atoms, bonds):
     """Populate a n_atoms x n_atoms matrix of ideal bond lengths. All non-bonded atoms will have zero ideal bond lengths
     
     Arguments:
-        xyzs (list(list)): standard xyzs list
-        bonds (list(tuple)): list of bonds defined by tuples of atom ids
+        atoms (list(autode.atoms.Atom)): List of atoms
+        bonds (list(tuple)): List of bonds defined by tuples of atom ids
     
     Returns:
-        np.array: matrix of idea bond lengths
+        np.array: Matrix of ideal bond lengths. n_atoms x n_atoms
     """
 
     logger.info('Getting ideal bond length matrix')
 
-    n_atoms = len(xyzs)
-    ideal_bondl_matrix = np.zeros((n_atoms, n_atoms))
+    n_atoms = len(atoms)
+    ideal_bond_length_matrix = np.zeros((n_atoms, n_atoms))
 
     for i in range(n_atoms):
         for j in range(n_atoms):
             if (i, j) in bonds or (j, i) in bonds:
-                ideal_bondl_matrix[i, j] = get_avg_bond_length(atom_i_label=xyzs[i][0], atom_j_label=xyzs[j][0])
-
-    return ideal_bondl_matrix
+                ideal_bond_length_matrix[i, j] = get_avg_bond_length(atom_i_label=atoms[i].label,
+                                                                     atom_j_label=atoms[j].label)
+    return ideal_bond_length_matrix
 
 
-def get_avg_bond_length(atom_i_label=None, atom_j_label=None, mol=None, bond=None):
+def get_avg_bond_length(atom_i_label, atom_j_label):
     """Get the average bond length from either a molecule and a bond or two atom labels (e.g. atom_i_label = 'C'
     atom_j_label = 'H')
     
     Keyword Arguments:
-        atom_i_label (str): atom label e.g 'C' (default: {None})
-        atom_j_label (str): atom label e.g 'C' (default: {None})
-        mol (molecule obj): molecule object
-        bond (tuple): bond defined by the two atom ids
+        atom_i_label (str): atom label e.g 'C'
+        atom_j_label (str): atom label e.g 'C'
     
     Returns:
         float: avg bond length of the bond
     """
-    if mol is not None and bond is not None:
-        atom_i_label = mol.get_atom_label(bond[0])
-        atom_j_label = mol.get_atom_label(bond[1])
-
     key1, key2 = atom_i_label + atom_j_label, atom_j_label + atom_i_label
 
     if key1 in avg_bond_lengths.keys():