two_spins.py

two_spins.py

from __future__ import print_function

import numpy as np


def make_spin_ops():
    ops = {}

    # S^z
    ops['Sz'] = np.diag((0.5, -0.5))

    # S^+
    sp = np.zeros((2,2))
    sp[0, 1] = 1
    ops['S+'] = sp

    # S^-
    sm = sp.transpose()
    ops['S-'] = sm

    # Sx, Sy
    ops['Sx'] = (sp + sm) / 2.0
    ops['Sy'] = (sp - sm) / 2.0j

    # S^0
    ops['I'] = np.identity(2)

    return ops


def make_matrix(s1, s2):
    return np.kron(s1, s2)


def main():
    sp_ops = make_spin_ops()
    sz = sp_ops['Sz']
    sp = sp_ops['S+']
    sm = sp_ops['S-']
    s0 = sp_ops['I']

    hamil = make_matrix(sz, sz) + (make_matrix(sp, sm) + make_matrix(sm, sp)) / 2.0
    print("H =\n", hamil)

    eigval, eigvec = np.linalg.eigh(hamil)

    print("\nEigenvalues =\n", eigval)
    print("\nEigenvectors =")
    for i in range(4):
        print(eigvec[:, i])


if __name__ == '__main__':
    main()