commit b4b2f30771e8d92c53fc691c7865a2aff9652d25
Author: Carles Marti <carles.marti-aliod@ens-lyon.fr>
Date:   Tue Jun 30 11:08:10 2020 +0200

    Implemented the possibility to dissociatie H atoms

diff --git a/modules/dos_input.py b/modules/dos_input.py
index 7ca22c4..34f7a42 100644
--- a/modules/dos_input.py
+++ b/modules/dos_input.py
@@ -536,12 +536,19 @@ def get_min_coll_height():
     return min_coll_height
 
 
-def get_try_disso():
+def get_disso_atoms():
     err_msg = "try_disso should be have a boolean value (True or False)"
-    try_disso = try_command(dos_inp.getboolean,
-                            [(ValueError, err_msg)],
-                            'Screening', 'try_disso', fallback=False)
-    return try_disso
+    disso_atoms_str = dos_inp.get('Screening', 'disso_atoms', fallback="False")
+    disso_atoms = []
+    if disso_atoms_str.lower() in turn_false_answers:
+        pass
+    else:
+        for el in disso_atoms_str.split():
+            try:
+                disso_atoms.append(int(el))
+            except ValueError:
+                disso_atoms.append(el)
+    return disso_atoms
 
 
 # Refinement
@@ -675,7 +682,7 @@ def read_input(in_file):
         return_vars['select_magns'] = get_select_magns()
         return_vars['confs_per_magn'] = get_confs_per_magn()
         return_vars['surf_norm_vect'] = get_surf_norm_vect()
-        return_vars['try_disso'] = get_try_disso()
+        return_vars['disso_atoms'] = get_disso_atoms()
         return_vars['ads_algo'] = get_ads_algo()
         return_vars['sample_points_per_angle'] = get_pts_per_angle()
         return_vars['collision_threshold'] = get_coll_thrsld()
diff --git a/modules/screening.py b/modules/screening.py
index aadc66e..1540bf0 100644
--- a/modules/screening.py
+++ b/modules/screening.py
@@ -220,13 +220,95 @@ def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
             return False
     else:
         slab_molec_cutoffs = natural_cutoffs(slab_molec, mult=coll_coeff)
-        slab_molec_nghbs = len(neighbor_list("i", slab_molec, slab_molec_cutoffs))
+        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
 
 
+def dissociate_h(slab_molec_orig, h_idx, num_atoms_slab, neigh_cutoff=1):
+    # TODO rethink
+    """Moves a H atom from a molecule to the surface if the distance is short
+    enough.
+
+    @param slab_molec_orig: The ase.Atoms object of the system adsorbate-slab.
+    @param h_idx: The index of the hydrogen atom to carry out adsorption of.
+    @param num_atoms_slab: The number of atoms of the slab without adsorbate.
+    @param neigh_cutoff: half the maximum distance between the surface and the
+        H for it to carry out dissociation.
+    @return: An ase.Atoms object of the system adsorbate-surface with H
+    """
+    from copy import deepcopy
+    from ase.neighborlist import NeighborList
+    slab_molec = deepcopy(slab_molec_orig)
+    cutoffs = len(slab_molec) * [neigh_cutoff]
+    nl = NeighborList(cutoffs, self_interaction=False, bothways=True)
+    nl.update(slab_molec)
+    surf_h_vect = np.array([np.infty] * 3)
+    for neigh_idx in nl.get_neighbors(h_idx)[0]:
+        if neigh_idx < num_atoms_slab:
+            dist = np.linalg.norm(slab_molec[neigh_idx].position -
+                                  slab_molec[h_idx].position)
+            if dist < np.linalg.norm(surf_h_vect):
+                surf_h_vect = slab_molec[neigh_idx].position \
+                              - slab_molec[h_idx].position
+    if np.linalg.norm(surf_h_vect) != np.infty:
+        trans_vect = surf_h_vect - surf_h_vect / np.linalg.norm(surf_h_vect)
+        slab_molec[h_idx].position = slab_molec[h_idx].position + trans_vect
+        return slab_molec
+
+
+def dissociation(slab_molec, disso_atoms, num_atoms_slab):
+    # TODO rethink
+    # TODO multiple dissociation
+    """Decides which H atoms to dissociate according to a list of atoms.
+
+    Given a list of chemical symbols or atom indices it checks for every atom
+    or any of its neighbor if it's a H and calls dissociate_h to try to carry
+    out dissociation of that H. For atom indices, it checks both whether
+    the atom index or its neighbors are H, for chemical symbols, it only checks
+    if there is a neighbor H.
+    @param slab_molec: The ase.Atoms object of the system adsorbate-slab.
+    @param disso_atoms: The indices or chemical symbols of the atoms
+    @param num_atoms_slab:
+    @return:
+    """
+    from ase.neighborlist import natural_cutoffs, NeighborList
+    molec = slab_molec[num_atoms_slab:]
+    cutoffs = natural_cutoffs(molec)
+    nl = NeighborList(cutoffs, self_interaction=False, bothways=True)
+    nl.update(molec)
+    disso_structs = []
+    for el in disso_atoms:
+        if isinstance(el, int):
+            if molec[el].symbol == 'H':
+                disso_struct = dissociate_h(slab_molec, el + num_atoms_slab,
+                                            num_atoms_slab)
+                if disso_struct is not None:
+                    disso_structs.append(disso_struct)
+            else:
+                for neigh_idx in nl.get_neighbors(el)[0]:
+                    if molec[neigh_idx].symbol == 'H':
+                        disso_struct = dissociate_h(slab_molec, neigh_idx +
+                                                    num_atoms_slab,
+                                                    num_atoms_slab)
+                        if disso_struct is not None:
+                            disso_structs.append(disso_struct)
+        else:
+            for atom in molec:
+                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 \
+                                                        + num_atoms_slab,
+                                                        num_atoms_slab)
+                            if disso_struct is not None:
+                                disso_structs.append(disso_struct)
+    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):
@@ -273,7 +355,8 @@ def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
 
 
 def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
-              min_coll_height, coll_coeff, norm_vect, slab_nghbs, molec_nghbs):
+              min_coll_height, coll_coeff, norm_vect, slab_nghbs, molec_nghbs,
+              disso_atoms):
     """Generates adsorbate-surface structures by sampling over Euler angles.
 
     This function generates a number of adsorbate-surface structures at
@@ -293,6 +376,7 @@ def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
     @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.
     @return: list of ase.Atoms object conatining all the orientations of a given
         conformer
     """
@@ -314,6 +398,8 @@ def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
                                                      coll_coeff)
                 if not collision:
                     slab_ads_list.append(slab_molec)
+                    slab_ads_list.extend(dissociation(slab_molec, disso_atoms,
+                                                      len(slab)))
 
     return slab_ads_list
 
@@ -323,7 +409,7 @@ def ads_chemcat(site, ctr, pts_angle):
 
 
 def adsorb_confs(conf_list, surf, molec_ctrs, sites, algo, num_pts, neigh_ctrs,
-                 norm_vect, min_coll_height, coll_coeff):
+                 norm_vect, min_coll_height, coll_coeff, disso_atoms):
     """Generates a number of adsorbate-surface structure coordinates.
 
     Given a list of conformers, a surface, a list of atom indices (or list of
@@ -343,6 +429,7 @@ def adsorb_confs(conf_list, surf, molec_ctrs, sites, algo, num_pts, neigh_ctrs,
     @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 disso_atoms: List of atom types or atom numbers to try to dissociate.
     @return: list of ase.Atoms for the adsorbate-surface structures
     """
     from ase.neighborlist import natural_cutoffs, neighbor_list
@@ -367,7 +454,8 @@ def adsorb_confs(conf_list, surf, molec_ctrs, sites, algo, num_pts, neigh_ctrs,
                     surf_ads_list.extend(ads_euler(conf, surf, ctr, site,
                                                    num_pts, min_coll_height,
                                                    coll_coeff, norm_vect,
-                                                   surf_nghbs, molec_nghbs))
+                                                   surf_nghbs, molec_nghbs,
+                                                   disso_atoms))
                 elif algo == 'chemcat':
                     surf_ads_list.extend(ads_chemcat(site, ctr, num_pts))
     return surf_ads_list
@@ -401,7 +489,8 @@ def run_screening(inp_vars):
                                  inp_vars['molec_neigh_ctrs'],
                                  inp_vars['surf_norm_vect'],
                                  inp_vars['min_coll_height'],
-                                 inp_vars['collision_threshold'])
+                                 inp_vars['collision_threshold'],
+                                 inp_vars['disso_atoms'])
     logger.info(f'Generated {len(surf_ads_list)} adsorbate-surface atomic '
                 f'configurations, to carry out a calculation of.')
     run_calc('screening', inp_vars, surf_ads_list)
diff --git a/tests/test_dos_input.py b/tests/test_dos_input.py
index 315d1a5..2da6551 100644
--- a/tests/test_dos_input.py
+++ b/tests/test_dos_input.py
@@ -123,7 +123,7 @@ class TestGoodInput(unittest.TestCase):
         self.assertEqual(get_surf_norm_vect(), np.array([1.0, 0.0, 0.0]))
 
     def test_try_disso(self):
-        self.assertEqual(get_try_disso(), True)
+        self.assertEqual(get_disso_atoms(), True)
 
     def test_pts_per_angle(self):
         self.assertEqual(get_pts_per_angle(), 4)
@@ -196,7 +196,7 @@ class TestBadInput(unittest.TestCase):
         self.assertRaises(ValueError, get_surf_norm_vect)
 
     def test_try_disso(self):
-        self.assertRaises(ValueError, get_try_disso)
+        self.assertRaises(ValueError, get_disso_atoms)
 
     def test_pts_per_angle(self):
         self.assertRaises(ValueError, get_pts_per_angle)
