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 arraylike 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 arraylike 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 coadsorptionprohibited 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