Passivate edges in carbon nanostructures

[decerio]

Hi,

I work with graphene nanostructures like graphene nanoribbons and nanoporous graphene, and I think there may be many other SISL users that very find useful the agnr, zgnr, graphene_nanoribbons. However if one wants to generate carbon geometries that can be used as input geometries for DFT codes, it is desirable to have the edges correctly passivated with hydrogen (H) atoms.

Even more, edge functionalization with chemical species different from H has been experimentally achieved:

• Ketone groups: "Magnetic Interactions Between Radical Pairs in Chiral Graphene Nanoribbons", Nano Lett. 22, 164–171 (2022). DOI: https://doi.org/10.1021/acs.nanolett.1c03578

• Fluorine: "Stabilizing Edge Fluorination in Graphene Nanoribbons", ACS Nano 14, 11120-11120 (2020). DOI: https://doi.org/10.1021/acsnano.0c01837

• Amino groups: "Band Depopulation of Graphene Nanoribbons Induced by Chemical Gating with Amino Groups", ACS Nano 14, 1895-1901 (2020). DOI: https://doi.org/10.1021/acsnano.9b08162

As far as I know, edge passivation is not implemented in SISL, is it? At least, I didn't find it. So it might be interested to add a passivate parameter to agnr, zgnr, graphene_nanoribbons that gives the passivated geometries by setting passivate = True. It could also be added as a method to passivate/functionalize geometries that have already been created (e.g. a porous ribbon) by doing something like geometry.passivate_edges().

My suggestion is something like this:

def passivate_edges(g,bond=1.1,atoms=None,species='H'):
    """
    Passivates atoms with dangling bonds.
    
    Parameters
    ----------
    g: sisl.Geometry
        Geometry
    bond: int, float
        Distance between the passivated edge atom and the new added atom.
        In the case of H passivation this is the C-H bond-length.
    atoms: array, list
        Atomic indices to be considered in the passivation. Only the atoms
        with a dangling bonds (less than 3 nearest neighbours) will be
        passivated. If None is passed, all atoms will be considered
    species: string, int
        Chemical label or atomic number of the atoms that will passivate the
        edge.
    """
    
    gtmp = g.copy()
    if atoms is None:
        atoms = np.arange(gtmp.na)
    for ia in atoms:
        nn = gtmp.close(ia,R=[0.1,1.5])[1]
        if len(nn) < 3:
            vec1 = gtmp.axyz(nn[0]) - gtmp.xyz[ia]
            vec2 = gtmp.axyz(nn[1]) - gtmp.xyz[ia]
            vec3 = -bond*(vec1+vec2)/np.linalg.norm(vec1+vec2)
            edge_atom = sisl.Geometry(xyz=[gtmp.xyz[ia]+vec3],atoms=sisl.Atom(species))
            gtmp = gtmp.add(edge_atom)
    return gtmp

I attach a jupyter notebook that serves as an example of how could this function be used. example_passivation.ipynb.zip

Do you think it could be interesting?