Example snippet for 0D example

Code reproducing the Hubbard dimer from https://github.com/joselado/SolidStatePhysics2024/blob/main/notebooks/session7.ipynb

import numpy as np
import matplotlib.pyplot as plt
import meanfi

s0 = np.identity(2)
sx = np.array([[0, 1], [1, 0]])
sy = np.array([[0, -1j], [1j, 0]])
sz = np.diag([1, -1])

hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
h_0 = {(): hopp + hopp.T.conj()}

filling = 2

s_list = [sx, sy, sz]
Us = np.linspace(0, 12, 201)
svals = []

for U in Us:
    h_int = {
    (): U * np.kron(np.eye(2), sx), }
    full_model = meanfi.Model(h_0, h_int, filling)
    guess = meanfi.guess_tb(frozenset(h_int), ndof=4)
    mf_sol = meanfi.solver(full_model, guess)

    h_mf = meanfi.add_tb(h_0, mf_sol)
    rho, _ = meanfi.density_matrix(h_mf, filling=filling, nk=0)

    s_value = 0
    for s_i in s_list:
        s_operator_i = {(): np.kron(np.diag([1,0]), s_i)}
        s_value += np.abs(meanfi.expectation_value(rho, s_operator_i)) ** 2
    svals.append(s_value)

plt.plot(Us, svals)
plt.xlabel(r"$U/t$")
plt.ylabel(r"$\langle S \rangle$")