three_spins_chiral.py

three_spins_chiral.py

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()