enh: finalized few missing details in bilayer method

Signed-off-by: Nick Papior <[email protected]>

Author: zerothi

sisl/geom/__init__.py

@@ -22,6 +22,7 @@
    honeycomb
    nanotube
    diamond
+   bilayer
 
 
 Basic

sisl/geom/bilayer.py

@@ -1,45 +1,52 @@
 import numpy as np
 
-from sisl import Atom, geom, Cuboid
+from sisl import geom, Atom, Cuboid
 
 __all__ = ['bilayer']
 
 
-def bilayer(bond=1.42, bottom_atom=None, top_atom=None, stacking='AB', twist=(0, 0), separation=3.35, ret_angle=False, layer='both'):
-    """ Commensurate unit cell of a hexagonal bilayer structure, possibly with a twist angle.
+def bilayer(bond=1.42, bottom_atom=None, top_atom=None, stacking='AB',
+            twist=(0, 0), separation=3.35, ret_angle=False, layer='both'):
+    r""" Commensurate unit cell of a hexagonal bilayer structure, possibly with a twist angle.
 
-    This routine follows the prescription of twisted bilayer graphene found in
-    Laissardiere et al., Nano Lett. 10, 804-808 (2010).
+    This routine follows the prescription of twisted bilayer graphene found in [1]_.
 
+    Notes
+    -----
+    This routine may change in the future to force some of the arguments.
 
     Parameters
     ----------
-    bond : 1.42
-       bond length (in Angstrom) between atoms in the honeycomb lattice
-    bottom_atom : Atom(6)
-       atom (or atoms) in the bottom layer
-    top_atom : Atom(6)
-       atom (or atoms) in the top layer
+    bond : float, optional
+       bond length between atoms in the honeycomb lattice
+    bottom_atom : Atom, optional
+       atom (or atoms) in the bottom layer. Defaults to ``Atom(6)``
+    top_atom : Atom, optional
+       atom (or atoms) in the top layer, defaults to `bottom_atom`
     stacking : {'AB', 'AA'}
        stacking sequence of the bilayer
-    twist : (0, 0)
+    twist : tuple of int, optional
        integer coordinates (m, n) defining a commensurate twist angle
-    separation : float
+    separation : float, optional
        distance between the two layers (in Angstrom)
-    ret_angle : bool
+    ret_angle : bool, optional
        return the twist angle (in degrees) in addition to the geometry instance
     layer : {'both', 'bottom', 'top'}
-       control which layer(s) to include in geometry
+       control which layer(s) to return
+
+    References
+    ----------
+    .. [1] G. Trambly de Laissardiere, D. Mayou, L. Magaud, "Localization of Dirac Electrons in Rotated Graphene Bilayers", Nano Letts. 10, 804-808 (2010)
     """
     if bottom_atom is None:
         bottom_atom = Atom(Z=6, R=bond * 1.01)
     if top_atom is None:
-        top_atom = Atom(Z=6, R=bond * 1.01)
+        top_atom = bottom_atom
 
     # Construct two layers
-    bottom = geom.honeycomb(bond=bond, atom=bottom_atom)
-    top = geom.honeycomb(bond=bond, atom=top_atom)
-    honeycomb = bottom.copy()
+    bottom = geom.honeycomb(bond, bottom_atom)
+    top = geom.honeycomb(bond, top_atom)
+    ref_cell = bottom.cell.copy()
 
     if stacking.lower() == 'aa':
         top = top.move([0, 0, separation])
@@ -54,7 +61,7 @@ def bilayer(bond=1.42, bottom_atom=None, top_atom=None, stacking='AB', twist=(0,
         m, n = n, m
 
     if not (isinstance(n, int) and isinstance(m, int)):
-        raise RuntimeError('twist coordinates need to be integers!')
+        raise ValueError("bilayer: twist coordinates need to be integers!")
 
     if m == n:
         # No twisting
@@ -63,22 +70,22 @@ def bilayer(bond=1.42, bottom_atom=None, top_atom=None, stacking='AB', twist=(0,
         natoms = 2
     else:
         # Twisting
-        cos_theta = (n ** 2 + 4 * n * m + m ** 2) / (2 * (n ** 2 + n * m + m ** 2) )
-        theta = np.arccos(cos_theta) * 360 / (2 * np.pi)
+        cos_theta = (n ** 2 + 4 * n * m + m ** 2) / (2 * (n ** 2 + n * m + m ** 2))
+        theta = np.arccos(cos_theta) * 180 / np.pi
         rep = 4 * (n + m)
         natoms = 2 * (n ** 2 + n * m + m ** 2)
 
     if rep > 1:
         # Set origo through an A atom near the middle of the geometry
         bottom = bottom.tile(rep, axis=0).tile(rep, axis=1)
         top = top.tile(rep, axis=0).tile(rep, axis=1)
-        tvec = rep * (honeycomb.cell[0] + honeycomb.cell[1]) / 2
+        tvec = rep * (ref_cell[0] + ref_cell[1]) / 2
         bottom = bottom.move(-tvec)
         top = top.move(-tvec)
 
         # Set new lattice vectors
-        bottom.cell[0] = n * honeycomb.cell[0] + m * honeycomb.cell[1]
-        bottom.cell[1] = -m * honeycomb.cell[0] + (n + m) * honeycomb.cell[1]
+        bottom.cell[0] = n * ref_cell[0] + m * ref_cell[1]
+        bottom.cell[1] = -m * ref_cell[0] + (n + m) * ref_cell[1]
 
         # Rotate top layer around A atom in bottom layer
         top = top.rotate(theta, [0, 0, 1])
@@ -96,16 +103,19 @@ def bilayer(bond=1.42, bottom_atom=None, top_atom=None, stacking='AB', twist=(0,
 
     if rep > 1:
         # Remove atoms outside cell
-        center = 1e-3 * np.array([1, 1, 1])
-        cell_box = Cuboid(1.0 * bilayer.cell[:], center=center)
+        cell_box = Cuboid(bilayer.cell)
+        cell_box.set_origo([-0.0001] * 3)
         inside_idx = cell_box.within_index(bilayer.xyz)
         bilayer = bilayer.sub(inside_idx)
 
         # Rotate whole cell
         vec = bilayer.cell[0] + bilayer.cell[1]
-        vec_costh = np.dot([1, 0, 0], vec) / np.dot(vec, vec) ** .5
-        vec_th = np.arccos(vec_costh) * 360 / (2 * np.pi)
+        vec_costh = vec[0] / vec.dot(vec) ** 0.5
+        vec_th = np.arccos(vec_costh) * 180 / np.pi
         bilayer = bilayer.rotate(vec_th, [0, 0, 1])
+    else:
+        # Shift back
+        bilayer = bilayer.move(tvec)
 
     # Sanity check
     assert len(bilayer) == natoms

sisl/geom/tests/test_geom.py

@@ -7,37 +7,39 @@
 import numpy as np
 
 
-class Test:
-
-    def test_basis(self):
-        a = sc(2.52, Atom['Fe'])
-        a = bcc(2.52, Atom['Fe'])
-        a = bcc(2.52, Atom['Fe'], orthogonal=True)
-        a = fcc(2.52, Atom['Au'])
-        a = fcc(2.52, Atom['Au'], orthogonal=True)
-        a = hcp(2.52, Atom['Au'])
-        a = hcp(2.52, Atom['Au'], orthogonal=True)
-
-    def test_flat(self):
-        a = graphene()
-        a = graphene(atom='C')
-        a = graphene(orthogonal=True)
-
-    def test_nanotube(self):
-        a = nanotube(1.42)
-        a = nanotube(1.42, chirality=(3, 5))
-        a = nanotube(1.42, chirality=(6, -3))
-
-    def test_diamond(self):
-        a = diamond()
-
-    def test_bilayer(self):
-        a = bilayer(1.42)
-        a = bilayer(1.42, stacking='AA')
-        for m in range(7):
-            a = bilayer(1.42, twist=(m, m + 1))
-        a = bilayer(1.42, twist=(6, 7), layer='bottom')
-        a = bilayer(1.42, twist=(6, 7), layer='TOP')
-        a = bilayer(1.42, bottom_atom=(Atom['B'], Atom['N']), twist=(6, 7))
-        a = bilayer(1.42, top_atom=(Atom(5), Atom(7)), twist=(6, 7))
-        a, th = bilayer(1.42, twist=(6, 7), ret_angle=True)
+def test_basis():
+    a = sc(2.52, Atom['Fe'])
+    a = bcc(2.52, Atom['Fe'])
+    a = bcc(2.52, Atom['Fe'], orthogonal=True)
+    a = fcc(2.52, Atom['Au'])
+    a = fcc(2.52, Atom['Au'], orthogonal=True)
+    a = hcp(2.52, Atom['Au'])
+    a = hcp(2.52, Atom['Au'], orthogonal=True)
+
+
+def test_flat():
+    a = graphene()
+    a = graphene(atom='C')
+    a = graphene(orthogonal=True)
+
+
+def test_nanotube():
+    a = nanotube(1.42)
+    a = nanotube(1.42, chirality=(3, 5))
+    a = nanotube(1.42, chirality=(6, -3))
+
+
+def test_diamond():
+    a = diamond()
+
+
+def test_bilayer():
+    a = bilayer(1.42)
+    a = bilayer(1.42, stacking='AA')
+    for m in range(7):
+        a = bilayer(1.42, twist=(m, m + 1))
+    a = bilayer(1.42, twist=(6, 7), layer='bottom')
+    a = bilayer(1.42, twist=(6, 7), layer='TOP')
+    a = bilayer(1.42, bottom_atom=(Atom['B'], Atom['N']), twist=(6, 7))
+    a = bilayer(1.42, top_atom=(Atom(5), Atom(7)), twist=(6, 7))
+    a, th = bilayer(1.42, twist=(6, 7), ret_angle=True)