Révision b4b2f307 modules/screening.py

b/modules/screening.py
220 220
            return False
221 221
    else:
222 222
        slab_molec_cutoffs = natural_cutoffs(slab_molec, mult=coll_coeff)
223
        slab_molec_nghbs = len(neighbor_list("i", slab_molec, slab_molec_cutoffs))
223
        slab_molec_nghbs = len(
224
            neighbor_list("i", slab_molec, slab_molec_cutoffs))
224 225
        if slab_molec_nghbs > nn_slab + nn_molec:
225 226
            return True
226 227
        else:
227 228
            return False
228 229

  
229 230

  
231
def dissociate_h(slab_molec_orig, h_idx, num_atoms_slab, neigh_cutoff=1):
232
    # TODO rethink
233
    """Moves a H atom from a molecule to the surface if the distance is short
234
    enough.
235

  
236
    @param slab_molec_orig: The ase.Atoms object of the system adsorbate-slab.
237
    @param h_idx: The index of the hydrogen atom to carry out adsorption of.
238
    @param num_atoms_slab: The number of atoms of the slab without adsorbate.
239
    @param neigh_cutoff: half the maximum distance between the surface and the
240
        H for it to carry out dissociation.
241
    @return: An ase.Atoms object of the system adsorbate-surface with H
242
    """
243
    from copy import deepcopy
244
    from ase.neighborlist import NeighborList
245
    slab_molec = deepcopy(slab_molec_orig)
246
    cutoffs = len(slab_molec) * [neigh_cutoff]
247
    nl = NeighborList(cutoffs, self_interaction=False, bothways=True)
248
    nl.update(slab_molec)
249
    surf_h_vect = np.array([np.infty] * 3)
250
    for neigh_idx in nl.get_neighbors(h_idx)[0]:
251
        if neigh_idx < num_atoms_slab:
252
            dist = np.linalg.norm(slab_molec[neigh_idx].position -
253
                                  slab_molec[h_idx].position)
254
            if dist < np.linalg.norm(surf_h_vect):
255
                surf_h_vect = slab_molec[neigh_idx].position \
256
                              - slab_molec[h_idx].position
257
    if np.linalg.norm(surf_h_vect) != np.infty:
258
        trans_vect = surf_h_vect - surf_h_vect / np.linalg.norm(surf_h_vect)
259
        slab_molec[h_idx].position = slab_molec[h_idx].position + trans_vect
260
        return slab_molec
261

  
262

  
263
def dissociation(slab_molec, disso_atoms, num_atoms_slab):
264
    # TODO rethink
265
    # TODO multiple dissociation
266
    """Decides which H atoms to dissociate according to a list of atoms.
267

  
268
    Given a list of chemical symbols or atom indices it checks for every atom
269
    or any of its neighbor if it's a H and calls dissociate_h to try to carry
270
    out dissociation of that H. For atom indices, it checks both whether
271
    the atom index or its neighbors are H, for chemical symbols, it only checks
272
    if there is a neighbor H.
273
    @param slab_molec: The ase.Atoms object of the system adsorbate-slab.
274
    @param disso_atoms: The indices or chemical symbols of the atoms
275
    @param num_atoms_slab:
276
    @return:
277
    """
278
    from ase.neighborlist import natural_cutoffs, NeighborList
279
    molec = slab_molec[num_atoms_slab:]
280
    cutoffs = natural_cutoffs(molec)
281
    nl = NeighborList(cutoffs, self_interaction=False, bothways=True)
282
    nl.update(molec)
283
    disso_structs = []
284
    for el in disso_atoms:
285
        if isinstance(el, int):
286
            if molec[el].symbol == 'H':
287
                disso_struct = dissociate_h(slab_molec, el + num_atoms_slab,
288
                                            num_atoms_slab)
289
                if disso_struct is not None:
290
                    disso_structs.append(disso_struct)
291
            else:
292
                for neigh_idx in nl.get_neighbors(el)[0]:
293
                    if molec[neigh_idx].symbol == 'H':
294
                        disso_struct = dissociate_h(slab_molec, neigh_idx +
295
                                                    num_atoms_slab,
296
                                                    num_atoms_slab)
297
                        if disso_struct is not None:
298
                            disso_structs.append(disso_struct)
299
        else:
300
            for atom in molec:
301
                if atom.symbol.lower() == el.lower():
302
                    for neigh_idx in nl.get_neighbors(atom.index)[0]:
303
                        if molec[neigh_idx].symbol == 'H':
304
                            disso_struct = dissociate_h(slab_molec, neigh_idx \
305
                                                        + num_atoms_slab,
306
                                                        num_atoms_slab)
307
                            if disso_struct is not None:
308
                                disso_structs.append(disso_struct)
309
    return disso_structs
310

  
311

  
230 312
def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
231 313
                 min_coll_height, norm_vect, slab_nghbs, molec_nghbs,
232 314
                 coll_coeff):
......
273 355

  
274 356

  
275 357
def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
276
              min_coll_height, coll_coeff, norm_vect, slab_nghbs, molec_nghbs):
358
              min_coll_height, coll_coeff, norm_vect, slab_nghbs, molec_nghbs,
359
              disso_atoms):
277 360
    """Generates adsorbate-surface structures by sampling over Euler angles.
278 361

  
279 362
    This function generates a number of adsorbate-surface structures at
......
293 376
    @param norm_vect: The vector perpendicular to the surface.
294 377
    @param slab_nghbs: Number of neigbors in the surface.
295 378
    @param molec_nghbs: Number of neighbors in the molecule.
379
    @param disso_atoms: List of atom types or atom numbers to try to dissociate.
296 380
    @return: list of ase.Atoms object conatining all the orientations of a given
297 381
        conformer
298 382
    """
......
314 398
                                                     coll_coeff)
315 399
                if not collision:
316 400
                    slab_ads_list.append(slab_molec)
401
                    slab_ads_list.extend(dissociation(slab_molec, disso_atoms,
402
                                                      len(slab)))
317 403

  
318 404
    return slab_ads_list
319 405

  
......
323 409

  
324 410

  
325 411
def adsorb_confs(conf_list, surf, molec_ctrs, sites, algo, num_pts, neigh_ctrs,
326
                 norm_vect, min_coll_height, coll_coeff):
412
                 norm_vect, min_coll_height, coll_coeff, disso_atoms):
327 413
    """Generates a number of adsorbate-surface structure coordinates.
328 414

  
329 415
    Given a list of conformers, a surface, a list of atom indices (or list of
......
343 429
    @param coll_coeff: The coefficient that multiplies the covalent radius of
344 430
        atoms resulting in a distance that two atoms being closer to that is
345 431
        considered as atomic collision.
432
    @param disso_atoms: List of atom types or atom numbers to try to dissociate.
346 433
    @return: list of ase.Atoms for the adsorbate-surface structures
347 434
    """
348 435
    from ase.neighborlist import natural_cutoffs, neighbor_list
......
367 454
                    surf_ads_list.extend(ads_euler(conf, surf, ctr, site,
368 455
                                                   num_pts, min_coll_height,
369 456
                                                   coll_coeff, norm_vect,
370
                                                   surf_nghbs, molec_nghbs))
457
                                                   surf_nghbs, molec_nghbs,
458
                                                   disso_atoms))
371 459
                elif algo == 'chemcat':
372 460
                    surf_ads_list.extend(ads_chemcat(site, ctr, num_pts))
373 461
    return surf_ads_list
......
401 489
                                 inp_vars['molec_neigh_ctrs'],
402 490
                                 inp_vars['surf_norm_vect'],
403 491
                                 inp_vars['min_coll_height'],
404
                                 inp_vars['collision_threshold'])
492
                                 inp_vars['collision_threshold'],
493
                                 inp_vars['disso_atoms'])
405 494
    logger.info(f'Generated {len(surf_ads_list)} adsorbate-surface atomic '
406 495
                f'configurations, to carry out a calculation of.')
407 496
    run_calc('screening', inp_vars, surf_ads_list)

Formats disponibles : Unified diff