Statistics
| Branch: | Tag: | Revision:

dockonsurf / modules / add_adsorbate_euler+dissociation.py @ 3c26e233

History | View | Annotate | Download (10.5 kB)

1
#!/usr/bin/env python
2
# -*- coding: utf-8 -*-
3

    
4
import sys
5
import os
6
import copy
7
import numpy as np
8
import ase
9
import ase.io
10
from ase import Atom
11
from ase import Atoms
12
from ase.neighborlist import NeighborList, neighbor_list
13
import ase_Fe1_Fe2
14

    
15
#############################################################
16
# Définition de la fonction permettant de calculer les angles
17
#############################################################
18

    
19
def get_proper_angle(v1, v2, degrees=True):
20
    norm_dot = np.dot(v1/np.linalg.norm(v1), v2/np.linalg.norm(v2))
21
    angle = np.arccos(np.sign(norm_dot) if np.abs(norm_dot) >= 1 else norm_dot)
22
    return(angle*180/np.pi if degrees else angle)
23

    
24

    
25
####################################################
26
# Re-definition of the add-adsorbate function of ase
27
####################################################
28

    
29
def new_add_adsorbate(slab, adsorbate, height, position=(0, 0, 0), offset=None, mol_index=0):
30
    # works as the add_adsorbate fonction of ase but takes into account not only
31
    # the x and y coordinates of the adsorption site but also the z in order to 
32
    # put the adsorbate at a given distance above the adsorption site (and not at
33
    # a given height above the highest point of the surface)
34
    # 
35
    #                                                   Sarah Blanck 22/01/2020
36

    
37
    info = slab.info.get('adsorbate_info', {})
38
    pos = np.array([0.0, 0.0]) # (x, y) part
39
    pos2 = np.array([0.0, 0.0, 0.0])  # (x, y, z) part
40
    spos = np.array([0.0, 0.0])  # part relative to unit cell
41
    if offset is not None:
42
        spos += np.asarray(offset, float)
43
    if isinstance(position, basestring):
44
        # A site-name:
45
        if 'sites' not in info:
46
            raise TypeError('If the atoms are not made by an ' + 'ase.build function, ' + 'position cannot be a name.')
47
        if position not in info['sites']:
48
            raise TypeError('Adsorption site %s not supported.' % position)
49
        spos += info['sites'][position]
50
    else:
51
        pos2 += position
52
        pos += [pos2[0], pos2[1]]
53
    if 'cell' in info:
54
        cell = info['cell']
55
    else:
56
        cell = slab.get_cell()[:2, :2]
57
    pos += np.dot(spos, cell)
58
    # Convert the adsorbate to an Atoms object
59
    if isinstance(adsorbate, Atoms):
60
        ads = adsorbate
61
    elif isinstance(adsorbate, Atom):
62
        ads = Atoms([adsorbate])
63
    else:
64
        # Assume it is a string representing a single Atom
65
        ads = Atoms([Atom(adsorbate)])
66
    # Get the z-coordinate:
67
    z = pos2[2] + height
68
    # Move adsorbate into position
69
    ads.translate([pos[0], pos[1], z] - ads.positions[mol_index])
70
    # Attach the adsorbate
71
    slab.extend(ads)
72

    
73
def adsorb_check_rotate(molec, surf, molec_atom_1, molec_atom_2, site, init_height):
74
    # This function tries to adsorb a molecule on top of a surface ensuring 
75
    # that the adsorbate doesn't overlap with the surface. To achieve that it
76
    # tries subsequently: rotation around one axe, rotation around a second axe
77
    # and increasing the height of adsorption. The check is made by comparing 
78
    # the sum of the total number of neighbors in the molecule alone and 
79
    # the surface alone with the total number of neighbors of the system once
80
    # the molecule has been adsorbed. The cutoff_coeff sets the tolerance for 
81
    # the check of overlaping structure by scaling the natural_cutoffs from ase.
82
    #
83
    #                                                  Carles Martí 14/02/2020
84

    
85
    cutoff_coeff = 0.85
86
    molec_cutoffs = natural_cutoffs(molec, mult=cutoff_coeff)
87
    molec_nghbs = len(neighbor_list("i",molec,molec_cutoffs))
88
    surf_cutoffs = natural_cutoffs(surf, mult=cutoff_coeff)
89
    surf_nghbs = len(neighbor_list("i",surf,surf_cutoffs))
90
    sys_cutoffs = natural_cutoffs(surf+molec, mult=cutoff_coeff)
91
    
92
    d_angle=20
93
    vec1 = molec[molec_atom_2].position - molec[molec_atom_1].position
94
    if np.linalg.norm(np.cross(vec1,(0,0,1))) != 0: # If vec1 is not parallel to z.
95
        vec2 = np.cross(vec1,(0,0,1))
96
    else:
97
        vec2 = (1,1,0)
98
    
99
    collision = True
100
    for d_height in np.arange(0, 2.1, 0.4):
101
        height = init_height + d_height
102
        # print "Height =", height # For debuggig purposes.
103
        for rot1 in range(0,360,d_angle):
104
            # print "Angle 1 =", rot1 # For debuggig purposes.
105
            new_molec = copy.deepcopy(molec)
106
            new_molec.rotate(rot1,vec1,center=new_molec[molec_atom_1].position)
107
            for rot2 in range(0,360,d_angle):
108
                # print "Angle 2 =", rot2 # For debuggig purposes.
109
                new_surf = copy.deepcopy(surf)
110
                new_add_adsorbate(new_surf, new_molec, height, position=site, mol_index=molec_atom_1)
111
                new_surf_nghbs = len(neighbor_list("i",new_surf, sys_cutoffs))
112
                # print new_surf_nghbs , surf_nghbs+molec_nghbs
113
                if new_surf_nghbs == (surf_nghbs + molec_nghbs):
114
                    collision = False
115
                    return new_surf, collision
116
                else:
117
                    new_molec.rotate(d_angle,vec2,center=new_molec[molec_atom_1].position)
118
    return new_surf, collision
119

    
120
def natural_cutoffs(atoms, mult=1, **kwargs):
121
    """Generate a radial cutoff for every atom based on covalent radii
122

123
    The covalent radii are a reasonable cutoff estimation for bonds in
124
    many applications such as neighborlists, so function generates an
125
    atoms length list of radii based on this idea.
126

127
    atoms: An atoms object
128
    mult: A multiplier for all cutoffs, useful for coarse grained adjustment
129
    kwargs: Symbol of the atom and its corresponding cutoff, used to override
130
            the covalent radii
131
    """
132
    return [kwargs.get(atom.symbol, ase.data.covalent_radii[atom.number] * mult)
133
            for atom in atoms]
134

    
135
##########################################
136
# Lecture des différents arguments de base
137
##########################################
138

    
139
surface_file = sys.argv[1]
140
atom_surface = int(sys.argv[2])
141
molecule_file = sys.argv[3]
142
atom_molecule = int(sys.argv[4])
143
distance = float(sys.argv[5])
144
atom_molecule_2 = int(sys.argv[6])
145
phi = int(sys.argv[7])
146
theta = int(sys.argv[8])
147
psi = int(sys.argv[9])
148
output = sys.argv[10]
149

    
150
surface = ase.io.read(surface_file)
151
molecule = ase.io.read(molecule_file)
152

    
153
#################################
154
# Lecture des positions de atomes
155
#################################
156

    
157
atom_surf_1 = surface[atom_surface].position
158
x_atom_surf_1 = atom_surf_1[0]
159
y_atom_surf_1 = atom_surf_1[1]
160
z_atom_surf_1 = atom_surf_1[2]
161

    
162
atom_mol_1 = molecule[atom_molecule].position
163
x_atom_mol_1 = atom_mol_1[0]
164
y_atom_mol_1 = atom_mol_1[1]
165
z_atom_mol_1 = atom_mol_1[2]
166

    
167
nb_atom_mol = int(len(molecule))
168

    
169
theta_min = theta - 0.01
170
theta_max = theta + 0.01
171
phi_min = phi - 0.01
172
phi_max = phi + 0.01
173
psi_min = psi - 0.01
174
psi_max = psi + 0.01
175

    
176
atom_mol_2 = molecule[atom_molecule_2].position
177
x_atom_mol_2 = atom_mol_2[0]
178
y_atom_mol_2 = atom_mol_2[1]
179
z_atom_mol_2 = atom_mol_2[2]
180

    
181
##################################
182
# Rotation avec les angles d'Euler
183
##################################
184

    
185
#Initialisation : placement du second atome de la molecule sur Oy
186
vect_y = [0, 1, 0]
187
vect_atom_mol_2 = molecule[atom_molecule_2].position - molecule[atom_molecule].position
188
angle_1 = get_proper_angle(vect_y, vect_atom_mol_2)
189
if (angle_1 != 0) :
190
    vect_normal = np.cross(vect_y, vect_atom_mol_2)
191
    molecule.rotate(-angle_1, vect_normal, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
192
    vect_atom_mol_2_verif = molecule[atom_molecule_2].position - molecule[atom_molecule].position
193
    angle_1_verif = get_proper_angle(vect_atom_mol_2_verif, vect_y)
194
    angle_max = 0.01
195
    angle_min = -0.01
196
    if (angle_1_verif < angle_min or angle_1_verif > angle_max) :
197
        print 'Error in initialisation'
198

    
199
#Rotation Euler
200
molecule.euler_rotate(phi, theta, psi, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
201

    
202
# Try to adsorb molecule
203
adsorbed_system, collision = adsorb_check_rotate(molecule, surface, atom_molecule, atom_molecule_2, atom_surf_1, distance)
204

    
205

    
206
if collision:
207
    print 'Error: the collision could not be corrected'
208
    pos1 = output.rfind('/')
209
    pos = pos1 + 1
210
    num_coll=output[pos:]
211
    j = output[:pos1]
212
    pos2 = j.rfind('/')
213
    molecule_directory = j[:pos2]
214
    fichier = open("%s/errors" % molecule_directory , "a")
215
    fichier.write("Adsorption %s : ERROR the collision could not be corrected\n" % num_coll )
216
    fichier.close()
217

    
218

    
219
##########################################
220
# Adsorption de la molecule sur la surface
221
##########################################
222

    
223
out=output+".xyz"
224
ase.io.write(out, adsorbed_system)
225

    
226

    
227
############################
228
# Dissociation si necessaire
229
############################
230

    
231
if not collision:
232
    molecule_reminder = molecule
233
    
234
    list_cutoffs = natural_cutoffs(molecule)
235
    nl=NeighborList(list_cutoffs, self_interaction=False, bothways=True)
236

    
237
    nl.update(molecule)
238
    neighbor_indices, offsets = nl.get_neighbors(atom_molecule)
239

    
240
    symbols = molecule.get_chemical_symbols()
241

    
242
    neighbors_symbols=[]
243
    for i in neighbor_indices:
244
            neighbors_symbols.append(symbols[i])
245

    
246
    nb_neighbors = len(neighbor_indices)
247

    
248
    diss=0
249
    for i in range(0,nb_neighbors):
250
        if neighbors_symbols[i] == "H":
251
            diss=1
252
            H_atom=neighbor_indices[i] #nb of the H atom
253
            print "Dissociation possible"
254

    
255
    if diss == 1:
256
        list_cutoffs_surface = natural_cutoffs(surface)
257
        nl2=NeighborList(list_cutoffs_surface, self_interaction=False, bothways=True)
258
        nl2.update(surface)
259
        neighbor_indices_surf, offsets_surf = nl2.get_neighbors(atom_surface)
260
        symbols_surf = surface.get_chemical_symbols()
261
        for i in neighbor_indices_surf:
262
            surface = ase.io.read(surface_file)
263
            molecule = molecule_reminder
264
            neighbor_surf_coord = surface[i].position
265
            coord_H_atom=molecule[H_atom].position
266
            if (neighbor_surf_coord[2] > z_atom_surf_1-1 and symbols_surf[i] == "O"):
267
                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]]
268
                translation_matrix_dissociation=[]
269
                for k in range (0, nb_atom_mol):
270
                    if k != H_atom:
271
                        translation_matrix_dissociation.append([0,0,0])
272
                    else:
273
                        translation_matrix_dissociation.append(vector_H_diss)
274
                molecule.translate(translation_matrix_dissociation)
275
                new_add_adsorbate(surface, molecule, distance, (x_atom_surf_1,y_atom_surf_1,z_atom_surf_1), mol_index=atom_molecule)
276
                output_diss = output + "_diss_" + str(i) + ".xyz"
277
                ase.io.write(output_diss, surface)
278