mnt: allowing h-shells all-around

I have not worked out the cartesian names for these, but perhaps
they will also be too long to have. In any case they are not that
important and can easily be added later.

Fixes #491

sisl/io/openmx/omx.py

@@ -311,11 +311,12 @@ def decompose_basis(l):
                 1: [1, -1, 0], # px, py, pz
                 2: [0, 2, -2, 1, -1], # d3z^2-r^2, dx^2-y^2, dxy, dxz, dyz
                 3: [0, 1, -1, 2, -2, 3, -3], # f5z^2-3r^2, f5xz^2-xr^2, f5yz^2-yr^2, fzx^2-zy^2, fxyz, fx^3-3*xy^2, f3yx^2-y^3
-                4: [0, 1, -1, 2, -2, 3, -3, 4, -4]
+                4: [0, 1, -1, 2, -2, 3, -3, 4, -4],
+                5: [0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5]
             }
             for i, c in enumerate(spec):
                 try:
-                    l = 'spdfg'.index(c)
+                    l = "spdfgh".index(c)
                     try:
                         nZ = int(spec[i+1])
                     except Exception:

sisl/orbital.py

@@ -33,13 +33,16 @@ def _rfact(l, m):
         return -msqrt(2*(2*l + 1)/pi4 * fact(l-m)/fact(l+m)) * (-1) ** m
     return msqrt(2*(2*l + 1)/pi4 * fact(l-m)/fact(l+m))
 
-# This is a list of dict
+# This is a tuple of dicts
 #  [0]{0} is l==0, m==0
 #  [1]{-1} is l==1, m==-1
 #  [1]{1} is l==1, m==1
 # and so on.
 # Calculate it up to l == 5 which is the h shell
-_rspher_harm_fact = [{m: _rfact(l, m) for m in range(-l, l+1)} for l in range(6)]
+_rspher_harm_fact = tuple(
+    {m: _rfact(l, m) for m in range(-l, l+1)}
+    for l in range(6)
+)
 # Clean-up
 del _rfact
 
@@ -884,7 +887,7 @@ def __init__(self, *args, **kwargs):
         self._zeta = zeta
         self._P = P
 
-        if self.l > 4:
+        if self.l > 5:
             raise ValueError(f"{self.__class__.__name__} does not implement shell h and above!")
         if abs(self.m) > self.l:
             raise ValueError(f"{self.__class__.__name__} requires |m| <= l.")
@@ -979,29 +982,33 @@ def equal(self, other, psi=False, radial=False):
     def name(self, tex=False):
         """ Return named specification of the atomic orbital """
         if tex:
-            name = "{}{}".format(self.n, "spdfg"[self.l])
+            name = "{}{}".format(self.n, "spdfgh"[self.l])
             if self.l == 1:
-                name += ["_y", "_z", "_x"][self.m+1]
+                name += ("_y", "_z", "_x")[self.m+1]
             elif self.l == 2:
-                name += ["_{xy}", "_{yz}", "_{z^2}", "_{xz}", "_{x^2-y^2}"][self.m+2]
+                name += ("_{xy}", "_{yz}", "_{z^2}", "_{xz}", "_{x^2-y^2}")[self.m+2]
             elif self.l == 3:
-                name += ["_{y(3x^2-y^2)}", "_{xyz}", "_{z^2y}", "_{z^3}", "_{z^2x}", "_{z(x^2-y^2)}", "_{x(x^2-3y^2)}"][self.m+3]
+                name += ("_{y(3x^2-y^2)}", "_{xyz}", "_{z^2y}", "_{z^3}", "_{z^2x}", "_{z(x^2-y^2)}", "_{x(x^2-3y^2)}")[self.m+3]
             elif self.l == 4:
-                name += ["_{_{xy(x^2-y^2)}}", "_{zy(3x^2-y^2)}", "_{z^2xy}", "_{z^3y}", "_{z^4}",
-                         "_{z^3x}", "_{z^2(x^2-y^2)}", "_{zx(x^2-3y^2)}", "_{x^4+y^4}"][self.m+4]
+                name += ("_{_{xy(x^2-y^2)}}", "_{zy(3x^2-y^2)}", "_{z^2xy}", "_{z^3y}", "_{z^4}",
+                         "_{z^3x}", "_{z^2(x^2-y^2)}", "_{zx(x^2-3y^2)}", "_{x^4+y^4}")[self.m+4]
+            elif self.l == 5:
+                name = f"{name}_{{m={self.m}}}"
             if self.P:
                 return name + fr"\zeta^{self.zeta}\mathrm{{P}}"
             return name + fr"\zeta^{self.zeta}"
-        name = "{}{}".format(self.n, "spdfg"[self.l])
+        name = "{}{}".format(self.n, "spdfgh"[self.l])
         if self.l == 1:
-            name += ["y", "z", "x"][self.m+1]
+            name += ("y", "z", "x")[self.m+1]
         elif self.l == 2:
-            name += ["xy", "yz", "z2", "xz", "x2-y2"][self.m+2]
+            name += ("xy", "yz", "z2", "xz", "x2-y2")[self.m+2]
         elif self.l == 3:
-            name += ["y(3x2-y2)", "xyz", "z2y", "z3", "z2x", "z(x2-y2)", "x(x2-3y2)"][self.m+3]
+            name += ("y(3x2-y2)", "xyz", "z2y", "z3", "z2x", "z(x2-y2)", "x(x2-3y2)")[self.m+3]
         elif self.l == 4:
-            name += ["xy(x2-y2)", "zy(3x2-y2)", "z2xy", "z3y", "z4",
-                     "z3x", "z2(x2-y2)", "zx(x2-3y2)", "x4+y4"][self.m+4]
+            name += ("xy(x2-y2)", "zy(3x2-y2)", "z2xy", "z3y", "z4",
+                     "z3x", "z2(x2-y2)", "zx(x2-3y2)", "x4+y4")[self.m+4]
+        elif self.l == 5:
+            name = f"{name}(m={self.m})"
         if self.P:
             return name + f"Z{self.zeta}P"
         return name + f"Z{self.zeta}"

sisl/tests/test_orbital.py

@@ -202,7 +202,7 @@ def test_toatomicorbital1(self):
         assert ao[0].m == 0
 
         # Check m and l
-        for l in range(1, 5):
+        for l in range(1, 6):
             orb = SphericalOrbital(l, rf)
             ao = orb.toAtomicOrbital()
             assert len(ao) == 2*l + 1
@@ -243,7 +243,7 @@ def test_toatomicorbital_q0(self):
         assert ao[0].q0 == 2.
 
         # Check m and l
-        for l in range(1, 5):
+        for l in range(1, 6):
             orb = SphericalOrbital(l, rf, 2.)
             ao = orb.toAtomicOrbital()
             assert ao[0].q0 == pytest.approx(2. / (2*l+1))
@@ -295,7 +295,7 @@ def test_init2(self):
 
     def test_init3(self):
         rf = r_f(6)
-        for l in range(5):
+        for l in range(6):
             a = AtomicOrbital(l=l, m=0, spherical=rf)
             a.name()
             a.name(True)
@@ -319,12 +319,12 @@ def test_init4(self):
 
     def test_init5(self):
         with pytest.raises(ValueError):
-            AtomicOrbital(5, 5, 0)
+            AtomicOrbital(5, 6, 0)
 
     def test_radial1(self):
         rf = r_f(6)
         r = np.linspace(0, 6, 100)
-        for l in range(5):
+        for l in range(6):
             so = SphericalOrbital(l, rf)
             sor = so.radial(r)
             for m in range(-l, l+1):
@@ -336,7 +336,7 @@ def test_radial1(self):
     def test_phi1(self):
         rf = r_f(6)
         r = np.linspace(0, 6, 999).reshape(-1, 3)
-        for l in range(5):
+        for l in range(6):
             so = SphericalOrbital(l, rf)
             for m in range(-l, l+1):
                 o = AtomicOrbital(l=l, m=m, spherical=rf)

toolbox/siesta/atom/_atom.py

@@ -73,7 +73,7 @@
     '4f', '5d', '6s', '6p', # [Rn]
     '7s', '5f', '6d', '7p' # [Og]
 ]
-_spdfg = 'spdfg'
+_spdfgh = 'spdfgh'
 
 
 class AtomInput:
@@ -168,7 +168,7 @@ def __init__(self, atom,
         #  _shell_order.index("2p") == 2
         # which has 1s and 2s occupied.
         try:
-            core = reduce(min, (_shell_order.index(f"{orb.n}{_spdfg[orb.l]}")
+            core = reduce(min, (_shell_order.index(f"{orb.n}{_spdfgh[orb.l]}")
                                 for orb in atom), len(_shell_order))
         except Exception:
             core = -1
@@ -472,10 +472,11 @@ def excite(self, *charge, **lq):
             charge = 0
         charge = lq.pop("charge", charge)
         # get indices of the orders
-        shell_idx = [_shell_order.index(f"{orb.n}{_spdfg[orb.l]}")
+        shell_idx = [_shell_order.index(f"{orb.n}{_spdfgh[orb.l]}")
                      for orb in self.atom]
         def _charge(idx):
-            return {'s': 2, 'p': 6, 'd': 10, 'f': 14}[_shell_order[idx][1]]
+            # while the g and h shells are never used, we just have them...
+            return {'s': 2, 'p': 6, 'd': 10, 'f': 14, 'g': 18, 'h': 22}[_shell_order[idx][1]]
         orig_charge = [_charge(idx) for idx in shell_idx]
         atom = self.atom.copy()
         # find order of orbitals (from highest index to lowest)
@@ -484,7 +485,7 @@ def _charge(idx):
         shell_sort.sort(reverse=charge > 0)
 
         for l, q in lq.items():
-            l = _spdfg.index(l)
+            l = _spdfgh.index(l)
             for orb in atom:
                 if orb.l == l:
                     orb._q0 -= q
@@ -575,13 +576,13 @@ def plot_wavefunction(ax):
 
                 r, w = get_xy(f"AEWFNR{il}")
                 if not r is None:
-                    p = ax.plot(r, w, label=f"AE {_spdfg[il]}")
+                    p = ax.plot(r, w, label=f"AE {_spdfgh[il]}")
                     color = p[0].get_color()
                     ax.axvline(orb.R * _Ang2Bohr, color=color, alpha=0.5)
 
                 r, w = get_xy(f"PSWFNR{il}")
                 if not r is None:
-                    ax.plot(r, w, '--', label=f"PS {_spdfg[il]}")
+                    ax.plot(r, w, '--', label=f"PS {_spdfgh[il]}")
 
             ax.set_xlim(0, atom_r * 5)
             ax.autoscale(enable=True, axis='y', tight=True)
@@ -663,7 +664,7 @@ def plot_log(ax):
                     emark.shape = (1, -1)
                 emark = emark[:, 0]
                 if not e is None:
-                    p = ax.plot(e, log, label=f"AE {_spdfg[il]}")
+                    p = ax.plot(e, log, label=f"AE {_spdfgh[il]}")
 
                     idx_mark = (np.fabs(e.reshape(-1, 1) - emark.reshape(1, -1))
                                 .argmin(axis=0))
@@ -676,7 +677,7 @@ def plot_log(ax):
                     emark.shape = (1, -1)
                 emark = emark[:, 0]
                 if not e is None:
-                    p = ax.plot(e, log, ':', label=f"PS {_spdfg[il]}")
+                    p = ax.plot(e, log, ':', label=f"PS {_spdfgh[il]}")
 
                     idx_mark = (np.fabs(e.reshape(-1, 1) - emark.reshape(1, -1))
                                 .argmin(axis=0))
@@ -695,7 +696,7 @@ def plot_potential(ax):
 
                 r, V = get_xy(f"PSPOTR{il}")
                 if not r is None:
-                    p = ax.plot(r, V, label=f"PS {_spdfg[il]}")
+                    p = ax.plot(r, V, label=f"PS {_spdfgh[il]}")
                     color = p[0].get_color()
                     ax.axvline(orb.R * _Ang2Bohr, color=color, alpha=0.5)
 

toolbox/siesta/minimizer/_atom_basis.py

@@ -94,7 +94,7 @@ def get_radius(orbs, zeta):
             if nl not in _shell_order:
                 continue
 
-            n, l = int(nl[0]), 'spdfg'.index(nl[1])
+            n, l = int(nl[0]), "spdfgh".index(nl[1])
             # Now we are sure we are dealing with valence shells
             basis = dic[nl].get("basis", {})
 
@@ -414,7 +414,7 @@ def add_variable(var):
                                     update_func=partial(update, d=self.opts, key="ion_charge")))
 
         # parse depending on shells in the atom
-        spdf = 'spdfg'
+        spdf = "spdfgh"
         for orb in self.atom:
             n, l = orb.n, orb.l
             nl = f"{n}{spdf[l]}"

toolbox/siesta/minimizer/_atom_pseudo.py

@@ -13,7 +13,7 @@
 
 
 _log = logging.getLogger("sisl_toolbox.siesta.minimize")
-_spdfg = 'spdfg'
+_spdfgh = "spdfgh"
 
 
 class AtomPseudo(AtomInput):
@@ -64,7 +64,7 @@ def add_variable(var):
 
         # parse depending on shells in the atom
         for orb in self.atom:
-            nl = f"{orb.n}{_spdfg[orb.l]}"
+            nl = f"{orb.n}{_spdfgh[orb.l]}"
             pseudo = dic.get(nl, {}).get("pseudo")
             if pseudo is None:
                 _log.info(f"{self.__class__.__name__} skipping node: {nl}.pseudo")