Adaptive path default (duartegroup#28)

* Improved benchmark set

* commit NCI FF example

* Simple adaptive k

* Remove gets

* CI-NEB

* NEB partitioning from initial path

* Various CI improvements

* Skip NEB without maximum

* Consider all rings in constraints

* Variable wait iteration

* Set final distance from product if available

* More accurate interpolation

* Reset bond in ring deletion

* Default step sizes

* Adaptive division of path to find peaks

* Replace linear interpolation with const opt

* Even less stringent imaginary frequency threshold

* Fix for in place bond modification

* NEB and doc improvements

* Started adaptive initial path

* Initial implementation

* XTB grad without constraints

* Heuristic for small bond ring TSs

* Adaptive path testing

* Early stopping criteria

* Documentation update

* Special case for 3 active bonds forming a ring

* Fix for Gaussian gradient extraction

* JS changes

* Remove unused import

* Fix for MOPAC gradients

autode/__init__.py

@@ -11,9 +11,12 @@
 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
 
-__version__ = '1.0.0b2'
+__version__ = '1.0.0b3'
 
 
 __all__ = [
@@ -30,7 +33,10 @@
     'KcalMol',
     'Calculation',
     'KjMol',
+    'pes',
+    'neb',
     'geom',
     'methods',
+    'mol_graphs',
     'utils'
 ]

autode/bond_rearrangement.py

@@ -1,16 +1,14 @@
 import itertools
 import os
 from autode.atoms import get_maximal_valance
-from autode.geom import get_neighbour_list
-from autode.geom import get_points_on_sphere
+from autode.geom import get_neighbour_list, get_points_on_sphere
 from autode.log import logger
-from autode.mol_graphs import get_bond_type_list
-from autode.mol_graphs import get_fbonds
-from autode.mol_graphs import is_isomorphic
-from autode.mol_graphs import connected_components
+from autode.mol_graphs import (get_bond_type_list, get_fbonds, is_isomorphic,
+                               connected_components, find_cycles)
+from autode.config import Config
 
 
-def get_bond_rearrangs(reactant, product, name):
+def get_bond_rearrangs(reactant, product, name, save=True):
     """For a reactant and product (complex) find the set of breaking and
     forming bonds that will turn reactants into products. This works by
     determining the types of bonds that have been made/broken (i.e CH) and
@@ -21,6 +19,9 @@ def get_bond_rearrangs(reactant, product, name):
         product (autode.complex.ProductComplex):
         name (str):
 
+    Keyword Arguments:
+        save (bool): Save bond rearrangements to a file for fast reloading
+
     Returns:
         (list(autode.bond_rearrangements.BondRearrangement)):
     """
@@ -110,10 +111,12 @@ def get_bond_rearrangs(reactant, product, name):
             if n_bond_rearrangs > 1:
                 logger.info(f'Multiple *{n_bond_rearrangs}* possible bond '
                             f'breaking/makings are possible')
-                possible_brs = strip_equiv_bond_rearrs(reactant, possible_brs)
+                possible_brs = strip_equiv_bond_rearrs(possible_brs, reactant)
+                prune_small_ring_rearrs(possible_brs, reactant)
 
-            save_bond_rearrangs_to_file(possible_brs,
-                                        filename=f'{name}_bond_rearrangs.txt')
+            if save:
+                save_bond_rearrangs_to_file(possible_brs,
+                                            filename=f'{name}_BRs.txt')
 
             logger.info(f'Found *{len(possible_brs)}* bond '
                         f'rearrangement(s) that lead to products')
@@ -164,22 +167,22 @@ def get_bond_rearrangs_from_file(filename='bond_rearrangs.txt'):
 
     with open(filename, 'r') as br_file:
         fbonds_block = False
-        bbonds_block = True
-        fbonds = []
-        bbonds = []
+        fbonds, bbonds = [], []
         for line in br_file:
             if 'fbonds' in line:
                 fbonds_block = True
-                bbonds_block = False
+
             if 'bbonds' in line:
                 fbonds_block = False
-                bbonds_block = True
-            if fbonds_block and len(line.split()) == 2:
-                atom_id_string = line.split()
-                fbonds.append((int(atom_id_string[0]), int(atom_id_string[1])))
-            if bbonds_block and len(line.split()) == 2:
-                atom_id_string = line.split()
-                bbonds.append((int(atom_id_string[0]), int(atom_id_string[1])))
+
+            if len(line.split()) == 2:
+                atom_idx0, atom_idx1 = (int(val) for val in line.split())
+
+                if fbonds_block:
+                    fbonds.append((atom_idx0, atom_idx1))
+                if not fbonds_block:
+                    bbonds.append((atom_idx0, atom_idx1))
+
             if 'end' in line:
                 bond_rearrangs.append(BondRearrangement(forming_bonds=fbonds,
                                                         breaking_bonds=bbonds))
@@ -482,13 +485,13 @@ def get_fbonds_bbonds_2b2f(reac, prod, possible_brs, all_possible_bbonds, all_po
     return possible_brs
 
 
-def strip_equiv_bond_rearrs(mol, possible_bond_rearrs, depth=6):
-    """Remove any bond rearrangement from possible_bond_rearrs for which
+def strip_equiv_bond_rearrs(possible_brs, mol, depth=6):
+    """Remove any bond rearrangement from possible_brs for which
     there is already an equivalent in the unique_bond_rearrangements list
 
     Arguments:
+        possible_brs (list(BondRearrangement)):
         mol (molecule object): reactant object
-        possible_bond_rearrs (list(object)): list of BondRearrangement objects
 
     Keyword Arguments:
         depth (int): Depth of neighbour list that must be identical for a set
@@ -502,7 +505,7 @@ def strip_equiv_bond_rearrs(mol, possible_bond_rearrs, depth=6):
 
     unique_bond_rearrs = []
 
-    for bond_rearr in possible_bond_rearrs:
+    for bond_rearr in possible_brs:
         bond_rearrang_is_unique = True
 
         # Compare bond_rearrang to all those already considered to be unique,
@@ -515,12 +518,72 @@ def strip_equiv_bond_rearrs(mol, possible_bond_rearrs, depth=6):
         if bond_rearrang_is_unique:
             unique_bond_rearrs.append(bond_rearr)
 
-    logger.info(
-        f'Stripped {len(possible_bond_rearrs)-len(unique_bond_rearrs)} '
-        f'bond rearrangements')
+    logger.info(f'Stripped {len(possible_brs) - len(unique_bond_rearrs)} '
+                'bond rearrangements')
     return unique_bond_rearrs
 
 
+def prune_small_ring_rearrs(possible_brs, mol):
+    """
+    Remove any bond rearrangements that go via small (3, 4) rings if there is
+    an alternative that goes vie
+
+    Arguments:
+        possible_brs (list(BondRearrangement)):
+        mol (molecule object): reactant object
+    """
+    small_ring_sizes = (3, 4)
+
+    if not Config.skip_small_ring_tss:
+        logger.info('Not pruning small ring TSs')
+        return None
+
+    # Membered-ness of rings in each bond rearrangement
+    n_mem_rings = [br.n_membered_rings(mol) for br in possible_brs]
+
+    # Unique elements involved in each bond rearrangement
+    elems = [set(mol.atoms[i].label for i in range(mol.n_atoms)
+                 if i in br.active_atoms) for br in possible_brs]
+
+    logger.info(f'Pruning {len(possible_brs)} to remove any '
+                f'{small_ring_sizes}-membered rings where others are possible')
+
+    excluded_idxs = []
+    for i, br in enumerate(possible_brs):
+        logger.info(f'Checking bond rearrangement {i} with rings:'
+                    f' {n_mem_rings[i]} and atom indexes: {br}')
+
+        # Only consider brs with at least one small ring
+        if not any(n_mem in small_ring_sizes for n_mem in n_mem_rings[i]):
+            continue
+
+        # Check against all other rearrangements
+        for j, other_br in enumerate(possible_brs):
+
+            # Only consider brs with the same set of elements
+            if elems[i] != elems[j]:
+                continue
+
+            # Needs to have the same number of rings
+            if len(n_mem_rings[i]) != len(n_mem_rings[j]):
+                continue
+
+            # Exclude i if j has a larger smallest ring size
+            if min(n_mem_rings[i]) < min(n_mem_rings[j]):
+                excluded_idxs.append(i)
+                break
+
+    logger.info(f'Excluding {len(excluded_idxs)} bond rearrangements based on '
+                f'small rings')
+
+    # Delete the excluded bond rearrangements (sorted high ->  low, so the
+    # idxs remain the same while deleting)
+    for idx in sorted(excluded_idxs, reverse=True):
+        del possible_brs[idx]
+
+    return None
+
+
 class BondRearrangement:
 
     def __eq__(self, other):
@@ -568,16 +631,53 @@ def nl(atom):
 
         return self.active_atom_nl
 
-    def _set_active_atom_list(self, bonds, ls):
+    def n_membered_rings(self, mol):
+        """
+        Find the membered-ness of the rings involved in this bond rearrangement
+        will add the forming bonds to the graph to determine
 
-        for bond in bonds:
-            for atom in bond:
-                if atom not in ls:
-                    ls.append(atom)
-                if atom not in self.active_atoms:
-                    self.active_atoms.append(atom)
+        Arguments:
+            (autode.species.Species):
 
-        return None
+        Returns:
+            (list(int)):
+        """
+        assert mol.graph is not None
+        graph = mol.graph.copy()
+
+        for fbond in self.fbonds:
+            if fbond not in graph.edges:
+                graph.add_edge(*fbond)
+
+        rings = find_cycles(graph)
+        n_mem_rings = []
+
+        # Full enumeration over all atoms and rings - could be faster..
+        for ring in rings:
+            for atom_idx in self.active_atoms:
+                if atom_idx in ring:
+                    # This ring has at least one active atom in
+                    n_mem_rings.append(len(ring))
+
+                    # don't add the same ring more than once
+                    break
+
+        return n_mem_rings
+
+    @property
+    def fatoms(self):
+        """Unique atoms indexes involved in forming bonds"""
+        return list(set([i for bond in self.fbonds for i in bond]))
+
+    @property
+    def batoms(self):
+        """Unique atoms indexes involved in breaking bonds"""
+        return list(set([i for bond in self.bbonds for i in bond]))
+
+    @property
+    def active_atoms(self):
+        """Unique atom indexes in forming or breaking bonds"""
+        return list(set(self.fatoms + self.batoms))
 
     @property
     def n_fbonds(self):
@@ -602,11 +702,5 @@ def __init__(self, forming_bonds=None, breaking_bonds=None):
         self.fbonds = forming_bonds if forming_bonds is not None else []
         self.bbonds = breaking_bonds if breaking_bonds is not None else []
 
-        self.active_atoms = []
-        self.fatoms = []
-        self.batoms = []
         self.active_atom_nl = None
         self.all = self.fbonds + self.bbonds
-
-        self._set_active_atom_list(self.fbonds, self.fatoms)
-        self._set_active_atom_list(self.bbonds, self.batoms)

autode/bond_lengths.py → autode/bonds.py

@@ -3,6 +3,83 @@
 from autode.log import logger
 
 
+class ScannedBond:
+
+    def __str__(self):
+        i, j = self.atom_indexes
+        return f'{i}-{j}'
+
+    def __getitem__(self, item):
+        return self.atom_indexes[item]
+
+    @property
+    def dr(self):
+        """Change in distance for this bond (∆r / Å)"""
+        if self.curr_dist is None or self.final_dist is None:
+            return 0
+
+        return self.final_dist - self.curr_dist
+
+    def __init__(self, atom_indexes):
+        """
+        Bond with a current and final distance which will be scanned over
+
+        Arguments:
+            atom_indexes (tuple(int)): Atom indexes that make this
+            'bond' e.g. (0, 1)
+        """
+        assert len(atom_indexes) == 2
+
+        self.atom_indexes = atom_indexes
+
+        self.curr_dist = None
+        self.final_dist = None
+
+        self.forming = False
+        self.breaking = False
+
+
+class FormingBond(ScannedBond):
+
+    def __init__(self, atom_indexes, species, final_species=None):
+        """
+        Forming bond with current and final distances
+
+        Arguments:
+            atom_indexes (tuple(int)):
+            species (autode.species.Species):
+        """
+        super().__init__(atom_indexes)
+        self.forming = True
+
+        i, j = self.atom_indexes
+        self.curr_dist = species.distance(i=i, j=j)
+
+        if final_species is None:
+            self.final_dist = get_avg_bond_length(species.atoms[i].label,
+                                                  species.atoms[j].label)
+        else:
+            self.final_dist = final_species.distance(*atom_indexes)
+
+
+class BreakingBond(ScannedBond):
+
+    def __init__(self, atom_indexes, species):
+        """
+        Form a breaking bond with current and final distances
+
+        Arguments:
+            atom_indexes (tuple(int)):
+            species (autode.species.Species):
+            reaction (autode.reaction.Reaction):
+        """
+        super().__init__(atom_indexes)
+        self.breaking = True
+
+        self.curr_dist = species.distance(*self.atom_indexes)
+        self.final_dist = 2.0 * self.curr_dist
+
+
 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