Revert change of Hamiltonian phase sign cf. #105

- while doing this, also correct the direction shown for the hopping in the tight-binding model visualization to align with actual electron hopping direction

Author: thchr

ext/SymmetricTightBindingMakieExt.jl

@@ -194,19 +194,21 @@ _cubic_basis(::Val{1}) = crystal(1)
 _cubic_basis(::Val{D}) where {D} = error("Unsupported dimension: $D")
 ## --------------------------------------------------------------------------------------- #
 
-# simple case: take no account of a coefficient vector, just plot all hoppings in `h`
+# Simple case: take no account of a coefficient vector, just plot all hoppings in `h`.
+# The hopping `δ` is such that `δ = b+R-a`, so the hopping direction is from `b+R`
+# (annihilated) to `a` (created).
 function _origins_and_destinations_from_hoppingorbit(
     h::HoppingOrbit{D},
     Rm::AbstractMatrix{<:Real}
 ) where D
     P = Point{D, Float32}
-    origins = mapreduce(vcat, h.hoppings) do hs
+    destinations = mapreduce(vcat, h.hoppings) do hs
         map(hs) do h
             a = Rm * constant(h[1])
             P(a)
         end
     end
-    destinations = mapreduce(vcat, h.hoppings) do hs
+    origins = mapreduce(vcat, h.hoppings) do hs
         map(hs) do h
             b_plus_R = Rm * (constant(h[2]) + constant(h[3]))
             P(b_plus_R)
@@ -244,8 +246,8 @@ function _origins_and_destinations_from_coefficients(
             has_hop || continue
             a = P(Rm * constant(h[1]))
             b_plus_R = P(Rm * (constant(h[2]) + constant(h[3])))
-            push!(origins, a)
-            push!(destinations, b_plus_R)
+            push!(destinations, a)
+            push!(origins, b_plus_R)
         end
     end
     return origins, destinations

src/gradients.jl

@@ -232,7 +232,7 @@ function evaluate_tight_binding_momentum_gradient_term!(
     # NB: ↓ one more case of assuming no free parameters in 
     #     `δs = constant.(orbit(block.h_orbit))` (we don't materialize `δs` though)
     δs = orbit(block.h_orbit)
-    v_conj = cispi.(dot.(Ref(2 .* k), constant.(δs)))
+    v_conj = cispi.(dot.(Ref(-2 .* k), constant.(δs)))
     δ_mult_v_conj = similar(v_conj) # preallocate for reuse below
     for (idx, component) in enumerate(components)
         ∇H = ∇Hs[idx]
@@ -243,7 +243,7 @@ function evaluate_tight_binding_momentum_gradient_term!(
         for (local_i, i) in enumerate(is)
             for (local_j, j) in enumerate(js)
                 ∇Hᵢⱼ = @inbounds dot(δ_mult_v_conj, @view MmtC[:, local_i, local_j])
-                ∇Hᵢⱼ *= -2im * π # (-2πi) factor from ∂/∂kᵢ of e^{-2πik·δ}
+                ∇Hᵢⱼ *= 2im * π # (-2πi) factor from ∂/∂kᵢ of e^{-2πik·δ}
                 isnothing(c) || (∇Hᵢⱼ *= c) # multiply by coefficient if provided
                 ∇H[i, j] += ∇Hᵢⱼ
                 i == j && continue # don't add diagonal elements twice

src/tightbinding.jl

@@ -11,7 +11,9 @@ WP of `brᵦ`, displaced by a set of primitive lattice-vector representatives `R
 
 ## Implementation
 1. Take a point `a` in the WP of `brₐ` and a point `b` in the WP of `brᵦ`. 
-    Compute the displacement vector `δ = b + R - a`, where `R ∈ Rs`.
+    Compute the displacement vector `δ = b + R - a`, where `R ∈ Rs`, giving the displacement
+    vector from `a` (created site) to `b + R` (annihilated site): i.e., `δ` points opposite
+    to the hopping direction.
 2. If `δ ∈ representatives`, add `δ => (a, b, R)` to the list of hoppings 
     for that representative and continue. Otherwise, search all representatives for one 
     whose list of hoppings contains `δ => (a, b, R)`. If none is found, add `δ` as a new 
@@ -218,14 +220,14 @@ function add_reversed_orbits!(h_orbits::Vector{HoppingOrbit{D}}) where {D}
             deleteat!(h_orbits, n′)
             continue
         end
-    # we now know that at least some δ doesn't have a -δ counterpart; we assume that
-    # this implies that none of the δs have a -δ counterpart, so we need to add them all
+        # we now know that at least some δ doesn't have a -δ counterpart; we assume that
+        # this implies that none of the δs have a -δ counterpart, so we need to add them all
         @assert all(δ -> !isapproxin(-δ, δs, nothing, false; atol = VEC_CMP_ATOL), δs)
 
         # first append the new δs into the orbit
         append!(δs, -δs)
 
-    # add the "reversed" hopping terms: i.e., for every a → b + R, add b + R → a 
+        # add the "reversed" hopping terms: i.e., for every a → b + R, add b + R → a 
         # => -δ = a - (b + R) = a - b - R = b - 2b - R -a + 2a = b + (2a - 2b - R) - a = b + R' - a
         hoppings = h_orbit.hoppings
         hoppings′ = map(hoppings) do hops
@@ -238,35 +240,6 @@ function add_reversed_orbits!(h_orbits::Vector{HoppingOrbit{D}}) where {D}
 end
 
 # ---------------------------------------------------------------------------- #
-# EBRs: (q|A), (w|B)
-# Wyckoff positions: q, w
-#   q: q1, ..., qN
-#   w: w1, ..., wM
-# Site symmetry irreps: A, B
-#   A: A1, ..., AJ
-#   B: B1, ..., BK
-# δs = [δ1, δ2, ..., δn]
-#   δ1: qi₁¹ -> wj₁¹, qi₁² -> wj₁², ...
-#   δ2: qi₂¹ -> wj₂¹, qi₂² -> wj₂², ...
-# v = [exp(-ik⋅δ1), exp(-ik⋅δ2), ..., exp(-ik⋅δn)]
-# t = [[t(δ1) ...], [t(δ2) ...], ..., [t(δn) ...]]
-#   t(δ1): [t(qi₁ᵅ -> wj₁ᵅ, A_f -> B_g) ...]
-
-# Current example: (1a|E), (2c|A)
-#   ___w2__
-#  |   x   |
-#  |q1 x   x w1
-#  |_______|
-#   δs = [1/2x, -1/2x, 1/2y, -1/2y]
-#      δ1: q1 -> w1 + G1
-#      δ2: q1 -> w1 + G2
-#      δ3: q1 -> w2 + G3
-#      δ4: q1 -> w2 + G4
-# t = [t(δ1)..., t(δ2)..., t(δ3)..., t(δ4)...]
-#   t(δ1): [t(q1 -> w1, G1, E1 -> A1), t(q1 -> w1, G1, E2 -> A1)]
-#   t(δ2): [t(q1 -> w1, G2, E1 -> A1), t(q1 -> w1, G2, E2 -> A1)]
-#   t(δ3): [t(q1 -> w2, G3, E1 -> A1), t(q1 -> w2, G3, E2 -> A1)]
-#   t(δ4): [t(q1 -> w2, G4, E1 -> A1), t(q1 -> w2, G4, E2 -> A1)]
 
 # NB: The ordering of `t` is a bit subtle: `t` is a kind of vector-flattened
 #        tensor `T`, with the following "indexing convention" for `T`:
@@ -420,8 +393,6 @@ function construct_M_matrix(
     return Mm
 end
 
-# H_{s,t} = v_i M_{i,j,s,t} t_j
-
 """
     representation_constraints_matrices(
         Mm::AbstractArray{Int,4}, 
@@ -749,7 +720,7 @@ Build the P matrix for a particular symmetry operation acting on k-space, which
 rows of the M matrix.
 
 To obtain the P matrix, we exploit that the action is on exponentials of the type
-``exp(-2πi𝐤⋅δ)``, and instead act on δ ∈ `h_orbit.orbit` rather than on k. Because of this,
+``exp(2πi𝐤⋅δ)``, and instead act on δ ∈ `h_orbit.orbit` rather than on k. Because of this,
 we need to use the inverse of the rotation part of the symmetry operation.
 
 !!! details "Sketch of proof"

src/types.jl

@@ -221,7 +221,7 @@ function _getindex(
         if !iszero(δ) # print the exponential: our notation is 𝕖(δ) ≡ exp(2πi𝐤⋅δ)
             print(io, "𝕖(")
             # TODO: could do a little better here, since we know that `-δ` is in the orbit?
-            conjugate || print(io, "-")
+            conjugate && print(io, "-")
             print(io, "δ", Crystalline.subscriptify(string(n)), ")")
         else
             print(io, "1")
@@ -416,14 +416,14 @@ function evaluate_tight_binding_term!(
     MmtC = block.MmtC # contracted product of `Mm` and (complexified) `t`
 
     # NB: ↓ one more case of assuming no free parameters in `δ`
-    v_conj = cispi.(dot.(Ref(2 .* k), constant.(orbit(block.h_orbit))))
+    v_conj = cispi.(dot.(Ref(-2 .* k), constant.(orbit(block.h_orbit))))
     # NB: ↑ this is `v` conjugated: we do this because the `dot`-product below conjugates
     #     its first argument; so by conjugating twice we get the unconjugated result.
     # NB: ↑ each term in the Hamiltonian is associated to an annihilation+creation operator
     #     pair `aᵢ† aⱼ`. Since we use Convention 1 for the Fourier transform, we have
-    #     aᵢ† = ∑ₖ e^{-ik·(tᵢ + rᵢ)} aₖ†, such that each term will be multiplied by a phase
-    #     e^{-ik·δᵢⱼ} with δᵢⱼ ≡ Rᵢⱼ + rᵢ - rⱼ, i.e., the hopping vector in the orbit; this
-    #     is what `orbit(block.h_orbit)` gives us above
+    #     aᵢ† = √N⁻¹ ∑ₖ e^{-ik·(tᵢ + rᵢ)} aₖ†, such that each term will be multiplied by a
+    #     phase e^{ik·δᵢⱼ} with δᵢⱼ ≡ Rᵢⱼ + rᵢ - rⱼ, i.e., the hopping vector in the orbit;
+    #     this is what `orbit(block.h_orbit)` gives us above
     for (local_i, i) in enumerate(is)
         for (local_j, j) in enumerate(js)
             Hᵢⱼ = @inbounds dot(v_conj, @view MmtC[:, local_i, local_j])

test/berry.jl

@@ -120,10 +120,10 @@ end
             kc = cartesianize(k, dualbasis(directbasis(13, Val(2))))
 
             H = (
-                t₁ * (cos(dot(kc, a₁)) + cos(dot(kc, a₂)) + cos(dot(kc, a₃))) .* [0 1; 1 0] +
-                t₁ * (sin(dot(kc, a₁)) + sin(dot(kc, a₂)) + sin(dot(kc, a₃))) .* [0 -1im; 1im 0] +
-                2 * t₂* cos(ϕ) * (cos(dot(kc, b₁)) + cos(dot(kc, b₂)) + cos(dot(kc, b₃))) .* [1 0; 0 1] +
-                2 * t₂* sin(ϕ) * (sin(dot(kc, b₁)) + sin(dot(kc, b₂)) + sin(dot(kc, b₃))) .* [1 0; 0 -1] +
+                t₁ * (cos(kc ⋅ a₁) + cos(kc ⋅ a₂) + cos(kc ⋅ a₃)) .* [0 1; 1 0] +
+                t₁ * (sin(kc ⋅ a₁) + sin(kc ⋅ a₂) + sin(kc ⋅ a₃)) .* [0 -1im; 1im 0] +
+                2t₂ * cos(ϕ) * (cos(kc ⋅ b₁) + cos(kc ⋅ b₂) + cos(kc ⋅ b₃)) .* [1 0; 0 1] +
+                2t₂ * sin(ϕ) * (sin(kc ⋅ b₁) + sin(kc ⋅ b₂) + sin(kc ⋅ b₃)) .* [1 0; 0 -1] +
                 m .* [1 0; 0 -1]
             )
             return H
@@ -136,9 +136,7 @@ end
             H_ptbm = haldane_model(t₁, m, t₂, ϕ)(k)
             H_original = manual_haldane_model(k, t₁, m, t₂, ϕ)
 
-            # TODO: the hardcoded coefficients in `manual_haldane_model` use the old
-            #       Hamiltonian phase convention; update after correcting (cf. PR #89)
-            @test_broken isapprox(H_ptbm, H_original)
+            @test isapprox(H_ptbm, H_original)
         end
 
         # test that we also get the right Chern numbers

test/show.jl

@@ -57,8 +57,8 @@ test_tp_show(v, expected::AbstractString) = test_show(repr(MIME"text/plain"(), v
 
         str = """
         2×2 TightBindingTerm{2} over [(2b|A₁)]:
-         0                     𝕖(-δ₄)+𝕖(-δ₅)+𝕖(-δ₆)
-         𝕖(-δ₁)+𝕖(-δ₂)+𝕖(-δ₃)  0                   
+         0                  𝕖(δ₄)+𝕖(δ₅)+𝕖(δ₆)
+         𝕖(δ₁)+𝕖(δ₂)+𝕖(δ₃)  0                
         δ₁=[-1/3,1/3], δ₂=[-1/3,-2/3], δ₃=[2/3,1/3], δ₄=-δ₁, δ₅=-δ₂, δ₆=-δ₃"""
         test_tp_show(tbm[2], str)
     end
@@ -71,8 +71,8 @@ test_tp_show(v, expected::AbstractString) = test_show(repr(MIME"text/plain"(), v
         │  ⎣ 0  1 ⎦
         └─ (2b|A₁) self-term
         ┌─
-        2. ⎡ 0                     𝕖(-δ₄)+𝕖(-δ₅)+𝕖(-δ₆) ⎤
-        │  ⎣ 𝕖(-δ₁)+𝕖(-δ₂)+𝕖(-δ₃)  0                    ⎦
+        2. ⎡ 0                  𝕖(δ₄)+𝕖(δ₅)+𝕖(δ₆) ⎤
+        │  ⎣ 𝕖(δ₁)+𝕖(δ₂)+𝕖(δ₃)  0                 ⎦
         └─ (2b|A₁) self-term:  δ₁=[-1/3,1/3], δ₂=[-1/3,-2/3], δ₃=[2/3,1/3], δ₄=-δ₁, δ₅=-δ₂, δ₆=-δ₃"""
         test_tp_show(tbm, str)
     end

test/symmetry_analysis.jl

@@ -41,16 +41,16 @@ end
 # topological EBRs): i.e., we don't know if the model is a single connected band or not,
 # only what the "sum" of symmetry vectors will be.
 
-function _test_symmetry_analysis(brs, i, αβγ; rng_seed=1234)
+function _test_symmetry_analysis(brs, i, αβγ; rng = seed!(copy(Random.default_rng()), 1234))
     coefs = zeros(length(brs))
     coefs[i] = 1
     cbr = CompositeBandRep(coefs, brs)
     iszero(free(position(brs[i]))) || pin_free!(brs, i => αβγ)
     tbm = tb_hamiltonian(cbr)
-    rng = seed!(copy(Random.default_rng()), rng_seed)
     ptbm = tbm(randn(rng, length(tbm)))
     ns = collect_compatible(ptbm)
-    return !isempty(ns) && sum(ns) == SymmetryVector(cbr)
+    @test !isempty(ns)
+    @test sum(ns) == SymmetryVector(cbr)
 end
 
 @testset "Symmetry analysis (full scan over EBRs)" begin
@@ -60,7 +60,7 @@ end
             @testset "Space group $sgnum in dimension $D" begin
                 brs = calc_bandreps(sgnum, Val(D))
                 for i in eachindex(brs)
-                    @test _test_symmetry_analysis(brs, i, αβγ)
+                    _test_symmetry_analysis(brs, i, αβγ)
                 end
             end
         end