Revision 36f722bc modules/add_adsorbate_euler+dissociation.py

b/modules/add_adsorbate_euler+dissociation.py
2 2
# -*- coding: utf-8 -*-
3 3

  
4 4
import sys
5
import os
5 6
import numpy as np
6 7
import ase
7 8
from ase import io
......
42 43
    angle = np.arccos(np.sign(norm_dot) if np.abs(norm_dot) >= 1 else norm_dot)
43 44
    return(angle*180/np.pi if degrees else angle)
44 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
45 112

  
113
    return collision
114
    
46 115
##########################################
47 116
# Lecture des différents arguments de base
48 117
##########################################
......
51 120
atom_surface = int(sys.argv[2])
52 121
molecule_file = sys.argv[3]
53 122
atom_molecule = int(sys.argv[4])
54
distance = int(sys.argv[5])
55
if len(sys.argv) < 8 :
56
    output = sys.argv[6]
123
distance = float(sys.argv[5])
124
if len(sys.argv) < 8 : output = sys.argv[6]
57 125

  
58 126
surface = read(surface_file)
59 127
molecule = read(molecule_file)
......
74 142

  
75 143
nb_atom_mol = int(len(molecule))
76 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]
77 147

  
78 148
#######################################################
79 149
# Cas avec un plus grand nombre d'arguments (rotations)
......
91 161
    theta = int(sys.argv[8])
92 162
    psi = int(sys.argv[9])
93 163
    output = sys.argv[10]
164
    coads_box = sys.argv[11]
165
    if coads_box == "true": first_atom_ads1 = int(sys.argv[12])
94 166

  
95 167
    theta_min = theta - 0.01
96 168
    theta_max = theta + 0.01
......
104 176
    y_atom_mol_2 = atom_mol_2[1]
105 177
    z_atom_mol_2 = atom_mol_2[2]
106 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

  
107 200
    ##################################
108 201
    # Rotation avec les angles d'Euler
109 202
    ##################################
......
128 221
    #Correction des collisions entre la molécule et la surface
129 222
    z_atom_surf_max = z_atom_surf_1 + 1
130 223

  
131
    collision = 0.5
224
    collision = 0
132 225
    rotation_tot = 0
133 226
    nb_atom_collision_min = 0
134 227
    rotation_opt = 0
135 228
    rotation_modif = 5
229
    corrected_collision=False
136 230
    
137
    for i in range(0,nb_atom_mol) :
231
    for i in range(0,nb_atom_mol):
138 232
        atom_test = molecule[i].position
139 233
        z_atom_test = atom_test[2]
140 234
        z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
141
        if (z_atom_test_final < z_atom_surf_max) :
235
        if (z_atom_test_final < z_atom_surf_max):
142 236
            collision = 1
143
            nb_atom_collision_min += 1
144
    if 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 :
145 242
        print 'Collision between the molecule and the surface - modification of theta angle'
146 243

  
147 244
    vect_z = [0, 0, 1]
148 245

  
149
    while (collision == 1 and rotation_tot < 354) :
246
    while (collision > 0 and rotation_tot < 354) :
150 247
        vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position
151 248
        
152 249
        if (vect_atom_mol_2_post_euler[0] != 0 and vect_atom_mol_2_post_euler[1] != 0) :
......
157 254
        rotation_tot += 5
158 255
        
159 256
        collision = 0
160
        nb_atom_collision = 0
161

  
162 257
        for i in range(0,nb_atom_mol) :
163 258
            atom_test = molecule[i].position
164 259
            z_atom_test = atom_test[2]
165 260
            z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
166 261
            if (z_atom_test_final < z_atom_surf_max) :
167
                nb_atom_collision += 1
168 262
                collision = 1
169
        if (nb_atom_collision < nb_atom_collision_min) :
170
            nb_atom_collision_min = nb_atom_collision
171
            rotation_opt = rotation_tot
172

  
173
        if nb_atom_collision_min != 0 :
174
            rotation_tot_2 = 0
175
            rotation_opt_2 = 0
176

  
177
            while (collision == 1 and rotation_tot_2 < 354) :
178
                vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position
179
                molecule.rotate(rotation_modif, vect_atom_mol_2_post_euler, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))
180
                rotation_tot_2 += 5
181
                
182
                collision = 0
183
                nb_atom_collision = 0
184

  
185
                for i in range(0,nb_atom_mol) :
186
                    atom_test = molecule[i].position
187
                    z_atom_test = atom_test[2]
188
                    z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance
189
                    if (z_atom_test_final < z_atom_surf_max) :
190
                        nb_atom_collision += 1
191
                        collision = 1
192
                if (nb_atom_collision < nb_atom_collision_min) :
193
                    nb_atom_collision_min = nb_atom_collision
194
                    rotation_opt_2 = rotation_tot_2
195

  
196
    if collision == 0 :
197
        print 'Collision corrected'
198
    elif 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:
199 289
        print 'Error: the collision could not be corrected'
200 290
        pos1 = output.rfind('/')
201 291
        pos = pos1 + 1
......
212 302
# Adsorption de la molecule sur la surface
213 303
##########################################
214 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

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

  
217 314
out=output+".xyz"

Also available in: Unified diff