Chimie Théorique » scripts_chimie4psmn » DockOnSurf

coads_box=false

    -b|--coadsorption-box)

      coads_box=true;;

  -b, --coadsorption-box \t When the surface file already contains an

\t \t \t \t adsorbed molecule, it creates a box around it to

\t \t \t \t prevent the second adsorbate to overlap the first.

	if [ "$coads_box" = true ] ; then

		read -a first_atom_ads1 -p 'Number of the adsorbate first atom in the surface file '

	fi

  ${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

from ase import io

from ase.io import read

from ase.io import write

from ase import build

from ase.build import add_adsorbate

from ase import constraints

from ase.constraints import FixAtoms

from ase import neighborlist

from ase.neighborlist import NeighborList

from ase import utils

if len(sys.argv) < 7 :

    print "This script has to be used with at least 6 arguments."

    print "1. Name of the file containing the surface"

    print "2. Atom number of the atom of the surface where the adsorbate will be adsorbed"

    print "3. Name of the file containing the molecule to be adsorbed"

    print "4. Atom number of the atom of the molecule that will be adsorbed of the surface"

    print "5. Distance between the molecule and the surface"

    print "6. Name of the output file"

    print "Other arguments can be added before the name of the output file :"

    print "  Second atom of the molecule"

    print "  phi angle (rotation around z)"

    print "  theta angle (rotation around new x)"

    print "  psi angle (rotation around new z) "

    sys.exit(1)

###########################################################

# Definition of PBC conditions for molecules on the surface

###########################################################

def point_in_cell(p,h_matrix):

    # This function checks wether a given point defined by [x,y,z] is inside a 

    # given cell defined by h_matrix=[a,b,c] where a=[xa,ya,za], b=[xb,yb,zb], 

    # c=[xc,yc,zc]. If not it translates the point inside the cell by adding 

    # integer times the cell vectors, this way preserving its relative position

    # respect the cell copy. It basically applies periodic boundary conditions.

    # The check is performed by considering that a point p lies inside the cell, 

    # when it can be expressed as:

    # p = alpha*a + beta*b + gamma*c,   with:  0 < alpha, beta, gamma < 1.

    # alpha, beta and gamma are calculated as the determinant of the matrix 

    # defined by substituting the cell component whose coefficient we want to 

    # calculate by p, everything divided by the cell volume (det[a,b,c]) 

    # e.g: alpha = det([p,b,c])/det([a,b,c]); beta = det([a,p,c])/det([a,b,c])

    # gamma = det([a,b,p])/det([a,b,c])

    #                                                   Carles Martí 12/12/2019

    try: assert len(p) == 3, "p has not 3 elements"

    except TypeError: raise TypeError("p is not an array-like variable")

    try: [float(coord) for coord in p]

    except ValueError: raise ValueError("The elements of p are not numbers")

    try: assert len(h_matrix) == 3, "h_matrix has not 3 elements"

    except ValueError: raise TypeError("h_matrix is not an array-like variable")

    try: [[float(coord) for coord in v] for v in h_matrix],

    except ValueError: raise ValueError("The elements of h_matrix are not numbers")

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

    h_matrix = np.array(h_matrix)

    p = np.array(p)

    volume = np.linalg.det(h_matrix)

    alpha = np.linalg.det(np.array([p, h_matrix[1], h_matrix[2]])) / volume

    beta = np.linalg.det(np.array([h_matrix[0], p, h_matrix[2]])) / volume

    gamma = np.linalg.det(np.array([h_matrix[0], h_matrix[1], p])) / volume

    p_prime = np.array(p)

    if 0 > alpha: p_prime = p_prime +  (1 + int(-alpha)) * h_matrix[0]

    if 1 < alpha: p_prime = p_prime + int(-alpha) * h_matrix[0]

    if 0 > beta: p_prime = p_prime +  (1 + int(-beta)) * h_matrix[1]

    if 1 < beta: p_prime = p_prime + int(-beta) * h_matrix[1]

    if 0 > gamma: p_prime = p_prime +  (1 + int(-gamma)) * h_matrix[2]

    if 1 < gamma: p_prime = p_prime + int(-gamma) * h_matrix[2]

    return p_prime,[alpha,beta,gamma]

def check_coads_collision(position, ads_site_pos, ads_atom_pos, box_low_lim, box_high_lim, h_matrix):

    # This routine checks wether an atom of the second adsorbate collides

    # with the box around the first molecule.           Carles Martí 12/12/2019

    collision=0

    # Define the position of the atom in the surface coordinates.

    x_i = position[0] - ads_atom_pos[0] + ads_site_pos[0]

    y_i = position[1] - ads_atom_pos[1] + ads_site_pos[1]

    z_i = position[2] - ads_atom_pos[2] + ads_site_pos[2] + 1

    x_i = point_in_cell([x_i,y_i,z_i],h_matrix)[0][0]

    y_i = point_in_cell([x_i,y_i,z_i],h_matrix)[0][1]

    z_i = point_in_cell([x_i, y_i,z_i],h_matrix)[0][2]

    # Check if the atom lies inside the coadsorption-prohibited box

    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]]):

        collision = 1

    return collision

if len(sys.argv) < 8 : output = sys.argv[6]

surface = read(surface_file)

molecule = read(molecule_file)

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")]]

h_matrix=[[float(coord) for coord in v] for v in h_matrix]

#######################################################

# Cas avec un plus grand nombre d'arguments (rotations)

#######################################################

if len(sys.argv) > 7 :

    #######################

    # Lecture des arguments

    #######################

    atom_molecule_2 = int(sys.argv[6])

    phi = int(sys.argv[7])

    theta = int(sys.argv[8])

    psi = int(sys.argv[9])

    output = sys.argv[10]

    coads_box = sys.argv[11]

    if coads_box == "true": first_atom_ads1 = int(sys.argv[12])

    theta_min = theta - 0.01

    theta_max = theta + 0.01

    phi_min = phi - 0.01

    phi_max = phi + 0.01

    psi_min = psi - 0.01

    psi_max = psi + 0.01

    atom_mol_2 = molecule[atom_molecule_2].position

    x_atom_mol_2 = atom_mol_2[0]

    y_atom_mol_2 = atom_mol_2[1]

    z_atom_mol_2 = atom_mol_2[2]

    # In the case of coadsorption, definition of a box around the first

    # adsorbate to prevent the second one to overlap with the former.

    if coads_box == "true":        

        max_float=sys.float_info.max

        box_low_lim=[max_float,max_float,max_float]

        box_high_lim=[-max_float,-max_float,-max_float]

        for adsorbate_atom in surface[first_atom_ads1:]:

            if adsorbate_atom.position[0] < box_low_lim[0]: box_low_lim[0] = adsorbate_atom.position[0]

            if adsorbate_atom.position[1] < box_low_lim[1]: box_low_lim[1] = adsorbate_atom.position[1]

            if adsorbate_atom.position[2] < box_low_lim[2]: box_low_lim[2] = adsorbate_atom.position[2]

            if adsorbate_atom.position[0] > box_high_lim[0]: box_high_lim[0] = adsorbate_atom.position[0]

            if adsorbate_atom.position[1] > box_high_lim[1]: box_high_lim[1] = adsorbate_atom.position[1]

            if adsorbate_atom.position[2] > box_high_lim[2]: box_high_lim[2] = adsorbate_atom.position[2]

        # Adding an offset to the limit

        z_max_ads_1 = box_high_lim[2]

        box_low_lim = [low_lim - 1 for low_lim in box_low_lim]

        box_high_lim = [ high_lim + 1 for high_lim in box_high_lim]

    ##################################

    # Rotation avec les angles d'Euler

    ##################################

    #Initialisation : placement du second atome de la molecule sur Oy

    vect_y = [0, 1, 0]

    vect_atom_mol_2 = molecule[atom_molecule_2].position - molecule[atom_molecule].position

    angle_1 = get_proper_angle(vect_y, vect_atom_mol_2)

    if (angle_1 != 0) :

        vect_normal = np.cross(vect_y, vect_atom_mol_2)

        molecule.rotate(-angle_1, vect_normal, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))

        vect_atom_mol_2_verif = molecule[atom_molecule_2].position - molecule[atom_molecule].position

        angle_1_verif = get_proper_angle(vect_atom_mol_2_verif, vect_y)

        angle_max = 0.01

        angle_min = -0.01

        if (angle_1_verif < angle_min or angle_1_verif > angle_max) :

            print 'Error in initialisation'

    #Rotation Euler

    molecule.euler_rotate(phi, theta, psi, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))

    #Correction des collisions entre la molécule et la surface

    z_atom_surf_max = z_atom_surf_1 + 1

    collision = 0

    rotation_tot = 0

    nb_atom_collision_min = 0

    rotation_opt = 0

    rotation_modif = 5

    corrected_collision=False

    for i in range(0,nb_atom_mol):

        atom_test = molecule[i].position

        z_atom_test = atom_test[2]

        z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance

        if (z_atom_test_final < z_atom_surf_max):

            collision = 1

        if coads_box == "true":

            collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)

            collision += collision_box

    if collision > 0 :

        print 'Collision between the molecule and the surface - modification of theta angle'

    vect_z = [0, 0, 1]

    while (collision > 0 and rotation_tot < 354) :

        vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position

        if (vect_atom_mol_2_post_euler[0] != 0 and vect_atom_mol_2_post_euler[1] != 0) :

            vect_rotation_theta = np.cross(vect_atom_mol_2_post_euler, vect_z)

        else :

            vect_rotation_theta = [0.0001, 0.0001 , 0]

        molecule.rotate(rotation_modif, vect_rotation_theta, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))

        rotation_tot += 5

        collision = 0

        for i in range(0,nb_atom_mol) :

            atom_test = molecule[i].position

            z_atom_test = atom_test[2]

            z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance

            if (z_atom_test_final < z_atom_surf_max) :

                collision = 1

            if coads_box == "true":

                collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)

                collision += collision_box

        rotation_tot_2 = 0

        while (collision > 0 and rotation_tot_2 < 354) :

            vect_atom_mol_2_post_euler = molecule[atom_molecule_2].position - molecule[atom_molecule].position

            molecule.rotate(rotation_modif, vect_atom_mol_2_post_euler, center=(x_atom_mol_1,y_atom_mol_1,z_atom_mol_1))

            rotation_tot_2 += 5

            collision = 0

            for i in range(0,nb_atom_mol) :

                atom_test = molecule[i].position

                z_atom_test = atom_test[2]

                z_atom_test_final = z_atom_test - z_atom_mol_1 + z_atom_surf_1 + distance

                if (z_atom_test_final < z_atom_surf_max):

                    collision = 1

                if coads_box == "true":

                    collision_box = check_coads_collision(molecule[i].position, atom_surf_1, atom_mol_1, box_low_lim, box_high_lim, h_matrix)

                    collision += collision_box

        if collision == 0: corrected_collision = True

    if corrected_collision == True:  print 'Collision corrected'

    if collision > 0:

        print 'Error: the collision could not be corrected'

        pos1 = output.rfind('/')

        pos = pos1 + 1

        num_coll=output[pos:]

        j = output[:pos1]

        pos2 = j.rfind('/')

        molecule_directory = j[:pos2]

        fichier = open("%s/errors" % molecule_directory , "a")

        fichier.write("Adsorption %s : ERROR the collision could not be corrected\n" % num_coll )

        fichier.close()

if coads_box == "true":

    surf_atom_z=-max_float

    for atom in surface[:first_atom_ads1]:

        if atom.position[2] > surf_atom_z: surf_atom_z = atom.position[2]        

    distance = distance - z_max_ads_1 + surf_atom_z

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

write(out, surface)

if collision == 0 :

    surface = read(surface_file)

    list_cutoffs = utils.natural_cutoffs(molecule)

        list_cutoffs_surface = utils.natural_cutoffs(surface)

            surface = read(surface_file)

                write(output_diss, surface)

distance_mol_surf=2

coads_box=$6

first_atom_ads1=$7

			${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