from __future__ import print_function
import numpy as np
from two_spins import make_spin_ops
from three_spins import make_matrix
def does_commute(m1, m2):
"""Return True if m1 and m2 commute"""
return np.allclose(m1 @ m2, m2 @ m1)
def main():
sp_ops = make_spin_ops()
sz = sp_ops['Sz']
sp = sp_ops['S+']
sm = sp_ops['S-']
s0 = sp_ops['I']
sx = sp_ops['Sx']
sy = sp_ops['Sy']
# Hamiltonian, H
hamil = make_matrix(sz, sz, s0) + (make_matrix(sp, sm, s0) + make_matrix(sm, sp, s0)) / 2.0 \
+ make_matrix(s0, sz, sz) + (make_matrix(s0, sp, sm) + make_matrix(s0, sm, sp)) / 2.0 \
+ make_matrix(sz, s0, sz) + (make_matrix(sm, s0, sp) + make_matrix(sp, s0, sm)) / 2.0
print("H =\n", hamil)
# total spin: Sx, Sy, Sz, S^2
tot_sx = make_matrix(sx, s0, s0) + make_matrix(s0, sx, s0) + make_matrix(s0, s0, sx)
tot_sy = make_matrix(sy, s0, s0) + make_matrix(s0, sy, s0) + make_matrix(s0, s0, sy)
tot_sz = make_matrix(sz, s0, s0) + make_matrix(s0, sz, s0) + make_matrix(s0, s0, sz)
tot_s = tot_sx @ tot_sx + tot_sy @ tot_sy + tot_sz @ tot_sz
print("S^x =\n", tot_sx)
print("S^2 =\n", tot_s)
# commutation relations
print("[H, Sz] :", does_commute(hamil, tot_sz))
print("[H, S^2] :", does_commute(hamil, tot_s))
print("[S^2, Sz]:", does_commute(tot_s, tot_sz))
print("[Sz, Sx] :", does_commute(tot_sz, tot_sx))
# scalar spin chirality, C
scalar_chiral_op = make_matrix(sx, sy, sz) - make_matrix(sx, sz, sy) \
+ make_matrix(sy, sz, sx) - make_matrix(sy, sx, sz) \
+ make_matrix(sz, sx, sy) - make_matrix(sz, sy, sx)
print("\nC =\n", scalar_chiral_op)
# commutation relations
print("[H, C] :", does_commute(hamil, scalar_chiral_op))
print("[C, Sz] :", does_commute(scalar_chiral_op, tot_sz))
print("[C, S^2] :", does_commute(scalar_chiral_op, tot_s))
# diagonalize C
scalar_chirality, _ = np.linalg.eigh(scalar_chiral_op)
print("\nEigenvalues of C\n", scalar_chirality)
if __name__ == '__main__':
main()