Statistics
| Branch: | Tag: | Revision:

dockonsurf / modules / add_adsorbate_euler+dissociation.py @ 36f722bc

History | View | Annotate | Download (14.9 kB)

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

    
4
import sys
5
import os
6
import numpy as np
7
import ase
8
from ase import io
9
from ase.io import read
10
from ase.io import write
11
from ase import Atom
12
from ase import build
13
from ase.build import add_adsorbate
14
from ase import constraints
15
from ase.constraints import FixAtoms
16
from ase import neighborlist
17
from ase.neighborlist import NeighborList
18
from ase import utils
19
import ase_Fe1_Fe2
20

    
21
if len(sys.argv) < 7 :
22
    print "This script has to be used with at least 6 arguments."
23
    print "1. Name of the file containing the surface"
24
    print "2. Atom number of the atom of the surface where the adsorbate will be adsorbed"
25
    print "3. Name of the file containing the molecule to be adsorbed"
26
    print "4. Atom number of the atom of the molecule that will be adsorbed of the surface"
27
    print "5. Distance between the molecule and the surface"
28
    print "6. Name of the output file"
29
    print "Other arguments can be added before the name of the output file :"
30
    print "  Second atom of the molecule"
31
    print "  phi angle (rotation around z)"
32
    print "  theta angle (rotation around new x)"
33
    print "  psi angle (rotation around new z) "
34
    sys.exit(1)
35

    
36

    
37
#############################################################
38
# Définition de la fonction permettant de calculer les angles
39
#############################################################
40

    
41
def get_proper_angle(v1, v2, degrees=True):
42
    norm_dot = np.dot(v1/np.linalg.norm(v1), v2/np.linalg.norm(v2))
43
    angle = np.arccos(np.sign(norm_dot) if np.abs(norm_dot) >= 1 else norm_dot)
44
    return(angle*180/np.pi if degrees else angle)
45

    
46
def point_in_cell(p,h_matrix):
47
    # This function checks wether a given point defined by [x,y,z] is inside a 
48
    # given cell defined by h_matrix=[a,b,c] where a=[xa,ya,za], b=[xb,yb,zb], 
49
    # c=[xc,yc,zc]. If not it translates the point inside the cell by adding 
50
    # integer times the cell vectors, this way preserving its relative position
51
    # respect the cell copy. It basically applies periodic boundary conditions.
52
    # The check is performed by considering that a point p lies inside the cell, 
53
    # when it can be expressed as:
54
    # p = alpha*a + beta*b + gamma*c,   with:  0 < alpha, beta, gamma < 1.
55
    # alpha, beta and gamma are calculated as the determinant of the matrix 
56
    # defined by substituting the cell component whose coefficient we want to 
57
    # calculate by p, everything divided by the cell volume (det[a,b,c]) 
58
    # e.g: alpha = det([p,b,c])/det([a,b,c]); beta = det([a,p,c])/det([a,b,c])
59
    # gamma = det([a,b,p])/det([a,b,c])
60
    #
61
    #                                                   Carles Martí 12/12/2019
62

    
63
    try: assert len(p) == 3, "p has not 3 elements"
64
    except TypeError: raise TypeError("p is not an array-like variable")
65
    try: [float(coord) for coord in p]
66
    except ValueError: raise ValueError("The elements of p are not numbers")
67

    
68
    try: assert len(h_matrix) == 3, "h_matrix has not 3 elements"
69
    except ValueError: raise TypeError("h_matrix is not an array-like variable")
70
    try: [[float(coord) for coord in v] for v in h_matrix],
71
    except ValueError: raise ValueError("The elements of h_matrix are not numbers")
72

    
73
    assert np.linalg.det(h_matrix) != 0, "The three vectors in h_matrix are not lineally independent"
74

    
75
    h_matrix = np.array(h_matrix)
76
    p = np.array(p)
77
    
78
    volume = np.linalg.det(h_matrix)
79
    alpha = np.linalg.det(np.array([p, h_matrix[1], h_matrix[2]])) / volume
80
    beta = np.linalg.det(np.array([h_matrix[0], p, h_matrix[2]])) / volume
81
    gamma = np.linalg.det(np.array([h_matrix[0], h_matrix[1], p])) / volume
82

    
83
    p_prime = np.array(p)
84

    
85
    if 0 > alpha: p_prime = p_prime +  (1 + int(-alpha)) * h_matrix[0]
86
    if 1 < alpha: p_prime = p_prime + int(-alpha) * h_matrix[0]
87
    if 0 > beta: p_prime = p_prime +  (1 + int(-beta)) * h_matrix[1]
88
    if 1 < beta: p_prime = p_prime + int(-beta) * h_matrix[1]
89
    if 0 > gamma: p_prime = p_prime +  (1 + int(-gamma)) * h_matrix[2]
90
    if 1 < gamma: p_prime = p_prime + int(-gamma) * h_matrix[2]
91

    
92
    return p_prime,[alpha,beta,gamma]
93

    
94
def check_coads_collision(position, ads_site_pos, ads_atom_pos, box_low_lim, box_high_lim, h_matrix):
95
    # This routine checks wether an atom of the second adsorbate collides
96
    # with the box around the first molecule.           Carles Martí 12/12/2019
97

    
98
    collision=0
99
    
100
    # Define the position of the atom in the surface coordinates.
101
    x_i = position[0] - ads_atom_pos[0] + ads_site_pos[0]
102
    y_i = position[1] - ads_atom_pos[1] + ads_site_pos[1]
103
    z_i = position[2] - ads_atom_pos[2] + ads_site_pos[2] + 1
104

    
105
    x_i = point_in_cell([x_i,y_i,z_i],h_matrix)[0][0]
106
    y_i = point_in_cell([x_i,y_i,z_i],h_matrix)[0][1]
107
    z_i = point_in_cell([x_i, y_i,z_i],h_matrix)[0][2]
108

    
109
    # Check if the atom lies inside the coadsorption-prohibited box
110
    if all([box_low_lim[0] < x_i < box_high_lim[0], box_low_lim[1] < y_i < box_high_lim[1], box_low_lim[2]< z_i < box_high_lim[2]]):
111
        collision = 1
112

    
113
    return collision
114
    
115
##########################################
116
# Lecture des différents arguments de base
117
##########################################
118

    
119
surface_file = sys.argv[1]
120
atom_surface = int(sys.argv[2])
121
molecule_file = sys.argv[3]
122
atom_molecule = int(sys.argv[4])
123
distance = float(sys.argv[5])
124
if len(sys.argv) < 8 : output = sys.argv[6]
125

    
126
surface = read(surface_file)
127
molecule = read(molecule_file)
128

    
129
#################################
130
# Lecture des positions de atomes
131
#################################
132

    
133
atom_surf_1 = surface[atom_surface].position
134
x_atom_surf_1 = atom_surf_1[0]
135
y_atom_surf_1 = atom_surf_1[1]
136
z_atom_surf_1 = atom_surf_1[2]
137

    
138
atom_mol_1 = molecule[atom_molecule].position
139
x_atom_mol_1 = atom_mol_1[0]
140
y_atom_mol_1 = atom_mol_1[1]
141
z_atom_mol_1 = atom_mol_1[2]
142

    
143
nb_atom_mol = int(len(molecule))
144

    
145
h_matrix=[[os.getenv("a1"),os.getenv("a2"),os.getenv("a3")],[os.getenv("b1"),os.getenv("b2"),os.getenv("b3")],[os.getenv("c1"),os.getenv("c2"),os.getenv("c3")]]
146
h_matrix=[[float(coord) for coord in v] for v in h_matrix]
147

    
148
#######################################################
149
# Cas avec un plus grand nombre d'arguments (rotations)
150
#######################################################
151

    
152

    
153
if len(sys.argv) > 7 :
154

    
155
    #######################
156
    # Lecture des arguments
157
    #######################
158

    
159
    atom_molecule_2 = int(sys.argv[6])
160
    phi = int(sys.argv[7])
161
    theta = int(sys.argv[8])
162
    psi = int(sys.argv[9])
163
    output = sys.argv[10]
164
    coads_box = sys.argv[11]
165
    if coads_box == "true": first_atom_ads1 = int(sys.argv[12])
166

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

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

    
179
    # In the case of coadsorption, definition of a box around the first
180
    # adsorbate to prevent the second one to overlap with the former.
181
    if coads_box == "true":        
182
        max_float=sys.float_info.max
183
        box_low_lim=[max_float,max_float,max_float]
184
        box_high_lim=[-max_float,-max_float,-max_float]
185
        
186
        for adsorbate_atom in surface[first_atom_ads1:]:
187
            if adsorbate_atom.position[0] < box_low_lim[0]: box_low_lim[0] = adsorbate_atom.position[0]
188
            if adsorbate_atom.position[1] < box_low_lim[1]: box_low_lim[1] = adsorbate_atom.position[1]
189
            if adsorbate_atom.position[2] < box_low_lim[2]: box_low_lim[2] = adsorbate_atom.position[2]
190

    
191
            if adsorbate_atom.position[0] > box_high_lim[0]: box_high_lim[0] = adsorbate_atom.position[0]
192
            if adsorbate_atom.position[1] > box_high_lim[1]: box_high_lim[1] = adsorbate_atom.position[1]
193
            if adsorbate_atom.position[2] > box_high_lim[2]: box_high_lim[2] = adsorbate_atom.position[2]
194

    
195
        # Adding an offset to the limit
196
        z_max_ads_1 = box_high_lim[2]
197
        box_low_lim = [low_lim - 1 for low_lim in box_low_lim]
198
        box_high_lim = [ high_lim + 1 for high_lim in box_high_lim]
199

    
200
    ##################################
201
    # Rotation avec les angles d'Euler
202
    ##################################
203

    
204
    #Initialisation : placement du second atome de la molecule sur Oy
205
    vect_y = [0, 1, 0]
206
    vect_atom_mol_2 = molecule[atom_molecule_2].position - molecule[atom_molecule].position
207
    angle_1 = get_proper_angle(vect_y, vect_atom_mol_2)
208
    if (angle_1 != 0) :
209
        vect_normal = np.cross(vect_y, vect_atom_mol_2)
210
        molecule.rotate(-angle_1, vect_normal, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
211
        vect_atom_mol_2_verif = molecule[atom_molecule_2].position - molecule[atom_molecule].position
212
        angle_1_verif = get_proper_angle(vect_atom_mol_2_verif, vect_y)
213
        angle_max = 0.01
214
        angle_min = -0.01
215
        if (angle_1_verif < angle_min or angle_1_verif > angle_max) :
216
            print 'Error in initialisation'
217

    
218
    #Rotation Euler
219
    molecule.euler_rotate(phi, theta, psi, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
220

    
221
    #Correction des collisions entre la molécule et la surface
222
    z_atom_surf_max = z_atom_surf_1 + 1
223

    
224
    collision = 0
225
    rotation_tot = 0
226
    nb_atom_collision_min = 0
227
    rotation_opt = 0
228
    rotation_modif = 5
229
    corrected_collision=False
230
    
231
    for i in range(0,nb_atom_mol):
232
        atom_test = molecule[i].position
233
        z_atom_test = atom_test[2]
234
        z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
235
        if (z_atom_test_final < z_atom_surf_max):
236
            collision = 1
237
        if coads_box == "true":
238
            collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)
239
            collision += collision_box
240

    
241
    if collision > 0 :
242
        print 'Collision between the molecule and the surface - modification of theta angle'
243

    
244
    vect_z = [0, 0, 1]
245

    
246
    while (collision > 0 and rotation_tot < 354) :
247
        vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position
248
        
249
        if (vect_atom_mol_2_post_euler[0] != 0 and vect_atom_mol_2_post_euler[1] != 0) :
250
            vect_rotation_theta = np.cross(vect_atom_mol_2_post_euler, vect_z)
251
        else :
252
            vect_rotation_theta = [0.0001, 0.0001 , 0]
253
        molecule.rotate(rotation_modif, vect_rotation_theta, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
254
        rotation_tot += 5
255
        
256
        collision = 0
257
        for i in range(0,nb_atom_mol) :
258
            atom_test = molecule[i].position
259
            z_atom_test = atom_test[2]
260
            z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
261
            if (z_atom_test_final < z_atom_surf_max) :
262
                collision = 1
263
            if coads_box == "true":
264
                collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)
265
                collision += collision_box
266

    
267
        rotation_tot_2 = 0
268

    
269
        while (collision > 0 and rotation_tot_2 < 354) :
270
            vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position
271
            molecule.rotate(rotation_modif, vect_atom_mol_2_post_euler, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
272
            rotation_tot_2 += 5
273
            
274
            collision = 0
275

    
276
            for i in range(0,nb_atom_mol) :
277
                atom_test = molecule[i].position
278
                z_atom_test = atom_test[2]
279
                z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
280
                if (z_atom_test_final < z_atom_surf_max):
281
                    collision = 1
282
                if coads_box == "true":
283
                    collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)
284
                    collision += collision_box
285
        if collision == 0: corrected_collision = True
286

    
287
    if corrected_collision == True:  print 'Collision corrected'
288
    if collision > 0:
289
        print 'Error: the collision could not be corrected'
290
        pos1 = output.rfind('/')
291
        pos = pos1 + 1
292
        num_coll=output[pos:]
293
        j = output[:pos1]
294
        pos2 = j.rfind('/')
295
        molecule_directory = j[:pos2]
296
        fichier = open("%s/errors" % molecule_directory , "a")
297
        fichier.write("Adsorption %s : ERROR the collision could not be corrected\n" % num_coll )
298
        fichier.close()
299

    
300

    
301
##########################################
302
# Adsorption de la molecule sur la surface
303
##########################################
304

    
305
if coads_box == "true":
306
    surf_atom_z=-max_float
307
    for atom in surface[:first_atom_ads1]:
308
        if atom.position[2] > surf_atom_z: surf_atom_z = atom.position[2]        
309
    
310
    distance = distance - z_max_ads_1 + surf_atom_z
311

    
312
add_adsorbate(surface, molecule, distance, (x_atom_surf_1,y_atom_surf_1), mol_index=atom_molecule)
313

    
314
out=output+".xyz"
315
write(out, surface)
316

    
317

    
318
############################
319
# Dissociation si necessaire
320
############################
321

    
322
if collision == 0 :
323
    surface = read(surface_file)
324
    molecule_reminder = molecule
325
    
326
    list_cutoffs = utils.natural_cutoffs(molecule)
327
    nl=NeighborList(list_cutoffs, self_interaction=False, bothways=True)
328

    
329
    nl.update(molecule)
330
    neighbor_indices, offsets = nl.get_neighbors(atom_molecule)
331

    
332
    symbols = molecule.get_chemical_symbols()
333

    
334
    neighbors_symbols=[]
335
    for i in neighbor_indices:
336
            neighbors_symbols.append(symbols[i])
337

    
338
    nb_neighbors = len(neighbor_indices)
339

    
340
    diss=0
341
    for i in range(0,nb_neighbors):
342
        if neighbors_symbols[i] == "H":
343
            diss=1
344
            H_atom=neighbor_indices[i] #nb of the H atom
345
            print "Dissociation possible"
346

    
347
    if diss == 1:
348
        list_cutoffs_surface = utils.natural_cutoffs(surface)
349
        nl2=NeighborList(list_cutoffs_surface, self_interaction=False, bothways=True)
350
        nl2.update(surface)
351
        neighbor_indices_surf, offsets_surf = nl2.get_neighbors(atom_surface)
352
        symbols_surf = surface.get_chemical_symbols()
353
        for i in neighbor_indices_surf:
354
            surface = read(surface_file)
355
            molecule = molecule_reminder
356
            neighbor_surf_coord = surface[i].position
357
            coord_H_atom=molecule[H_atom].position
358
            if (neighbor_surf_coord[2] > z_atom_surf_1-1 and symbols_surf[i] == "O"):
359
                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]]
360
                translation_matrix_dissociation=[]
361
                for k in range (0, nb_atom_mol):
362
                    if k != H_atom:
363
                        translation_matrix_dissociation.append([0,0,0])
364
                    else:
365
                        translation_matrix_dissociation.append(vector_H_diss)
366
                molecule.translate(translation_matrix_dissociation)
367
                add_adsorbate(surface, molecule, distance, (x_atom_surf_1,y_atom_surf_1), mol_index=atom_molecule)
368
                output_diss = output + "_diss_" + str(i) + ".xyz"
369
                write(output_diss, surface)
370
    
371

    
372

    
373

    
374

    
375