commit 5261a07fc5a73321162f95bc63d82f3045c7a21d
Author: Carles Marti <carles.marti-aliod@ens-lyon.fr>
Date:   Wed Sep 9 10:27:17 2020 +0200

    Adapted error messages, reformat code for readability, removed TODO, refactor variables, updated docstrings.

diff --git a/modules/screening.py b/modules/screening.py
index 44f0899..39b2ce8 100644
--- a/modules/screening.py
+++ b/modules/screening.py
@@ -226,7 +226,7 @@ def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
                     if atom.position[i] * coord < min_height:
                         return True
             return False
-    else:  # TODO Implement PBC for sphere collision
+    else:
         slab_molec_cutoffs = natural_cutoffs(slab_molec, mult=coll_coeff)
         slab_molec_nghbs = len(
             neighbor_list("i", slab_molec, slab_molec_cutoffs))
@@ -310,7 +310,7 @@ def dissociation(slab_molec, disso_atoms, num_atoms_slab):
                 if atom.symbol.lower() == el.lower():
                     for neigh_idx in nl.get_neighbors(atom.index)[0]:
                         if molec[neigh_idx].symbol == 'H':
-                            disso_struct = dissociate_h(slab_molec, neigh_idx \
+                            disso_struct = dissociate_h(slab_molec, neigh_idx
                                                         + num_atoms_slab,
                                                         num_atoms_slab)
                             if disso_struct is not None:
@@ -346,6 +346,7 @@ def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
     """
     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)
@@ -372,14 +373,15 @@ def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
     This function generates a number of adsorbate-surface structures at
     different orientations of the adsorbate sampled at multiple Euler (zxz)
     angles.
-    @param orig_molec: ase.Atoms object of the molecule to adsorb
-    @param slab: ase.Atoms object of the surface on which to adsorb the molecule
+    @param orig_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
+        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 min_coll_height: The lowest height for which to detect a collision.
     @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.
@@ -388,7 +390,7 @@ def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
     @param molec_nghbs: Number of neighbors in the molecule.
     @param disso_atoms: List of atom types or atom numbers to try to dissociate.
     @return: list of ase.Atoms object conatining all the orientations of a given
-        conformer
+        conformer.
     """
     from copy import deepcopy
     slab_ads_list = []
@@ -406,9 +408,9 @@ def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
                                                      norm_vect,
                                                      slab_nghbs, molec_nghbs,
                                                      coll_coeff)
-                if not collision and \
-                        not any([(slab_molec.positions == conf.positions).all()
-                                 for conf in slab_ads_list]):
+                if not collision and not any([np.allclose(slab_molec.positions,
+                                                          conf.positions)
+                                              for conf in slab_ads_list]):
                     slab_ads_list.append(slab_molec)
                     slab_ads_list.extend(dissociation(slab_molec, disso_atoms,
                                                       len(slab)))
@@ -619,9 +621,10 @@ def chemcat_rotate(molecule, surf, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
         # Check bond rotation is unmodified
         bond_angle = get_vect_angle(-vect_inter, vect_mol)
         if np.isclose((dihedral_angle - dihedral_angle_target + 90) % 180 - 90,
-                      0, atol=1e-3) and np.isclose(
-            (bond_angle - bond_angle_target + 90) %
-            180 - 90, 0, atol=1e-5) and np.allclose(get_atom_coords(molecule,
+                      0, atol=1e-3) \
+                and np.isclose((bond_angle - bond_angle_target + 90) %
+                               180 - 90, 0, atol=1e-5) \
+                and np.allclose(get_atom_coords(molecule,
                                                                     ctr1_mol)
                                                     - get_atom_coords(surf,
                                                                       ctr1_surf),
@@ -645,9 +648,9 @@ def chemcat_rotate(molecule, surf, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
         # vect_mol)
         vect_mol_inner = np.dot(vect_mol, vect_mol)
         bond_inter_reject = -vect_inter - np.dot(-vect_inter, vect_mol) / \
-                            vect_mol_inner * vect_mol
+            vect_mol_inner * vect_mol
         bond2_mol_reject = vect2_mol - np.dot(vect2_mol, vect_mol) / \
-                           vect_mol_inner * vect_mol
+            vect_mol_inner * vect_mol
         dihedral_angle_ini = get_vect_angle(bond_inter_reject,
                                             bond2_mol_reject, vect_mol)
 
@@ -657,23 +660,23 @@ def chemcat_rotate(molecule, surf, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
 
         # Update molecular bonds
         vect_mol = get_atom_coords(molecule, ctr2_mol) \
-                   - get_atom_coords(molecule, ctr1_mol)
+            - get_atom_coords(molecule, ctr1_mol)
         vect2_mol = get_atom_coords(molecule, ctr3_mol) \
-                    - get_atom_coords(molecule, ctr2_mol)
+            - get_atom_coords(molecule, ctr2_mol)
 
         # Check if rotation was successful
         # Check adsorbate dihedral rotation
         vect_mol_inner = np.dot(vect_mol, vect_mol)
         bond2_mol_reject = vect2_mol - np.dot(vect2_mol, vect_mol) / \
-                           vect_mol_inner * vect_mol
+            vect_mol_inner * vect_mol
         mol_dihedral_angle = get_vect_angle(bond_inter_reject,
                                             bond2_mol_reject, vect_mol)
         # Check dihedral rotation
         vect_inter_inner = np.dot(vect_inter, vect_inter)
         vect_surf_reject = vect_surf - np.dot(vect_surf, vect_inter) / \
-                           vect_inter_inner * vect_inter
+            vect_inter_inner * vect_inter
         vect_mol_reject = vect_mol - np.dot(vect_mol, vect_inter) / \
-                          vect_inter_inner * vect_inter
+            vect_inter_inner * vect_inter
         dihedral_angle = get_vect_angle(vect_surf_reject, vect_mol_reject,
                                         vect_inter)
         # Check bond rotation is unmodified
@@ -702,12 +705,41 @@ def chemcat_rotate(molecule, surf, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
 def ads_chemcat(orig_molec, slab, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
                 ctr2_surf, num_pts, min_coll_height, coll_coeff, norm_vect,
                 slab_nghbs, molec_nghbs, disso_atoms, max_h):
+    """Generates adsorbate-surface structures by sampling over chemcat angles.
+
+    @param orig_molec: ase.Atoms object of the molecule to adsorb (adsorbate).
+    @param slab: ase.Atoms object of the surface on which to adsorb the molecule
+    @param ctr1_mol: The index/es of the center in the adsorbate to use as
+        adsorption center.
+    @param ctr2_mol: The index/es of the center in the adsorbate to use for the
+        definition of the surf-adsorbate angle, surf-adsorbate dihedral angle
+        and adsorbate dihedral angle.
+    @param ctr3_mol: The index/es of the center in the adsorbate to use for the
+        definition of the adsorbate dihedral angle.
+    @param ctr1_surf: The index/es of the center in the surface to use as
+        adsorption center.
+    @param ctr2_surf: The index/es of the center in the surface to use for the
+        definition of the surf-adsorbate dihedral angle.
+    @param num_pts: Number on which to sample Euler angles.
+    @param min_coll_height: The lowest height for which to detect a collision
+    @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 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 disso_atoms: List of atom types or atom numbers to try to dissociate.
+    @param max_h: Maximum value for sampling the helicopter
+        (surf-adsorbate dihedral) angle.
+    @return: list of ase.Atoms object conatining all the orientations of a given
+        conformer.
+    """
     from copy import deepcopy
     slab_ads_list = []
     # Rotation over bond angle
     for alpha in np.arange(90, 180, 90 / min(num_pts, 2)):
         # Rotation over surf-adsorbate dihedral angle
-        for beta in np.arange(0, max_h, max_h / num_pts):  # Maybe 360
+        for beta in np.arange(0, max_h, max_h / num_pts):
             # Rotation over adsorbate bond dihedral angle
             for gamma in np.arange(0, 360, 360 / num_pts):
                 new_molec = deepcopy(orig_molec)
@@ -723,7 +755,8 @@ def ads_chemcat(orig_molec, slab, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
                                                      molec_nghbs, coll_coeff,
                                                      2.5)
                 if not collision and \
-                        not any([(slab_molec.positions == conf.positions).all()
+                        not any([np.allclose(slab_molec.positions,
+                                             conf.positions)
                                  for conf in slab_ads_list]):
                     slab_ads_list.append(slab_molec)
                     slab_ads_list.extend(dissociation(slab_molec, disso_atoms,
