- #1
dieon
- 5
- 0
I’ve been working on solving gap equations related to superconductivity using Python, in the context of RMFT using the t-J model of the Hubbard Hamiltonian. My goal is to calculate and plot the variables (\Delta_x), (\Delta_y), (\xi_x), and (\xi_y) against the doping concentration (x). However, despite several attempts, I’m encountering errors in my code.
$$\Delta_r = \frac{1}{L} \sum_k \cos(kr) \frac{\Delta_k}{E_k} $$
&
$$\xi_r = -\frac{1}{L} \sum_k \cos(kr) \frac{\xi_k}{E_k}$$
Where, ##E_k = \sqrt{(\xi_k^2 + \Delta_k^2)}##
Hole doping concentration
$$x = \frac{1}{L} \sum_k \frac{xi_k} { E_k}$$
Code Overview: I’ve written a Python script that includes a loop over a range of values for (x). Within the loop, I calculate the chemical potential ((\mu)) and update the gap parameters. The main components of my code involve wave vectors, gap equations, and iterative refinement.
Specific Questions:
Actual Output: The graph I’m currently getting doesn’t match my expectations.
$$\Delta_r = \frac{1}{L} \sum_k \cos(kr) \frac{\Delta_k}{E_k} $$
&
$$\xi_r = -\frac{1}{L} \sum_k \cos(kr) \frac{\xi_k}{E_k}$$
Where, ##E_k = \sqrt{(\xi_k^2 + \Delta_k^2)}##
Hole doping concentration
$$x = \frac{1}{L} \sum_k \frac{xi_k} { E_k}$$
Code Overview: I’ve written a Python script that includes a loop over a range of values for (x). Within the loop, I calculate the chemical potential ((\mu)) and update the gap parameters. The main components of my code involve wave vectors, gap equations, and iterative refinement.
Python:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import root
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, int(N_1 // 2)):
for n_2 in range(-N_2 // 2, int(N_2 // 2)):
k = (n_1 / N_1) * b1 + (n_2 / N_2) * b2
K.append(k)
return np.array(K)
def compute_Delta_k_xi_k(Delta_x, Delta_y, xi_x, xi_y, k, mu, J, g_s, g_t, t):
cos_kx = np.cos(k[0])
cos_ky = np.cos(k[1])
Delta_k = (3 * g_s * J / 4) * (Delta_x * cos_kx + Delta_y * cos_ky)
xi_k = -(2 * g_t * t + 3 * g_s * J * xi_x / 4) * cos_kx - (2 * g_t * t + 3 * g_s * J * xi_y / 4) * cos_ky - mu
return Delta_k, xi_k
def solve_for_mu(Delta_x, Delta_y, xi_x, xi_y, K, x, L, J, g_s, g_t, t, mu_initial_guess):
def equation(mu):
mu_sum = 0.0
for k in K:
Delta_k, xi_k = compute_Delta_k_xi_k(Delta_x, Delta_y, xi_x, xi_y, k, mu, J, g_s, g_t, t)
E_k = np.sqrt(xi_k**2 + Delta_k**2)
mu_sum += xi_k / E_k
return (mu_sum / L) - x
result = root(equation, mu_initial_guess, method='hybr', tol=1e-16, options={'maxfev': 10000})
return float(result.x[0])
def write_gap_equations(Delta_x, Delta_y, xi_x, xi_y, K, x, L, J, g_s, g_t, t, mu_initial_guess):
mu = solve_for_mu(Delta_x, Delta_y, xi_x, xi_y, K, x, L, J, g_s, g_t, t, mu_initial_guess)
Delta_x_sum = 0.0
Delta_y_sum = 0.0
xi_x_sum = 0.0
xi_y_sum = 0.0
for k in K:
Delta_k, xi_k = compute_Delta_k_xi_k(Delta_x, Delta_y, xi_x, xi_y, k, mu, J, g_s, g_t, t)
E_k = np.sqrt(xi_k**2 + Delta_k**2)
cos_kx = np.cos(k[0])
cos_ky = np.cos(k[1])
Delta_x_sum += (Delta_k / E_k) * cos_kx
Delta_y_sum += (Delta_k / E_k) * cos_ky
xi_x_sum += (xi_k / E_k) * cos_kx
xi_y_sum += (xi_k / E_k) * cos_ky
new_Delta_x = Delta_x_sum / L
new_Delta_y = Delta_y_sum / L
new_xi_x = -xi_x_sum / L
new_xi_y = -xi_y_sum / L
return new_Delta_x, new_Delta_y, new_xi_x, new_xi_y, mu
def iterate_gap_equations(Delta_x, Delta_y, xi_x, xi_y, K, x, L, J, g_s, g_t, t, mu_initial_guess, max_iter=1000, tol=1e-6):
mu = solve_for_mu(Delta_x, Delta_y, xi_x, xi_y, K, x, L, J, g_s, g_t, t, mu_initial_guess)
for iteration in range(max_iter):
new_Delta_x, new_Delta_y, new_xi_x, new_xi_y, mu = write_gap_equations(Delta_x, Delta_y, xi_x, xi_y, K, x, L, J, g_s, g_t, t, mu_initial_guess)
# Print intermediate values for diagnostics
#print(f"Iteration {iteration}: Δx={new_Delta_x:.6f}, Δy={new_Delta_y:.6f}, ξx={new_xi_x:.6f}, ξy={new_xi_y:.6f}, μ={mu_value:.6f}")
if (np.abs(new_Delta_x - Delta_x) < tol and
np.abs(new_Delta_y - Delta_y) < tol and
np.abs(new_xi_x - xi_x) < tol and
np.abs(new_xi_y - xi_y) < tol):
break
Delta_x, Delta_y, xi_x, xi_y = new_Delta_x, new_Delta_y, new_xi_x, new_xi_y
return Delta_x, Delta_y, xi_x, xi_y, mu
def points(a, N_1, N_2, t, J, g_s, g_t, initial_values, x_range):
L = N_1 * N_2
K = generate_wave_vectors(a, N_1, N_2)
x_values = []
delta_values = []
xi_values = []
mu_values = []
mu_initial_guess = 1.0 # Initial guess for the first iteration
for x in x_range:
Delta_x, Delta_y, xi_x, xi_y = initial_values
Delta_x_new, Delta_y_new, xi_x_new, xi_y_new, mu = iterate_gap_equations(Delta_x, Delta_y, xi_x, xi_y, K, x, L, J, g_s, g_t, t, mu_initial_guess)
Delta_avg = np.sqrt(Delta_x_new**2 + Delta_y_new**2)
xi_avg = (xi_x_new + xi_y_new) / 2
x_values.append(x)
delta_values.append(Delta_avg)
xi_values.append(xi_avg)
mu_values.append(mu)
# Update initial guess for mu for the next iteration
mu_initial_guess = mu
plt.figure(figsize=(10, 6))
plt.plot(x_values, delta_values, label='Average Delta', marker='o')
plt.plot(x_values, xi_values, label='Average xi', marker='+')
plt.xlabel('x')
plt.ylabel('Values')
plt.title('Plot of x vs. Average Delta and xi')
plt.legend()
plt.grid(True)
plt.show()
plt.figure(figsize=(10, 6))
plt.plot(x_values, mu_values, label='mu', marker='*')
plt.xlabel('x')
plt.ylabel('mu')
plt.title('Plot of x vs. mu')
plt.legend()
plt.grid(True)
plt.show()
# Constants
a = 1.0
N_1 = 10
N_2 = 10
t = 1.0
J = t / 3
g_s = 1.0
g_t = 1.0
initial_values = (1.0, 1.0, 0.35, 1.0)
x_range = np.linspace(0.0, 1.0, 60)
# Run the function
points(a, N_1, N_2, t, J, g_s, g_t, initial_values, x_range)Specific Questions:
- Are there any mistakes in my loops or the reassignment of values?
- What could be causing the discrepancies between the expected graph and the actual graph I’m receiving?
Actual Output: The graph I’m currently getting doesn’t match my expectations.
Last edited by a moderator: