Chimie Théorique » scripts_chimie4psmn » DockOnSurf

import logging

import numpy as np

import ase

logger = logging.getLogger('DockOnSurf')

def assign_prop(atoms: ase.Atoms, prop_name: str, prop_val):  # TODO Needed?

    atoms.info[prop_name] = prop_val

def select_confs(orig_conf_list: list, calc_dirs: list, magns: list,

                 num_sel: int, code: str):

    """Takes a list ase.Atoms and selects the most different magnitude-wise.

    Given a list of ase.Atoms objects and a list of magnitudes, it selects a

    number of the most different conformers according to every magnitude

    specified.

    @param orig_conf_list: list of ase.Atoms objects to select among.

    @param calc_dirs: List of directories where to read the energies from.

    @param magns: list of str with the names of the magnitudes to use for the

        conformer selection. Supported magnitudes: 'energy', 'moi'.

    @param num_sel: number of conformers to select for every of the magnitudes.

    @param code: The code that generated the magnitude information.

         Supported codes: See formats.py

    @return: list of the selected ase.Atoms objects.

    """

    from copy import deepcopy

    from modules.formats import collect_energies

    conf_list = deepcopy(orig_conf_list)

    conf_enrgs, mois, selected_ids = [], [], []

    if num_sel >= len(conf_list):

        logger.warning('Number of conformers per magnitude is equal or larger '

                       'than the total number of conformers. Using all '

                       f'available conformers: {len(conf_list)}.')

        return conf_list

    # Read properties

    if 'energy' in magns:

        if code == 'cp2k':

            conf_enrgs = collect_energies(calc_dirs, code, 'isolated')

        elif code == 'vasp':

            conf_enrgs = np.array([conf.get_total_energy()

                                   for conf in orig_conf_list])

    if 'moi' in magns:

        mois = np.array([conf.get_moments_of_inertia() for conf in conf_list])

    # Assign values

    for i, conf in enumerate(conf_list):

        assign_prop(conf, 'idx', i)

        assign_prop(conf, 'iso', calc_dirs[i])

        if 'energy' in magns:

            assign_prop(conf, 'energy', conf_enrgs[i])

        if 'moi' in magns:

            assign_prop(conf, 'moi', mois[i, 2])

    # pick ids

    for magn in magns:

        sorted_list = sorted(conf_list, key=lambda conf: abs(conf.info[magn]))

        if sorted_list[-1].info['idx'] not in selected_ids:

            selected_ids.append(sorted_list[-1].info['idx'])

        if num_sel > 1:

            for i in range(0, len(sorted_list) - 1,

                           len(conf_list) // (num_sel - 1)):

                if sorted_list[i].info['idx'] not in selected_ids:

                    selected_ids.append(sorted_list[i].info['idx'])

    logger.info(f'Selected {len(selected_ids)} conformers for adsorption.')

    return [conf_list[idx] for idx in selected_ids]

def get_vect_angle(v1: list, v2: list, ref=None, degrees=True):

    """Computes the angle between two vectors.

    @param v1: The first vector.

    @param v2: The second vector.

    @param ref: Orthogonal vector to both v1 and v2,

        along which the sign of the rotation is defined (i.e. positive if

        counterclockwise angle when facing ref)

    @param degrees: Whether the result should be in radians (True) or in

        degrees (False).

    @return: The angle in radians if degrees = False, or in degrees if

        degrees =True

    """

    v1_u = v1 / np.linalg.norm(v1)

    v2_u = v2 / np.linalg.norm(v2)

    angle = np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

    if ref is not None:

        # Give sign according to ref direction

        angle *= (1 if np.dot(np.cross(v1, v2), ref) >= 0 else -1)

    return angle if not degrees else angle * 180 / np.pi

def vect_avg(vects):

    """Computes the element-wise mean of a set of vectors.

    @param vects: list of lists-like: containing the vectors (num_vectors,

        length_vector).

    @return: vector average computed doing the element-wise mean.

    """

    from modules.utilities import try_command

    err = "vect_avg parameter vects must be a list-like, able to be converted" \

          " np.array"

    array = try_command(np.array, [(ValueError, err)], vects)

    if len(array.shape) == 1:

        return array

    else:

        num_vects = array.shape[1]

        return np.array([np.average(array[:, i]) for i in range(num_vects)])

def get_atom_coords(atoms: ase.Atoms, center=None):

    """Gets the coordinates of the specified center for an ase.Atoms object.

    If center is not an index but a list of indices, it computes the

    element-wise mean of the coordinates of the atoms specified in the inner

    list.

    @param atoms: ase.Atoms object for which to obtain the coordinates of.

    @param center: index/list of indices of the atoms for which the coordinates

                   should be extracted.

    @return: np.ndarray of atomic coordinates.

    """

    err_msg = "Argument 'ctr' must be an integer or a list of integers. "\

              "Every integer must be in the range [0, num_atoms)"

    if center is None:

        center = list(range(len(atoms)))

    if isinstance(center, int):

        if center not in list(range(len(atoms))):

            logger.error(err_msg)

            raise ValueError(err_msg)

        return atoms[center].position

    elif isinstance(center, list):

        for elm in center:

            if elm not in list(range(len(atoms))):

                logger.error(err_msg)

                raise ValueError(err_msg)

        return vect_avg([atoms[idx].position for idx in center])

    else:

        logger.error(err_msg)

        raise ValueError(err_msg)

def compute_norm_vect(atoms, idxs, cell):

    """Computes the local normal vector of a surface at a given site.

    Given an ase.Atoms object and a site defined as a linear combination of

    atoms it computes the vector perpendicular to the surface, considering the

    local environment of the site.

    @param atoms: ase.Atoms object of the surface.

    @param idxs: list or int: Index or list of indices of the atom/s that define

        the site

    @param cell: Unit cell. A 3x3 matrix (the three unit cell vectors)

    @return: numpy.ndarray of the coordinates of the vector locally

    perpendicular to the surface.

    """

    from modules.ASANN import coordination_numbers as coord_nums

    if isinstance(idxs, list):

        atm_vect = [-np.round(coord_nums(atoms.get_scaled_positions(),

                                         pbc=np.any(cell),

                                         cell_vectors=cell)[3][i], 2)

                    for i in idxs]

        norm_vec = vect_avg([vect / np.linalg.norm(vect) for vect in atm_vect])

    elif isinstance(idxs, int):

        norm_vec = -coord_nums(atoms.get_scaled_positions(),

                               pbc=np.any(cell),

                               cell_vectors=cell)[3][idxs]

    else:

        err = "'idxs' must be either an int or a list"

        logger.error(err)

        raise ValueError(err)

    norm_vec = np.round(norm_vec, 2) / np.linalg.norm(np.round(norm_vec, 2))

    logger.info(f"The perpendicular vector to the surface at site '{idxs}' is "

                f"{norm_vec}")

    return norm_vec

def align_molec(orig_molec, ctr_coord, ref_vect):

    """Align a molecule to a vector by a center.

    Given a reference vector to be aligned to and some coordinates acting as

    alignment center, it first averages the vectors pointing to neighboring

    atoms and then tries to align this average vector to the target. If the

    average vector is 0 it takes the vector to the nearest neighbor.

    @param orig_molec: The molecule to align.

    @param ctr_coord: The coordinates to use ase alignment center.

    @param ref_vect: The vector to be aligned with.

    @return: ase.Atoms of the aligned molecule.

    """

    from copy import deepcopy

    from ase import Atom

    from ase.neighborlist import natural_cutoffs, neighbor_list

    molec = deepcopy(orig_molec)

    if len(molec) == 1:

        err_msg = "Cannot align a single atom"

        logger.error(err_msg)

        ValueError(err_msg)

    cutoffs = natural_cutoffs(molec, mult=1.2)

    # Check if ctr_coord are the coordinates of an atom and if not creates a

    # dummy one to extract the neighboring atoms.

    ctr_idx = None

    dummy_atom = False

    for atom in molec:

        if np.allclose(ctr_coord, atom.position, rtol=1e-2):

            ctr_idx = atom.index

            break

    if ctr_idx is None:

        molec.append(Atom('X', position=ctr_coord))

        cutoffs.append(0.2)

        ctr_idx = len(molec) - 1

        dummy_atom = True

    # Builds the neighbors and computes the average vector

    refs, vects = neighbor_list("iD", molec, cutoffs, self_interaction=False)

    neigh_vects = [vects[i] for i, atm in enumerate(refs) if atm == ctr_idx]

    # If no neighbors are present, the cutoff of the alignment center is

    # set to a value where at least one atom is a neighbor and neighbors are

    # recalculated.

    if len(neigh_vects) == 0:

        min_dist, min_idx = (np.inf, np.inf)

        for atom in molec:

            if atom.index == ctr_idx:

                continue

            if molec.get_distance(ctr_idx, atom.index) < min_dist:

                min_dist = molec.get_distance(ctr_idx, atom.index)

                min_idx = atom.index

        cutoffs[ctr_idx] = min_dist - cutoffs[min_idx] + 0.05

        refs, vects = neighbor_list("iD", molec, cutoffs,

                                    self_interaction=False)

        neigh_vects = [vects[i] for i, atm in enumerate(refs) if atm == ctr_idx]

    target_vect = vect_avg(neigh_vects)

    # If the target vector is 0 (the center is at the baricenter of its

    # neighbors). Assuming the adsorption center is coplanar or colinear to its

    # neighbors (it would not make a lot of sense to chose a center which is

    # the baricenter of neighbors distributed in 3D), the target_vector is

    # chosen perpendicular to the nearest neighbor.

    if np.allclose(target_vect, 0, rtol=1e-3):

        nn_vect = np.array([np.inf] * 3)

        for vect in neigh_vects:

            if np.linalg.norm(vect) < np.linalg.norm(nn_vect):

                nn_vect = vect

        cart_axes = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

        axis = cart_axes[int(np.argmax([np.linalg.norm(np.cross(ax, nn_vect))

                                        for ax in cart_axes]))]

        target_vect = np.cross(axis, nn_vect)

    rot_vect = np.cross(target_vect, ref_vect)

    if np.allclose(rot_vect, 0):

        cart_axes = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

        axis = cart_axes[int(np.argmax([np.linalg.norm(np.cross(ax, ref_vect))

                                        for ax in cart_axes]))]

        rot_vect = np.cross(ref_vect, axis)

    rot_angle = -get_vect_angle(ref_vect, target_vect, rot_vect)

    molec.rotate(rot_angle, rot_vect, ctr_coord)

    if dummy_atom:

        del molec[-1]

    return molec

def add_adsorbate(slab, adsorbate, site_coord, ctr_coord, height, offset=None,

                  norm_vect=(0, 0, 1)):

    """Add an adsorbate to a surface.

    This function extends the functionality of ase.build.add_adsorbate

    (https://wiki.fysik.dtu.dk/ase/ase/build/surface.html#ase.build.add_adsorbate)

    by enabling to change the z coordinate and the axis perpendicular to the

    surface.

    @param slab: ase.Atoms object containing the coordinates of the surface

    @param adsorbate: ase.Atoms object containing the coordinates of the

        adsorbate.

    @param site_coord: The coordinates of the adsorption site on the surface.

    @param ctr_coord: The coordinates of the adsorption center in the molecule.

    @param height: The height above the surface where to adsorb.

    @param offset: Offsets the adsorbate by a number of unit cells. Mostly

        useful when adding more than one adsorbate.

    @param norm_vect: The vector perpendicular to the surface.

    """

    from copy import deepcopy

    info = slab.info.get('adsorbate_info', {})

    pos = np.array([0.0, 0.0, 0.0])  # part of absolute coordinates

    spos = np.array([0.0, 0.0, 0.0])  # part relative to unit cell

    norm_vect_u = np.array(norm_vect) / np.linalg.norm(norm_vect)

    if offset is not None:

        spos += np.asarray(offset, float)

    if isinstance(site_coord, str):

        # A site-name:

        if 'sites' not in info:

            raise TypeError('If the atoms are not made by an ase.build '

                            'function, position cannot be a name.')

        if site_coord not in info['sites']:

            raise TypeError('Adsorption site %s not supported.' % site_coord)

        spos += info['sites'][site_coord]

    else:

        pos += site_coord

    if 'cell' in info:

        cell = info['cell']

    else:

        cell = slab.get_cell()

    pos += np.dot(spos, cell)

    # Convert the adsorbate to an Atoms object

    if isinstance(adsorbate, ase.Atoms):

        ads = deepcopy(adsorbate)

    elif isinstance(adsorbate, ase.Atom):

        ads = ase.Atoms([adsorbate])

    else:

        # Assume it is a string representing a single Atom

        ads = ase.Atoms([ase.Atom(adsorbate)])

    pos += height * norm_vect_u

    # Move adsorbate into position

    ads.translate(pos - ctr_coord)

    # Attach the adsorbate

    slab.extend(ads)

def check_collision(slab_molec, slab_num_atoms, min_height, vect, nn_slab=0,

                    nn_molec=0, coll_coeff=1.2, exclude_atom=False):

    """Checks whether a slab and a molecule collide or not.

    @param slab_molec: The system of adsorbate-slab for which to detect if there

        are collisions.

    @param nn_slab: Number of neigbors in the surface.

    @param nn_molec: Number of neighbors in the molecule.

    @param coll_coeff: The coefficient that multiplies the covalent radius of

        atoms resulting in a distance that two atoms being closer to that is

        considered as atomic collision.

    @param slab_num_atoms: Number of atoms of the bare slab.

    @param min_height: The minimum height atoms can have to not be considered as

        colliding.

    @param vect: The vector perpendicular to the slab.

    @param exclude_atom: Whether to exclude the adsorption center in the

        molecule in the collision detection.

    @return: bool, whether the surface and the molecule collide.

    """

    from copy import deepcopy

    from ase.neighborlist import natural_cutoffs, neighbor_list

    # Check structure overlap by height

    if min_height is not False:

        cart_axes = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0],

                     [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]]

        if vect.tolist() not in cart_axes:

            err_msg = "'min_coll_height' option is only implemented for " \

                      "'surf_norm_vect' to be one of the x, y or z axes. "

            logger.error(err_msg)

            raise ValueError(err_msg)

        for atom in slab_molec[slab_num_atoms:]:

            if exclude_atom is not False \

                    and atom.index == exclude_atom:

                continue

            for i, coord in enumerate(vect):

                if coord == 0:

                    continue

                if atom.position[i] * coord < min_height * coord:

                    return True

    # Check structure overlap by sphere collision

    if coll_coeff is not False:

        if exclude_atom is not False:

            slab_molec_wo_ctr = deepcopy(slab_molec)

            del slab_molec_wo_ctr[exclude_atom + slab_num_atoms]

            slab_molec_cutoffs = natural_cutoffs(slab_molec_wo_ctr,

                                                 mult=coll_coeff)

            slab_molec_nghbs = len(neighbor_list("i", slab_molec_wo_ctr,

                                                 slab_molec_cutoffs))

        else:

            slab_molec_cutoffs = natural_cutoffs(slab_molec, mult=coll_coeff)

            slab_molec_nghbs = len(neighbor_list("i", slab_molec,

                                                 slab_molec_cutoffs))

        if slab_molec_nghbs > nn_slab + nn_molec:

            return True

    return False

def correct_coll(molec, slab, ctr_coord, site_coord, num_pts,

                 min_coll_height, norm_vect, slab_nghbs, molec_nghbs,

                 coll_coeff, height=2.5, excl_atom=False):

    # TODO Rename this function

    """Tries to adsorb a molecule on a slab trying to avoid collisions by doing

    small rotations.

    @param molec: ase.Atoms object of the molecule to adsorb

    @param slab: ase.Atoms object of the surface on which to adsorb the

        molecule

    @param ctr_coord: The coordinates of the molecule to use as adsorption

        center.

    @param site_coord: The coordinates of the surface on which to adsorb the

        molecule

    @param num_pts: Number on which to sample Euler angles.

    @param min_coll_height: The lowermost height for which to detect a collision

    @param norm_vect: The vector perpendicular to the surface.

    @param slab_nghbs: Number of neigbors in the surface.

    @param molec_nghbs: Number of neighbors in the molecule.

    @param coll_coeff: The coefficient that multiplies the covalent radius of

        atoms resulting in a distance that two atoms being closer to that is

        considered as atomic collision.

    @param height: Height on which to try adsorption.

    @param excl_atom: Whether to exclude the adsorption center in the

        molecule in the collision detection.

    @return collision: bool, whether the structure generated has collisions

        between slab and adsorbate.

    """

    from copy import deepcopy

    slab_num_atoms = len(slab)

    slab_molec = []

    collision = True

    max_corr = 6  # Should be an even number

    d_angle = 180 / ((max_corr / 2.0) * num_pts)

    num_corr = 0

    while collision and num_corr <= max_corr:

        k = num_corr * (-1) ** num_corr

        slab_molec = deepcopy(slab)

        molec.euler_rotate(k * d_angle, k * d_angle / 2, k * d_angle,

                           center=ctr_coord)

        add_adsorbate(slab_molec, molec, site_coord, ctr_coord, height,

                      norm_vect=norm_vect)

        collision = check_collision(slab_molec, slab_num_atoms, min_coll_height,

                                    norm_vect, slab_nghbs, molec_nghbs,

                                    coll_coeff, excl_atom)

        num_corr += 1

    return slab_molec, collision

def dissociate_h(slab_molec_orig, h_idx, num_atoms_slab, h_acceptor,

                 neigh_cutoff=1):

    # TODO rethink

    """Tries to dissociate a H from the molecule and adsorbs it on the slab.

    Tries to dissociate a H atom from the molecule and adsorb in on top of the

    surface if the distance is shorter than two times the neigh_cutoff value.

    @param slab_molec_orig: The ase.Atoms object of the system adsorbate-slab.

    @param h_idx: The index of the hydrogen atom to carry out adsorption of.

    @param num_atoms_slab: The number of atoms of the slab without adsorbate.

    @param h_acceptor: List of atom types or atom numbers that are H-acceptors.

    @param neigh_cutoff: half the maximum distance between the surface and the

        H for it to carry out dissociation.

    @return: An ase.Atoms object of the system adsorbate-surface with H

    """

    from copy import deepcopy

    from ase.neighborlist import NeighborList

    slab_molec = deepcopy(slab_molec_orig)

    cutoffs = len(slab_molec) * [neigh_cutoff]

    nl = NeighborList(cutoffs, self_interaction=False, bothways=True, skin=0.0)

    nl.update(slab_molec)

    surf_h_vect = np.array([np.infty] * 3)

    if h_acceptor == 'all':

        h_acceptor = list(range(num_atoms_slab))

    for neigh_idx in nl.get_neighbors(h_idx)[0]:

        for elm in h_acceptor:

            if isinstance(elm, int):

                if neigh_idx == elm and neigh_idx < num_atoms_slab:

                    dist = np.linalg.norm(slab_molec[neigh_idx].position -

                                          slab_molec[h_idx].position)

                    if dist < np.linalg.norm(surf_h_vect):

                        surf_h_vect = slab_molec[neigh_idx].position \

                                      - slab_molec[h_idx].position

            else:

                if slab_molec[neigh_idx].symbol == elm \

                        and neigh_idx < num_atoms_slab:

                    dist = np.linalg.norm(slab_molec[neigh_idx].position -

                                          slab_molec[h_idx].position)

                    if dist < np.linalg.norm(surf_h_vect):

                        surf_h_vect = slab_molec[neigh_idx].position \

                                      - slab_molec[h_idx].position

    if np.linalg.norm(surf_h_vect) != np.infty:

        trans_vect = surf_h_vect - surf_h_vect / np.linalg.norm(surf_h_vect)

        slab_molec[h_idx].position = slab_molec[h_idx].position + trans_vect

        return slab_molec

def dissociation(slab_molec, h_donor, h_acceptor, num_atoms_slab):

    # TODO multiple dissociation

    """Decides which H atoms to dissociate according to a list of atoms.

    Given a list of chemical symbols or atom indices it checks for every atom

    or any of its neighbor if it's a H and calls dissociate_h to try to carry

    out dissociation of that H. For atom indices, it checks both whether

    the atom index or its neighbors are H, for chemical symbols, it only checks

    if there is a neighbor H.

    @param slab_molec: The ase.Atoms object of the system adsorbate-slab.

    @param h_donor: List of atom types or atom numbers that are H-donors.

    @param h_acceptor: List of atom types or atom numbers that are H-acceptors.

    @param num_atoms_slab: Number of atoms of the bare slab.

    @return:

    """

    from ase.neighborlist import natural_cutoffs, NeighborList

    molec = slab_molec[num_atoms_slab:]

    cutoffs = natural_cutoffs(molec)

    nl = NeighborList(cutoffs, self_interaction=False, bothways=True)

    nl.update(molec)

    disso_structs = []

    for el in h_donor:

        if isinstance(el, int):

            if molec[el].symbol == 'H':

                disso_struct = dissociate_h(slab_molec, el + num_atoms_slab,

                                            num_atoms_slab, h_acceptor)

                if disso_struct is not None:

                    disso_structs.append(disso_struct)

            else:

                for neigh_idx in nl.get_neighbors(el)[0]:

                    if molec[neigh_idx].symbol == 'H':

                        disso_struct = dissociate_h(slab_molec, neigh_idx +

                                                    num_atoms_slab,

                                                    num_atoms_slab, h_acceptor)

                        if disso_struct is not None:

                            disso_structs.append(disso_struct)

        else:

            for atom in molec:

                if atom.symbol.lower() == el.lower():

                    for neigh_idx in nl.get_neighbors(atom.index)[0]:

                        if molec[neigh_idx].symbol == 'H':

                            disso_struct = dissociate_h(slab_molec, neigh_idx

                                                        + num_atoms_slab,

                                                        num_atoms_slab,

                                                        h_acceptor)

                            if disso_struct is not None:

                                disso_structs.append(disso_struct)

    return disso_structs

def ads_euler(orig_molec, slab, ctr_coord, site_coord, num_pts,

              min_coll_height, coll_coeff, norm_vect, slab_nghbs, molec_nghbs,

              h_donor, h_acceptor, height, excl_atom):

    """Generates adsorbate-surface structures by sampling over Euler angles.

    This function generates a number of adsorbate-surface structures at

    different orientations of the adsorbate sampled at multiple Euler (zxz)

    angles.

    @param orig_molec: ase.Atoms object of the molecule to adsorb.

    @param slab: ase.Atoms object of the surface on which to adsorb the

        molecule.

    @param ctr_coord: The coordinates of the molecule to use as adsorption

        center.

    @param site_coord: The coordinates of the surface on which to adsorb the

        molecule.

    @param num_pts: Number on which to sample Euler angles.

    @param min_coll_height: The lowest height for which to detect a collision.

    @param coll_coeff: The coefficient that multiplies the covalent radius of

        atoms resulting in a distance that two atoms being closer to that is

        considered as atomic collision.

    @param norm_vect: The vector perpendicular to the surface.

    @param slab_nghbs: Number of neigbors in the surface.

    @param molec_nghbs: Number of neighbors in the molecule.

    @param h_donor: List of atom types or atom numbers that are H-donors.

    @param h_acceptor: List of atom types or atom numbers that are H-acceptors.

    @param height: Height on which to try adsorption.

    @param excl_atom: Whether to exclude the adsorption center in the

        molecule in the collision detection.

    @return: list of ase.Atoms object conatining all the orientations of a given

        conformer.

    """

    from copy import deepcopy

    slab_ads_list = []

    prealigned_molec = align_molec(orig_molec, ctr_coord, norm_vect)

    # rotation around z

    for alpha in np.arange(0, 360, 360 / num_pts):

        # rotation around x'

        for beta in np.arange(0, 180, 180 / num_pts):

            # rotation around z'

            for gamma in np.arange(0, 360, 360 / num_pts):

                if beta == 0 and gamma > 0:

                    continue

                molec = deepcopy(prealigned_molec)

                molec.euler_rotate(alpha, beta, gamma, center=ctr_coord)

                slab_molec, collision = correct_coll(molec, slab, ctr_coord,

                                                     site_coord, num_pts,

                                                     min_coll_height, norm_vect,

                                                     slab_nghbs, molec_nghbs,

                                                     coll_coeff, height,

                                                     excl_atom)

                if not collision and not any([np.allclose(slab_molec.positions,

                                                          conf.positions)

                                              for conf in slab_ads_list]):

                    slab_ads_list.append(slab_molec)

                    if h_donor is not False:

                        slab_ads_list.extend(dissociation(slab_molec, h_donor,

                                                          h_acceptor,

                                                          len(slab)))

    return slab_ads_list

def internal_rotate(molecule, surf, ctr1_mol, ctr2_mol, ctr3_mol, ctr1_surf,

                    ctr2_surf, bond_vector, bond_angle_target,

                    dihedral_angle_target=None, mol_dihedral_angle_target=None):

    """Performs translation and rotation of an adsorbate defined by an

    adsorption bond length, direction, angle and dihedral angle

    Carles modification of chemcat's transform_adsorbate to work with

    coordinates instead of ase.Atom

    Parameters:

        molecule (ase.Atoms): The molecule to adsorb.

        surf (ase.Atoms): The surface ontop of which to adsorb.

        ctr1_mol (int/list): The position of the adsorption center in the

        molecule that will be bound to the surface.

        ctr2_mol (int/list): The position of a second center of the

        adsorbate used to define the adsorption bond angle, and the dihedral

        adsorption angle.

        ctr3_mol (int/list): The position of a third center in the

        adsorbate used to define the adsorbate dihedral angle.

        ctr1_surf (int/list): The position of the site on the surface that

        will be bound to the molecule.

        ctr2_surf (int/list): The position of a second center of the

        surface used to define the dihedral adsorption angle.

        bond_vector (numpy.ndarray): The adsorption bond desired.

            Details: offset = vect(atom1_surf, atom1_mol)

        bond_angle_target (float or int): The adsorption bond angle desired (in

            degrees).

            Details: bond_angle_target = angle(atom1_surf, atom1_mol, atom2_mol)

        dihedral_angle_target (float or int): The dihedral adsorption angle

            desired (in degrees).

            Details: dihedral_angle_target = dihedral(atom2_surf, atom1_surf,

            atom1_mol, atom2_mol)

                The rotation vector is facing the adsorbate from the surface

                (i.e. counterclockwise rotation when facing the surface (i.e.

                view from top))

        mol_dihedral_angle_target (float or int): The adsorbate dihedral angle

            desired (in degrees).

            Details: mol_dihedral_angle_target = dihedral(atom1_surf, atom1_mol,

            atom2_mol, atom3_mol)

                The rotation vector is facing atom2_mol from atom1_mol

    Returns:

        None (the `molecule` object is modified in-place)

    """

    vect_surf = get_atom_coords(surf, ctr2_surf) - get_atom_coords(surf,

                                                                   ctr1_surf)

    vect_inter = get_atom_coords(molecule, ctr1_mol) \

        - get_atom_coords(surf, ctr1_surf)

    vect_mol = get_atom_coords(molecule, ctr2_mol) - get_atom_coords(molecule,

                                                                     ctr1_mol)

    vect2_mol = get_atom_coords(molecule, ctr3_mol) - get_atom_coords(molecule,

                                                                      ctr2_mol)

    # Check if dihedral angles can be defined

    do_dihedral = dihedral_angle_target is not None

    do_mol_dihedral = mol_dihedral_angle_target is not None

    dihedral_use_mol2 = False

    # Check if vect_surf and bond_vector are aligned

    if np.allclose(np.cross(vect_surf, bond_vector), 0):

        do_dihedral = False

    # Check if requested bond angle is not flat

    if np.isclose((bond_angle_target + 90) % 180 - 90, 0):

        do_mol_dihedral = False

        dihedral_use_mol2 = True

    # Check if vect_mol and vect2_mol are not aligned

    if np.allclose(np.cross(vect_mol, vect2_mol), 0):

        do_mol_dihedral = False

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

    #       Translation       #

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

    # Compute and apply translation of adsorbate

    translation = bond_vector - vect_inter

    molecule.translate(translation)

    # Update adsorption bond

    vect_inter = get_atom_coords(molecule, ctr1_mol) - \

        get_atom_coords(surf, ctr1_surf)

    # Check if translation was successful

    if np.allclose(vect_inter, bond_vector):

        pass  # print("Translation successfully applied (error: ~ {:.5g} unit "

        # "length)".format(np.linalg.norm(vect_inter - bond_vector)))

    else:

        err = 'An unknown error occured during the translation'

        logger.error(err)

        raise AssertionError(err)

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

    #   Bond angle rotation   #

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

    # Compute rotation vector

    rotation_vector = np.cross(-vect_inter, vect_mol)

    if np.allclose(rotation_vector, 0, atol=1e-5):

        # If molecular bonds are aligned, any vector orthogonal to vect_inter

        # can be used Such vector can be found as the orthogonal rejection of

        # either X-axis, Y-axis or Z-axis onto vect_inter (since they cannot

        # be all aligned)

        non_aligned_vector = np.zeros(3)

        # Select the most orthogonal axis (lowest dot product):

        non_aligned_vector[np.argmin(np.abs(vect_inter))] = 1

        rotation_vector = non_aligned_vector - np.dot(non_aligned_vector,

                                                      vect_inter) / np.dot(

            vect_inter, vect_inter) * vect_inter

    # Retrieve current bond angle

    bond_angle_ini = get_vect_angle(-vect_inter, vect_mol, rotation_vector)

    # Apply rotation to reach desired bond_angle

    molecule.rotate(bond_angle_target - bond_angle_ini, v=rotation_vector,

                    center=get_atom_coords(molecule, ctr1_mol))

    # Update molecular bonds

    vect_mol = get_atom_coords(molecule, ctr2_mol) - get_atom_coords(molecule,

                                                                     ctr1_mol)

    vect2_mol = get_atom_coords(molecule, ctr3_mol) - get_atom_coords(molecule,

                                                                      ctr2_mol)

    # Check if rotation was successful

    bond_angle = get_vect_angle(-vect_inter, vect_mol)

    if np.isclose((bond_angle - bond_angle_target + 90) % 180 - 90, 0,

                  atol=1e-3) and np.allclose(get_atom_coords(molecule, ctr1_mol)

                                             - get_atom_coords(surf,

                                                               ctr1_surf),

                                             vect_inter):

        pass  # print("Rotation successfully applied (error: {:.5f}°)".format(

        # (bond_angle - bond_angle_target + 90) % 180 - 90))

    else:

        err = 'An unknown error occured during the rotation'

        logger.error(err)

        raise AssertionError(err)

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

    # Dihedral angle rotation #

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

    # Perform dihedral rotation if possible

    if do_dihedral:

        # Retrieve current dihedral angle (by computing the angle between the

        # vector rejection of vect_surf and vect_mol onto vect_inter)

        vect_inter_inner = np.dot(vect_inter, vect_inter)

        vect_surf_reject = vect_surf - np.dot(vect_surf, vect_inter) / \

            vect_inter_inner * vect_inter

        if dihedral_use_mol2:

            vect_mol_reject = vect2_mol - np.dot(vect2_mol, vect_inter) / \

                              vect_inter_inner * vect_inter

        else:

            vect_mol_reject = vect_mol - np.dot(vect_mol, vect_inter) / \

                              vect_inter_inner * vect_inter

        dihedral_angle_ini = get_vect_angle(vect_surf_reject, vect_mol_reject,

                                            vect_inter)

        # Apply dihedral rotation along vect_inter

        molecule.rotate(dihedral_angle_target - dihedral_angle_ini,

                        v=vect_inter, center=get_atom_coords(molecule,

                                                             ctr1_mol))

        # Update molecular bonds

        vect_mol = get_atom_coords(molecule, ctr2_mol) - \

            get_atom_coords(molecule, ctr1_mol)

        vect2_mol = get_atom_coords(molecule, ctr3_mol) - \

            get_atom_coords(molecule, ctr2_mol)

        # Check if rotation was successful

        # Check dihedral rotation

        if dihedral_use_mol2:

            vect_mol_reject = vect2_mol - np.dot(vect2_mol, vect_inter) / \

                              vect_inter_inner * vect_inter

        else:

            vect_mol_reject = vect_mol - np.dot(vect_mol, vect_inter) / \

                              vect_inter_inner * vect_inter

        dihedral_angle = get_vect_angle(vect_surf_reject, vect_mol_reject,

                                        vect_inter)

        # Check bond rotation is unmodified

        bond_angle = get_vect_angle(-vect_inter, vect_mol)

        if np.isclose((dihedral_angle - dihedral_angle_target + 90) % 180 - 90,

                      0, atol=1e-3) \

                and np.isclose((bond_angle - bond_angle_target + 90) %