Coverage for tbkit/slater_koster.py: 100%
128 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 06:13 +0000
1r"""
2Slater-Koster two-centre integrals (J. C. Slater and G. F. Koster, Phys.
3Rev. 94, 1498 (1954)) between s, p and d orbitals.
5The matrix element between orbital :math:`\alpha` on one atom and orbital
6:math:`\beta` on a neighbour at :math:`\mathbf{d} = d\,(l, m, n)` depends
7only on a handful of bond integrals :math:`V_{ll'\mu}`
8(:math:`\mu = \sigma, \pi, \delta`: the angular momentum about the bond),
9through the direction cosines :math:`(l, m, n)` -- e.g.
11.. math::
13 E_{s,x} = l\,V_{sp\sigma}\, ,\qquad
14 E_{x,x} = l^2 V_{pp\sigma} + (1-l^2) V_{pp\pi}\, ,\qquad
15 E_{x,y} = lm\,(V_{pp\sigma} - V_{pp\pi})\, .
17Rather than tabulating all 81 such formulas, they are generated: in a
18frame whose z axis lies along the bond, the matrix is diagonal in the
19angular momentum about the bond, with entries :math:`V_{ll'\mu}`; rotating
20the frame onto :math:`\hat{\mathbf{d}}` gives every entry,
21:math:`E(\hat{\mathbf{d}}) = D_{l}(R)\,E(\hat z)\,D_{l'}(R)^T`, with
22:math:`D_1 = R` for the p orbitals and :math:`D_2` its action on the
23quadratic forms of the d orbitals. This is exactly the Slater-Koster table
24(checked against it in the tests), in the usual real-orbital basis:
26``'s'``, ``'px'``, ``'py'``, ``'pz'``, ``'dxy'``, ``'dyz'``, ``'dzx'``,
27``'dx2-y2'``, ``'dz2'``.
29Bond integrals are given as a dictionary, with keys ``'ss_sigma'``,
30``'sp_sigma'``, ``'pp_sigma'``, ``'pp_pi'``, ``'sd_sigma'``,
31``'pd_sigma'``, ``'pd_pi'``, ``'dd_sigma'``, ``'dd_pi'``, ``'dd_delta'``
32(missing keys are zero). Between different atoms, :math:`V_{sp\sigma}`
33(s on the first atom, p on the second) and :math:`V_{ps\sigma}` (p on the
34first, s on the second) differ: give ``'ps_sigma'``, ``'ds_sigma'``,
35``'dp_sigma'``, ``'dp_pi'`` for the latter (they default to their
36homonuclear counterparts).
37"""
38from __future__ import annotations
40import numpy as np
41from numpy.typing import ArrayLike, NDArray
43import tbkit.error_handling as error_handling
46#: Real orbitals, with their angular momentum l.
47ORBITALS = {'s': 0, 'px': 1, 'py': 1, 'pz': 1, 'dxy': 2, 'dyz': 2, 'dzx': 2,
48 'dx2-y2': 2, 'dz2': 2}
50#: Bond-integral keys (the heteronuclear ps/ds/dp keys are optional).
51PARAMS = ('ss_sigma', 'sp_sigma', 'ps_sigma', 'pp_sigma', 'pp_pi', 'sd_sigma',
52 'ds_sigma', 'pd_sigma', 'pd_pi', 'dp_sigma', 'dp_pi', 'dd_sigma',
53 'dd_pi', 'dd_delta')
55_P = ['px', 'py', 'pz']
56_D = ['dxy', 'dyz', 'dzx', 'dx2-y2', 'dz2']
57_S3 = np.sqrt(3.)
58# the d orbitals as quadratic forms r^T M r, all of Frobenius norm^2 = 3/2
59_M = {'dxy': np.array([[0., _S3/2, 0.], [_S3/2, 0., 0.], [0., 0., 0.]]),
60 'dyz': np.array([[0., 0., 0.], [0., 0., _S3/2], [0., _S3/2, 0.]]),
61 'dzx': np.array([[0., 0., _S3/2], [0., 0., 0.], [_S3/2, 0., 0.]]),
62 'dx2-y2': np.diag([_S3/2, -_S3/2, 0.]),
63 'dz2': np.diag([-0.5, -0.5, 1.])}
66def _rotation(d_hat: NDArray[np.float64]) -> NDArray[np.float64]:
67 '''
68 Private function. A rotation matrix taking z onto d_hat.
69 '''
70 z = np.array([0., 0., 1.])
71 c = float(np.dot(z, d_hat))
72 if c > 1. - 1e-14:
73 return np.eye(3)
74 if c < -1. + 1e-14:
75 return np.diag([1., -1., -1.]) # a half turn about x
76 axis = np.cross(z, d_hat)
77 s = np.linalg.norm(axis)
78 k = axis / s
79 kx = np.array([[0., -k[2], k[1]], [k[2], 0., -k[0]], [-k[1], k[0], 0.]])
80 return np.eye(3) + s * kx + (1. - c) * kx @ kx
83def _rep(l: int, rot: NDArray[np.float64]) -> NDArray[np.float64]:
84 '''
85 Private function. Real representation D_l(R) of a rotation on the real
86 orbitals of angular momentum l (columns: rotated orbitals).
87 '''
88 if l == 0:
89 return np.ones((1, 1))
90 if l == 1:
91 return rot
92 return np.array([[np.sum(_M[a] * (rot @ _M[b] @ rot.T)) / 1.5 for b in _D] for a in _D])
95def _bond_frame(l1: int, l2: int, params: dict) -> NDArray[np.float64]:
96 '''
97 Private function. Matrix between the orbitals of angular momentum l1 (first
98 atom) and l2 (second atom, at +z), in the bond frame.
99 '''
100 def get(key, fallback=None):
101 return params.get(key, params.get(fallback, 0.) if fallback else 0.)
103 if (l1, l2) == (0, 0):
104 return np.array([[get('ss_sigma')]])
105 if (l1, l2) == (0, 1):
106 return np.array([[0., 0., get('sp_sigma')]])
107 if (l1, l2) == (1, 0):
108 return np.array([[0.], [0.], [-get('ps_sigma', 'sp_sigma')]])
109 if (l1, l2) == (1, 1):
110 return np.diag([get('pp_pi'), get('pp_pi'), get('pp_sigma')])
111 if (l1, l2) == (0, 2):
112 return np.array([[0., 0., 0., 0., get('sd_sigma')]])
113 if (l1, l2) == (2, 0):
114 return np.array([[0.], [0.], [0.], [0.], [get('ds_sigma', 'sd_sigma')]])
115 if (l1, l2) == (1, 2):
116 # (px, dzx), (py, dyz): pi; (pz, dz2): sigma
117 out = np.zeros((3, 5))
118 out[0, 2] = out[1, 1] = get('pd_pi')
119 out[2, 4] = get('pd_sigma')
120 return out
121 if (l1, l2) == (2, 1):
122 out = np.zeros((5, 3))
123 out[2, 0] = out[1, 1] = -get('dp_pi', 'pd_pi')
124 out[4, 2] = -get('dp_sigma', 'pd_sigma')
125 return out
126 return np.diag([get('dd_delta'), get('dd_pi'), get('dd_pi'), get('dd_delta'), get('dd_sigma')])
129def sk_block(orbs1: list[str], orbs2: list[str], d: ArrayLike, params: dict) -> NDArray[np.float64]:
130 r'''
131 Get the Slater-Koster hopping matrix between the orbitals *orbs1* of an
132 atom and the orbitals *orbs2* of a neighbour at displacement :math:`\mathbf{d}`.
134 :param orbs1: List of orbital names (see :data:`ORBITALS`).
135 :param orbs2: List of orbital names.
136 :param d: Tuple of 2 or 3 reals. Displacement from the first atom to the
137 second (a 2D vector lies in the (x, y) plane).
138 :param params: Dictionary of bond integrals (see the module docstring).
140 :returns:
141 * **block** -- Real ndarray, shape (len(orbs1), len(orbs2)).
142 '''
143 error_handling.orbitals(orbs1)
144 error_handling.orbitals(orbs2)
145 error_handling.sk_params(params)
146 d = np.asarray(d, dtype='f8')
147 error_handling.bond_vector(d)
148 d = np.concatenate([d, np.zeros(3 - len(d))])
149 rot = _rotation(d / np.linalg.norm(d))
150 names = {0: ['s'], 1: _P, 2: _D}
151 block = np.zeros((len(orbs1), len(orbs2)))
152 for l1 in {ORBITALS[o] for o in orbs1}:
153 for l2 in {ORBITALS[o] for o in orbs2}:
154 full = _rep(l1, rot) @ _bond_frame(l1, l2, params) @ _rep(l2, rot).T
155 rows = [(i, names[l1].index(o)) for i, o in enumerate(orbs1) if ORBITALS[o] == l1]
156 cols = [(j, names[l2].index(o)) for j, o in enumerate(orbs2) if ORBITALS[o] == l2]
157 for i, a in rows:
158 for j, b in cols:
159 block[i, j] = full[a, b]
160 return block
163def orbital_angular_momentum() -> list[NDArray[np.complex128]]:
164 r'''
165 Get the orbital angular momentum :math:`(L_x, L_y, L_z)` (:math:`\hbar = 1`)
166 on the real p orbitals ``(px, py, pz)``: :math:`(L_k)_{ij} = -i\epsilon_{kij}`.
168 :returns:
169 * **L** -- List of three complex ndarrays, shape (3, 3).
170 '''
171 eps = np.zeros((3, 3, 3))
172 eps[0, 1, 2] = eps[1, 2, 0] = eps[2, 0, 1] = 1.
173 eps[0, 2, 1] = eps[2, 1, 0] = eps[1, 0, 2] = -1.
174 return [-1j * eps[k] for k in range(3)]
177def sk_kspace(lat, orbitals: dict[str, list[str]], bonds: dict[int, dict],
178 onsite: dict | None = None, overlap: dict[int, dict] | None = None):
179 r'''
180 Get the Slater-Koster Bloch Hamiltonian of a periodic lattice: a
181 **KSpace** with one orbital per (site, orbital) of the unit cell, and
182 the hoppings of *sk_block* between the *n*-th neighbours listed in
183 *bonds* -- the reciprocal-space twin of
184 :meth:`tbkit.orbital.OrbitalSystem.set_slater_koster`.
186 :param lat: **Lattice** instance (unit cell and primitive vectors).
187 :param orbitals: Dictionary {tag: list of orbital names}.
188 :param bonds: Dictionary {n: bond integrals}: the neighbour orders to
189 connect (1 for nearest neighbours), each with its Slater-Koster
190 parameters (see the module docstring).
191 :param onsite: Dictionary. Default value None. {tag: energy} (every
192 orbital of the site) or {tag: {orbital: energy}}.
193 :param overlap: Dictionary. Default value None. {n: bond integrals} of
194 the overlap matrix (see *KSpace.set_overlap*).
196 :returns:
197 * **ks** -- **KSpace** instance; *ks.sk_orbitals* lists the (site,
198 orbital) of each of its orbitals.
199 '''
200 from tbkit.kspace import KSpace
201 from tbkit.lattice import Lattice
202 error_handling.lat(lat)
203 error_handling.orbital_dict(orbitals, lat.tags)
204 error_handling.sk_bonds(bonds)
205 if overlap is not None:
206 error_handling.sk_bonds(overlap)
207 cell = lat.unit_cell
208 table = [(i, o) for i, dic in enumerate(cell) for o in orbitals[dic['tag']]]
209 index = {key: n for n, key in enumerate(table)}
210 new_cell = [{'tag': cell[i]['tag'], 'r0': cell[i]['r0']} for i, _ in table]
211 ks = KSpace(Lattice(unit_cell=new_cell, prim_vec=lat.prim_vec))
212 setattr(ks, 'sk_orbitals', table) # an attribute of this KSpace only
213 if onsite is not None:
214 error_handling.set_onsite_orb(onsite, lat.tags, orbitals)
215 for n, (i, o) in enumerate(table):
216 val = onsite.get(cell[i]['tag'], 0.)
217 ks.onsite[n] = val.get(o, 0.) if isinstance(val, dict) else val
218 # candidate bonds: every site of the home cell to every site of the
219 # neighbouring cells, each unordered bond once
220 a = np.array(lat.prim_vec, dtype='f8')
221 r0 = np.array([dic['r0'] for dic in cell], dtype='f8')
222 shifts = np.array(np.meshgrid(*[range(-3, 4)] * len(a), indexing='ij')).reshape(len(a), -1).T
223 cand = []
224 for shift in shifts:
225 for i in range(len(cell)):
226 for j in range(len(cell)):
227 d = r0[j] + shift @ a - r0[i]
228 dist = np.linalg.norm(d)
229 if dist < 1e-9:
230 continue
231 if (i, j) > (j, i) or (i == j and tuple(-shift) > tuple(shift)):
232 continue
233 cand.append((round(dist, 4), i, j, tuple(int(c) for c in shift), d))
234 distances = sorted({c[0] for c in cand})
235 for store, spec in ((ks.set_hopping, bonds), (ks.set_overlap, overlap or {})):
236 hops = []
237 for n, params in spec.items():
238 error_handling.positive_int_lim(n, 'n', len(distances))
239 for dist, i, j, R, d in cand:
240 if dist != distances[n - 1]:
241 continue
242 orbs_i, orbs_j = orbitals[cell[i]['tag']], orbitals[cell[j]['tag']]
243 block = sk_block(orbs_i, orbs_j, d, params)
244 for p, oi in enumerate(orbs_i):
245 for q, oj in enumerate(orbs_j):
246 if block[p, q]:
247 hops.append({'i': index[(i, oi)], 'j': index[(j, oj)], 'R': R, 't': float(block[p, q])})
248 if hops:
249 store(hops)
250 return ks