@@ -8,6 +8,7 @@ import SPG.Algebra.Group
88import SPG.Geometry.SpatialOps
99import SPG.Geometry.SpinOps
1010import SPG.Interface.Notation
11+ import SPG.Physics.Hamiltonian
1112import Mathlib.LinearAlgebra.Matrix.Determinant.Basic
1213
1314namespace Demo.Altermagnet
@@ -17,6 +18,7 @@ open SPG.Geometry.SpatialOps
1718open SPG.Geometry.SpinOps
1819open SPG.Interface
1920open SPG.Algebra
21+ open SPG.Physics.Hamiltonian
2022
2123-- 1. Generators for D4h altermagnet
2224def mat_4_z : Matrix (Fin 3 ) (Fin 3 ) ℚ := ![![0 , -1 , 0 ], ![1 , 0 , 0 ], ![0 , 0 , 1 ]]
@@ -34,136 +36,7 @@ def gen_Inv : SPGElement := Op[mat_inv, ^1]
3436
3537def Altermagnet_Group : List SPGElement := generate_group [gen_C4z_TR, gen_C2xy, gen_Inv]
3638
37- -- 2. Hamiltonian Analysis
38- -- We want to find H(k) ~ c_ij k_i k_j terms allowed by symmetry.
39- -- General form: H(k) = sum_{i,j} A_{ij} k_i k_j
40- -- Symmetry constraint: g H(k) g^{-1} = H(g k)
41- -- For spin-independent hopping (kinetic energy), H is a scalar in spin space (Identity).
42- -- But altermagnetism involves spin-dependent terms.
43- -- Altermagnet Hamiltonian usually looks like: H = (k^2 terms) * I_spin + (k-dependent field) * sigma
44- -- Let's analyze terms of form: k_a k_b * sigma_c
45-
46- -- Basis for quadratic k terms (k_x^2, k_y^2, k_z^2, k_x k_y, k_y k_z, k_z k_x)
47- inductive QuadTerm
48- | xx | yy | zz | xy | yz | zx
49- deriving Repr, DecidableEq, Inhabited
50-
51- def eval_quad (q : QuadTerm) (k : Vec3) : ℚ :=
52- match q with
53- | .xx => k 0 * k 0
54- | .yy => k 1 * k 1
55- | .zz => k 2 * k 2
56- | .xy => k 0 * k 1
57- | .yz => k 1 * k 2
58- | .zx => k 2 * k 0
59-
60- def all_quads : List QuadTerm := [.xx, .yy, .zz, .xy, .yz, .zx]
61-
62- -- Basis for spin matrices (sigma_x, sigma_y, sigma_z)
63- inductive SpinComp
64- | I | x | y | z
65- deriving Repr, DecidableEq, Inhabited
66-
67- -- Action on k-vector (spatial part)
68- def act_on_k (g : SPGElement) (k : Vec3) : Vec3 :=
69- Matrix.mulVec g.spatial k
70-
71- -- Action on spin component (spin part)
72- -- sigma_prime = U sigma U^dagger
73- -- Since our spin matrices are +/- I, the conjugation is trivial?
74- -- WAIT. The user wants "d-wave altermagnet".
75- -- In altermagnets, the spin symmetry is usually non-trivial (e.g. spin flip).
76- -- But our SPGElement definition only supports spin matrices as "numbers" (Matrix 3x3).
77- -- Currently `spin` is just +/- I.
78- -- If spin part is just +/- I (Time Reversal), then:
79- -- T sigma T^{-1} = - sigma (Time reversal flips spin)
80- -- So if g.spin = -I (Time Reversal), the spin vector flips.
81- -- If g.spin = I, spin vector stays.
82-
83- def act_on_spin (g : SPGElement) (s : SpinComp) : Vec3 :=
84- let s_vec : Vec3 := match s with
85- | .x => ![1 , 0 , 0 ]
86- | .y => ![0 , 1 , 0 ]
87- | .z => ![0 , 0 , 1 ]
88- | .I => ![0 , 0 , 0 ] -- Handle I separately
89-
90- if s == .I then ![0 , 0 , 0 ] -- I transforms to I (scalar)
91- else
92- -- Apply spatial rotation R to the axial vector sigma
93- -- sigma' = (det R) * R * sigma
94- -- If spin part has T (-I), then sigma' = - sigma'
95- let detR := Matrix.det g.spatial
96- let rotated := Matrix.mulVec g.spatial s_vec
97- let axial_rotated := detR • rotated -- Scalar multiplication
98-
99- if g.spin == spin_neg_I then
100- -axial_rotated
101- else
102- axial_rotated
103-
104- -- Helper to check if a transformed spin vector matches a target basis component (with sign)
105- -- Returns 0 if orthogonal, 1 or -1 if parallel/antiparallel
106- def project_spin (v : Vec3) (target : SpinComp) : ℚ :=
107- match target with
108- | .x => v 0
109- | .y => v 1
110- | .z => v 2
111- | .I => 0
112-
113- def check_invariant (q : QuadTerm) (s : SpinComp) : Bool :=
114- Altermagnet_Group.all fun g =>
115- -- Symmetry constraint: H(k) = U H(g^-1 k) U^dagger
116- -- H(k) = f(k) * sigma_s
117- -- Transform: f(g^-1 k) * (g sigma_s g^-1)
118- -- We check: f(g k) * (transformed sigma) == f(k) * sigma_s
119- -- Note: using g k instead of g^-1 k for group average equivalence.
120-
121- let test_ks : List Vec3 := [![1 , 0 , 0 ], ![0 , 1 , 0 ], ![0 , 0 , 1 ], ![1 , 1 , 0 ], ![1 , 0 , 1 ], ![0 , 1 , 1 ], ![1 , 2 , 3 ]]
122-
123- test_ks.all fun k =>
124- let val_gk := eval_quad q (act_on_k g k)
125- let val_k := eval_quad q k
126-
127- if s == .I then
128- -- Scalar term: val_gk == val_k
129- val_gk == val_k
130- else
131- -- Vector term: val_gk * (transformed sigma) == val_k * sigma_s
132- -- We only support cases where transformed sigma is parallel to sigma_s
133- -- (i.e. no mixing like x -> y).
134- -- Let's project transformed sigma onto s.
135- let s_prime := act_on_spin g s
136- let coeff := project_spin s_prime s
137-
138- -- Check if it stays in the same direction (possibly with sign change)
139- -- And check if orthogonal components are zero (no mixing)
140- let is_eigen :=
141- (s == .x && s_prime 1 == 0 && s_prime 2 == 0 ) ||
142- (s == .y && s_prime 0 == 0 && s_prime 2 == 0 ) ||
143- (s == .z && s_prime 0 == 0 && s_prime 1 == 0 )
144-
145- if is_eigen then
146- val_gk * coeff == val_k
147- else
148- false -- If mixing occurs, this single term is not an invariant by itself.
149-
150- def find_invariants : List (QuadTerm × SpinComp) :=
151- let terms := (all_quads.product [.I, .x, .y, .z])
152- -- Also consider symmetric combinations like (kx^2 + ky^2) and antisymmetric (kx^2 - ky^2)
153- -- Currently we only check pure basis terms.
154- -- But (kx^2 + ky^2) * I should be invariant.
155- -- Let's manually check (kx^2 + ky^2) * I and (kx^2 - ky^2) * sigma_z?
156- -- Or better, construct a linear combination solver?
157- -- For now, let's just stick to pure basis.
158- -- If kx^2 and ky^2 transform into each other, neither is invariant alone.
159- -- But their sum might be.
160- -- To properly find Hamiltonian, we need representation theory (projectors).
161- -- Simplification: Check "Is kx^2 + ky^2 invariant?" explicitly.
162-
163- let simple_invariants := terms.filter fun (q, s) => check_invariant q s
164- simple_invariants
165-
166- -- 3. ICE Symbol Output (Simplistic)
39+ -- 2. ICE Symbol Output (Simplistic)
16740-- We just print generators for now, as full ICE symbol logic is complex.
16841def ice_symbol : String :=
16942 "4'22 (D4h magnetic)" -- Placeholder, deriving from generators
@@ -173,33 +46,35 @@ end Demo.Altermagnet
17346def main : IO Unit := do
17447 IO.println s! "Generated Group Size: { Demo.Altermagnet.Altermagnet_Group.length} "
17548 IO.println s! "ICE Symbol (Approx): { Demo.Altermagnet.ice_symbol} "
176- IO.println "Allowed Quadratic Hamiltonian Terms H(k):"
49+ IO.println "Allowed Hamiltonian Terms H(k) (up to quadratic ):"
17750
178- let invariants := Demo.Altermagnet.find_invariants
179- for (q, s) in invariants do
180- let q_str := match q with
51+ let invariants := SPG.Physics.Hamiltonian.find_invariants Demo.Altermagnet.Altermagnet_Group
52+ for (p, s) in invariants do
53+ let p_str := match p with
54+ | .const => "1"
55+ | .x => "kx" | .y => "ky" | .z => "kz"
18156 | .xx => "kx^2" | .yy => "ky^2" | .zz => "kz^2"
18257 | .xy => "kx ky" | .yz => "ky kz" | .zx => "kz kx"
18358 let s_str := match s with
18459 | .I => "I" | .x => "σx" | .y => "σy" | .z => "σz"
185- IO.println s! " { q_str } * { s_str} "
60+ IO.println s! " { p_str } * { s_str} "
18661
18762 -- Manually check d-wave altermagnet term: (kx^2 - ky^2) * sigma_z?
18863 -- Or kx ky * sigma_z ?
18964 -- Let's define a custom checker for kx^2 - ky^2
190- let check_custom (f : SPG.Vec3 → ℚ) (s : Demo.Altermagnet .SpinComp) : Bool :=
65+ let check_custom (f : SPG.Vec3 → ℚ) (s : SPG.Physics.Hamiltonian .SpinComp) : Bool :=
19166 let test_ks : List SPG.Vec3 := [![1 , 0 , 0 ], ![0 , 1 , 0 ], ![0 , 0 , 1 ], ![1 , 1 , 0 ], ![1 , 0 , 1 ], ![0 , 1 , 1 ], ![1 , 2 , 3 ]]
19267 test_ks.all fun k =>
19368 -- Re-implement loop
19469 Demo.Altermagnet.Altermagnet_Group.all fun g =>
195- let val_gk := f (Demo.Altermagnet .act_on_k g k)
70+ let val_gk := f (SPG.Physics.Hamiltonian .act_on_k g k)
19671 let val_k := f k
19772 -- Re-use projection logic from check_invariant
19873 if s == .I then
19974 val_gk == val_k
20075 else
201- let s_prime := Demo.Altermagnet .act_on_spin g s
202- let coeff := Demo.Altermagnet .project_spin s_prime s
76+ let s_prime := SPG.Physics.Hamiltonian .act_on_spin g s
77+ let coeff := SPG.Physics.Hamiltonian .project_spin s_prime s
20378
20479 -- Check eigenstate property (no mixing)
20580 let is_eigen :=
@@ -218,9 +93,19 @@ def main : IO Unit := do
21893
21994 if check_custom kx2_minus_ky2 .z then
22095 IO.println " (kx^2 - ky^2) * σz [d-wave altermagnetism (x^2-y^2 type)]"
96+ else
97+ IO.println " (kx^2 - ky^2) * σz [FORBIDDEN]"
98+ -- Analyze why it is forbidden (or allowed)
99+ let _ ← SPG.Physics.Hamiltonian.analyze_term_symmetry Demo.Altermagnet.Altermagnet_Group kx2_minus_ky2 .z "(kx^2 - ky^2)" "σz"
221100
222101 if check_custom kx_ky .z then
223102 IO.println " kx ky * σz [d-wave altermagnetism (xy type)]"
103+ else
104+ IO.println " kx ky * σz [FORBIDDEN]"
105+ let _ ← SPG.Physics.Hamiltonian.analyze_term_symmetry Demo.Altermagnet.Altermagnet_Group kx_ky .z "kx ky" "σz"
224106
225107 if check_custom kx2_plus_ky2 .I then
226108 IO.println " (kx^2 + ky^2) * I [Standard kinetic term]"
109+ else
110+ IO.println " (kx^2 + ky^2) * I [FORBIDDEN]"
111+ let _ ← SPG.Physics.Hamiltonian.analyze_term_symmetry Demo.Altermagnet.Altermagnet_Group kx2_plus_ky2 .I "(kx^2 + ky^2)" "I"
0 commit comments