二维伊辛模型Monte Carlo模拟代码报错求助：ValueError问题排查

Hey there! Let’s work through your 2D Ising Model code issues one by one—starting with that frustrating ValueError and then fixing the bugs that are stopping your plots from generating correctly.

1. Fixing the ValueError: setting an array element with a sequence #

This error is a classic name collision issue:

• You defined two functions: Energy(Q) and Magnetization(Q) to calculate system energy and magnetization.
• Later, you initialized numpy arrays with the exact same names:

  Magnetization=np.zeros(temp_points)
  Energy=np.zeros(temp_points)
  

• When you tried to call Energy[Q] and Magnetization[Q] in your simulation loop, Python thought you were trying to index these arrays with a 2D lattice (Q) instead of calling your functions. That’s why you get the "sequence as array element" error.

Fix: Rename the arrays to avoid conflict. For example, use energy_arr, magnetization_arr, specific_heat_arr, and susceptibility_arr.

2. Other Critical Bugs to Fix for Working Simulations & Plots #

Let’s go through the rest of the issues that are breaking your code or giving incorrect results:

a. Temperature Variable Mismatch #

You generated your temperature array as Temp, but tried to plot using T—this would throw a NameError. Stick to one name (I’ll use Temp in the corrected code).

b. Wrong Equilibration Step Order #

Right now, you run your data-collecting Monte Carlo steps before equilibrating the system. That means you’re collecting data from an unbalanced system, which will give garbage results. You need to run the equilibration steps first, then start collecting data.

c. Incorrect Indexing in Physical Quantity Calculation #

Lines like Energy[j]=num_a*E_a[i][j] are completely wrong: E_a is a scalar that accumulates energy values over Monte Carlo steps, not a 2D array. You should just use E_a (not E_a[i][j]) since it’s the total sum over all steps.

d. Function Call Syntax Error #

You were using array indexing syntax (Energy[Q]) instead of function call syntax (Energy(Q)). This is tied to the name collision issue above, but even after renaming, you need to fix the parentheses.

Corrected Full Code #

Here’s the fixed version with all the above issues addressed, plus some minor readability improvements:

import numpy as np
import matplotlib.pyplot as plt
import numpy.random as rnd

# Simulation parameters
num = 16
steps_monte = 1000  # Increased for better statistics (100 was too low)
temp_points = 50
steps_equil = 500   # Increased for proper equilibration
num_a = 1 / (steps_monte * (num ** 2))
num_b = 1 / ((steps_monte ** 2) * (num ** 2))

# Generate initial spin state
def startspin(num):
    spin = rnd.randint(2, size=(num, num)) - 1
    return spin

# Calculate total energy of the spin lattice
def calculate_energy(Q):
    starting_energy = 0
    for i in range(num):
        for j in range(num):
            g = Q[i, j]
            n_y = Q[(i+1)%num, j] + Q[i, (j+1)%num] + Q[(i-1)%num, j] + Q[i, (j-1)%num]
            starting_energy += g * (-n_y)
    return starting_energy / 4  # Divide by 4 to avoid double-counting bonds

# Metropolis Monte Carlo step
def monte_carlo_step(Q, beta):
    for _ in range(num * num):  # Try flipping every spin once per MC step
        x = rnd.randint(0, num)
        y = rnd.randint(0, num)
        spin = Q[x, y]
        # Calculate energy change if we flip this spin
        neighbor_sum = Q[(x+1)%num, y] + Q[x, (y+1)%num] + Q[(x-1)%num, y] + Q[x, (y-1)%num]
        delta_E = 2 * spin * neighbor_sum
        
        # Metropolis acceptance criterion
        if delta_E < 0:
            Q[x, y] *= -1
        elif rnd.rand() < np.exp(-beta * delta_E):
            Q[x, y] *= -1
    return Q

# Calculate total magnetization of the lattice
def calculate_magnetization(Q):
    return np.sum(Q)

# Generate temperature points (1 < T < 5, centered at Tc ~ 2.269)
temp_mid = 2.25
Temp = rnd.normal(temp_mid, 0.5, temp_points)
Temp = Temp[(Temp > 1) & (Temp < 5)]
temp_points = np.size(Temp)
Temp = np.sort(Temp)  # Sort for cleaner plots

# Initialize arrays to store physical quantities (renamed to avoid function conflicts)
energy_arr = np.zeros(temp_points)
magnetization_arr = np.zeros(temp_points)
specific_heat_arr = np.zeros(temp_points)
susceptibility_arr = np.zeros(temp_points)

# Run simulation for each temperature
for j in range(temp_points):
    E_total = M_total = 0
    E_sq_total = M_sq_total = 0
    Q = startspin(num)
    beta = 1 / Temp[j]
    beta_sq = beta ** 2
    
    # First: Equilibrate the system
    for _ in range(steps_equil):
        monte_carlo_step(Q, beta)
    
    # Second: Collect data over Monte Carlo steps
    for _ in range(steps_monte):
        monte_carlo_step(Q, beta)
        E = calculate_energy(Q)
        M = calculate_magnetization(Q)
        
        E_total += E
        M_total += M
        E_sq_total += E ** 2
        M_sq_total += M ** 2
    
    # Calculate averaged physical quantities
    energy_arr[j] = num_a * E_total
    magnetization_arr[j] = num_a * M_total
    # Specific heat formula: C = (1/(kT²))(<E²> - <E>²)/N (k=1 here)
    specific_heat_arr[j] = (num_a * E_sq_total - num_b * (E_total ** 2)) * beta_sq
    # Susceptibility formula: χ = (1/(kT))(<M²> - <M>²)/N (k=1 here)
    susceptibility_arr[j] = (num_a * M_sq_total - num_b * (M_total ** 2)) * beta

# Plot results
plt.figure(figsize=(10, 8))

plt.subplot(2, 2, 1)
plt.plot(Temp, energy_arr, 'd', color="#8A2BE2")
plt.xlabel("Temperature (Kelvin)")
plt.ylabel("Energy (Arbitrary Units)")
plt.title("Energy vs Temperature")

plt.subplot(2, 2, 2)
plt.plot(Temp, np.abs(magnetization_arr), 'x', color="#7FFF00")
plt.xlabel("Temperature (Kelvin)")
plt.ylabel("Magnetization (Arbitrary Units)")
plt.title("Magnetization vs Temperature")

plt.subplot(2, 2, 3)
plt.plot(Temp, specific_heat_arr, 'd', color="#00FFFF")
plt.xlabel("Temperature (Kelvin)")
plt.ylabel("Specific Heat (Arbitrary Units)")
plt.title("Specific Heat vs Temperature")

plt.subplot(2, 2, 4)
plt.plot(Temp, susceptibility_arr, 'x', color="#0000FF")
plt.xlabel("Temperature (Kelvin)")
plt.ylabel("Susceptibility (Arbitrary Units)")
plt.title("Susceptibility vs Temperature")

plt.tight_layout()
plt.show()

Key Improvements in the Corrected Code: #

• Renamed functions/arrays to avoid name collisions
• Moved equilibration steps before data collection
• Fixed function call syntax (used calculate_energy(Q) instead of Energy[Q])
• Corrected physical quantity indexing
• Sorted temperature points for cleaner plots
• Increased steps_monte and steps_equil for better statistical accuracy (100 steps was too few to get smooth curves)
• Combined all plots into a single figure for easier comparison

内容的提问来源于stack exchange，提问作者M00N KNIGHT