Revision 3c26e233 modules/add_adsorbate_euler+dissociation.py

b/modules/add_adsorbate_euler+dissociation.py
3 3

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

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

  
37

  
38 15
#############################################################
39 16
# Définition de la fonction permettant de calculer les angles
40 17
#############################################################
......
93 70
    # Attach the adsorbate
94 71
    slab.extend(ads)
95 72

  
96

  
97
###########################################################
98
# Definition of PBC conditions for molecules on the surface
99
###########################################################
100

  
101
def point_in_cell(p,h_matrix):
102
    # This function checks wether a given point defined by [x,y,z] is inside a 
103
    # given cell defined by h_matrix=[a,b,c] where a=[xa,ya,za], b=[xb,yb,zb], 
104
    # c=[xc,yc,zc]. If not it translates the point inside the cell by adding 
105
    # integer times the cell vectors, this way preserving its relative position
106
    # respect the cell copy. It basically applies periodic boundary conditions.
107
    # The check is performed by considering that a point p lies inside the cell, 
108
    # when it can be expressed as:
109
    # p = alpha*a + beta*b + gamma*c,   with:  0 < alpha, beta, gamma < 1.
110
    # alpha, beta and gamma are calculated as the determinant of the matrix 
111
    # defined by substituting the cell component whose coefficient we want to 
112
    # calculate by p, everything divided by the cell volume (det[a,b,c]) 
113
    # e.g: alpha = det([p,b,c])/det([a,b,c]); beta = det([a,p,c])/det([a,b,c])
114
    # gamma = det([a,b,p])/det([a,b,c])
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.
115 82
    #
116
    #                                                   Carles Martí 12/12/2019
117

  
118
    try: assert len(p) == 3, "p has not 3 elements"
119
    except TypeError: raise TypeError("p is not an array-like variable")
120
    try: [float(coord) for coord in p]
121
    except ValueError: raise ValueError("The elements of p are not numbers")
122

  
123
    try: assert len(h_matrix) == 3, "h_matrix has not 3 elements"
124
    except ValueError: raise TypeError("h_matrix is not an array-like variable")
125
    try: [[float(coord) for coord in v] for v in h_matrix],
126
    except ValueError: raise ValueError("The elements of h_matrix are not numbers")
127

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

  
130
    h_matrix = np.array(h_matrix)
131
    p = np.array(p)
132
    
133
    volume = np.linalg.det(h_matrix)
134
    alpha = np.linalg.det(np.array([p, h_matrix[1], h_matrix[2]])) / volume
135
    beta = np.linalg.det(np.array([h_matrix[0], p, h_matrix[2]])) / volume
136
    gamma = np.linalg.det(np.array([h_matrix[0], h_matrix[1], p])) / volume
137

  
138
    p_prime = np.array(p)
139

  
140
    if 0 > alpha: p_prime = p_prime +  (1 + int(-alpha)) * h_matrix[0]
141
    if 1 < alpha: p_prime = p_prime + int(-alpha) * h_matrix[0]
142
    if 0 > beta: p_prime = p_prime +  (1 + int(-beta)) * h_matrix[1]
143
    if 1 < beta: p_prime = p_prime + int(-beta) * h_matrix[1]
144
    if 0 > gamma: p_prime = p_prime +  (1 + int(-gamma)) * h_matrix[2]
145
    if 1 < gamma: p_prime = p_prime + int(-gamma) * h_matrix[2]
146

  
147
    return p_prime,[alpha,beta,gamma]
148

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

  
153
    collision=0
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)
154 91
    
155
    # Define the position of the atom in the surface coordinates.
156
    x_i = position[0] - ads_atom_pos[0] + ads_site_pos[0]
157
    y_i = position[1] - ads_atom_pos[1] + ads_site_pos[1]
158
    z_i = position[2] - ads_atom_pos[2] + ads_site_pos[2] + 1
159

  
160
    x_i = point_in_cell([x_i,y_i,z_i],h_matrix)[0][0]
161
    y_i = point_in_cell([x_i,y_i,z_i],h_matrix)[0][1]
162
    z_i = point_in_cell([x_i, y_i,z_i],h_matrix)[0][2]
163

  
164
    # Check if the atom lies inside the coadsorption-prohibited box
165
    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]]):
166
        collision = 1
167

  
168
    return collision
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)
169 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]
170 134

  
171 135
##########################################
172 136
# Lecture des différents arguments de base
......
177 141
molecule_file = sys.argv[3]
178 142
atom_molecule = int(sys.argv[4])
179 143
distance = float(sys.argv[5])
180
if len(sys.argv) < 8 : output = sys.argv[6]
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]
181 149

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

  
185 153
#################################
186 154
# Lecture des positions de atomes
......
198 166

  
199 167
nb_atom_mol = int(len(molecule))
200 168

  
201
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")]]
202
h_matrix=[[float(coord) for coord in v] for v in h_matrix]
203

  
204
#######################################################
205
# Cas avec un plus grand nombre d'arguments (rotations)
206
#######################################################
207

  
208

  
209
if len(sys.argv) > 7 :
210

  
211
    #######################
212
    # Lecture des arguments
213
    #######################
214

  
215
    atom_molecule_2 = int(sys.argv[6])
216
    phi = int(sys.argv[7])
217
    theta = int(sys.argv[8])
218
    psi = int(sys.argv[9])
219
    output = sys.argv[10]
220
    coads_box = sys.argv[11]
221
    if coads_box == "true": first_atom_ads1 = int(sys.argv[12])
222

  
223
    theta_min = theta - 0.01
224
    theta_max = theta + 0.01
225
    phi_min = phi - 0.01
226
    phi_max = phi + 0.01
227
    psi_min = psi - 0.01
228
    psi_max = psi + 0.01
229

  
230
    atom_mol_2 = molecule[atom_molecule_2].position
231
    x_atom_mol_2 = atom_mol_2[0]
232
    y_atom_mol_2 = atom_mol_2[1]
233
    z_atom_mol_2 = atom_mol_2[2]
234

  
235
    # In the case of coadsorption, definition of a box around the first
236
    # adsorbate to prevent the second one to overlap with the former.
237
    if coads_box == "true":        
238
        max_float=sys.float_info.max
239
        box_low_lim=[max_float,max_float,max_float]
240
        box_high_lim=[-max_float,-max_float,-max_float]
241
        
242
        for adsorbate_atom in surface[first_atom_ads1:]:
243
            if adsorbate_atom.position[0] < box_low_lim[0]: box_low_lim[0] = adsorbate_atom.position[0]
244
            if adsorbate_atom.position[1] < box_low_lim[1]: box_low_lim[1] = adsorbate_atom.position[1]
245
            if adsorbate_atom.position[2] < box_low_lim[2]: box_low_lim[2] = adsorbate_atom.position[2]
246

  
247
            if adsorbate_atom.position[0] > box_high_lim[0]: box_high_lim[0] = adsorbate_atom.position[0]
248
            if adsorbate_atom.position[1] > box_high_lim[1]: box_high_lim[1] = adsorbate_atom.position[1]
249
            if adsorbate_atom.position[2] > box_high_lim[2]: box_high_lim[2] = adsorbate_atom.position[2]
250

  
251
        # Adding an offset to the limit
252
        z_max_ads_1 = box_high_lim[2]
253
        box_low_lim = [low_lim - 1 for low_lim in box_low_lim]
254
        box_high_lim = [ high_lim + 1 for high_lim in box_high_lim]
255

  
256
    ##################################
257
    # Rotation avec les angles d'Euler
258
    ##################################
259

  
260
    #Initialisation : placement du second atome de la molecule sur Oy
261
    vect_y = [0, 1, 0]
262
    vect_atom_mol_2 = molecule[atom_molecule_2].position - molecule[atom_molecule].position
263
    angle_1 = get_proper_angle(vect_y, vect_atom_mol_2)
264
    if (angle_1 != 0) :
265
        vect_normal = np.cross(vect_y, vect_atom_mol_2)
266
        molecule.rotate(-angle_1, vect_normal, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
267
        vect_atom_mol_2_verif = molecule[atom_molecule_2].position - molecule[atom_molecule].position
268
        angle_1_verif = get_proper_angle(vect_atom_mol_2_verif, vect_y)
269
        angle_max = 0.01
270
        angle_min = -0.01
271
        if (angle_1_verif < angle_min or angle_1_verif > angle_max) :
272
            print 'Error in initialisation'
273

  
274
    #Rotation Euler
275
    molecule.euler_rotate(phi, theta, psi, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
276

  
277
    #Correction des collisions entre la molécule et la surface
278
    z_atom_surf_max = z_atom_surf_1 + 1
279

  
280
    collision = 0
281
    rotation_tot = 0
282
    nb_atom_collision_min = 0
283
    rotation_opt = 0
284
    rotation_modif = 5
285
    corrected_collision=False
286
    
287
    for i in range(0,nb_atom_mol):
288
        atom_test = molecule[i].position
289
        z_atom_test = atom_test[2]
290
        z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
291
        if (z_atom_test_final < z_atom_surf_max):
292
            collision = 1
293
        if coads_box == "true":
294
            collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)
295
            collision += collision_box
296

  
297
    if collision > 0 :
298
        print 'Collision between the molecule and the surface - modification of theta angle'
299

  
300
    vect_z = [0, 0, 1]
301

  
302
    while (collision > 0 and rotation_tot < 354) :
303
        vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position
304
        
305
        if (vect_atom_mol_2_post_euler[0] != 0 and vect_atom_mol_2_post_euler[1] != 0) :
306
            vect_rotation_theta = np.cross(vect_atom_mol_2_post_euler, vect_z)
307
        else :
308
            vect_rotation_theta = [0.0001, 0.0001 , 0]
309
        molecule.rotate(rotation_modif, vect_rotation_theta, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
310
        rotation_tot += 5
311
        
312
        collision = 0
313
        for i in range(0,nb_atom_mol) :
314
            atom_test = molecule[i].position
315
            z_atom_test = atom_test[2]
316
            z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
317
            if (z_atom_test_final < z_atom_surf_max) :
318
                collision = 1
319
            if coads_box == "true":
320
                collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)
321
                collision += collision_box
322

  
323
        rotation_tot_2 = 0
324

  
325
        while (collision > 0 and rotation_tot_2 < 354) :
326
            vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position
327
            molecule.rotate(rotation_modif, vect_atom_mol_2_post_euler, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
328
            rotation_tot_2 += 5
329
            
330
            collision = 0
331

  
332
            for i in range(0,nb_atom_mol) :
333
                atom_test = molecule[i].position
334
                z_atom_test = atom_test[2]
335
                z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
336
                if (z_atom_test_final < z_atom_surf_max):
337
                    collision = 1
338
                if coads_box == "true":
339
                    collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)
340
                    collision += collision_box
341
        if collision == 0: corrected_collision = True
342

  
343
    if corrected_collision == True:  print 'Collision corrected'
344
    if collision > 0:
345
        print 'Error: the collision could not be corrected'
346
        pos1 = output.rfind('/')
347
        pos = pos1 + 1
348
        num_coll=output[pos:]
349
        j = output[:pos1]
350
        pos2 = j.rfind('/')
351
        molecule_directory = j[:pos2]
352
        fichier = open("%s/errors" % molecule_directory , "a")
353
        fichier.write("Adsorption %s : ERROR the collision could not be corrected\n" % num_coll )
354
        fichier.close()
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()
355 217

  
356 218

  
357 219
##########################################
358 220
# Adsorption de la molecule sur la surface
359 221
##########################################
360 222

  
361
if coads_box == "true":
362
    surf_atom_z=-max_float
363
    for atom in surface[:first_atom_ads1]:
364
        if atom.position[2] > surf_atom_z: surf_atom_z = atom.position[2]        
365
    
366
    distance = distance - z_max_ads_1 + surf_atom_z
367

  
368
new_add_adsorbate(surface, molecule, distance, (x_atom_surf_1,y_atom_surf_1,z_atom_surf_1), mol_index=atom_molecule)
369

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

  
373 226

  
374 227
############################
375 228
# Dissociation si necessaire
376 229
############################
377 230

  
378
if collision == 0 :
379
    surface = read(surface_file)
231
if not collision:
380 232
    molecule_reminder = molecule
381 233
    
382
    list_cutoffs = utils.natural_cutoffs(molecule)
234
    list_cutoffs = natural_cutoffs(molecule)
383 235
    nl=NeighborList(list_cutoffs, self_interaction=False, bothways=True)
384 236

  
385 237
    nl.update(molecule)
......
401 253
            print "Dissociation possible"
402 254

  
403 255
    if diss == 1:
404
        list_cutoffs_surface = utils.natural_cutoffs(surface)
256
        list_cutoffs_surface = natural_cutoffs(surface)
405 257
        nl2=NeighborList(list_cutoffs_surface, self_interaction=False, bothways=True)
406 258
        nl2.update(surface)
407 259
        neighbor_indices_surf, offsets_surf = nl2.get_neighbors(atom_surface)
408 260
        symbols_surf = surface.get_chemical_symbols()
409 261
        for i in neighbor_indices_surf:
410
            surface = read(surface_file)
262
            surface = ase.io.read(surface_file)
411 263
            molecule = molecule_reminder
412 264
            neighbor_surf_coord = surface[i].position
413 265
            coord_H_atom=molecule[H_atom].position
......
422 274
                molecule.translate(translation_matrix_dissociation)
423 275
                new_add_adsorbate(surface, molecule, distance, (x_atom_surf_1,y_atom_surf_1,z_atom_surf_1), mol_index=atom_molecule)
424 276
                output_diss = output + "_diss_" + str(i) + ".xyz"
425
                write(output_diss, surface)
426
    
427

  
428

  
429

  
430

  
277
                ase.io.write(output_diss, surface)
431 278

  

Also available in: Unified diff