commit e8bebccafd75428059af6e3eea78e4394022a5bd
Author: Carles Marti <carles.marti-aliod@ens-lyon.fr>
Date:   Thu Sep 10 11:25:41 2020 +0200

    Implemented the possibility to check adsorbate-surface overlap with both methods at the same time, height and sphere collision. Changed function order.

diff --git a/modules/screening.py b/modules/screening.py
index ab1c2f2..da82d10 100644
--- a/modules/screening.py
+++ b/modules/screening.py
@@ -198,7 +198,7 @@ def add_adsorbate(slab, adsorbate, site_coord, ctr_coord, height, offset=None,
 
 
 def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
-                    nn_molec=0, coll_coeff=0.9):
+                    nn_molec=0, coll_coeff=1.2):
     """Checks whether a slab and a molecule collide or not.
 
     @param slab_molec: The system of adsorbate-slab for which to detect if there
@@ -213,27 +213,81 @@ def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
         colliding.
     @param vect: The vector perpendicular to the slab.
     @return: bool, whether the surface and the molecule collide.
-    """  # TODO Enable both checks: 1st height, 2nd sphere collision
+    """
     from ase.neighborlist import natural_cutoffs, neighbor_list
+
+    # Check structure overlap by height
     if min_height is not False:
         cart_axes = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0],
                      [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]]
-        if vect.tolist() in cart_axes:
-            for atom in slab_molec[slab_num_atoms:]:
-                for i, coord in enumerate(vect):
-                    if coord == 0:
-                        continue
-                    if atom.position[i] * coord < min_height:
-                        return True
-            return False
-    else:
+        if vect.tolist() not in cart_axes:
+            err_msg = "'min_coll_height' option is only implemented for " \
+                      "'surf_norm_vect' to be one of the x, y or z axes. "
+            logger.error(err_msg)
+            raise ValueError(err_msg)
+        for atom in slab_molec[slab_num_atoms:]:
+            for i, coord in enumerate(vect):
+                if coord == 0:
+                    continue
+                if atom.position[i] * coord < min_height * coord:
+                    return True
+
+    # Check structure overlap by sphere collision
+    if coll_coeff is not False:
         slab_molec_cutoffs = natural_cutoffs(slab_molec, mult=coll_coeff)
         slab_molec_nghbs = len(
             neighbor_list("i", slab_molec, slab_molec_cutoffs))
         if slab_molec_nghbs > nn_slab + nn_molec:
             return True
-        else:
-            return False
+
+    return False
+
+
+def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
+                 min_coll_height, norm_vect, slab_nghbs, molec_nghbs,
+                 coll_coeff, height=2.5):
+    # TODO Rethink this function
+    """Tries to adsorb a molecule on a slab trying to avoid collisions by doing
+    small rotations.
+
+    @param molec: ase.Atoms object of the molecule to adsorb
+    @param slab: ase.Atoms object of the surface on which to adsorb the
+        molecule
+    @param ctr_coord: The coordinates of the molecule to use as adsorption
+        center.
+    @param site_coord: The coordinates of the surface on which to adsorb the
+        molecule
+    @param num_pts: Number on which to sample Euler angles.
+    @param min_coll_height: The lowermost height for which to detect a collision
+    @param norm_vect: The vector perpendicular to the surface.
+    @param slab_nghbs: Number of neigbors in the surface.
+    @param molec_nghbs: Number of neighbors in the molecule.
+    @param coll_coeff: The coefficient that multiplies the covalent radius of
+        atoms resulting in a distance that two atoms being closer to that is
+        considered as atomic collision.
+    @param height: Height on which to try adsorption
+    @return collision: bool, whether the structure generated has collisions
+        between slab and adsorbate.
+    """
+    from copy import deepcopy
+    slab_num_atoms = len(slab)
+    slab_molec = []
+    collision = True
+    max_corr = 6  # Should be an even number
+    d_angle = 180 / ((max_corr / 2.0) * num_pts)
+    num_corr = 0
+    while collision and num_corr <= max_corr:
+        k = num_corr * (-1) ** num_corr
+        slab_molec = deepcopy(slab)
+        molec.euler_rotate(k * d_angle, k * d_angle / 2, k * d_angle,
+                           center=ctr_coord)
+        add_adsorbate(slab_molec, molec, site_coord, ctr_coord, height,
+                      norm_vect=norm_vect)
+        collision = check_collision(slab_molec, slab_num_atoms, min_coll_height,
+                                    norm_vect, slab_nghbs, molec_nghbs,
+                                    coll_coeff)
+        num_corr += 1
+    return slab_molec, collision
 
 
 def dissociate_h(slab_molec_orig, h_idx, num_atoms_slab, h_acceptor,
@@ -334,53 +388,6 @@ def dissociation(slab_molec, h_donor, h_acceptor, num_atoms_slab):
     return disso_structs
 
 
-def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
-                 min_coll_height, norm_vect, slab_nghbs, molec_nghbs,
-                 coll_coeff, height=2.5):
-    # TODO Rethink this function
-    """Tries to adsorb a molecule on a slab trying to avoid collisions by doing
-    small rotations.
-
-    @param molec: ase.Atoms object of the molecule to adsorb
-    @param slab: ase.Atoms object of the surface on which to adsorb the
-        molecule
-    @param ctr_coord: The coordinates of the molecule to use as adsorption
-        center.
-    @param site_coord: The coordinates of the surface on which to adsorb the
-        molecule
-    @param num_pts: Number on which to sample Euler angles.
-    @param min_coll_height: The lowermost height for which to detect a collision
-    @param norm_vect: The vector perpendicular to the surface.
-    @param slab_nghbs: Number of neigbors in the surface.
-    @param molec_nghbs: Number of neighbors in the molecule.
-    @param coll_coeff: The coefficient that multiplies the covalent radius of
-        atoms resulting in a distance that two atoms being closer to that is
-        considered as atomic collision.
-    @param height: Height on which to try adsorption
-    @return collision: bool, whether the structure generated has collisions
-        between slab and adsorbate.
-    """
-    from copy import deepcopy
-    slab_num_atoms = len(slab)
-    slab_molec = []
-    collision = True
-    max_corr = 6  # Should be an even number
-    d_angle = 180 / ((max_corr / 2.0) * num_pts)
-    num_corr = 0
-    while collision and num_corr <= max_corr:
-        k = num_corr * (-1) ** num_corr
-        slab_molec = deepcopy(slab)
-        molec.euler_rotate(k * d_angle, k * d_angle / 2, k * d_angle,
-                           center=ctr_coord)
-        add_adsorbate(slab_molec, molec, site_coord, ctr_coord, height,
-                      norm_vect=norm_vect)
-        collision = check_collision(slab_molec, slab_num_atoms, min_coll_height,
-                                    norm_vect, slab_nghbs, molec_nghbs,
-                                    coll_coeff)
-        num_corr += 1
-    return slab_molec, collision
-
-
 def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
               min_coll_height, coll_coeff, norm_vect, slab_nghbs, molec_nghbs,
               h_donor, h_acceptor):
@@ -770,8 +777,7 @@ def ads_chemcat(orig_molec, slab, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
                                                      ctr_coords, site_coords,
                                                      num_pts, min_coll_height,
                                                      norm_vect, slab_nghbs,
-                                                     molec_nghbs, coll_coeff,
-                                                     2.5)
+                                                     molec_nghbs, coll_coeff)
                 if not collision and \
                         not any([np.allclose(slab_molec.positions,
                                              conf.positions)
@@ -814,7 +820,7 @@ def adsorb_confs(conf_list, surf, inp_vars):
 
     surf_ads_list = []
     sites_coords = get_atom_coords(surf, sites)
-    if min_coll_height is False:
+    if coll_coeff is not False:
         surf_cutoffs = natural_cutoffs(surf, mult=coll_coeff)
         surf_nghbs = len(neighbor_list("i", surf, surf_cutoffs))
     else:
@@ -824,7 +830,7 @@ def adsorb_confs(conf_list, surf, inp_vars):
         if inp_vars['pbc_cell'] is not False:
             conf.set_pbc(True)
             conf.set_cell(inp_vars['pbc_cell'])
-        if min_coll_height is False:
+        if coll_coeff is not False:
             conf_cutoffs = natural_cutoffs(conf, mult=coll_coeff)
             molec_nghbs = len(neighbor_list("i", conf, conf_cutoffs))
         else:
