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