用c++写一个简单的分子动力学程序,并且写上注释
时间: 2023-03-29 13:04:29 浏览: 94
好的,以下是一个简单的分子动力学程序,用 C 语言编写,注释已经加上:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 100 // 分子数
#define L 10. // 盒子边长
#define T 1. // 温度
#define dt .01 // 时间步长
#define steps 100 // 总步数
double x[N], y[N], z[N]; // 分子的位置
double vx[N], vy[N], vz[N]; // 分子的速度
double fx[N], fy[N], fz[N]; // 分子的受力
double rand_normal(double mean, double stddev) {
double u1 = rand() / (double)RAND_MAX;
double u2 = rand() / (double)RAND_MAX;
double z = sqrt(-2. * log(u1)) * cos(2. * M_PI * u2);
return mean + stddev * z;
}
int main() {
int i, j, k;
double r2, r6, f, e, v, ke, pe;
// 初始化分子的位置和速度
for (i = ; i < N; i++) {
x[i] = rand() / (double)RAND_MAX * L;
y[i] = rand() / (double)RAND_MAX * L;
z[i] = rand() / (double)RAND_MAX * L;
vx[i] = rand_normal(., sqrt(T));
vy[i] = rand_normal(., sqrt(T));
vz[i] = rand_normal(., sqrt(T));
}
// 开始模拟
for (k = ; k < steps; k++) {
// 计算分子之间的相互作用力
for (i = ; i < N; i++) {
fx[i] = fy[i] = fz[i] = .;
for (j = ; j < N; j++) {
if (i != j) {
r2 = (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]) + (z[i] - z[j]) * (z[i] - z[j]);
r6 = r2 * r2 * r2;
f = 24. * (2. / r6 - 1.) / r6 / r2;
fx[i] += f * (x[i] - x[j]);
fy[i] += f * (y[i] - y[j]);
fz[i] += f * (z[i] - z[j]);
}
}
}
// 计算分子的速度和位置
ke = .;
for (i = ; i < N; i++) {
vx[i] += fx[i] * dt;
vy[i] += fy[i] * dt;
vz[i] += fz[i] * dt;
x[i] += vx[i] * dt;
y[i] += vy[i] * dt;
z[i] += vz[i] * dt;
if (x[i] < .) x[i] += L;
if (x[i] > L) x[i] -= L;
if (y[i] < .) y[i] += L;
if (y[i] > L) y[i] -= L;
if (z[i] < .) z[i] += L;
if (z[i] > L) z[i] -= L;
ke += .5 * (vx[i] * vx[i] + vy[i] * vy[i] + vz[i] * vz[i]);
}
// 计算分子之间的势能
pe = .;
for (i = ; i < N; i++) {
for (j = i + 1; j < N; j++) {
r2 = (x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]) + (z[i] - z[j]) * (z[i] - z[j]);
r6 = r2 * r2 * r2;
e = 4. * (1. / r6 - 1.) / r6;
pe += e;
}
}
// 计算总能量和温度
v = ke + pe;
T = 2. * ke / (3. * N);
printf("%d %lf %lf %lf\n", k, v, ke, pe);
}
return ;
}