tb_harrison.py ダウンロード/コピー
tb_harrison.py
tb_harrison.py
1"""
2Slater-Kosterタイトバインディングモデルを実装するモジュール。
3
4s, p, d軌道を持つ有限クラスター（分子）用の最小限のSlater-Kosterタイトバインディング（TB）モデルを提供します。
5直交基底（S = I）と実対称ハミルトニアン（H）を仮定します。ハリスン則による積分パラメータのスケーリングに対応しています。
6
7:doc:`tb_harrison_usage`
8"""
9import numpy as np
10from scipy.linalg import eigh
11
12
13class SlaterKosterTB:
14 """Slater-Kosterタイトバインディング法を実装するクラス。
15
16 s, p, d軌道を持つ有限クラスター（分子）用の最小限のSlater-Kosterタイトバインディング（TB）モデルを提供します。
17 直交基底（S = I）と実対称ハミルトニアン（H）を仮定します。
18
19 ハリスン則によるスケーリングをサポートしており、`sk_params` の "harrison" エントリを通じて積分キーごとの指数を指定できます。
20
21 `harrison` パラメータの例:
22 ::
23
24 "harrison": {
25 "d0": 1.0,
26 "power_default": 2.0,
27 "power_map": {
28 "sd_sigma": 3.5,
29 "pd_sigma": 3.5,
30 "pd_pi": 3.5,
31 # "dd_sigma": 5.0, "dd_pi": 5.0, "dd_delta": 5.0, # if you implement dd later
32 },
33 # オプション: キーごとのd0
34 # "d0_map": {"pd_sigma": 2.0}
35 }
36
37 `harrison` が省略された場合、スケーリングは行われません（スケールは1）。
38 後方互換性のため、`{"d0":..., "power":...}` のみが与えられた場合、それがデフォルトとして使用されます。
39 """
40
41 def __init__(self):
42 """SlaterKosterTBクラスの新しいインスタンスを初期化します。
43
44 既知の軌道タイプ（s, p, d）とそれに対応する軌道のリストを設定し、
45 原子リストとハミルトニアン行列を空の状態に初期化します。
46 """
47 self.orbitals = {
48 "s": ["s"],
49 "p": ["px", "py", "pz"],
50 "d": ["dxy", "dyz", "dzx", "dx2-y2", "dz2"],
51 }
52 self.atoms = []
53 self.hamiltonian = None
54
55 def add_atom(self, symbol, orb_type, x, y, z):
56 """モデルに原子を追加します。
57
58 指定された記号、軌道タイプ、および空間座標を持つ原子を内部リストに追加します。
59 軌道タイプが`orbitals`辞書に存在しない場合はエラーを発生させます。
60
61 :param symbol: str: 原子の元素記号。
62 :param orb_type: str: 原子の軌道タイプ（例: "s", "p", "d"）。
63 :param x: float: 原子のX座標。
64 :param y: float: 原子のY座標。
65 :param z: float: 原子のZ座標。
66 :returns: None
67 :raises ValueError: 未知の軌道タイプが指定された場合。
68 """
69 if orb_type not in self.orbitals:
70 raise ValueError(f"Unknown orb_type='{orb_type}'. Choose from {list(self.orbitals.keys())}")
71 self.atoms.append(
72 {
73 "symbol": symbol,
74 "orb_type": orb_type,
75 "pos": np.array([x, y, z], dtype=float),
76 "orbs": self.orbitals[orb_type],
77 }
78 )
79
80 @staticmethod
81 def _need_param(p, key):
82 """指定されたSKパラメータが存在するか確認し、その値を返します。
83
84 `sk_params` 辞書から `key` に対応する値を取得します。
85 `key` が存在しない場合は `KeyError` を発生させます。
86
87 :param p: dict: Slater-Kosterパラメータを格納する辞書。
88 :param key: str: 取得するパラメータのキー。
89 :returns: float or int: 指定されたキーに対応するパラメータ値。
90 :raises KeyError: 指定されたキーが `sk_params` に存在しない場合。
91 """
92 if key not in p:
93 raise KeyError(f"Missing SK parameter '{key}' in sk_params.")
94 return p[key]
95
96 @staticmethod
97 def _harrison_scale(d, sk_params, key=None):
98 """ハリソン則に基づいて、距離によるスケーリング因子を計算します。
99
100 サイト間距離 `d` と `sk_params` の "harrison" セクションに定義されたルールに従って、
101 スケーリング因子を計算します。
102 `d0` (基準距離) と `power` (指数) は、グローバルまたはキーごとに指定できます。
103 `harrison` パラメータが指定されていない場合、または `d` が0以下の場合、スケーリングは行われず1.0を返します。
104 指数の選択順序は、1) `power_map[key]` -> 2) `power_default` -> 3) 従来の `power` です。
105
106 :param d: float: サイト間の距離。
107 :param sk_params: dict: Slater-Kosterパラメータを含む辞書。
108 :param key: str, optional: スケーリング対象の積分キー（例: "ss_sigma"）。キー固有のスケーリング設定に使用されます。デフォルトはNone。
109 :returns: float: 計算されたハリスン則によるスケーリング因子。
110 """
111 hs = sk_params.get("harrison", None)
112 if hs is None:
113 return 1.0
114 if d <= 0:
115 return 1.0
116
117 # d0 (キーごとに上書き可能)
118 d0 = float(hs.get("d0", 1.0))
119 d0_map = hs.get("d0_map", None)
120 if key is not None and isinstance(d0_map, dict) and key in d0_map:
121 d0 = float(d0_map[key])
122
123 # 指数 (power) の選択順序:
124 # 1) power_map[key]
125 # 2) power_default
126 # 3) 従来の power
127 power_map = hs.get("power_map", None)
128 if key is not None and isinstance(power_map, dict) and key in power_map:
129 power = float(power_map[key])
130 elif "power_default" in hs:
131 power = float(hs["power_default"])
132 else:
133 power = float(hs.get("power", 2.0)) # 従来の形式
134
135 return (d0 / d) ** power
136
137 def _scaled_param(self, d, sk_params, key):
138 """ハリスン則でスケーリングされたSlater-Kosterパラメータを返します。
139
140 まず `_need_param` を使用して元のパラメータ値を取得し、
141 次に `_harrison_scale` を使用してスケーリング因子を計算します。
142 最後に、元のパラメータ値にスケーリング因子を乗じて結果を返します。
143
144 :param d: float: サイト間の距離。
145 :param sk_params: dict: Slater-Kosterパラメータを含む辞書。
146 :param key: str: スケーリング対象のSlater-Koster積分キー。
147 :returns: float: ハリスン則でスケーリングされたパラメータ値。
148 """
149 v = self._need_param(sk_params, key)
150 s = self._harrison_scale(d, sk_params, key=key)
151 return v * s
152
153 def _get_sk_elements(self, vec, orb1, orb2, p):
154 """2つの軌道間のSlater-Koster行列要素を計算します。
155
156 2つの原子間の相対位置ベクトル `vec` とそれぞれの軌道タイプ `orb1`, `orb2` を基に、
157 ハミルトニアンのオフサイト要素を計算します。
158 距離 `d`、方向余弦 `l, m, n` を用いて、s-s, s-p, s-d, p-p, p-d の各相互作用を処理します。
159 d-d 相互作用は現在未実装です。
160 軌道間の対称性に基づいて、特定の相互作用では符号が反転することがあります。
161
162 :param vec: numpy.ndarray: 原子間の相対位置ベクトル `(x, y, z)`。
163 :param orb1: str: 1番目の原子の軌道タイプ（例: "s", "px", "dxy"）。
164 :param orb2: str: 2番目の原子の軌道タイプ（例: "s", "px", "dxy"）。
165 :param p: dict: Slater-Kosterパラメータを含む辞書。
166 :returns: float: 計算されたSlater-Koster行列要素の値。
167 :raises NotImplementedError: d-d ホッピングまたはサポートされていない軌道ペアが指定された場合。
168 """
169 d = np.linalg.norm(vec)
170 if d < 1e-12:
171 return 0.0
172
173 l, m, n = vec / d
174
175 # quick type flags
176 is_s1, is_s2 = (orb1 == "s"), (orb2 == "s")
177 is_p1, is_p2 = orb1.startswith("p"), orb2.startswith("p")
178 is_d1, is_d2 = orb1.startswith("d"), orb2.startswith("d")
179
180 # 未実装チャネルでの高速失敗
181 if is_d1 and is_d2:
182 raise NotImplementedError("d-d hopping (dd_sigma, dd_pi, dd_delta) is not implemented.")
183
184 # --- 1. s-s ---
185 if is_s1 and is_s2:
186 return self._scaled_param(d, p, "ss_sigma")
187
188 # --- 2. s-p ---
189 if is_s1 and is_p2:
190 vsp = self._scaled_param(d, p, "sp_sigma")
191 if orb2 == "px":
192 return l * vsp
193 if orb2 == "py":
194 return m * vsp
195 if orb2 == "pz":
196 return n * vsp
197
198 if is_p1 and is_s2:
199 return -self._get_sk_elements(vec, orb2, orb1, p)
200
201 # --- 3. s-d ---
202 if is_s1 and is_d2:
203 vsd = self._scaled_param(d, p, "sd_sigma")
204 if orb2 == "dz2":
205 return (n**2 - 0.5 * (l**2 + m**2)) * vsd
206 if orb2 == "dx2-y2":
207 return 0.5 * np.sqrt(3.0) * (l**2 - m**2) * vsd
208 if orb2 == "dxy":
209 return np.sqrt(3.0) * l * m * vsd
210 if orb2 == "dyz":
211 return np.sqrt(3.0) * m * n * vsd
212 if orb2 == "dzx":
213 return np.sqrt(3.0) * l * n * vsd
214
215 if is_d1 and is_s2:
216 return self._get_sk_elements(vec, orb2, orb1, p)
217
218 # --- 4. p-p ---
219 if is_p1 and is_p2:
220 vpps = self._scaled_param(d, p, "pp_sigma")
221 vppp = self._scaled_param(d, p, "pp_pi")
222 coord = {"px": l, "py": m, "pz": n}
223 c1, c2 = coord[orb1], coord[orb2]
224 if orb1 == orb2:
225 return c1**2 * vpps + (1.0 - c1**2) * vppp
226 else:
227 return c1 * c2 * (vpps - vppp)
228
229 # --- 5. p-d ---
230 if is_p1 and is_d2:
231 vpds = self._scaled_param(d, p, "pd_sigma")
232 vpdp = self._scaled_param(d, p, "pd_pi")
233
234 if orb1 == "px":
235 if orb2 == "dxy":
236 return m * (l**2 * vpds + (1 - 2 * l**2) * vpdp)
237 if orb2 == "dyz":
238 return l * m * n * (vpds - 2 * vpdp)
239 if orb2 == "dzx":
240 return n * (l**2 * vpds + (1 - 2 * l**2) * vpdp)
241 if orb2 == "dx2-y2":
242 return 0.5 * l * (l**2 - m**2) * vpds + l * (1 - l**2 + m**2) * vpdp
243 if orb2 == "dz2":
244 return l * (n**2 - 0.5 * (l**2 + m**2)) * vpds - np.sqrt(3.0) * l * n**2 * vpdp
245
246 if orb1 == "py":
247 if orb2 == "dxy":
248 return l * (m**2 * vpds + (1 - 2 * m**2) * vpdp)
249 if orb2 == "dyz":
250 return n * (m**2 * vpds + (1 - 2 * m**2) * vpdp)
251 if orb2 == "dzx":
252 return l * m * n * (vpds - 2 * vpdp)
253 if orb2 == "dx2-y2":
254 return 0.5 * m * (l**2 - m**2) * vpds - m * (1 + l**2 - m**2) * vpdp
255 if orb2 == "dz2":
256 return m * (n**2 - 0.5 * (l**2 + m**2)) * vpds - np.sqrt(3.0) * m * n**2 * vpdp
257
258 if orb1 == "pz":
259 if orb2 == "dxy":
260 return l * m * n * (vpds - 2 * vpdp)
261 if orb2 == "dyz":
262 return m * (n**2 * vpds + (1 - 2 * n**2) * vpdp)
263 if orb2 == "dzx":
264 return l * (n**2 * vpds + (1 - 2 * n**2) * vpdp)
265 if orb2 == "dx2-y2":
266 return 0.5 * n * (l**2 - m**2) * vpds - n * (l**2 - m**2) * vpdp
267 if orb2 == "dz2":
268 return n * (n**2 - 0.5 * (l**2 + m**2)) * vpds + np.sqrt(3.0) * n * (l**2 + m**2) * vpdp
269
270 if is_d1 and is_p2:
271 return -self._get_sk_elements(vec, orb2, orb1, p)
272
273 raise NotImplementedError(f"Unsupported orbital pair: ({orb1}, {orb2})")
274
275 def build_hamiltonian(self, sk_params, onsite_energies):
276 """システムのハミルトニアン行列を構築します。
277
278 追加された原子とその軌道情報、Slater-Kosterパラメータ、オンサイトエネルギーに基づいて、
279 全体のハミルトニアン行列を構築します。
280 まず、オンサイト項を対角要素に設定します。次に、異なる原子間のオフサイト項を
281 `_get_sk_elements` を使用して計算し、行列を埋めます。
282 ハミルトニアンは実対称行列として構築されます。
283
284 :param sk_params: dict: Slater-Kosterパラメータを含む辞書。
285 :param onsite_energies: dict: 原子の記号と軌道タイプごとのオンサイトエネルギーを含む辞書。
286 :returns: numpy.ndarray: 構築されたハミルトニアン行列。
287 """
288 total_orbs = sum(len(a["orbs"]) for a in self.atoms)
289 H = np.zeros((total_orbs, total_orbs), dtype=float)
290
291 # Onsite terms
292 idx = 0
293 for a in self.atoms:
294 n_orb = len(a["orbs"])
295 e_on = onsite_energies.get(a["symbol"], {}).get(a["orb_type"], 0.0)
296 for k in range(n_orb):
297 H[idx + k, idx + k] = e_on
298 idx += n_orb
299
300 # Offsite terms
301 idx_i = 0
302 for i, atom_i in enumerate(self.atoms):
303 idx_j = 0
304 for j, atom_j in enumerate(self.atoms):
305 if i < j:
306 vec = atom_j["pos"] - atom_i["pos"]
307 for io, orb_i in enumerate(atom_i["orbs"]):
308 for jo, orb_j in enumerate(atom_j["orbs"]):
309 val = self._get_sk_elements(vec, orb_i, orb_j, sk_params)
310 H[idx_i + io, idx_j + jo] = val
311 H[idx_j + jo, idx_i + io] = val
312 idx_j += len(atom_j["orbs"])
313 idx_i += len(atom_i["orbs"])
314
315 self.hamiltonian = H
316 return H
317
318 def orbital_labels(self):
319 """各原子の軌道に対応するラベルのリストを生成します。
320
321 各原子の位置と軌道タイプを組み合わせて、ハミルトニアン行列の各要素に対応する
322 識別のための文字列ラベルを生成します。
323 例: "Ba(0.000,0.000,0.000)_s", "N(0.000,0.000,0.800)_px"
324
325 :returns: list[str]: 各軌道に対応するラベルのリスト。
326 """
327 labels = []
328 for a in self.atoms:
329 x, y, z = a["pos"]
330 for o in a["orbs"]:
331 labels.append(f"{a['symbol']}({x:.3f},{y:.3f},{z:.3f})_{o}")
332 return labels
333
334
335def main():
336 """`SlaterKosterTB` クラスの使用例を示します。
337
338 `SlaterKosterTB` クラスを初期化し、複数の原子をモデルに追加します。
339 定義済みのSlater-Kosterパラメータとオンサイトエネルギーを使用してハミルトニアンを構築し、
340 `scipy.linalg.eigh` を使用してハミルトニアンを対角化し、固有値（軌道エネルギー）と固有ベクトルを計算します。
341 結果として得られたエネルギーと、それぞれの固有状態における支配的な軌道成分を表示します。
342 """
343 tb = SlaterKosterTB()
344
345 tb.add_atom("Ba", "d", 0, 0, 0)
346 tb.add_atom("N", "p", 0, 0, 0.8)
347 tb.add_atom("N", "p", 1, 0, 0)
348 tb.add_atom("N", "p", -1, 0, 0)
349 tb.add_atom("N", "p", 0, 1, 0)
350 tb.add_atom("N", "p", 0, -1, 0)
351
352 params = {
353 "ss_sigma": -1.40,
354 "sp_sigma": 1.84,
355 "pp_sigma": 3.24,
356 "pp_pi": -0.81,
357 "sd_sigma": -2.0,
358 "pd_sigma": -2.2,
359 "pd_pi": 1.0,
360 # キーに依存するハリスン則スケーリング:
361 "harrison": {
362 "d0": 1.0,
363 "power_default": 2.0, # ss/sp/pp はデフォルトで 2.0
364 "power_map": { # 選択された積分でのオーバーライド
365 "sd_sigma": 3.5,
366 "pd_sigma": 3.5,
367 "pd_pi": 3.5,
368 },
369 },
370 }
371
372 onsite = {"Ba": {"d": 0.0}, "N": {"p": -3.0}}
373
374 H = tb.build_hamiltonian(params, onsite)
375 energies, vectors = eigh(H)
376
377 idx_sort = energies.argsort()
378 evals = energies[idx_sort]
379 evecs = vectors[:, idx_sort]
380 labels = tb.orbital_labels()
381
382 print("Eigenvalues (Orbital Energies) [eV]:")
383 print(evals)
384
385 print("\nNo. | Energy (eV) | Dominant components")
386 print("-" * 72)
387 for i in range(len(evals)):
388 coeffs = evecs[:, i]
389 w = coeffs**2
390 top = np.argsort(w)[::-1][:5]
391# desc = ", ".join([f"{labels[k]}:{coeffs[k]:+.2f}" for k in top])
392# print(f"{i+1:2d} | {evals[i]:8.4f} | {desc}")
393 print(f"{i+1:2d} | {evals[i]:8.4f}")
394 for l, c in zip(labels, coeffs):
395 print(f" {l}:{c:+.3f}")
396
397
398if __name__ == "__main__":
399 main()