analytic.py

from math import sin, cos, sqrt, pi

import numpy as np
from scipy.optimize import root
import time
import matplotlib.pyplot as plt


def ground_energy_OBC(N, J, g):
    lamb = J/(2*g)

    def fun(k_vec):
        f_vec = np.sin(k_vec*(N+1))
        f_vec = np.divide(f_vec, np.sin(k_vec*N))
        return f_vec + lamb

    k_guess = np. arange (0.001 , np. pi ,0.01)
    k_sol = []

    k_sol = root(fun, k_guess, tol=10**(-16)).x

    k_sol_r = np.round(k_sol, 4)

    _, idx = np.unique(k_sol_r, True)
    k_sol = k_sol[idx]


    k_dist = []
    for k in k_sol:
        if 0 <= k <= pi:
            k_dist.append(k)
    k_sol = k_dist

    Nk = len(k_sol)

    if Nk != N:
        print(Nk)
        print("WARNING: could not find ground-state energy for N = ", N)
        return 0

    E0 = 0
    for k in k_dist:
        E0 += sqrt(1+lamb**2 + 2*lamb*cos(k))
    E0 *= -g/2

    return E0


def save_energies(J, g, N_max, file_name):
    E0_list = []
    with open(file_name, "w") as file:
        for N in range(2, N_max+1):
            E0 = ground_energy_OBC(N, J, g)
            if E0:
                file.write(str(N) + ", " + str(E0) + "\n")
                E0_list.append([N, E0])
    return np.array(E0_list)

def plot_energies(E0_list):
    plt.plot(E0_list[:,0], E0_list[:,1], 'ro')
    plt.xlabel("Chain length N")
    plt.ylabel("Ground state energy [J/hbar]")
    plt.show()


J = 0.5
g = 1
N = 3

print(ground_energy_OBC(N, J, g))
E0_list = save_energies(J, g, 50, "E0_analytic.txt")
plot_energies(E0_list)