Coverage for tbkit/bdg.py: 100%
164 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-30 11:04 +0000
1r"""
2Superconductivity in the Bogoliubov-de Gennes (BdG) formalism. With a
3pairing term :math:`\sum_{ij}\Delta_{ij}c^\dagger_ic^\dagger_j + h.c.`,
4the Hamiltonian is :math:`\frac12\Psi^\dagger H_{BdG}\Psi` in the Nambu
5basis :math:`\Psi = (c_1, \dots, c_N, c_1^\dagger, \dots, c_N^\dagger)^T`,
7.. math::
9 H_{BdG} = \begin{pmatrix} h - \mu & \Delta \\ \Delta^\dagger & -(h-\mu)^T \end{pmatrix}\, ,
11with :math:`\Delta^T = -\Delta` (Fermi statistics). Its spectrum is
12symmetric, :math:`\pm E_n`, by the particle-hole symmetry
13:math:`\tau_x H_{BdG}^*\tau_x = -H_{BdG}`; the quasiparticle energies are
14the :math:`E_n \ge 0`, and a zero mode is a Majorana mode. Real-space
15models come from *bdg_ham*, periodic ones from *bdg_kspace* (a
16:class:`tbkit.kspace.KSpace` with the holes as extra orbitals, so its bands,
17Berry phases, *finite_ham*, ... all apply).
18"""
19from __future__ import annotations
21import numpy as np
22from numpy.typing import NDArray
23import scipy.sparse as sparse
25import tbkit.error_handling as error_handling
26from tbkit.lattice import Lattice
29def bdg_ham(ham, pairing, mu: float = 0.) -> sparse.csr_matrix:
30 r'''
31 Get the real-space BdG Hamiltonian (see the module docstring).
33 Parameters
34 ----------
35 ham : Square Hermitian matrix, shape (N, N) (e.g. *System.ham*)
36 pairing
37 Square antisymmetric matrix, shape (N, N): :math:`\Delta`
38 (see *pairing_bonds*, *pairing_s_wave*).
39 mu : Real number, default 0
40 Chemical potential.
42 Returns
43 -------
44 ham_bdg : Sparse CSR matrix, shape (2N, 2N)
45 '''
46 h = sparse.csr_matrix(ham, dtype='c16')
47 delta = sparse.csr_matrix(pairing, dtype='c16')
48 error_handling.square_matrix(h, 'ham')
49 error_handling.hermitian(h)
50 error_handling.pairing(delta, h.shape[0])
51 error_handling.real_number(mu, 'mu')
52 h = h - mu * sparse.identity(h.shape[0], dtype='c16', format='csr')
53 return sparse.csr_matrix(sparse.bmat([[h, delta], [delta.conj().T, -h.T]]))
56def pairing_bonds(sys, delta: complex) -> sparse.csr_matrix:
57 r'''
58 Get a p-wave (odd) pairing on the bonds of *sys.hop*:
59 :math:`\Delta_{ij} = \Delta`, :math:`\Delta_{ji} = -\Delta` for every
60 stored hopping :math:`i\to j` (oriented as *System.set_hopping* stores
61 them: angle in :math:`[0, 180)`). On a chain it is Kitaev's
62 :math:`\Delta\,c^\dagger_jc^\dagger_{j+1}`.
64 Parameters
65 ----------
66 sys : **System** instance, with hoppings
67 delta : Complex number
68 :math:`\Delta`.
70 Returns
71 -------
72 pairing : Sparse CSR matrix, shape (sites, sites)
73 '''
74 error_handling.empty_hop(sys.hop)
75 error_handling.number(delta, 'delta')
76 i, j = sys.hop['i'].astype(int), sys.hop['j'].astype(int)
77 n = sys.lat.sites
78 vals = np.full(len(i), delta, dtype='c16')
79 return sparse.csr_matrix((np.concatenate([vals, -vals]), (np.concatenate([i, j]), np.concatenate([j, i]))),
80 shape=(n, n))
83def pairing_s_wave(n_sites: int, delta: complex) -> sparse.csr_matrix:
84 r'''
85 Get an onsite singlet (s-wave) pairing
86 :math:`\Delta\sum_i c^\dagger_{i\uparrow}c^\dagger_{i\downarrow} + h.c.`
87 for a spinful model with rows ordered site by site, spin up then down
88 (as *KSpace* with ``spin=True`` and *OrbitalSystem* with one orbital).
90 Parameters
91 ----------
92 n_sites : Positive integer
93 Number of sites.
94 delta : Complex number
95 :math:`\Delta`.
97 Returns
98 -------
99 pairing : Sparse CSR matrix, shape (2 n_sites, 2 n_sites)
100 '''
101 error_handling.positive_int(n_sites, 'n_sites')
102 error_handling.number(delta, 'delta')
103 block = np.array([[0., delta], [-delta, 0.]], dtype='c16')
104 return sparse.csr_matrix(sparse.kron(sparse.identity(n_sites), block))
107def bdg_kspace(ks, pairing: list[dict], mu: float = 0.):
108 r'''
109 Get the BdG Bloch Hamiltonian of a periodic model,
111 .. math::
113 H_{BdG}(\mathbf{k}) = \begin{pmatrix} h(\mathbf{k})-\mu & \Delta(\mathbf{k}) \\
114 \Delta(\mathbf{k})^\dagger & -[h(-\mathbf{k})-\mu]^T\end{pmatrix}\, ,\qquad
115 \Delta(\mathbf{k}) = \sum_{\mathbf{R}}\Delta_{ij}(\mathbf{R})\,e^{i\mathbf{k}\cdot\mathbf{R}}\, ,
117 as a new **KSpace** whose orbitals are those of *ks* (particles,
118 ``0 ... norb-1``) followed by their holes (``norb ... 2 norb-1``, at the
119 same positions). Each pairing :math:`\Delta\,c^\dagger_{i,\mathbf{0}}c^\dagger_{j,\mathbf{R}} + h.c.`
120 enters antisymmetrized, :math:`\Delta_{ij}(\mathbf{R}) = -\Delta_{ji}(-\mathbf{R})`.
122 Parameters
123 ----------
124 ks : **KSpace** instance (Hermitian, orthogonal basis)
125 pairing
126 List of dictionaries with keys ('i', 'j', 'R', 'delta'):
127 orbital indices of *ks* (``0 ... norb-1``; with ``spin=True``,
128 ``2*site + spin``), lattice vector, amplitude.
129 mu : Real number, default 0
130 Chemical potential.
132 Returns
133 -------
134 ks_bdg : **KSpace** instance, 2 norb orbitals
136 Examples
137 --------
138 Example usage::
140 # Kitaev chain: E(k) = +-sqrt((2t cos k - mu)^2 + 4 Delta^2 sin^2 k)
141 chain = KSpace(lattices.chain())
142 chain.set_hopping([{'i': 0, 'j': 0, 'R': (1,), 't': t}])
143 kit = bdg_kspace(chain, [{'i': 0, 'j': 0, 'R': (1,), 'delta': delta}], mu=mu)
144 '''
145 from tbkit.kspace import KSpace
146 error_handling.hermitian_kspace(ks.is_hermitian(), ks._overlap_hop)
147 error_handling.pairing_kspace(pairing, ks.norb, ks.dim)
148 error_handling.real_number(mu, 'mu')
149 norb = ks.norb
150 tau = ks.orbital_positions()
151 tags = np.repeat(ks.tags, 2) if ks.spin else ks.tags
152 cell = [{'tag': str(t), 'r0': tuple(float(c) for c in r)} for t, r in zip(tags, tau)]
153 new = KSpace(Lattice(unit_cell=cell + [dict(c) for c in cell], prim_vec=ks.lat.prim_vec))
154 h0 = np.diag(ks.onsite) + ks._onsite_offdiag - mu * np.eye(norb)
155 new.onsite = np.concatenate([np.diag(h0), -np.diag(h0)]).astype('c16')
156 off = h0 - np.diag(np.diag(h0))
157 new._onsite_offdiag = np.block([[off, np.zeros((norb, norb))], [np.zeros((norb, norb)), -off.T]]).astype('c16')
158 hops = [(i, j, R, t) for i, j, R, t in ks._hop]
159 hops += [(j + norb, i + norb, -R, -t) for i, j, R, t in ks._hop]
160 for dic in pairing:
161 R = np.zeros(ks.space_dim)
162 for n, a in zip(dic['R'], ks.lat.prim_vec):
163 R += n * np.array(a)
164 i, j, d = dic['i'], dic['j'], complex(dic['delta'])
165 hops += [(i, j + norb, R, d), (j, i + norb, -R, -d),
166 (j + norb, i, -R, np.conj(d)), (i + norb, j, R, -np.conj(d))]
167 new._hop = hops
168 return new
171def particle_hole(n: int) -> NDArray[np.float64]:
172 r'''
173 Get the unitary part :math:`\tau_x` of the BdG particle-hole symmetry
174 :math:`C = \tau_x\mathcal{K}` (:math:`C^2 = +1`), for *n* particle
175 orbitals: pass it to *KSpace.symmetry_error* (``antiunitary=True``,
176 ``anti=True``) or *KSpace.tenfold_class* (``particle_hole=``).
178 Parameters
179 ----------
180 n : Positive integer
181 Number of particle orbitals.
183 Returns
184 -------
185 tau_x : Real ndarray, shape (2n, 2n)
186 '''
187 error_handling.positive_int(n, 'n')
188 return np.kron(np.array([[0., 1.], [1., 0.]]), np.eye(n))
191def pfaffian(a) -> complex:
192 r'''
193 Get the Pfaffian of an antisymmetric matrix,
194 :math:`\mathrm{Pf}(A)^2 = \det A`, by the Parlett-Reid reduction to
195 tridiagonal form with pivoting (Wimmer, ACM Trans. Math. Softw. 38, 30
196 (2012)). For a 2x2 block :math:`\begin{pmatrix}0 & a\\ -a & 0\end{pmatrix}`
197 it is :math:`a`; it is 0 for an odd dimension.
199 Parameters
200 ----------
201 a : Square antisymmetric matrix, shape (n, n)
203 Returns
204 -------
205 pf : Complex number (a float for a real matrix)
206 '''
207 real = not np.iscomplexobj(a)
208 a = np.array(a, dtype='c16')
209 error_handling.square_matrix(a, 'a')
210 error_handling.antisymmetric(a)
211 n = len(a)
212 if n % 2:
213 return 0.
214 pf = 1. + 0j
215 for k in range(0, n - 1, 2):
216 # bring the largest entry of column k below the diagonal to row k + 1
217 p = k + 1 + int(np.argmax(np.abs(a[k + 1:, k])))
218 if p != k + 1:
219 a[[k + 1, p], :] = a[[p, k + 1], :]
220 a[:, [k + 1, p]] = a[:, [p, k + 1]]
221 pf = -pf
222 if a[k + 1, k] == 0.:
223 return 0.
224 pf *= a[k, k + 1]
225 if k + 2 < n:
226 tau = a[k, k + 2:] / a[k, k + 1]
227 a[k + 2:, k + 2:] += np.outer(tau, a[k + 2:, k + 1]) - np.outer(a[k + 2:, k + 1], tau)
228 return float(pf.real) if real else complex(pf)
231def majorana_number(ks) -> int:
232 r'''
233 Get Kitaev's Majorana number of a BdG model (Kitaev, Physics-Uspekhi
234 44, 131 (2001)): the product, over the time-reversal-invariant momenta
235 :math:`\Gamma_i`, of the signs of the Pfaffians of the Hamiltonian in
236 the Majorana basis,
238 .. math::
240 \mathcal{M} = \prod_i \mathrm{sign}\,\mathrm{Pf}\,A(\Gamma_i)\, ,\qquad
241 A = -i\,U H_{BdG} U^\dagger\, ,\qquad
242 U = \frac{1}{\sqrt2}\begin{pmatrix}1 & 1\\ -i & i\end{pmatrix}\, ,
244 :math:`U` taking the Nambu spinor :math:`(c, c^\dagger)` to the Majorana
245 operators :math:`(c + c^\dagger, -i(c - c^\dagger))`. Particle-hole
246 symmetry makes :math:`A(\Gamma_i)` real and antisymmetric. In 1D,
247 :math:`\mathcal{M} = \mathrm{sign}\,\mathrm{Pf}A(0)\,\mathrm{Pf}A(\pi) = -1`
248 in the topological phase (a Majorana zero mode at each end of an open
249 chain), +1 in the trivial one: the :math:`\mathbb{Z}_2` invariant of
250 class D. In 2D the product over the four momenta is
251 :math:`(-1)^C`, the parity of the Chern number (Ghosh et al., Phys. Rev.
252 B 82, 184525 (2010)). The invariant changes only where the gap closes
253 at one of the :math:`\Gamma_i`.
255 Parameters
256 ----------
257 ks
258 **KSpace** instance of a BdG model, 2n orbitals: n particles
259 then their n holes (see *bdg_kspace*).
261 Returns
262 -------
263 m : Integer, -1 (topological) or +1 (trivial)
264 '''
265 from tbkit.kspace import KSpace
266 error_handling.kspace(ks, KSpace)
267 error_handling.bdg_orbitals(ks.norb)
268 n = ks.norb // 2
269 u = np.kron(np.array([[1., 1.], [-1j, 1j]]) / np.sqrt(2), np.eye(n))
270 sign = 1
271 for k in ks._trims():
272 a = -1j * u @ ks.get_ham(k) @ u.conj().T
273 error_handling.majorana_form(a)
274 pf = pfaffian(a.real)
275 error_handling.gapped_trim(pf)
276 sign *= int(np.sign(pf))
277 return sign
280class GapResult():
281 r'''
282 Self-consistent s-wave pairing of *s_wave_gap*.
284 Attributes
285 ----------
286 delta : Complex ndarray
287 Gap :math:`\Delta_i` on every site.
288 energies : Real ndarray, shape (2N,)
289 Spectrum of the reduced BdG
290 Hamiltonian (symmetric, :math:`\pm E_n`).
291 states : Complex ndarray, shape (2N, 2N)
292 Its eigenvectors
293 :math:`(u_n, v_n)` (columns).
294 iterations : Integer
295 Number of iterations.
296 '''
298 def __init__(self, delta, energies, states, iterations) -> None:
299 self.delta = delta
300 self.energies = energies
301 self.states = states
302 self.iterations = iterations
304 @property
305 def mean_gap(self) -> float:
306 r'''
307 Average :math:`|\Delta_i|` over the sites.
308 '''
309 return float(np.mean(np.abs(self.delta)))
312def s_wave_gap(
313 ham, V: float | NDArray, mu: float = 0., temperature: float = 0., delta0: complex | NDArray = 0.1,
314 mixing: float = 1., tol: float = 1e-8, max_iter: int = 10000,
315) -> GapResult:
316 r'''
317 Solve the real-space gap equation of an s-wave superconductor
318 self-consistently (the attractive Hubbard model
319 :math:`-\sum_iV_in_{i\uparrow}n_{i\downarrow}` in the BCS mean-field
320 approximation, de Gennes 1966). For a spin-independent *ham* the singlet
321 BdG Hamiltonian reduces to the Nambu basis
322 :math:`(c_{i\uparrow}, c^\dagger_{i\downarrow})`,
324 .. math::
326 H = \begin{pmatrix} h-\mu & \Delta\\ \Delta^* & -(h-\mu)^*\end{pmatrix}\, ,
327 \qquad \Delta = \mathrm{diag}(\Delta_i)\, ,
329 whose eigenvectors :math:`(u_n, v_n)`, :math:`E_n > 0`, give the gap
330 equation
332 .. math::
334 \Delta_i = V_i\,\langle c_{i\downarrow}c_{i\uparrow}\rangle
335 = V_i\sum_{E_n>0}u_n(i)\,v_n^*(i)\tanh\frac{E_n}{2T}\, ,
337 iterated (with linear mixing) until self-consistent. On a clean lattice
338 it is the BCS gap equation :math:`1 = \frac{V}{N}\sum_{\mathbf{k}}
339 \tanh(E_{\mathbf{k}}/2T)/2E_{\mathbf{k}}`, with
340 :math:`\Delta(T_c) = 0` at :math:`T_c \approx \Delta(0)/1.764` at weak
341 coupling. The chemical potential is fixed (the Hartree shift
342 :math:`-V_in_i/2` is left out, i.e. absorbed in :math:`\mu`). Near
343 :math:`T_c` the iteration slows down (the linearized map has an
344 eigenvalue close to 1): raise *max_iter* there.
346 Parameters
347 ----------
348 ham : Square Hermitian matrix, (N, N), spin independent (e.g. *System.ham*)
349 V : Positive real, or array of N non-negative reals
350 Attraction.
351 mu : Real number, default 0
352 Chemical potential.
353 temperature : Positive real or zero, default 0
354 delta0 : Number, or array of N numbers, default 0.1
355 Initial gap.
356 mixing : Real in (0, 1]
357 Default value 1 (plain iteration).
358 tol : Positive real, default 1e-8
359 Convergence threshold
360 on the largest change of :math:`\Delta_i`.
361 max_iter : Positive integer, default 10000
363 Returns
364 -------
365 result : :class:`GapResult`
366 '''
367 h = ham.toarray() if hasattr(ham, 'toarray') else np.asarray(ham)
368 error_handling.square_matrix(h, 'ham')
369 error_handling.hermitian_dense(h)
370 n = len(h)
371 error_handling.interaction(V, n)
372 error_handling.real_number(mu, 'mu')
373 error_handling.positive_real_zero(temperature, 'temperature')
374 error_handling.pairing_amplitudes(delta0, n)
375 error_handling.mixing(mixing)
376 error_handling.positive_real(tol, 'tol')
377 error_handling.positive_int(max_iter, 'max_iter')
378 V = np.broadcast_to(np.asarray(V, dtype='f8'), (n,))
379 delta = np.broadcast_to(np.asarray(delta0, dtype='c16'), (n,)).copy()
380 xi = h.astype('c16') - mu * np.eye(n)
381 big = np.zeros((2 * n, 2 * n), 'c16')
382 big[:n, :n], big[n:, n:] = xi, -xi.conj()
383 diag = np.arange(n)
384 for it in range(1, max_iter + 1):
385 big[diag, diag + n] = delta
386 big[diag + n, diag] = delta.conj()
387 en, vec = np.linalg.eigh(big)
388 pos = en > 0
389 weight = np.tanh(en[pos] / (2 * temperature)) if temperature > 0 else np.ones(pos.sum())
390 new = V * np.einsum('in,in,n->i', vec[:n, pos], vec[n:, pos].conj(), weight)
391 change = np.max(np.abs(new - delta))
392 delta = (1 - mixing) * delta + mixing * new
393 if change < tol:
394 break
395 else:
396 error_handling.converged(False, 's_wave_gap')
397 return GapResult(new, en, vec, it)
400def superfluid_weight(
401 ks, delta: complex | NDArray, mu: float = 0., nk: int | tuple[int, ...] = 60,
402 temperature: float = 0., dq: float = 1e-3,
403) -> NDArray[np.float64]:
404 r'''
405 Get the superfluid weight of an s-wave superconductor,
406 :math:`D_{\mu\nu} = \partial^2\Omega/\partial A_\mu\partial A_\nu` at
407 :math:`\mathbf{A} = 0` (:math:`\hbar = e = 1`): the curvature of the
408 grand potential per unit volume in a uniform vector potential, which
409 gives the supercurrent :math:`j_\mu = -D_{\mu\nu}A_\nu` (London). *ks* is
410 the spin-up model; the spin-down electrons are its time-reversed
411 partners, :math:`h_\downarrow(\mathbf{k}) = h(-\mathbf{k})^*`, so the
412 Kramers pairs :math:`(\mathbf{k}\uparrow, -\mathbf{k}\downarrow)` pair
413 even when *ks* is a Chern band (for real hoppings, the spin-independent
414 model). In the Nambu basis
415 :math:`(c_{\mathbf{k}\uparrow}, c^\dagger_{-\mathbf{k}\downarrow})`,
417 .. math::
419 H_{\mathbf{A}}(\mathbf{k}) = \begin{pmatrix} h_{\mathbf{A}}(\mathbf{k})-\mu & \Delta \\
420 \Delta^* & -[h_{-\mathbf{A}}(\mathbf{k})-\mu]\end{pmatrix}\, ,\qquad
421 \Omega(\mathbf{A}) = \frac{1}{N_kV}\sum_{\mathbf{k}}\Big[\mathrm{Tr}\,
422 [h_{-\mathbf{A}}(\mathbf{k})-\mu] - T\sum_n\ln\big(1 + e^{-E_n/T}\big)\Big]\, ,
424 with :math:`h_{\mathbf{A}}` the Peierls Hamiltonian of
425 *get_ham_peierls* (phases on the bond vectors, orbital positions
426 included), :math:`\Delta = \mathrm{diag}(\Delta_i)` the onsite pairing,
427 :math:`E_n` the eigenvalues and :math:`V` the cell volume.
428 :math:`\Delta` is held fixed (not re-solved in the field), which is
429 exact when symmetry fixes it, e.g. uniform pairing. The derivatives are
430 central differences in :math:`\mathbf{A}`.
432 For a dispersive band it is the conventional :math:`n/m`. For an
433 isolated flat band at filling :math:`\nu`, with uniform :math:`\Delta`
434 small compared to the band gaps, there is no dispersion, and
435 :math:`D` is purely geometric (Peotta and Törmä, Nat. Commun. 6, 8944
436 (2015)),
438 .. math::
440 D_{\mu\nu} = 8\Delta\sqrt{\nu(1-\nu)}\int\frac{d^d\mathbf{k}}{(2\pi)^d}\,g_{\mu\nu}(\mathbf{k})
441 \, ,\qquad \mathrm{Tr}\,D \ge \frac{4\Delta}{\pi}\sqrt{\nu(1-\nu)}\,|C|\ \ (2D)\, ,
443 with :math:`g` the quantum metric of the band (see
444 *KSpace.quantum_geometric_tensor*) and :math:`C` its Chern number.
446 Parameters
447 ----------
448 ks : **KSpace** instance (Hermitian, orthogonal basis, spinless)
449 Spin-up model.
450 delta : Number, or array of norb numbers
451 Onsite pairing :math:`\Delta_i`.
452 mu : Real number, default 0
453 Chemical potential.
454 nk : Positive integer, or tuple of *dim* positive integers, default 60
455 k-mesh (see *mesh_grid*).
456 temperature : Positive real or zero, default 0
457 dq : Positive real, default 1e-3
458 Finite-difference step in :math:`\mathbf{A}`.
460 Returns
461 -------
462 D : Real ndarray, shape (dim, dim), in the k coordinates of *get_ham*
464 Examples
465 --------
466 Example usage::
468 # BCS chain: D = (1/2 pi) int dk 2t cos k (1 - xi/E)
469 chain = KSpace(lattices.chain())
470 chain.set_hopping([{'i': 0, 'j': 0, 'R': (1,), 't': -1.}])
471 d = superfluid_weight(chain, 0.2, mu=-0.5)
472 '''
473 from tbkit.kspace import KSpace
474 error_handling.kspace(ks, KSpace)
475 error_handling.hermitian_kspace(ks.is_hermitian(), ks._overlap_hop)
476 error_handling.spin_up_model(ks.spin)
477 error_handling.pairing_amplitudes(delta, ks.norb, 'delta')
478 error_handling.real_number(mu, 'mu')
479 error_handling.nk(nk, ks.dim)
480 error_handling.positive_real_zero(temperature, 'temperature')
481 error_handling.positive_real(dq, 'dq')
482 n, dim = ks.norb, ks.dim
483 k_cart = ks.mesh_grid(nk)[1] @ ks.k_basis.T
484 diag = np.arange(n)
485 pair = np.broadcast_to(np.asarray(delta, dtype='c16'), (n,))
486 volume = (2 * np.pi) ** dim / abs(np.linalg.det(ks.rec_vec_k))
488 def omega(q):
489 A = ks.k_basis @ q
490 hole = ks._bloch_ham(k_cart, A=-A) - mu * np.eye(n)
491 big = np.zeros((len(k_cart), 2 * n, 2 * n), 'c16')
492 big[:, :n, :n] = ks._bloch_ham(k_cart, A=A) - mu * np.eye(n)
493 big[:, n:, n:] = -hole
494 big[:, diag, diag + n] = pair
495 big[:, diag + n, diag] = pair.conj()
496 en = np.linalg.eigvalsh(big)
497 free = (-temperature * np.logaddexp(0., -en / temperature) if temperature > 0
498 else np.minimum(en, 0.))
499 return (np.trace(hole, axis1=1, axis2=2).real.sum() + free.sum()) / (len(k_cart) * volume)
501 step = dq * np.eye(dim)
502 w0 = omega(np.zeros(dim))
503 d = np.zeros((dim, dim))
504 for a in range(dim):
505 d[a, a] = (omega(step[a]) - 2 * w0 + omega(-step[a])) / dq ** 2
506 for b in range(a):
507 d[a, b] = d[b, a] = (omega(step[a] + step[b]) - omega(step[a] - step[b])
508 - omega(step[b] - step[a]) + omega(-step[a] - step[b])) / (4 * dq ** 2)
509 return d