Statistics
| Branch: | Tag: | Revision:

 1 #!/usr/bin/env python  # -*- coding: utf-8 -*-  import sys  import os  import copy  import numpy as np  import ase  import ase.io  from ase import Atom  from ase import Atoms  from ase.neighborlist import NeighborList, neighbor_list  import ase_Fe1_Fe2  #############################################################  # Définition de la fonction permettant de calculer les angles  #############################################################  def get_proper_angle(v1, v2, degrees=True):   norm_dot = np.dot(v1/np.linalg.norm(v1), v2/np.linalg.norm(v2))   angle = np.arccos(np.sign(norm_dot) if np.abs(norm_dot) >= 1 else norm_dot)   return(angle*180/np.pi if degrees else angle)  ####################################################  # Re-definition of the add-adsorbate function of ase  ####################################################  def new_add_adsorbate(slab, adsorbate, height, position=(0, 0, 0), offset=None, mol_index=0):   # works as the add_adsorbate fonction of ase but takes into account not only   # the x and y coordinates of the adsorption site but also the z in order to   # put the adsorbate at a given distance above the adsorption site (and not at   # a given height above the highest point of the surface)   #   # Sarah Blanck 22/01/2020   info = slab.info.get('adsorbate_info', {})   pos = np.array([0.0, 0.0]) # (x, y) part   pos2 = np.array([0.0, 0.0, 0.0]) # (x, y, z) part   spos = np.array([0.0, 0.0]) # part relative to unit cell   if offset is not None:   spos += np.asarray(offset, float)   if isinstance(position, basestring):   # A site-name:   if 'sites' not in info:   raise TypeError('If the atoms are not made by an ' + 'ase.build function, ' + 'position cannot be a name.')   if position not in info['sites']:   raise TypeError('Adsorption site %s not supported.' % position)   spos += info['sites'][position]   else:   pos2 += position   pos += [pos2[0], pos2[1]]   if 'cell' in info:   cell = info['cell']   else:   cell = slab.get_cell()[:2, :2]   pos += np.dot(spos, cell)   # Convert the adsorbate to an Atoms object   if isinstance(adsorbate, Atoms):   ads = adsorbate   elif isinstance(adsorbate, Atom):   ads = Atoms([adsorbate])   else:   # Assume it is a string representing a single Atom   ads = Atoms([Atom(adsorbate)])   # Get the z-coordinate:   z = pos2[2] + height   # Move adsorbate into position   ads.translate([pos[0], pos[1], z] - ads.positions[mol_index])   # Attach the adsorbate   slab.extend(ads)  def adsorb_check_rotate(molec, surf, molec_atom_1, molec_atom_2, site, init_height):   # This function tries to adsorb a molecule on top of a surface ensuring   # that the adsorbate doesn't overlap with the surface. To achieve that it   # tries subsequently: rotation around one axe, rotation around a second axe   # and increasing the height of adsorption. The check is made by comparing   # the sum of the total number of neighbors in the molecule alone and   # the surface alone with the total number of neighbors of the system once   # the molecule has been adsorbed. The cutoff_coeff sets the tolerance for   # the check of overlaping structure by scaling the natural_cutoffs from ase.   #   # Carles Martí 14/02/2020   cutoff_coeff = 0.85   molec_cutoffs = natural_cutoffs(molec, mult=cutoff_coeff)   molec_nghbs = len(neighbor_list("i",molec,molec_cutoffs))   surf_cutoffs = natural_cutoffs(surf, mult=cutoff_coeff)   surf_nghbs = len(neighbor_list("i",surf,surf_cutoffs))   sys_cutoffs = natural_cutoffs(surf+molec, mult=cutoff_coeff)     d_angle=20   vec1 = molec[molec_atom_2].position - molec[molec_atom_1].position   if np.linalg.norm(np.cross(vec1,(0,0,1))) != 0: # If vec1 is not parallel to z.   vec2 = np.cross(vec1,(0,0,1))   else:   vec2 = (1,1,0)     collision = True   for d_height in np.arange(0, 2.1, 0.4):   height = init_height + d_height   # print "Height =", height # For debuggig purposes.   for rot1 in range(0,360,d_angle):   # print "Angle 1 =", rot1 # For debuggig purposes.   new_molec = copy.deepcopy(molec)   new_molec.rotate(rot1,vec1,center=new_molec[molec_atom_1].position)   for rot2 in range(0,360,d_angle):   # print "Angle 2 =", rot2 # For debuggig purposes.   new_surf = copy.deepcopy(surf)   new_add_adsorbate(new_surf, new_molec, height, position=site, mol_index=molec_atom_1)   new_surf_nghbs = len(neighbor_list("i",new_surf, sys_cutoffs))   # print new_surf_nghbs , surf_nghbs+molec_nghbs   if new_surf_nghbs == (surf_nghbs + molec_nghbs):   collision = False   return new_surf, collision   else:   new_molec.rotate(d_angle,vec2,center=new_molec[molec_atom_1].position)   return new_surf, collision  def natural_cutoffs(atoms, mult=1, **kwargs):   """Generate a radial cutoff for every atom based on covalent radii     The covalent radii are a reasonable cutoff estimation for bonds in   many applications such as neighborlists, so function generates an   atoms length list of radii based on this idea.     atoms: An atoms object   mult: A multiplier for all cutoffs, useful for coarse grained adjustment   kwargs: Symbol of the atom and its corresponding cutoff, used to override   the covalent radii   """   return [kwargs.get(atom.symbol, ase.data.covalent_radii[atom.number] * mult)   for atom in atoms]  ##########################################  # Lecture des différents arguments de base  ##########################################  surface_file = sys.argv[1]  atom_surface = int(sys.argv[2])  molecule_file = sys.argv[3]  atom_molecule = int(sys.argv[4])  distance = float(sys.argv[5])  atom_molecule_2 = int(sys.argv[6])  phi = int(sys.argv[7])  theta = int(sys.argv[8])  psi = int(sys.argv[9])  output = sys.argv[10]  surface = ase.io.read(surface_file)  molecule = ase.io.read(molecule_file)  #################################  # Lecture des positions de atomes  #################################  atom_surf_1 = surface[atom_surface].position  x_atom_surf_1 = atom_surf_1[0]  y_atom_surf_1 = atom_surf_1[1]  z_atom_surf_1 = atom_surf_1[2]  atom_mol_1 = molecule[atom_molecule].position  x_atom_mol_1 = atom_mol_1[0]  y_atom_mol_1 = atom_mol_1[1]  z_atom_mol_1 = atom_mol_1[2]  nb_atom_mol = int(len(molecule))  theta_min = theta - 0.01  theta_max = theta + 0.01  phi_min = phi - 0.01  phi_max = phi + 0.01  psi_min = psi - 0.01  psi_max = psi + 0.01  atom_mol_2 = molecule[atom_molecule_2].position  x_atom_mol_2 = atom_mol_2[0]  y_atom_mol_2 = atom_mol_2[1]  z_atom_mol_2 = atom_mol_2[2]  ##################################  # Rotation avec les angles d'Euler  ##################################  #Initialisation : placement du second atome de la molecule sur Oy  vect_y = [0, 1, 0]  vect_atom_mol_2 = molecule[atom_molecule_2].position - molecule[atom_molecule].position  angle_1 = get_proper_angle(vect_y, vect_atom_mol_2)  if (angle_1 != 0) :   vect_normal = np.cross(vect_y, vect_atom_mol_2)   molecule.rotate(-angle_1, vect_normal, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))   vect_atom_mol_2_verif = molecule[atom_molecule_2].position - molecule[atom_molecule].position   angle_1_verif = get_proper_angle(vect_atom_mol_2_verif, vect_y)   angle_max = 0.01   angle_min = -0.01   if (angle_1_verif < angle_min or angle_1_verif > angle_max) :   print 'Error in initialisation'  #Rotation Euler  molecule.euler_rotate(phi, theta, psi, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))  # Try to adsorb molecule  adsorbed_system, collision = adsorb_check_rotate(molecule, surface, atom_molecule, atom_molecule_2, atom_surf_1, distance)  if collision:   print 'Error: the collision could not be corrected'   pos1 = output.rfind('/')   pos = pos1 + 1   num_coll=output[pos:]   j = output[:pos1]   pos2 = j.rfind('/')   molecule_directory = j[:pos2]   fichier = open("%s/errors" % molecule_directory , "a")   fichier.write("Adsorption %s : ERROR the collision could not be corrected\n" % num_coll )   fichier.close()  ##########################################  # Adsorption de la molecule sur la surface  ##########################################  out=output+".xyz"  ase.io.write(out, adsorbed_system)  ############################  # Dissociation si necessaire  ############################  if not collision:   molecule_reminder = molecule     list_cutoffs = natural_cutoffs(molecule)   nl=NeighborList(list_cutoffs, self_interaction=False, bothways=True)   nl.update(molecule)   neighbor_indices, offsets = nl.get_neighbors(atom_molecule)   symbols = molecule.get_chemical_symbols()   neighbors_symbols=[]   for i in neighbor_indices:   neighbors_symbols.append(symbols[i])   nb_neighbors = len(neighbor_indices)   diss=0   for i in range(0,nb_neighbors):   if neighbors_symbols[i] == "H":   diss=1   H_atom=neighbor_indices[i] #nb of the H atom   print "Dissociation possible"   if diss == 1:   list_cutoffs_surface = natural_cutoffs(surface)   nl2=NeighborList(list_cutoffs_surface, self_interaction=False, bothways=True)   nl2.update(surface)   neighbor_indices_surf, offsets_surf = nl2.get_neighbors(atom_surface)   symbols_surf = surface.get_chemical_symbols()   for i in neighbor_indices_surf:   surface = ase.io.read(surface_file)   molecule = molecule_reminder   neighbor_surf_coord = surface[i].position   coord_H_atom=molecule[H_atom].position   if (neighbor_surf_coord[2] > z_atom_surf_1-1 and symbols_surf[i] == "O"):   vector_H_diss=[neighbor_surf_coord[0]-coord_H_atom[0], neighbor_surf_coord[1]-coord_H_atom[1], neighbor_surf_coord[2]+1-coord_H_atom[2]]   translation_matrix_dissociation=[]   for k in range (0, nb_atom_mol):   if k != H_atom:   translation_matrix_dissociation.append([0,0,0])   else:   translation_matrix_dissociation.append(vector_H_diss)   molecule.translate(translation_matrix_dissociation)   new_add_adsorbate(surface, molecule, distance, (x_atom_surf_1,y_atom_surf_1,z_atom_surf_1), mol_index=atom_molecule)   output_diss = output + "_diss_" + str(i) + ".xyz"   ase.io.write(output_diss, surface)