Révision e8bebcca

b/modules/screening.py
198 198

  
199 199

  
200 200
def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,
201
                    nn_molec=0, coll_coeff=0.9):
201
                    nn_molec=0, coll_coeff=1.2):
202 202
    """Checks whether a slab and a molecule collide or not.
203 203

  
204 204
    @param slab_molec: The system of adsorbate-slab for which to detect if there
......
213 213
        colliding.
214 214
    @param vect: The vector perpendicular to the slab.
215 215
    @return: bool, whether the surface and the molecule collide.
216
    """  # TODO Enable both checks: 1st height, 2nd sphere collision
216
    """
217 217
    from ase.neighborlist import natural_cutoffs, neighbor_list
218

  
219
    # Check structure overlap by height
218 220
    if min_height is not False:
219 221
        cart_axes = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0],
220 222
                     [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]]
221
        if vect.tolist() in cart_axes:
222
            for atom in slab_molec[slab_num_atoms:]:
223
                for i, coord in enumerate(vect):
224
                    if coord == 0:
225
                        continue
226
                    if atom.position[i] * coord < min_height:
227
                        return True
228
            return False
229
    else:
223
        if vect.tolist() not in cart_axes:
224
            err_msg = "'min_coll_height' option is only implemented for " \
225
                      "'surf_norm_vect' to be one of the x, y or z axes. "
226
            logger.error(err_msg)
227
            raise ValueError(err_msg)
228
        for atom in slab_molec[slab_num_atoms:]:
229
            for i, coord in enumerate(vect):
230
                if coord == 0:
231
                    continue
232
                if atom.position[i] * coord < min_height * coord:
233
                    return True
234

  
235
    # Check structure overlap by sphere collision
236
    if coll_coeff is not False:
230 237
        slab_molec_cutoffs = natural_cutoffs(slab_molec, mult=coll_coeff)
231 238
        slab_molec_nghbs = len(
232 239
            neighbor_list("i", slab_molec, slab_molec_cutoffs))
233 240
        if slab_molec_nghbs > nn_slab + nn_molec:
234 241
            return True
235
        else:
236
            return False
242

  
243
    return False
244

  
245

  
246
def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
247
                 min_coll_height, norm_vect, slab_nghbs, molec_nghbs,
248
                 coll_coeff, height=2.5):
249
    # TODO Rethink this function
250
    """Tries to adsorb a molecule on a slab trying to avoid collisions by doing
251
    small rotations.
252

  
253
    @param molec: ase.Atoms object of the molecule to adsorb
254
    @param slab: ase.Atoms object of the surface on which to adsorb the
255
        molecule
256
    @param ctr_coord: The coordinates of the molecule to use as adsorption
257
        center.
258
    @param site_coord: The coordinates of the surface on which to adsorb the
259
        molecule
260
    @param num_pts: Number on which to sample Euler angles.
261
    @param min_coll_height: The lowermost height for which to detect a collision
262
    @param norm_vect: The vector perpendicular to the surface.
263
    @param slab_nghbs: Number of neigbors in the surface.
264
    @param molec_nghbs: Number of neighbors in the molecule.
265
    @param coll_coeff: The coefficient that multiplies the covalent radius of
266
        atoms resulting in a distance that two atoms being closer to that is
267
        considered as atomic collision.
268
    @param height: Height on which to try adsorption
269
    @return collision: bool, whether the structure generated has collisions
270
        between slab and adsorbate.
271
    """
272
    from copy import deepcopy
273
    slab_num_atoms = len(slab)
274
    slab_molec = []
275
    collision = True
276
    max_corr = 6  # Should be an even number
277
    d_angle = 180 / ((max_corr / 2.0) * num_pts)
278
    num_corr = 0
279
    while collision and num_corr <= max_corr:
280
        k = num_corr * (-1) ** num_corr
281
        slab_molec = deepcopy(slab)
282
        molec.euler_rotate(k * d_angle, k * d_angle / 2, k * d_angle,
283
                           center=ctr_coord)
284
        add_adsorbate(slab_molec, molec, site_coord, ctr_coord, height,
285
                      norm_vect=norm_vect)
286
        collision = check_collision(slab_molec, slab_num_atoms, min_coll_height,
287
                                    norm_vect, slab_nghbs, molec_nghbs,
288
                                    coll_coeff)
289
        num_corr += 1
290
    return slab_molec, collision
237 291

  
238 292

  
239 293
def dissociate_h(slab_molec_orig, h_idx, num_atoms_slab, h_acceptor,
......
334 388
    return disso_structs
335 389

  
336 390

  
337
def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,
338
                 min_coll_height, norm_vect, slab_nghbs, molec_nghbs,
339
                 coll_coeff, height=2.5):
340
    # TODO Rethink this function
341
    """Tries to adsorb a molecule on a slab trying to avoid collisions by doing
342
    small rotations.
343

  
344
    @param molec: ase.Atoms object of the molecule to adsorb
345
    @param slab: ase.Atoms object of the surface on which to adsorb the
346
        molecule
347
    @param ctr_coord: The coordinates of the molecule to use as adsorption
348
        center.
349
    @param site_coord: The coordinates of the surface on which to adsorb the
350
        molecule
351
    @param num_pts: Number on which to sample Euler angles.
352
    @param min_coll_height: The lowermost height for which to detect a collision
353
    @param norm_vect: The vector perpendicular to the surface.
354
    @param slab_nghbs: Number of neigbors in the surface.
355
    @param molec_nghbs: Number of neighbors in the molecule.
356
    @param coll_coeff: The coefficient that multiplies the covalent radius of
357
        atoms resulting in a distance that two atoms being closer to that is
358
        considered as atomic collision.
359
    @param height: Height on which to try adsorption
360
    @return collision: bool, whether the structure generated has collisions
361
        between slab and adsorbate.
362
    """
363
    from copy import deepcopy
364
    slab_num_atoms = len(slab)
365
    slab_molec = []
366
    collision = True
367
    max_corr = 6  # Should be an even number
368
    d_angle = 180 / ((max_corr / 2.0) * num_pts)
369
    num_corr = 0
370
    while collision and num_corr <= max_corr:
371
        k = num_corr * (-1) ** num_corr
372
        slab_molec = deepcopy(slab)
373
        molec.euler_rotate(k * d_angle, k * d_angle / 2, k * d_angle,
374
                           center=ctr_coord)
375
        add_adsorbate(slab_molec, molec, site_coord, ctr_coord, height,
376
                      norm_vect=norm_vect)
377
        collision = check_collision(slab_molec, slab_num_atoms, min_coll_height,
378
                                    norm_vect, slab_nghbs, molec_nghbs,
379
                                    coll_coeff)
380
        num_corr += 1
381
    return slab_molec, collision
382

  
383

  
384 391
def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,
385 392
              min_coll_height, coll_coeff, norm_vect, slab_nghbs, molec_nghbs,
386 393
              h_donor, h_acceptor):
......
770 777
                                                     ctr_coords, site_coords,
771 778
                                                     num_pts, min_coll_height,
772 779
                                                     norm_vect, slab_nghbs,
773
                                                     molec_nghbs, coll_coeff,
774
                                                     2.5)
780
                                                     molec_nghbs, coll_coeff)
775 781
                if not collision and \
776 782
                        not any([np.allclose(slab_molec.positions,
777 783
                                             conf.positions)
......
814 820

  
815 821
    surf_ads_list = []
816 822
    sites_coords = get_atom_coords(surf, sites)
817
    if min_coll_height is False:
823
    if coll_coeff is not False:
818 824
        surf_cutoffs = natural_cutoffs(surf, mult=coll_coeff)
819 825
        surf_nghbs = len(neighbor_list("i", surf, surf_cutoffs))
820 826
    else:
......
824 830
        if inp_vars['pbc_cell'] is not False:
825 831
            conf.set_pbc(True)
826 832
            conf.set_cell(inp_vars['pbc_cell'])
827
        if min_coll_height is False:
833
        if coll_coeff is not False:
828 834
            conf_cutoffs = natural_cutoffs(conf, mult=coll_coeff)
829 835
            molec_nghbs = len(neighbor_list("i", conf, conf_cutoffs))
830 836
        else:

Formats disponibles : Unified diff