用C语言写模拟分子热运动的代码

以下是用C语言编写的简单模拟分子热运动的代码：

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

#define NUM_PARTICLES 1000 // 粒子数量
#define BOX_SIZE 10.0 // 盒子大小
#define TIME_STEPS 1000 // 时间步长

double rand_range(double min, double max) {
    return ((double)rand() / RAND_MAX) * (max - min) + min;
}

double distance(double x1, double y1, double z1, double x2, double y2, double z2) {
    double dx = x2 - x1;
    double dy = y2 - y1;
    double dz = z2 - z1;
    return sqrt(dx*dx + dy*dy + dz*dz);
}

int main() {
    srand(time(NULL)); // 设置随机数种子
    double positions[NUM_PARTICLES][3]; // 粒子位置
    double velocities[NUM_PARTICLES][3]; // 粒子速度
    double forces[NUM_PARTICLES][3]; // 粒子受力
    double masses[NUM_PARTICLES]; // 粒子质量
    double temperature = 300.0; // 温度
    double dt = 0.001; // 时间步长
    double kb = 1.38064852e-23; // 玻尔兹曼常数

    // 初始化粒子
    for (int i = 0; i < NUM_PARTICLES; i++) {
        positions[i][0] = rand_range(0.0, BOX_SIZE);
        positions[i][1] = rand_range(0.0, BOX_SIZE);
        positions[i][2] = rand_range(0.0, BOX_SIZE);
        velocities[i][0] = rand_range(-1.0, 1.0);
        velocities[i][1] = rand_range(-1.0, 1.0);
        velocities[i][2] = rand_range(-1.0, 1.0);
        masses[i] = rand_range(1.0, 10.0);
    }

    // 模拟运动
    for (int t = 0; t < TIME_STEPS; t++) {
        // 计算受力
        for (int i = 0; i < NUM_PARTICLES; i++) {
            forces[i][0] = 0.0;
            forces[i][1] = 0.0;
            forces[i][2] = 0.0;
            for (int j = 0; j < NUM_PARTICLES; j++) {
                if (i != j) {
                    double d = distance(positions[i][0], positions[i][1], positions[i][2],
                                        positions[j][0], positions[j][1], positions[j][2]);
                    double f = (masses[i] * masses[j]) / (d*d);
                    forces[i][0] += f * (positions[j][0] - positions[i][0]) / d;
                    forces[i][1] += f * (positions[j][1] - positions[i][1]) / d;
                    forces[i][2] += f * (positions[j][2] - positions[i][2]) / d;
                }
            }
        }

        // 更新位置和速度
        for (int i = 0; i < NUM_PARTICLES; i++) {
            velocities[i][0] += forces[i][0] * dt / masses[i];
            velocities[i][1] += forces[i][1] * dt / masses[i];
            velocities[i][2] += forces[i][2] * dt / masses[i];
            positions[i][0] += velocities[i][0] * dt;
            positions[i][1] += velocities[i][1] * dt;
            positions[i][2] += velocities[i][2] * dt;
        }

        // 控制温度
        double kinetic_energy = 0.0;
        for (int i = 0; i < NUM_PARTICLES; i++) {
            kinetic_energy += 0.5 * masses[i] * (velocities[i][0]*velocities[i][0] + 
                                                 velocities[i][1]*velocities[i][1] + 
                                                 velocities[i][2]*velocities[i][2]);
        }
        double current_temperature = (2.0 / 3.0) * kinetic_energy / (kb * NUM_PARTICLES);
        double temperature_scale = sqrt(temperature / current_temperature);
        for (int i = 0; i < NUM_PARTICLES; i++) {
            velocities[i][0] *= temperature_scale;
            velocities[i][1] *= temperature_scale;
            velocities[i][2] *= temperature_scale;
        }

        // 输出粒子位置
        for (int i = 0; i < NUM_PARTICLES; i++) {
            printf("%f,%f,%f\n", positions[i][0], positions[i][1], positions[i][2]);
        }
    }

    return 0;
}

该代码使用随机数生成粒子的位置、速度和质量，并通过计算受力和更新位置和速度来模拟粒子的运动。此外，该代码还通过控制温度来使粒子保持一定的热运动状态。最后，该代码输出粒子的位置信息，可以用于可视化粒子的运动轨迹。