commit 9cd032cfa1fda70e3593947c5dbf47bc43bd7229
Author: Carles Marti <carles.marti-aliod@ens-lyon.fr>
Date:   Wed Oct 7 09:52:51 2020 +0200

    Implemented the possibility to exclude the molecule adsorption center from the collision detection.

diff --git a/modules/dos_input.py b/modules/dos_input.py
index 461479a..ef6a07e 100644
--- a/modules/dos_input.py
+++ b/modules/dos_input.py
@@ -58,9 +58,11 @@ dos_inp = ConfigParser(inline_comment_prefixes='#',
 new_answers = {'n': False, 'none': False, 'nay': False,
                'y': True, 'sí': True, 'aye': True, 'sure': True}
 for answer, val in new_answers.items():
-    dos_inp.BOOLEAN_STATES[answer] = val
+    dos_inp.BOOLEAN_STATES[answer] = val  # TODO Check value 0
 turn_false_answers = [answer for answer in dos_inp.BOOLEAN_STATES
                       if dos_inp.BOOLEAN_STATES[answer] is False]
+turn_true_answers = [answer for answer in dos_inp.BOOLEAN_STATES
+                     if dos_inp.BOOLEAN_STATES[answer]]
 
 no_sect_err = "Section '%s' not found on input file"
 no_opt_err = "Option '%s' not found on section '%s'"
@@ -637,6 +639,14 @@ def get_min_coll_height(norm_vect):
     return min_coll_height
 
 
+def get_exclude_ads_ctr():
+    err_msg = "exclude_ads_ctr must have a boolean value."
+    exclude_ads_ctr = try_command(dos_inp.getboolean, [(ValueError, err_msg)],
+                                  "Screening", "exclude_ads_ctr",
+                                  fallback=False)
+    return exclude_ads_ctr
+
+
 def get_H_donor(spec_atoms):
     from ase.data import chemical_symbols
     err_msg = "The value of 'h_donor' must be either False, a chemical symbol " \
@@ -836,6 +846,7 @@ def read_input(in_file):
                 and inp_vars['collision_threshold'] is False:
             logger.warning("Collisions are deactivated: Overlapping of "
                            "adsorbate and surface is possible")
+        inp_vars['exclude_ads_ctr'] = get_exclude_ads_ctr()
         inp_vars['h_donor'] = get_H_donor(inp_vars['special_atoms'])
         inp_vars['max_structures'] = get_max_structures()
         inp_vars['use_molec_file'] = get_use_molec_file()
diff --git a/modules/screening.py b/modules/screening.py
index 500c621..21a17de 100644
--- a/modules/screening.py
+++ b/modules/screening.py
@@ -314,7 +314,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=1.2):
+                    nn_molec=0, coll_coeff=1.2, exclude_atom=False):
     """Checks whether a slab and a molecule collide or not.
 
     @param slab_molec: The system of adsorbate-slab for which to detect if there
@@ -328,8 +328,11 @@ def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
     @param min_height: The minimum height atoms can have to not be considered as
         colliding.
     @param vect: The vector perpendicular to the slab.
+    @param exclude_atom: Whether to exclude the adsorption center in the
+        molecule in the collision detection.
     @return: bool, whether the surface and the molecule collide.
     """
+    from copy import deepcopy
     from ase.neighborlist import natural_cutoffs, neighbor_list
 
     # Check structure overlap by height
@@ -342,6 +345,9 @@ def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
             logger.error(err_msg)
             raise ValueError(err_msg)
         for atom in slab_molec[slab_num_atoms:]:
+            if exclude_atom is not False \
+                    and atom.index == exclude_atom:
+                continue
             for i, coord in enumerate(vect):
                 if coord == 0:
                     continue
@@ -350,9 +356,17 @@ def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
 
     # 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 exclude_atom is not False:
+            slab_molec_wo_ctr = deepcopy(slab_molec)
+            del slab_molec_wo_ctr[exclude_atom + slab_num_atoms]
+            slab_molec_cutoffs = natural_cutoffs(slab_molec_wo_ctr,
+                                                 mult=coll_coeff)
+            slab_molec_nghbs = len(neighbor_list("i", slab_molec_wo_ctr,
+                                                 slab_molec_cutoffs))
+        else:
+            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
 
@@ -361,8 +375,8 @@ def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
 
 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
+                 coll_coeff, height=2.5, excl_atom=False):
+    # TODO Rename this function
     """Tries to adsorb a molecule on a slab trying to avoid collisions by doing
     small rotations.
 
@@ -382,6 +396,8 @@ def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
         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.
+    @param excl_atom: Whether to exclude the adsorption center in the
+        molecule in the collision detection.
     @return collision: bool, whether the structure generated has collisions
         between slab and adsorbate.
     """
@@ -401,7 +417,7 @@ def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
                       norm_vect=norm_vect)
         collision = check_collision(slab_molec, slab_num_atoms, min_coll_height,
                                     norm_vect, slab_nghbs, molec_nghbs,
-                                    coll_coeff)
+                                    coll_coeff, excl_atom)
         num_corr += 1
     return slab_molec, collision
 
@@ -506,7 +522,7 @@ def dissociation(slab_molec, h_donor, h_acceptor, num_atoms_slab):
 
 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, height):
+              h_donor, h_acceptor, height, excl_atom):
     """Generates adsorbate-surface structures by sampling over Euler angles.
 
     This function generates a number of adsorbate-surface structures at
@@ -530,6 +546,8 @@ def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
     @param h_donor: List of atom types or atom numbers that are H-donors.
     @param h_acceptor: List of atom types or atom numbers that are H-acceptors.
     @param height: Height on which to try adsorption.
+    @param excl_atom: Whether to exclude the adsorption center in the
+        molecule in the collision detection.
     @return: list of ase.Atoms object conatining all the orientations of a given
         conformer.
     """
@@ -546,12 +564,12 @@ def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
                     continue
                 molec = deepcopy(prealigned_molec)
                 molec.euler_rotate(alpha, beta, gamma, center=ctr_coord)
-                slab_molec, collision = correct_coll(molec, slab,
-                                                     ctr_coord, site_coord,
-                                                     num_pts, min_coll_height,
-                                                     norm_vect,
+                slab_molec, collision = correct_coll(molec, slab, ctr_coord,
+                                                     site_coord, num_pts,
+                                                     min_coll_height, norm_vect,
                                                      slab_nghbs, molec_nghbs,
-                                                     coll_coeff, height)
+                                                     coll_coeff, height,
+                                                     excl_atom)
                 if not collision and not any([np.allclose(slab_molec.positions,
                                                           conf.positions)
                                               for conf in slab_ads_list]):
@@ -832,7 +850,8 @@ def internal_rotate(molecule, surf, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
 
 def ads_internal(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, h_donor, h_acceptor, max_hel, height):
+                 slab_nghbs, molec_nghbs, h_donor, h_acceptor, max_hel, height,
+                 excl_atom):
     """Generates adsorbate-surface structures by sampling over internal angles.
 
     @param orig_molec: ase.Atoms object of the molecule to adsorb (adsorbate).
@@ -861,6 +880,8 @@ def ads_internal(orig_molec, slab, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
     @param max_hel: Maximum value for sampling the helicopter
         (surf-adsorbate dihedral) angle.
     @param height: Height on which to try adsorption.
+    @param excl_atom: Whether to exclude the adsorption center in the
+        molecule in the collision detection.
     @return: list of ase.Atoms object conatining all the orientations of a given
         conformer.
     """
@@ -887,7 +908,7 @@ def ads_internal(orig_molec, slab, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,
                                                      num_pts, min_coll_height,
                                                      norm_vect, slab_nghbs,
                                                      molec_nghbs, coll_coeff,
-                                                     height)
+                                                     height, excl_atom)
                 if not collision and \
                         not any([np.allclose(slab_molec.positions,
                                              conf.positions)
@@ -913,6 +934,7 @@ def adsorb_confs(conf_list, surf, inp_vars):
     @param inp_vars: Calculation parameters from input file.
     @return: list of ase.Atoms for the adsorbate-surface structures
     """
+    from copy import deepcopy
     from ase.neighborlist import natural_cutoffs, neighbor_list
     molec_ctrs = inp_vars['molec_ctrs']
     sites = inp_vars['sites']
@@ -921,6 +943,7 @@ def adsorb_confs(conf_list, surf, inp_vars):
     inp_norm_vect = inp_vars['surf_norm_vect']
     min_coll_height = inp_vars['min_coll_height']
     coll_coeff = inp_vars['collision_threshold']
+    exclude_ads_ctr = inp_vars['exclude_ads_ctr']
     h_donor = inp_vars['h_donor']
     h_acceptor = inp_vars['h_acceptor']
     height = inp_vars['adsorption_height']
@@ -941,11 +964,6 @@ 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 coll_coeff is not False:
-            conf_cutoffs = natural_cutoffs(conf, mult=coll_coeff)
-            molec_nghbs = len(neighbor_list("i", conf, conf_cutoffs))
-        else:
-            molec_nghbs = 0
         for s, site in enumerate(sites_coords):
             if isinstance(inp_norm_vect, str) and inp_norm_vect == 'auto':
                 norm_vect = compute_norm_vect(surf, sites[s],
@@ -953,12 +971,33 @@ def adsorb_confs(conf_list, surf, inp_vars):
             else:
                 norm_vect = inp_norm_vect
             for c, ctr in enumerate(molec_ctr_coords):
+                if exclude_ads_ctr and isinstance(molec_ctrs[c], int):
+                    exclude_atom = molec_ctrs[c]
+                else:
+                    exclude_atom = False
+                    if exclude_ads_ctr and not isinstance(molec_ctrs[c], int):
+                        logger.warning("'exclude_ads_ctr' only works for atomic"
+                                       "centers and not for many-atoms "
+                                       f"barycenters. {molec_ctrs[c]} are not "
+                                       f"going to be excluded from collison.")
+                if coll_coeff and exclude_atom is not False:
+                    conf_wo_ctr = deepcopy(conf)
+                    del conf_wo_ctr[exclude_atom]
+                    conf_cutoffs = natural_cutoffs(conf_wo_ctr, mult=coll_coeff)
+                    molec_nghbs = len(neighbor_list("i", conf_wo_ctr,
+                                                    conf_cutoffs))
+                elif coll_coeff and exclude_atom is False:
+                    conf_cutoffs = natural_cutoffs(conf, mult=coll_coeff)
+                    molec_nghbs = len(neighbor_list("i", conf, conf_cutoffs))
+                else:
+                    molec_nghbs = 0
                 if angles == 'euler':
                     surf_ads_list.extend(ads_euler(conf, surf, ctr, site,
                                                    num_pts, min_coll_height,
                                                    coll_coeff, norm_vect,
                                                    surf_nghbs, molec_nghbs,
-                                                   h_donor, h_acceptor, height))
+                                                   h_donor, h_acceptor, height,
+                                                   exclude_atom))
                 elif angles == 'internal':
                     mol_ctr1 = molec_ctrs[c]
                     mol_ctr2 = inp_vars["molec_ctrs2"][c]
@@ -973,7 +1012,8 @@ def adsorb_confs(conf_list, surf, inp_vars):
                                                       coll_coeff, norm_vect,
                                                       surf_nghbs, molec_nghbs,
                                                       h_donor, h_acceptor,
-                                                      max_h, height))
+                                                      max_h, height,
+                                                      exclude_atom))
     return surf_ads_list
 
 
