Revision 36f722bc

b/dockOnSurf.sh
5 5
only_ads=false
6 6
only_iso=false
7 7
only_ref=false
8
coads_box=false
8 9

  
9 10
for i in $@; do
10 11
  case $i in
......
18 19
      only_iso=true;;
19 20
    -r|--only-refinement)
20 21
      only_ref=true;;
22
    -b|--coadsorption-box)
23
      coads_box=true;;
21 24
  esac
22 25
done
23 26

  
......
27 30
  -q=NUM, --max-qw=NUM \t \t allows a maximum of NUM jobs on 'qw' status.
28 31
  -a, --only-ads \t \t carry out only adsorption part.
29 32
  -i, --only-isolated \t \t carry out only isolated molecules part.
30
  -r, --only-refinement \t carry out only refinement of adsorbed structure."
33
  -r, --only-refinement \t carry out only refinement of adsorbed structure.
34
  -b, --coadsorption-box \t when the surface file already contains an
35
\t \t \t \t adsorbed molecule, it creates a box around it to
36
\t \t \t \t prevent the second adsorbate to overlap the first."
31 37
  exit
32 38
fi
33 39

  
......
129 135
  	exit
130 136
  fi
131 137
  read -a list_atom_surf -p 'List of the atoms of the surface that can be adsorption sites '
138
	if [ "$coads_box" = true ] ; then
139
		read -a first_atom_ads1 -p 'Number of the adsorbate first atom in the surface file '
140
	fi
132 141
  echo '' 
133 142
  
134 143
  ### Adsorption de la molecule sur la surface
135 144
  
136
  ${DockOnSurf_path}/modules/script_add_adsorbate+diss.sh $molecule $surface "${list_atom_1_mol[*]}" "${list_atom_2_mol[*]}" "${list_atom_surf[*]}"
145
  ${DockOnSurf_path}/modules/script_add_adsorbate+diss.sh $molecule $surface "${list_atom_1_mol[*]}" "${list_atom_2_mol[*]}" "${list_atom_surf[*]}" $coads_box $first_atom_ads1
137 146
  
138 147
  echo ' -- ADSORPTION SUTRUCTURES HAVE BEEN CALCULATED --'
139 148
  
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"
b/modules/script_add_adsorbate+diss.sh
9 9
list_atom_1_mol=($3)
10 10
list_atom_2_mol=($4)
11 11
list_atom_surf=($5)
12
coads_box=$6
13
first_atom_ads1=$7
12 14

  
13 15
nombre_atom_mol=${#list_atom_1_mol[@]}
14 16

  
......
68 70
			angle_phi=${angles_phi[$k]}
69 71
			angle_theta=${angles_theta[$k]}
70 72
			angle_psi=${angles_psi[$k]}
71
			${DockOnSurf_path}/modules/add_adsorbate_euler+dissociation.py ${Surface_path}/${surface_d_adsorption}.xyz ${atom_surf} ${Molecule_results_path}/${nom_de_la_molecule}/${molecule_a_adsorber}/last_geo.xyz ${atom_1_mol} ${distance_mol_surf} ${atom_2_mol} ${angle_phi} ${angle_theta} ${angle_psi} ${MolOnSurf_results_path}/${nom_de_la_molecule}/coords/${nom_de_la_molecule}_$b
73
			${DockOnSurf_path}/modules/add_adsorbate_euler+dissociation.py ${Surface_path}/${surface_d_adsorption}.xyz ${atom_surf} ${Molecule_results_path}/${nom_de_la_molecule}/${molecule_a_adsorber}/last_geo.xyz ${atom_1_mol} ${distance_mol_surf} ${atom_2_mol} ${angle_phi} ${angle_theta} ${angle_psi} ${MolOnSurf_results_path}/${nom_de_la_molecule}/coords/${nom_de_la_molecule}_$b $coads_box $first_atom_ads1
72 74
			b=$((b+1))
73 75
		done
74 76
	done

Also available in: Unified diff