Function to produce armchair heteroribbons.

[pfebrer]

I was tired of having to think about how to build each heteroribbon and I've written a function to produce only valid heteroribbons (with no dangling bonds) specifying the width, units and shift of each section. Initially, I thought that it might be a useful thing to put in sisl_toolbox, but seeing that these structures are getting popular I thought that maybe it's worth it to add it to sisl.geom.

Here's the (pretty long) code.

```python import sisl import math import numpy as np def graphene_heteroribbon(Ws, units, shifts=None, **kwargs): """Builds a graphene nanoribbon with multiple widths. The function makes sure that only valid ribbons, that is ribbons with no dangling bonds, are being produced. Explanatory errors are raised if that's not possible. Only the connection at the periodic boundary remains unchecked. It only works for armchair ribbons for now. Parameters ---------- Ws: array-like of int The width of each section of the heteroribbon. units: array-like of int The number of units of each section. Note that a "unit" is not a unit cell, but half of it. I.e. a zigzag string. There are two situations in which you can choose the starting configuration of your border ("open", "closed"). At the beggining of the heteroribbon and on an even section after an odd one. In those cases, an integer imaginary number here will indicate that you want the open configuration (e.g. 2j indicates that you want 2 units, starting with an open configuration). If you pass an imaginary number on an invalid situation, an error will be raised. shifts: array-like of int, optional How each section must be shifted from the previous one. The exact meaning of the value depends on the sections that are being joined so that only valid ribbons are produced: - If both ribbons are odd, ribbons are aligned on the center and then shifted: ^ (shift) zig zag atoms if the incoming section is thinner than the previous one. ^ (2 * shift) zig zag atoms if they differ in a multiple of 4 atoms and both configurations are closed or if they differ in (2 + a multiple of 4 atoms) and configurations are different (one open and one closed.) ^ (1 + 2 * shift) zig zag atoms otherwise. - If both ribbons are even, they are aligned by the open end of the last section, not the incoming one, and then they are shifted: ^ (0) zig zag atoms if the shift is 0. In this case both sections are forced to have the same edge open. ^ (-1 + 2 * abs(shift)) zig zag atoms in the appropiate direction if abs(shift) is bigger than 1. Notice that you can only shift the open edges toward the center. - If the ribbons have different parities, they are aligned on the open edge of the even ribbon and then shifted: ^ (0) zig zag atoms if shift is 0 in the particular cases that (i) the incoming odd section is thinner than the last even section or (ii) the last section, being odd, has a border in the open configuration. Only in those cases does it make sense to exactly match the edges. ^ (bias + 2 * abs(shift)) zig zag atoms in the appropiate direction, where `bias` is either -1 or 1 depending on the exact situation and always ensures that the minimum shift is always `abs(shift) == 1`. If not provided, it defaults to an array full of zeros. **kwargs Keyword arguments are directly passed to `sisl.geom.graphene_nanoribbon`. Returns ---------- sisl.Geometry: The final structure of the heteroribbon. """ if kwargs.get('kind', 'armchair') == "zigzag": raise ValueError("Zigzag heteroribbons are not implemented yet.") # Initialize all variables that are going to be updated through the # iterations geom = None last_addition = None last_open = False last_W = -1 # Get the distance of an atom shift. atom_shift = kwargs.get('bond', 1.42) * math.sin(math.radians(60)) # If shifts were not provided it means the user wants no shifts, # so just initialize the array to all zeros. if shifts is None: shifts = np.zeros(len(Ws)) # Helper function to manipulate ribbons. def _open_start(geometry): """Changes the border used for a ribbon. It does so by shifting half unit cell. This must be done before any tiling of the geometry. """ geometry = geometry.move(geometry.cell[0] / 2) geometry.xyz = (geometry.fxyz % 1).dot(geometry.cell) return geometry # Loop through all the sections that compose the heteroribbon. for i, (W, W_units, shift) in enumerate(zip(Ws, units, shifts)): # Check that we are not joining an open odd ribbon with # a smaller ribbon, since no matter what the shift is there will # always be dangling bonds. if last_W % 2 == 1 and W < last_W: assert not last_open, (f"At the interface between {last_W} and {W}, {last_W} has an open end." "If the wider ribbon is odd, it must always have a closed end.") # Get the unit cell of the ribbon that we are going to add to the # heteroribbon. new_addition = sisl.geom.graphene_nanoribbon(W, **kwargs) # Unit cells returned by sisl are always in the "closed" configuration. # In even width ribbons, there will always be an open edge and a closed # edge. We use "open" in the situation where the top edge is open. open_start = False # If the number of units has been provided as an imaginary number, this means that # the user wants the open configuration of the ribbon. if W_units.imag > 0: # However this is only possible in very specific cases, since in all other # cases the border of the incoming section is fully determined by the border # of the previous section and the shift. if i == 0 or W % 2 == 0 and last_W % 2 == 1: new_addition = _open_start(new_addition) open_start = not open_start W_units = int(W_units.imag) else: raise ValueError("This section does not allow choosing the configuration, but an imaginary", f" number {W_units} was passed as the number of units." ) # If this is not the first ribbon section in the heteroribbon, we are going to # calculate the shifts. if not i == 0: # Get the difference in width between the previous and this ribbon section W_diff = W - last_W # And also the mod(4) because there are certain differences if the width differs # on 1, 2, 3 or 4 atoms. After that, the cycle just repeats (e.g. 5 == 1, etc). diff_mod = W_diff % 4 # Now, we need to calculate the offset that we have to apply to the incoming # section depending on several factors. if diff_mod % 2 == 0 and W % 2 == 1: # Both sections are odd # In this case, ribbons are aligned based on their centers. offset = (last_addition.center() - new_addition.center())[1] different_configs = last_open != open_start if W < last_W: # The incoming section is thinner than the last one. Note that at this point # we are sure that the last section has a closed border, otherwise we # would have raised an error. At this point, centers can differ by any # integer number of atoms, but we might need to open the start depending # on the shift so that the sections match. if shift % 2 == 0 and diff_mod == 2 or shift % 2 == 1 and diff_mod == 0: # We must open the start new_addition = _open_start(new_addition) open_start = not open_start # Shift limits are a bit complicated and are different for even and odd shifts. # This is because a closed incoming section can shift until there is no connection # between ribbons, while an open one needs to stop before its edge goes outside # the previous section. shift_lims = { "closed": last_W // 2 + W // 2 - 2, "open": last_W // 2 - W // 2 - 1 } shift_pars = { lim % 2: lim for k, lim in shift_lims.items() } assert (shift % 2 == 0 and abs(shift) <= shift_pars[0])\ or (shift % 2 == 1 and abs(shift) <= shift_pars[1]),\ (f"Shift must be an even number between {-shift_pars[0]} and" f" {shift_pars[0]} or an odd number between {-shift_pars[1]}" f" and {shift_pars[1]}") offset += shift * atom_shift else: # At this point, we already know that the incoming section is wider and # therefore it MUST have a closed start, otherwise there will be dangling bonds. if diff_mod == 2 and last_open or diff_mod == 0 and not last_open: # In these cases, centers must match or differ by an even number of atoms. offset += (2*shift) * atom_shift # And these are the limits for the shift if last_open: # To avoid the current open section leaving dangling bonds. shift_lim = (W - last_W) // 4 else: # To avoid sections being disconnected. shift_lim = max((last_W - 3) // 2, 0) + (W - last_W) // 4 min_shift, max_shift = - shift_lim, shift_lim else: # Otherwise, centers must differ by an odd number of atoms. offset += (1 + 2*shift) * atom_shift # And these are the limits for the shift if last_open: # To avoid the current open section leaving dangling bonds. min_shift = - max(0, (W - last_W) // 4) max_shift = max(0, (W - last_W) // 4 - 1) else: # To avoid sections being disconnected. shift_lim = max((last_W - 3) // 2, 0) + (W - last_W) // 4 min_shift, max_shift = -shift_lim - 1, shift_lim assert min_shift <= shift <= max_shift, (f"Shift must be between {min_shift}" f" and {max_shift}, but {shift} was provided." ) else: # There is at least one even section. # In this case ribbons are aligned based on one of their edges. # We choose the edge of the even ribbon that is open (i.e. has # dangling bonds). If both ribbons are even, we align it based # on the last section, not the incoming one. if last_W % 2 == 0 and last_open or W % 2 == 0 and open_start: # Align them on the top edge offset = last_addition.xyz[:, 1].max() - new_addition.xyz[:, 1].max() else: # Align them on the bottom edge offset = last_addition.xyz[:, 1].min() - new_addition.xyz[:, 1].min() # We have to make sure that the open edge of the even ribbons (however # many there are) is always shifted towards the center. Shifting in the # opposite direction would result in dangling bonds. no_shift = False if diff_mod % 2 == 0: min_shift = max(0, (W - last_W) - 1) max_shift = W // 2 # Both ribbons are even if (last_open and shift == 0) or (not last_open and shift != 0): # If there is no shift, we must make sure that both ribbons are open or close. # If there is shift, both ribbons must have different configurations to avoid # dangling bonds. new_addition = _open_start(new_addition) open_start = not open_start if shift == 0: no_shift = True else: bias, shift_sign = { True: (-1, 1), False: (1, -1) }[last_open] elif W % 2 == 1: # Last section was even, incoming section is odd. min_shift = 0 max_shift = last_W // 2 if W < last_W and shift == 0: # Incoming odd section is thinner than current section. # We can make the ends match if we ensure that the odd section is open.Make ends match (shift = 0): The odd ribbon must be open. new_addition = _open_start(new_addition) open_start = not open_start no_shift = True else: bias, shift_sign = { (True, True): (-1, 1), (True, False): (1, -1), (False, True): (1, 1), (False, False): (-1, -1) }[(W < last_W, last_open)] else: # Last section was odd, incoming section is even. min_shift = 0 max_shift = W // 2 if last_open: assert shift == 0, ("The shift between an open odd ribbon and an incoming" "even ribbon must be always 0, otherwise the odd section will have a dangling bond" ) no_shift = True else: bias, shift_sign = { True: (-1, -1), False: (1, 1) }[open_start] if not no_shift: # Check that the shift value is correct. There are two things to check here: # Provided shift is in the right direction right_shift_sign = shift*shift_sign - abs(shift) == 0 # Provided shift falls within the bounds. shift_in_bounds = min_shift <= abs(shift) <= max_shift assert (right_shift_sign and shift_in_bounds), ( f" Shift must be between {min_shift*shift_sign} and {max_shift*shift_sign}, but {shift} was provided." ) offset += (bias + 2*shift) * atom_shift # Apply the offset that we have calculated. new_addition = new_addition.move([0, offset, 0]) # Check how many times we have to tile the unit cell (tiles) and whether # we have to cut the last string of atoms (cut_last) tiles, res = divmod(W_units + 1, 2) cut_last = res == 0 # Tile the current section unit cell new_addition = new_addition.tile(tiles, 0) # Cut the last string of atoms. if cut_last: new_addition.cell[0,0] *= W_units / (W_units + 1) new_addition = new_addition.remove({"x": (new_addition.cell[0,0] - 0.01, None)}) # Initialize the heterorribon or append to it if it's already initialized. if i == 0: geom = new_addition else: geom = geom.append(new_addition, 0) # Update the state for the next iteration last_addition = new_addition last_open = open_start != cut_last last_W = W return geom ```

And some examples of structures you can get:

graphene_heteroribbon([7,9,6], units=[2,4,2]).plot(axes="xy")

graphene_heteroribbon([7,8,9,8,7], [3,4,2,2,3], shifts=[0,0,-1,1,0]).plot(axes="xy")

graphene_heteroribbon([7,8,9,8,7], [2j,1,2,1j,2]).plot(axes="xy")

And this is just trying random numbers, you always get the valid ribbons or an error if that's not possible.

Also, this might be useful to design TranSIESTA setups with graphene electrodes? Something like:

graphene_heteroribbon([31, 9, 7, 5, 7, 9, 31], [4,9,5,3,6,9,4]).plot(axes="xy")

and then you just cut the cell to get the graphene electrodes, for example.

[zerothi]

Great job @pfebrer.

An idea, could we make the interface a bit more easy to use, something along the lines:

# create a class with only static members
class heteroribbon:
     # I don't know if `bond` is necessary, it may be rather difficult to match them up, but one
    # could for instance shift it according to the previous bond, and then track and *accummulated* shift?
    segment = namedtuple("Segment", ["W", "L", "shift", "bond", "atoms"], (0, None, None))
    def __call__(self, segments, **kwargs):
        for s in segments:
            atoms = s.atoms or kwargs["atoms"]
            bond = s.bond or kwargs["bond"]
            ... continue with the creation

class graphene_heteroribbon(heteroribbon):
    def __call__(self, segments, *, atoms=sisl.Atom(6, 1.42), **kwargs):
         super()(segments, atoms=atoms, **kwargs)

It think this would make it a bit more readable?

S = heteroribbon.segment
heteroribbon(
    (S(W=10, L=3, atoms=sisl.Atom(5),
     S(W=11, L=5) # default atom
     S(W=12, L=3, shift=1, atoms=sisl.Atom(7)),
    atoms=sisl.Atom(6))

# It could still be done using 3 different lists
widths = [10, 11, 12]
lengths = [3, 5, 3]
shifts = [0, 0, 1]
heteroribbon(zip(widths, lengths, shifts), atoms=sisl.Atom(6))

Would this make sense?

I think this fits nicely in the geom.nanoribbon module?

[pfebrer]

Yes, sure. I'm not sure how easy it is to explain to new python users the named tuple thing though. So they would probably just go with lists. But

heteroribbon((
    [10, 3, 0, sisl.Atom(5)],
    [11, 5],
    [12, 3, 1, sisl.Atom(7)]
))

Seems fine to me. I think that this approach makes more sense but the other one is easier to use interactively, so I have divided opinions. But yeah since you can go with zip it probably doesn't make it much more tedious. Also, to design a graphical interface your approach is much better :)

Probably we could also assign atoms to names in the geometry. geom["section 1"] would be the atoms in section 1, etc... Or you don't want to use this names thing? :sweat_smile:

Specifying different bonds could be done. There's no need to track accumulated shift, since the shift is calculated always based only on the last section. However this would even add more complexity to the code. Would you need it currently? Otherwise I would say we could incorporate it later if someone needs it.

Also, I'm not sure that using classes makes the code clean/clear. Perhaps we could use plain functions and provide the heteroribbon_segment construct in the geom.nanoribbon module?

[pfebrer]

Probably we could also assign atoms to names in the geometry. geom["section 1"] would be the atoms in section 1, etc... Or you don't want to use this names thing? :sweat_smile:

Or return the sections breakpoints (in Ang).

[zerothi]

Yes, sure. I'm not sure how easy it is to explain to new python users the named tuple thing though. So they would probably just go with lists. But

heteroribbon((
    [10, 3, 0, sisl.Atom(5)],
    [11, 5],
    [12, 3, 1, sisl.Atom(7)]
))

Seems fine to me. I think that this approach makes more sense but the other one is easier to use interactively, so I have divided opinions. But yeah since you can go with zip it probably doesn't make it much more tedious. Also, to design a graphical interface your approach is much better :)

Ok, and also, the heteroribbon.segment isn't really needed, it is just a convenience that is possible. So definitely emphasis should be put on using it in the simple case. (as you have written it). Then additions may be added by using the namedtuple stuff.

Probably we could also assign atoms to names in the geometry. geom["section 1"] would be the atoms in section 1, etc... Or you don't want to use this names thing? sweat_smile

Hehe, yeah, we really need these names to get going. Again we have the problems of copying, subbing, etc... It would be fine if you add them (if you want) then we can always fix it later. Just don't mention in the docs.

Specifying different bonds could be done. There's no need to track accumulated shift, since the shift is calculated always based only on the last section. However this would even add more complexity to the code. Would you need it currently? Otherwise I would say we could incorporate it later if someone needs it.

No, not needed, just a thought when you start mixing C with N and B etc.

Also, I'm not sure that using classes makes the code clean/clear. Perhaps we could use plain functions and provide the heteroribbon_segment construct in the geom.nanoribbon module?

The end user does not notice. And the fact that extensions belong to the class makes it pretty clear, I think? In this case the class construct only serves as a namespace variable.

[zerothi]

Also, why do you think your (old) approach is easier interactively? Having to count indices where to put in new items seems awfully hard to me?

[pfebrer]

Also, why do you think your (old) approach is easier interactively? Having to count indices where to put in new items seems awfully hard to me?

Yeah it is a tradeoff between that and having to write new lists. But since you can go with zip you can choose, it's fine. Also I like to see the widths/units in the same place because then I can quickly identify which ribbon it is. Like: "Oh yes this is the [7,9,7] ribbon. May we try the [7,11,7]?" But this might be just me :laughing:

The end user does not notice. And the fact that extensions belong to the class makes it pretty clear, I think? In this case the class construct only serves as a namespace variable.

Some end user trying to understand how it works may look at it. Functions are more universal and easier to read than a class with a __call__ method. The chances they might just give up understanding it are pretty high. Also probably documentation would be more confusing. Unless you are planning to add more methods, which seems unlikely since __call__ wouldn't actually return an instance of the class, I would vote for functions. I think geom.nanoribbon is a good enough namespace for heteroribbon_section. And anyway it's not really an essential thing to use, it's more of an "advanced usage".

[zerothi]

You are absolutely right.

Lets change it to then:

def heteroribbon(...):
    ...
heteroribbon.Segment = namedtuple

;) Still a function. :)

As for the readability, zip is really easy :) And with the single list it enables adding options without disturbing backwards compatibility.

[pfebrer]

Good! Giving attributes to functions is accepted but nobody does it. This seems like a good case for it, I hope nobody gets mad :laughing: .

I'll submit a PR later in the day then.

[zerothi]

Good! Giving attributes to functions is accepted but nobody does it. This seems like a good case for it, I hope nobody gets mad laughing .

Exactly, it has been there since 2001, PEP232!

I'll submit a PR later in the day then.

:rocket:

[tfrederiksen]

Looks great!!! I agree that it would be very useful to have this in sisl.geom.

I have also been playing with heterostructures, but my function is not very sophisticated (can't do "open" configurations, for instance). I leave it here in case it could be of some inspiration:

def composite_agnr(width, bond=1.42, atoms=sisl.Atom(6, R=1.43)):
    if isinstance(width, Integral):
        return sisl.geom.agnr(width, bond=bond, atoms=atoms)
    hc = sisl.geom.honeycomb(bond, atoms, orthogonal=True)
    vx = hc.cell[0]
    vy = hc.cell[1]
    xyz = []
    for i, w in enumerate(width):
        if isinstance(w, Integral):
            rng = np.arange(0, w)
        else:
            rng = np.arange(w[0], w[1])
        a = 2 * (i % 2)
        b = a + 1
        for j in rng:
            if j % 2 == i % 2:
                xyz.append(hc.xyz[b] - j // 2 * vy + i // 2 * vx)
            else:
                xyz.append(hc.xyz[a] - j // 2 * vy + i // 2 * vx)
    ribbon = sisl.Geometry(xyz, atoms=atoms, sc=hc.cell)
    ribbon.xyz[:, 1] *= -1
    ribbon = ribbon.sort(axis=(0, 1))
    ribbon.cell[0, 0] *= len(width) / 2
    ribbon.cell[1, 1] += 20.
    xyz = ribbon.xyz.min(axis=0) * [1, 1, 0]
    ribbon.optimize_nsc()
    return ribbon.move(-xyz + [0, 10, 0])

It works by specifying first/last atoms in each "column":

composite_agnr([5, 6, 7, 7, 6, 5]).plot(axes="xy")

composite_agnr([5, [-1, 5], [-2, 5], [-2, 5], [-1, 5], 5]).plot(axes="xy")

It would be cool if one could do also heterostructures in the zigzag direction!

[tfrederiksen]

@pfebrer , since one is typically concerned with periodic structures, wouldn't it be simpler to create the "open" structures from the "closed" ones by shuffling a single "column" from one side to the other at the end of the construction?

[pfebrer]

Thanks @tfrederiksen!

I totally should use geom.optimize_nsc(), I didn't know it existed :)

@pfebrer , since one is typically concerned with periodic structures, wouldn't it be simpler to create the "open" structures from the "closed" ones by shuffling a single "column" from one side to the other at the end of the construction?

I don't know exactly what you mean. Do you mean doing so when the heteroribbon is completely build? There are cases in which both ends do not have the same width and there's only one string of atoms. E.g.:

graphene_heteroribbon([7, 8, 5], [1, 2j, 1]).plot(axes="xy", nsc=[2,1,1])

Also when an even section comes after an odd section, both "open" and "closed" starts of the even section are valid. This section is not at the borders of the ribbon and there must be some way of generating the different ribbons. That's why I introduced this 2j to mean "2 units in the open config". This was the most usable interface I could think of, instead of having a special argument for the configuration of each section :sweat_smile:.

I'm thinking about what would feel more natural for people, a combination of width and shift (my approach) or first atom/last atom (your approach). To me it seems more natural to specify the width and then (maybe) shift it. But that may be because I've become used to it. Specifying the ribbon by sections of length L instead of individual strings seems more natural as well and easier for interfaces to interact.

In your case the shift is always referenced to a fixed position and a given shift value always corresponds to the same shift. In my case the reference is the last section and the shifts do not always mean the same thing. This is so that only valid ribbons are produced and still valid shift values always feel natural (0 means align with previous, 1 means first valid shift up, etc...). This is completely deterministic but not trivial and I attempted to explain all cases in the best way on the docs. Do you think this makes sense or shift=1 should always mean 1 atom shift and then check if that's valid?

[pfebrer]

It would be cool if one could do also heterostructures in the zigzag direction!

Yes it would :sweat_smile:. Depending on when @zerothi is thinking on dropping the release I can give it a shot. Otherwise we can leave it for the future haha

(I only thought about them for a moment, but zig zag heterostructures seem to be less complex to treat :thinking: )

[zerothi]

Now I just really want to clarify the units in #418, and then I really want it out. But after that, I should be able to do releases more often (:crossed_fingers: )

[tfrederiksen]

I don't know exactly what you mean. Do you mean doing so when the heteroribbon is completely build? There are cases in which both ends do not have the same width and there's only one string of atoms.

I don't see this. In your example there are four vertical "strings" of atoms. In my function it would just be (apart from mirror):

composite_agnr([7, 8, 8, 5]).plot(axes="xy", nsc=[2, 1, 1])

Also when an even section comes after an odd section, both "open" and "closed" starts of the even section are valid. This section is not at the borders of the ribbon and there must be some way of generating the different ribbons. That's why I introduced this 2j to mean "2 units in the open config". This was the most usable interface I could think of, instead of having a special argument for the configuration of each section sweat_smile.

I'm thinking about what would feel more natural for people, a combination of width and shift (my approach) or first atom/last atom (your approach). To me it seems more natural to specify the width and then (maybe) shift it. But that may be because I've become used to it. Specifying the ribbon by sections of length L instead of individual strings seems more natural as well and easier for interfaces to interact.

In your case the shift is always referenced to a fixed position and a given shift value always corresponds to the same shift. In my case the reference is the last section and the shifts do not always mean the same thing. This is so that only valid ribbons are produced and still valid shift values always feel natural (0 means align with previous, 1 means first valid shift up, etc...). This is completely deterministic but not trivial and I attempted to explain all cases in the best way on the docs. Do you think this makes sense or shift=1 should always mean 1 atom shift and then check if that's valid?

"open" vs "closed" both refer to the same periodic structure, right? What I meant is that we can build the periodic structure first (in the simplest possible way) and then decide subsequently which bond should be crossing the unit cell borders.

[tfrederiksen]

"open" vs "closed" both refer to the same periodic structure, right? What I meant is that we can build the periodic structure first (in the simplest possible way) and then decide subsequently which bond should be crossing the unit cell borders.

Something like this:

g_closed = composite_agnr([7, 8, 8, 5])
g_open = g_closed.move(g_closed.cell[0], atoms=range(7)).move([-1.42, 0, 0])
g_open.plot(axes="xy", nsc=[2, 1, 1])

[pfebrer]

Indeed, this is something that can be done afterwards. But what I mean is that each particular section can start "open" or "closed". In our example, the first 7 atom string could start "open", for instance, and it would be a different ribbon.

[pfebrer]

See this ribbon, for example:

graphene_heteroribbon([7,10,9,8], [1j,2j,2,1], [0,0,1,0]).plot(axes="xy",nsc=[2,1,1])

[pfebrer]

I mean that you can not just build this ribbon with the "closed" configuration and then convert it.

Ok, certainly you could build it starting with the 9 atom section and then shift it. But I don't know, is this always the case? And it seems easier to be able to specify where your unit cell starts anyway in general, isn't it?

[tfrederiksen]

I mean that you can not just build this ribbon with the "closed" configuration and then convert it.

Ok, certainly you could build it starting with the 9 atom section and then shift it. But I don't know, is this the case always?

I think that's the case, yes.

[pfebrer]

Here both 7 and 9 start in the open configuration:

graphene_heteroribbon([7,10,9,8], [1j,2j,1,2], [0,0,0,1]).plot(axes="xy",nsc=[2,1,1])

I don't know, to me it seems easier to specify how you want to start than having the constraint of starting with a "closed" string and then calculate later how much you need to shift the ribbon.

[tfrederiksen]

Up to a shift of the cell, this structure is equivalent to

composite_agnr([8,7,[-4,6],[-4,6],[-3,6],8]).move([0,-5,0]).plot(axes='xy', nsc=[2,1,1])

Of course, personal preferences may differ :)

Note that the "open/closed" notation is not clear for even-numbered widths. In my example above the first column is of width 8 - should this be considered one or the other case?

[pfebrer]

Yes, but the point is that you had to find in your ribbon where can you start to build it and then understand how much you need to shift. Anyway, if you have the possibility to start "open" you can choose to shift it yourself if you want as well :)

Note that the "open/closed" notation is not clear for even-numbered widths. In my example above the first column is of width 8 - should this be considered one or the other case?

Yes, I know. I consider "closed" the one that sisl produces. Just to be consistent with the odd widths. That is, I consider "closed" the one with the bottom edge open.