Problem using Gutzwiller-RVB theory for high-temperature superconductivity

dieon · Jul 4, 2024

I'm currently working on reproducing a graph from the Gutzwiller-RVB theory of high-temperature superconductivity using Python. The graph plots χ_BBxy versus Δ, where Δ ranges from 16 to 26. I've implemented the theory's iterative parameter update scheme but encountered discrepancies in the plotted results compared to the expected graph.

Here's a simplified version of my Python code:

code:

import numpy as np
import matplotlib.pyplot as plt

# Constants and parameters (initial guesses)
U = 20.0
t = 1.0
a = 1.0
N_1 = 20
N_2 = 20
N = N_1 * N_2
max_iterations = 500
tolerance = 1e-5

# Initial guesses for the parameters
m_s = 0.8
delta = 0.8
chi_AB_up = 0.8
chi_AB_down = 0.8
chi_BB_up = 0.8
chi_BB_down = 0.8
chi_BBxy_up = 0.8
chi_BBxy_down = 0.8
Delta_AB = 0.8

# Generate wave vectors
def generate_wave_vectors(a, N_1, N_2):
    b1 = (2 * np.pi / a) * np.array([1, 0])
    b2 = (2 * np.pi / a) * np.array([0, 1])
    K = []
    for n_1 in range(-N_1 // 2, N_1 // 2):
        for n_2 in range(-N_2 // 2, N_2 // 2):
            k = (n_1 / N_1) * b1 + (n_2 / N_2) * b2
            K.append(k)
    return np.array(K)

K = generate_wave_vectors(a, N_1, N_2)

def gamma_k(k):
    return 2 * (np.cos(k[0]) + np.cos(k[1]))

def gamma_sc(k):
    return np.cos(k[0]) - np.cos(k[1])

def gamma_k_prime(k):
    return 2 * (np.cos(2*k[0]) + np.cos(2*k[1])) + 4 * (np.cos(k[0] + k[1]) + np.cos(k[0] - k[1]))

def h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U):
    g1 = 1
    gs = 4 / ((1 + delta)**2 - m_s**2)
    h1_up_val = (U - Delta) / 2 * (1 + delta - m_s) / 2 \
        - t**2 / Delta * (4 * (1 - 2 * delta) + g1 * (2 * chi_BB_up + 4 * chi_BBxy_up) + g1 * (1 - delta + m_s) / 2 * gamma_k_prime(k)) \
        - 2 * t**2 * gs * m_s / (U + Delta) + 2 * t**2 / (U + Delta) * (1 - delta)
    return h1_up_val

def h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U):
    g1 = 1
    gs = 4 / ((1 + delta)**2 - m_s**2)
    h1_down_val = (U - Delta) / 2 * (1 + delta + m_s) / 2 \
        - t**2 / Delta * (4 * (1 - 2 * delta) + g1 * (2 * chi_BB_down + 4 * chi_BBxy_down) + g1 * (1 - delta - m_s) / 2 * gamma_k_prime(k)) \
        + 2 * t**2 * gs * m_s / (U + Delta) + 2 * t**2 / (U + Delta) * (1 - delta)
    return h1_down_val

def h2_up(k, delta, m_s, chi_AB_up, t, Delta, U):
    g1 = 1
    g2 = 4 * delta / ((1 + delta)**2 - m_s**2)
    gs = 4 / ((1 + delta)**2 - m_s**2)
    h2_up_val = -t * g1 - t**2 / Delta * (-2 * chi_AB_up + 6 * g2 * chi_AB_up) \
        - t**2 / (U + Delta) * (gs * (0.5 * chi_AB_up + chi_AB_up) + 0.5 * chi_AB_up)
    return h2_up_val * gamma_k(k)

def h2_down(k, delta, m_s, chi_AB_down, t, Delta, U):
    g1 = 1
    g2 = 4 * delta / ((1 + delta)**2 - m_s**2)
    gs = 4 / ((1 + delta)**2 - m_s**2)
    h2_down_val = -t * g1 - t**2 / Delta * (-2 * chi_AB_down + 6 * g2 * chi_AB_down) \
        - t**2 / (U + Delta) * (gs * (0.5 * chi_AB_down + chi_AB_down) + 0.5 * chi_AB_down)
    return h2_down_val * gamma_k(k)

def h3(k, delta, m_s, Delta_AB, t, Delta, U):
    g2 = 4 * delta / ((1 + delta)**2 - m_s**2)
    gt_up = 2 * delta / (1 + delta + m_s)
    gt_down = 2 * delta / (1 + delta - m_s)
    gs = 4 / ((1 + delta)**2 - m_s**2)
    h3_val = 4 * t**2 / Delta * (1 - g2) + 4 * t**2 / (U + Delta) * (3 * gs / 4 - 1 / 4) - 2 * t**2 / Delta * (gt_up + gt_down)
    return h3_val * Delta_AB / 2 * gamma_sc(k)

def omega_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U):
    alpha_up, beta_up, _, _ = alpha_beta(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
    h1_up_k = h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U)
    h2_up_k = h2_up(k, delta, m_s, chi_AB_up, t, Delta, U)
    return h1_up_k * (alpha_up**2 - beta_up**2) - 2 * h2_up_k * alpha_up * beta_up

def omega_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U):
    _, _, alpha_down, beta_down = alpha_beta(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U)
    h1_down_k = h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U)
    h2_down_k = h2_down(k, delta, m_s, chi_AB_down, t, Delta, U)
    return h1_down_k * (alpha_down**2 - beta_down**2) - 2 * h2_down_k * alpha_down * beta_down

def nu(k, delta, m_s, Delta_AB, t, Delta, U):
    alpha_up, beta_up, alpha_down, beta_down = alpha_beta(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
    return -h3(k, delta, m_s, Delta_AB, t, Delta, U) * (alpha_up * beta_down + alpha_down * beta_up)

def alpha_beta(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U):
    h1_up_k = h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U)
    h1_down_k = h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U)
    E_up_k_val = E_up_k(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
    E_down_k_val = E_down_k(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U)

    alpha_up = np.sqrt(0.5 * (1 - h1_up_k / E_up_k_val))
    beta_up = np.sqrt(0.5 * (1 + h1_up_k / E_up_k_val))
    alpha_down = np.sqrt(0.5 * (1 - h1_down_k / E_down_k_val))
    beta_down = np.sqrt(0.5 * (1 + h1_down_k / E_down_k_val))
    return alpha_up, beta_up, alpha_down, beta_down

def E_up_k(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U):
    h1_up_k = h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U)
    h2_up_k = h2_up(k, delta, m_s, chi_AB_up, t, Delta, U)
    return np.sqrt(h1_up_k**2 + h2_up_k**2)

def E_down_k(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U):
    h1_down_k = h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U)
    h2_down_k = h2_down(k, delta, m_s, chi_AB_down, t, Delta, U)
    return np.sqrt(h1_down_k**2 + h2_down_k**2)

def u_v_squared(k, m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BBxy_up, chi_BB_down, chi_BBxy_down, t, Delta, U):
    omega_up_val = omega_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
    omega_down_val = omega_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U)
    nu_k = nu(k, delta, m_s, Delta_AB, t, Delta, U)

    denominator = np.sqrt((omega_up_val + omega_down_val)**2 + 4 * nu_k**2)
    u1_squared = 0.5 * (1 + (omega_up_val + omega_down_val) / denominator)
    u2_squared = 0.5 * (1 - (omega_up_val + omega_down_val) / denominator)
    v1_squared = u2_squared
    v2_squared = u1_squared
    return u1_squared, u2_squared, v1_squared, v2_squared


# Function to update parameters
def update_parameters(k_points, m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down, Delta_AB, t, Delta):
    sum_m_s = 0
    sum_delta = 0
    sum_chi_AB_up = 0
    sum_chi_AB_down = 0
    sum_chi_BB_up = 0
    sum_chi_BB_down = 0
    sum_chi_BBxy_up = 0
    sum_chi_BBxy_down = 0
    sum_Delta_AB = 0

    for k in k_points:
        N_1 = 20
        N_2 = 20
        N = N_1 * N_2

        alpha_up, beta_up, alpha_down, beta_down = alpha_beta(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
        h1_up_k = h1_up(k, m_s, delta, chi_BB_up, chi_BBxy_up, t, Delta, U)
        h1_down_k = h1_down(k, m_s, delta, chi_BB_down, chi_BBxy_down, t, Delta, U)
        h2_up_k = h2_up(k, delta, m_s, chi_AB_up, t, Delta, U)
        h2_down_k = h2_down(k, delta, m_s, chi_AB_down, t, Delta, U)
        E_up_k_val = E_up_k(k, m_s, delta, chi_BB_up, chi_BBxy_up, chi_AB_up, t, Delta, U)
        E_down_k_val = E_down_k(k, m_s, delta, chi_BB_down, chi_BBxy_down, chi_AB_down, t, Delta, U)

        u1_squared, u2_squared, v1_squared, v2_squared = u_v_squared(k, m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BBxy_up, chi_BB_down, chi_BBxy_down, t, Delta, U)

        u1 = np.sqrt(u1_squared)
        u2 = np.sqrt(u2_squared)
        v1 = np.sqrt(v1_squared)
        v2 = np.sqrt(v2_squared)

        # Update m_s and delta
        sum_m_s += ((alpha_up**2 - alpha_down**2) * v1_squared + (beta_up**2 - beta_down**2) * v2_squared)/N
        sum_delta += (0.5 * (alpha_up**2 * (v1_squared - v2_squared) + beta_up**2 * (v1_squared - v2_squared) + alpha_down**2 * (v1_squared - v2_squared) + beta_down**2 * (v1_squared - v2_squared)))/N

        # Update chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down
        sum_chi_AB_up += 0.25 * alpha_up * beta_up * (v2_squared - v1_squared)/N
        sum_chi_AB_down += 0.25 * alpha_down * beta_down * (v2_squared - v1_squared)/N
        sum_chi_BB_up += (np.cos(2*k[0]) + np.cos(2*k[1])) * (alpha_up**2 * v2_squared + beta_up**2 * v1_squared)/N
        sum_chi_BB_down += (np.cos(2*k[0]) + np.cos(2*k[1])) * (alpha_down**2 * v2_squared + beta_down**2 * v1_squared)/N
        sum_chi_BBxy_up += 2 * np.cos(k[0]) * np.cos(k[1]) * (alpha_up**2 * v2_squared + beta_up**2 * v1_squared)/N
        sum_chi_BBxy_down += 2 * np.cos(k[0]) * np.cos(k[1]) * (alpha_down**2 * v2_squared + beta_down**2 * v1_squared)/N
        sum_Delta_AB += (alpha_down * beta_up * u2 * v2 - alpha_up * beta_down * u1 * v1) * gamma_sc(k)

    new_m_s = sum_m_s
    new_delta = sum_delta
    new_chi_AB_up = sum_chi_AB_up
    new_chi_AB_down = sum_chi_AB_down
    new_chi_BB_up = sum_chi_BB_up
    new_chi_BB_down = sum_chi_BB_down
    new_chi_BBxy_up = sum_chi_BBxy_up
    new_chi_BBxy_down = sum_chi_BBxy_down
    new_Delta_AB = sum_Delta_AB

    return new_m_s, new_delta, new_chi_AB_up, new_chi_AB_down, new_chi_BB_up, new_chi_BB_down, new_chi_BBxy_up, new_chi_BBxy_down, new_Delta_AB

# Initializing the plot data lists
Delta_values = np.linspace(16, 26, 50)
chi_BBxy_up_values = []
chi_BBxy_down_values = []

# Loop over Delta values
for Delta in Delta_values:
    m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down, Delta_AB = 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8
    for iteration in range(max_iterations):
        new_m_s, new_delta, new_chi_AB_up, new_chi_AB_down, new_chi_BB_up, new_chi_BB_down, new_chi_BBxy_up, new_chi_BBxy_down, new_Delta_AB = update_parameters(
            K, m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down, Delta_AB, t, Delta
        )

        if np.abs(new_Delta_AB - Delta_AB) < tolerance:
            break

        m_s, delta, chi_AB_up, chi_AB_down, chi_BB_up, chi_BB_down, chi_BBxy_up, chi_BBxy_down, Delta_AB = new_m_s, new_delta, new_chi_AB_up, new_chi_AB_down, new_chi_BB_up, new_chi_BB_down, new_chi_BBxy_up, new_chi_BBxy_down, new_Delta_AB

    chi_BBxy_up_values.append(chi_BBxy_up)
    chi_BBxy_down_values.append(chi_BBxy_down)

# Plotting the results
plt.plot(Delta_values, chi_BBxy_down_values, label='χ_BBxy_down', marker='o')
plt.plot(Delta_values, chi_BBxy_up_values, label='χ_BBxy_up', marker='+')
plt.xlabel('Δ')
plt.ylabel('χ_BBxy')
plt.title('Plot of χ_BBxy vs Δ')
plt.grid(True)
plt.show()

Key Issues:

Graph Accuracy: The plotted graph of χ_BBxy against Δ does not closely match the expected results from the literature. Even small changes in initial parameters or the size of the wave vector grid (N_1 and N_2) significantly alter the output.
Parameter Sensitivity: I observe that varying initial guesses (m_s, delta, etc.) and the grid size (N_1, N_2) affect the final plotted values of χ_BBxy. Is this expected behavior, or could there be an error in how I'm updating parameters or generating wave vectors?

Request for Assistance:

Clarification: Is it typical for small changes in initial parameters or the grid size to affect the output significantly in numerical simulations of the Gutzwiller-RVB theory?
Code Review: Could you please review the provided Python code snippet to identify potential errors or improvements that might lead to more accurate results?

(I've also attached a file with the formulas I've used)

Any guidance or suggestions would be greatly appreciated!

Problem using Gutzwiller-RVB theory for high-temperature superconductivity

Attachments

Similar threads

Hot Threads

Recent Insights