二维周期性边界Ising模型模拟故障：无数据写入与段错误

问题现象 #

• data.txt无数据写入
• 循环仅执行一次就终止
• 触发段错误（core dumped）
• 无法定位错误来源

预期功能 #

• 接收N（格点大小）、J（耦合常数）、T（输入温度）作为命令行参数
• 输出自旋翻转次数与最终系统能量变化
• 在1.0~7.0温度范围内绘制磁化强度-温度（M-T）曲线，曲线应呈下降趋势，在临界温度（约2.269）后趋于平缓

错误分析与修复点 #

1. 未初始化变量导致段错误
   IsingM函数中计算邻域自旋和son时，i、j是未初始化的局部变量，直接访问数组会触发非法内存访问。修复：每次选择随机格点后重新计算邻域自旋和。

2. 磁化强度计算逻辑错误
   原代码计算总磁化强度时，循环内始终使用p[i][j]而非当前循环的p[l][m]，导致结果完全错误。修复：改为累加p[l][m]。

3. 二维数组内存分配错误
   原代码中int **p = (int **)malloc(N * sizeof(int));分配的是N个int的空间，而非N个int*的空间，导致后续行指针分配越界。修复：改为malloc(N * sizeof(int*))。

4. 参数检查后未终止程序
   命令行参数个数错误时，打印用法后未终止程序，会继续执行后续代码导致非法访问。修复：添加return 1;终止程序。

5. 能量差计算逻辑错误
   邻域自旋和son仅在函数开头计算一次，未随每次选择的格点更新，完全不符合Metropolis算法逻辑。修复：每次选择随机格点后重新计算son。

6. 文件操作错误处理缺失
   打开data.txt时未检查是否成功，若文件无法打开会导致无数据写入。修复：添加文件打开成功的判断。

7. 磁化强度归一化错误
   总磁化强度应除以总格点数N*N，原代码仅除以N，导致结果数量级错误。修复：改为M /= (N*N);。

8. 随机数未初始化
   未调用srand(time(NULL))初始化随机数生成器，每次运行得到相同随机序列。修复：在main函数开头添加随机数初始化。

9. 能量计算缺失
   原代码未实现最终能量变化的输出，新增能量计算逻辑。

修复后的完整代码 #

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define MCSTEPS 100000 // Monte Carlo steps

// 计算磁化强度或翻转次数，同时可返回最终能量
double IsingM(int **p, int N, double J, double T, int c, double *final_energy) {
    int i, j, k, iter;
    double M, dE, r, f, son;
    M = 0.0;
    iter = 0;

    for (k = 0; k < MCSTEPS; k++) {
        // 随机选择一个格点
        i = rand() % N;
        j = rand() % N;

        // 计算周期性边界条件下的邻域自旋和
        son = p[(i + 1) % N][j] + p[(i - 1 + N) % N][j] +
              p[i][(j + 1) % N] + p[i][(j - 1 + N) % N];

        // 计算翻转能量差
        dE = 2 * J * p[i][j] * son; // 修正符号：dE = E_new - E_old = 2J*s*sum(s_neighbors)
        if (dE <= 0) {
            p[i][j] *= -1;
            iter++;
        } else {
            // Metropolis准则
            f = exp(-dE / T);
            r = (double)rand() / RAND_MAX;
            if (r < f) {
                p[i][j] *= -1;
                iter++;
            }
        }
    }

    // 计算最终磁化强度
    for (int l = 0; l < N; l++) {
        for (int m = 0; m < N; m++) {
            M += p[l][m];
        }
    }
    M /= (N * N); // 归一化到单位格点

    // 计算最终系统能量
    if (final_energy != NULL) {
        *final_energy = 0.0;
        for (int l = 0; l < N; l++) {
            for (int m = 0; m < N; m++) {
                son = p[(l + 1) % N][m] + p[(l - 1 + N) % N][m] +
                      p[l][(m + 1) % N] + p[l][(m - 1 + N) % N];
                *final_energy += -J * p[l][m] * son;
            }
        }
        *final_energy /= 2; // 避免重复计算相互作用
    }

    if (c == 0) {
        return M;
    } else if (c == 1) {
        return (double)iter;
    }
    return 0.0;
}

int main(int argc, char **argv) {
    // 初始化随机数生成器
    srand(time(NULL));

    if (argc != 4) {
        printf("Usage : %s <N> <J> <T>\n", argv[0]);
        return 1;
    }

    int N = atoi(argv[1]);
    double J = atof(argv[2]);
    double T_input = atof(argv[3]);

    // 分配二维数组内存
    int **p = (int **)malloc(N * sizeof(int*));
    if (!p) {
        printf("Memory allocation failed\n");
        return 1;
    }
    for (int i = 0; i < N; i++) {
        p[i] = (int *)malloc(N * sizeof(int));
        if (!p[i]) {
            printf("Memory allocation failed for row %d\n", i);
            // 释放已分配的内存
            for (int k = 0; k < i; k++) {
                free(p[k]);
            }
            free(p);
            return 1;
        }
    }

    // 初始化自旋为随机值（1或-1）
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            p[i][j] = rand() % 2 ? 1 : -1;
        }
    }

    // 生成M-T曲线数据
    FILE *fp = fopen("data.txt", "w");
    if (!fp) {
        fprintf(stderr, "Error opening data.txt for writing\n");
        // 释放内存
        for (int i = 0; i < N; i++) {
            free(p[i]);
        }
        free(p);
        return 1;
    }

    double M, t;
    for (t = 1.0; t <= 7.0; t += 0.6) {
        // 每次温度循环前重新初始化自旋，避免前序温度的影响
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                p[i][j] = rand() % 2 ? 1 : -1;
            }
        }
        M = IsingM(p, N, J, t, 0, NULL);
        fprintf(fp, "%f %f\n", t, M);
    }
    fclose(fp);

    // 计算输入温度下的翻转次数和最终能量
    double final_energy;
    int iter = (int)IsingM(p, N, J, T_input, 1, &final_energy);
    printf("N: %d, J: %f, T: %f\n", N, J, T_input);
    printf("Spin flip count: %d\n", iter);
    printf("Final system energy: %f\n", final_energy);

    // 生成Gnuplot脚本并绘制图像
    FILE *gp = fopen("plot.gp", "w");
    if (!gp) {
        printf("Error creating Gnuplot script.\n");
        // 释放内存
        for (int i = 0; i < N; i++) {
            free(p[i]);
        }
        free(p);
        return 1;
    }

    fprintf(gp, "set terminal pngcairo enhanced\n");
    fprintf(gp, "set output 'plot.png'\n");
    fprintf(gp, "set title 'Magnetization vs Temperature (2D Ising Model)'\n");
    fprintf(gp, "set xlabel 'Temperature T'\n");
    fprintf(gp, "set ylabel 'Magnetization M'\n");
    fprintf(gp, "plot 'data.txt' using 1:2 with linespoints lw 2 title 'M(T)'\n");

    fclose(gp);
    system("gnuplot plot.gp");
    printf("Plot saved as plot.png\n");

    // 释放内存
    for (int i = 0; i < N; i++) {
        free(p[i]);
    }
    free(p);

    return 0;
}

使用说明 #

1. 编译：gcc ising.c -o ising -lm（需要链接数学库）
2. 运行：./ising <N> <J> <T>，例如./ising 50 1.0 2.0
3. 输出：
   • 控制台打印输入参数、自旋翻转次数、最终系统能量
   • 生成data.txt存储M-T数据
   • 生成plot.png为M-T曲线图像

内容的提问来源于stack exchange，提问作者UNKNOWN LEGEND