/** * engines/c/main.c * ----------------- * C 语言动力学模拟引擎。 * 与 Python 版 (compute.py) 算法保持一致。 * * 输入: param.json 数值参数（Python 从 YAML 转换得来） * /coord.txt * /connection.txt * /bond.txt * 输出: /trajectory.txt (JSON 格式，与 Python 版兼容) * * 编译: cmake --build build --target dynamics_c * 用法: ./build/dynamics_c */ #include #include #include #include #include /* ======================================================================== * 配置参数（从 param.json 读取） * ======================================================================== */ typedef struct { double box_a; /* 盒子半边长 */ int NT; /* 总步数 */ double DT; /* 时间步长 */ int NSTEP; /* 抽帧间隔 */ int warmup_steps; /* 预热步数 */ char method[32]; /* 算法名称 */ double G[3]; /* 重力分量 */ double B[3]; /* 阻尼分量 */ int gravity_field; /* 均匀重力场开关 */ int gravity_interaction; /* 原子间万有引力开关 */ int elastic_force; /* 弹簧键力开关 */ int damping_force; /* 阻尼开关 */ double gravity_strength; /* 万有引力强度 */ } SimParams; /* ======================================================================== * 原子数据 * ======================================================================== */ typedef struct { int n_atoms; int *atom_ids; double *masses; double *radii; double *pos_0; /* 初始位置 (n_atoms*3) */ double *vel_0; /* 初始速度 (n_atoms*3) */ int *fixed; /* 固定约束 (n_atoms*3) */ } AtomData; /* ======================================================================== * 成键数据 * ======================================================================== */ typedef struct { int n_bonds; int *pairs; /* (n_bonds*2) */ double *stiffness; double *rest_lengths; } BondData; /* ======================================================================== * 轨迹缓冲区 * ======================================================================== */ typedef struct { int n_steps; int n_atoms; double *x, *y, *z; double *vx, *vy, *vz; } Trajectory; /* ======================================================================== * 辅助函数 * ======================================================================== */ static void die(const char *msg) { fprintf(stderr, "[C-engine] 错误: %s\n", msg); exit(1); } static void *xmalloc(size_t sz) { void *p = malloc(sz); if (!p) die("内存分配失败"); return p; } /* 从 JSON 中读取 double 值 */ static double json_read_double(const char *json, const char *key) { char search[256]; snprintf(search, sizeof(search), "\"%s\"", key); const char *p = strstr(json, search); if (!p) return 0.0; p = strchr(p, ':'); if (!p) return 0.0; p++; while (*p == ' ' || *p == '\t' || *p == '\n') p++; return strtod(p, NULL); } static int json_read_int(const char *json, const char *key) { return (int)json_read_double(json, key); } /* 从 JSON 中读取字符串值（写入 dst，最多 dst_sz 字节） */ static void json_read_string(const char *json, const char *key, char *dst, int dst_sz) { char search[256]; snprintf(search, sizeof(search), "\"%s\"", key); const char *p = strstr(json, search); if (!p) { dst[0] = '\0'; return; } p = strchr(p, ':'); if (!p) { dst[0] = '\0'; return; } p++; while (*p == ' ' || *p == '\t' || *p == '\n') p++; if (*p != '"') { dst[0] = '\0'; return; } p++; int i = 0; while (*p && *p != '"' && i < dst_sz - 1) { dst[i++] = *p++; } dst[i] = '\0'; } /* 读取 JSON 数组 (如 "G": [0, 0, -9.8]) 到 double[3] */ static void json_read_double3(const char *json, const char *key, double out[3]) { char search[256]; snprintf(search, sizeof(search), "\"%s\"", key); const char *p = strstr(json, search); if (!p) { out[0]=out[1]=out[2]=0; return; } p = strchr(p, '['); if (!p) { out[0]=out[1]=out[2]=0; return; } p++; for (int i = 0; i < 3; i++) { while (*p == ' ' || *p == '\t' || *p == '\n' || *p == ',') p++; out[i] = strtod(p, (char**)&p); } } /* 读取 param.json */ static int g_gravity_field = 1; static int g_gravity_interaction = 0; static int g_elastic_force = 1; static int g_damping_force = 0; static double g_gravity_strength = 1.0; static SimParams read_params(const char *path) { FILE *f = fopen(path, "rb"); if (!f) die("无法打开 param.json"); fseek(f, 0, SEEK_END); long sz = ftell(f); fseek(f, 0, SEEK_SET); char *buf = (char*)xmalloc(sz + 1); fread(buf, 1, sz, f); buf[sz] = '\0'; fclose(f); SimParams p; p.box_a = json_read_double(buf, "box_a"); p.NT = json_read_int(buf, "NT"); p.DT = json_read_double(buf, "DT"); p.NSTEP = json_read_int(buf, "NSTEP"); p.warmup_steps = json_read_int(buf, "warmup_steps"); strcpy(p.method, "leapfrog"); /* 默认 */ json_read_string(buf, "method", p.method, sizeof(p.method)); json_read_double3(buf, "G", p.G); json_read_double3(buf, "B", p.B); p.gravity_field = json_read_int(buf, "gravity_field"); p.gravity_interaction = json_read_int(buf, "gravity_interaction"); p.elastic_force = json_read_int(buf, "elastic_force"); p.damping_force = json_read_int(buf, "damping_force"); p.gravity_strength = json_read_double(buf, "gravity_strength"); g_gravity_field = p.gravity_field; g_gravity_interaction = p.gravity_interaction; g_elastic_force = p.elastic_force; g_damping_force = p.damping_force; g_gravity_strength = p.gravity_strength; free(buf); return p; } /* 读取 coord.txt */ static AtomData read_coord(const char *input_dir) { char path[512]; snprintf(path, sizeof(path), "%s/coord.txt", input_dir); FILE *f = fopen(path, "r"); if (!f) die("无法打开 coord.txt"); /* 跳过第一行表头 */ char line[1024]; if (!fgets(line, sizeof(line), f)) die("coord.txt 为空"); int capacity = 16; AtomData a; a.n_atoms = 0; a.atom_ids = (int*)xmalloc(capacity * sizeof(int)); a.masses = (double*)xmalloc(capacity * sizeof(double)); a.radii = (double*)xmalloc(capacity * sizeof(double)); a.pos_0 = (double*)xmalloc(capacity * 3 * sizeof(double)); a.vel_0 = (double*)xmalloc(capacity * 3 * sizeof(double)); a.fixed = (int*)xmalloc(capacity * 3 * sizeof(int)); while (fgets(line, sizeof(line), f)) { if (a.n_atoms >= capacity) { capacity *= 2; a.atom_ids = realloc(a.atom_ids, capacity * sizeof(int)); a.masses = realloc(a.masses, capacity * sizeof(double)); a.radii = realloc(a.radii, capacity * sizeof(double)); a.pos_0 = realloc(a.pos_0, capacity * 3 * sizeof(double)); a.vel_0 = realloc(a.vel_0, capacity * 3 * sizeof(double)); a.fixed = realloc(a.fixed, capacity * 3 * sizeof(int)); } int id, fx, fy, fz; double mass, rad, px, py, pz, vx, vy, vz; int n_parsed = sscanf(line, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %d", &id, &mass, &rad, &px, &py, &pz, &vx, &vy, &vz, &fx, &fy, &fz); if (n_parsed == 9) { fx = fy = fz = 0; } else if (n_parsed != 12) { continue; } int i = a.n_atoms; a.atom_ids[i] = id; a.masses[i] = mass; a.radii[i] = rad; a.pos_0[i*3+0] = px; a.pos_0[i*3+1] = py; a.pos_0[i*3+2] = pz; a.vel_0[i*3+0] = vx; a.vel_0[i*3+1] = vy; a.vel_0[i*3+2] = vz; a.fixed[i*3+0] = fx; a.fixed[i*3+1] = fy; a.fixed[i*3+2] = fz; a.n_atoms++; } fclose(f); if (a.n_atoms <= 0) die("coord.txt 原子数无效"); return a; } /* 读取 connection.txt */ static BondData read_bonds(const char *input_dir, const AtomData *atoms) { char path[512]; BondData b; b.n_bonds = 0; b.pairs = NULL; b.stiffness = NULL; b.rest_lengths = NULL; snprintf(path, sizeof(path), "%s/connection.txt", input_dir); FILE *f = fopen(path, "r"); if (!f) return b; char line[256]; if (!fgets(line, sizeof(line), f)) { fclose(f); return b; } int n_lines = 0, tmp_a, tmp_b; char bond_name[256]; while (fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name) == 3) n_lines++; rewind(f); if (n_lines == 0) { fclose(f); return b; } b.n_bonds = n_lines; b.pairs = (int*)xmalloc(n_lines * 2 * sizeof(int)); b.stiffness = (double*)xmalloc(n_lines * sizeof(double)); b.rest_lengths = (double*)xmalloc(n_lines * sizeof(double)); char bond_path[512]; snprintf(bond_path, sizeof(bond_path), "%s/bond.txt", input_dir); FILE *fb = fopen(bond_path, "r"); for (int i = 0; i < n_lines; i++) { fscanf(f, "%d %d %s", &tmp_a, &tmp_b, bond_name); b.pairs[i*2+0] = tmp_a - 1; b.pairs[i*2+1] = tmp_b - 1; b.stiffness[i] = 1.0; b.rest_lengths[i] = 2.0; if (fb) { char name[256], header[256]; double k, r0; rewind(fb); fgets(header, sizeof(header), fb); // 跳过表头行 while (fscanf(fb, "%s %lf %lf", name, &k, &r0) == 3) { if (strcmp(name, bond_name) == 0) { b.stiffness[i] = k; b.rest_lengths[i] = r0; break; } } } } fclose(f); if (fb) fclose(fb); return b; } /* ======================================================================== * 物理核心（与 Python compute.py 对应） * ======================================================================== */ /* 加速度计算（各力独立开关控制） */ static void compute_acceleration( int n, const double *x, const double *y, const double *z, const double *vx, const double *vy, const double *vz, const double *m, const double G[3], const double B[3], const BondData *bonds, double *ax, double *ay, double *az) { /* 先清零 */ for (int i = 0; i < n; i++) { ax[i] = 0.0; ay[i] = 0.0; az[i] = 0.0; } /* 均匀重力场 */ if (g_gravity_field) { for (int i = 0; i < n; i++) { ax[i] += G[0]; ay[i] += G[1]; az[i] += G[2]; } } /* 阻尼 */ if (g_damping_force) { for (int i = 0; i < n; i++) { ax[i] -= B[0] * vx[i] / m[i]; ay[i] -= B[1] * vy[i] / m[i]; az[i] -= B[2] * vz[i] / m[i]; } } /* 弹簧键力 */ if (g_elastic_force) { for (int b = 0; b < bonds->n_bonds; b++) { int i = bonds->pairs[b*2+0]; int j = bonds->pairs[b*2+1]; double dx = x[j] - x[i]; double dy = y[j] - y[i]; double dz = z[j] - z[i]; double dist = sqrt(dx*dx + dy*dy + dz*dz); if (dist < 1e-12) continue; double stretch = dist - bonds->rest_lengths[b]; double fmag = bonds->stiffness[b] * stretch; double ux = dx / dist, uy = dy / dist, uz = dz / dist; double fx = fmag * ux, fy = fmag * uy, fz = fmag * uz; ax[i] += fx / m[i]; ay[i] += fy / m[i]; az[i] += fz / m[i]; ax[j] -= fx / m[j]; ay[j] -= fy / m[j]; az[j] -= fz / m[j]; } } /* 万有引力（所有原子对之间） */ if (g_gravity_interaction) { for (int i = 0; i < n; i++) { for (int j = i + 1; j < n; j++) { double dx = x[j] - x[i]; double dy = y[j] - y[i]; double dz = z[j] - z[i]; double r2 = dx*dx + dy*dy + dz*dz; if (r2 <= 1e-12) continue; double r = sqrt(r2); double f_mag = g_gravity_strength * m[i] * m[j] / r2; double fx_g = f_mag * dx / r; double fy_g = f_mag * dy / r; double fz_g = f_mag * dz / r; ax[i] += fx_g / m[i]; ay[i] += fy_g / m[i]; az[i] += fz_g / m[i]; ax[j] -= fx_g / m[j]; ay[j] -= fy_g / m[j]; az[j] -= fz_g / m[j]; } } } } /* 边界条件：clamp 位置 + 速度反转 ——与 Python Limit_in_box 一致 */ static void limit_in_box(double *pos, double *vel, double lo, double hi) { if (*pos > hi) { *pos = hi; *vel = -*vel; } if (*pos < lo) { *pos = lo; *vel = -*vel; } } /* ── 显式欧拉法 ──────────── */ static void explicit_euler_step( int n, double *x, double *y, double *z, double *vx, double *vy, double *vz, const double *m, const double G[3], const double B[3], const BondData *bonds, const int *fixed, double dt) { double *ax = (double*)alloca(n * sizeof(double)); double *ay = (double*)alloca(n * sizeof(double)); double *az = (double*)alloca(n * sizeof(double)); compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az); for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; x[i] += vx[i] * dt; y[i] += vy[i] * dt; z[i] += vz[i] * dt; vx[i] += ax[i] * dt; vy[i] += ay[i] * dt; vz[i] += az[i] * dt; } } /* ── 隐式欧拉法 ──────────── */ static void implicit_euler_step( int n, double *x, double *y, double *z, double *vx, double *vy, double *vz, const double *m, const double G[3], const double B[3], const BondData *bonds, const int *fixed, double dt) { double *vxn = (double*)alloca(n * sizeof(double)); double *vyn = (double*)alloca(n * sizeof(double)); double *vzn = (double*)alloca(n * sizeof(double)); for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) { vxn[i] = 0; vyn[i] = 0; vzn[i] = 0; continue; } double gx = B[0] / m[i], gy = B[1] / m[i], gz = B[2] / m[i]; vxn[i] = (vx[i] + G[0] * dt) / (1.0 + gx * dt); vyn[i] = (vy[i] + G[1] * dt) / (1.0 + gy * dt); vzn[i] = (vz[i] + G[2] * dt) / (1.0 + gz * dt); } double *ax = (double*)alloca(n * sizeof(double)); double *ay = (double*)alloca(n * sizeof(double)); double *az = (double*)alloca(n * sizeof(double)); compute_acceleration(n, x, y, z, vxn, vyn, vzn, m, G, B, bonds, ax, ay, az); for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; vx[i] += ax[i] * dt; vy[i] += ay[i] * dt; vz[i] += az[i] * dt; x[i] += vx[i] * dt; y[i] += vy[i] * dt; z[i] += vz[i] * dt; } } /* ── 中点法 ──────────── */ static void midpoint_step( int n, double *x, double *y, double *z, double *vx, double *vy, double *vz, const double *m, const double G[3], const double B[3], const BondData *bonds, const int *fixed, double dt) { double *ax = (double*)alloca(n * sizeof(double)); double *ay = (double*)alloca(n * sizeof(double)); double *az = (double*)alloca(n * sizeof(double)); compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az); double *xm = (double*)alloca(n * sizeof(double)); double *ym = (double*)alloca(n * sizeof(double)); double *zm = (double*)alloca(n * sizeof(double)); double *vxm = (double*)alloca(n * sizeof(double)); double *vym = (double*)alloca(n * sizeof(double)); double *vzm = (double*)alloca(n * sizeof(double)); for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) { xm[i]=ym[i]=zm[i]=vxm[i]=vym[i]=vzm[i]=0; continue; } xm[i] = x[i] + 0.5 * vx[i] * dt; ym[i] = y[i] + 0.5 * vy[i] * dt; zm[i] = z[i] + 0.5 * vz[i] * dt; vxm[i] = vx[i] + 0.5 * ax[i] * dt; vym[i] = vy[i] + 0.5 * ay[i] * dt; vzm[i] = vz[i] + 0.5 * az[i] * dt; x[i] = x[i] + vxm[i] * dt; y[i] = y[i] + vym[i] * dt; z[i] = z[i] + vzm[i] * dt; } compute_acceleration(n, xm, ym, zm, vxm, vym, vzm, m, G, B, bonds, ax, ay, az); for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; vx[i] += ax[i] * dt; vy[i] += ay[i] * dt; vz[i] += az[i] * dt; } } /* ── 蛙跳法（Velocity-Verlet）── */ static void leapfrog_step( int n, double *x, double *y, double *z, double *vx, double *vy, double *vz, const double *m, const double G[3], const double B[3], const BondData *bonds, const int *fixed, double dt) { double *ax = (double*)alloca(n * sizeof(double)); double *ay = (double*)alloca(n * sizeof(double)); double *az = (double*)alloca(n * sizeof(double)); compute_acceleration(n, x, y, z, vx, vy, vz, m, G, B, bonds, ax, ay, az); /* 半推速度：v_half = v + 0.5*a*dt */ for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; vx[i] += ax[i] * dt * 0.5; vy[i] += ay[i] * dt * 0.5; vz[i] += az[i] * dt * 0.5; } /* 全推位置（不含边界）*/ for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; x[i] += vx[i] * dt; /* vx 此时是 v_half */ y[i] += vy[i] * dt; z[i] += vz[i] * dt; } /* 显式预测器：v_pred = v_half + 0.5*a_old*dt，用第一次加速度外推半步 包含重力+阻尼+弹簧的所有贡献（标准 Velocity-Verlet 预测步）*/ double *pred_vx = (double*)alloca(n * sizeof(double)); double *pred_vy = (double*)alloca(n * sizeof(double)); double *pred_vz = (double*)alloca(n * sizeof(double)); for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; pred_vx[i] = vx[i] + 0.5 * ax[i] * dt; pred_vy[i] = vy[i] + 0.5 * ay[i] * dt; pred_vz[i] = vz[i] + 0.5 * az[i] * dt; } /* 用新位置 + 预测速度重算加速度 */ compute_acceleration(n, x, y, z, pred_vx, pred_vy, pred_vz, m, G, B, bonds, ax, ay, az); /* 速度后半步：v = v_half + 0.5*a_next*dt vx 仍为 v_half（未被覆盖）*/ for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; vx[i] += ax[i] * dt * 0.5; vy[i] += ay[i] * dt * 0.5; vz[i] += az[i] * dt * 0.5; } } /* ── 分发器：调用对应积分方法 + 边界条件 + 自由度约束（与 Python 一致）── */ static void apply_step( const char *method, int n, double *x, double *y, double *z, double *vx, double *vy, double *vz, const double *m, const double G[3], const double B[3], const BondData *bonds, const int *fixed, const double *pos_0, double box_a, double dt) { if (strcmp(method, "explicit_euler") == 0) { explicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt); } else if (strcmp(method, "implicit_euler") == 0) { implicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt); } else if (strcmp(method, "midpoint") == 0) { midpoint_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt); } else if (strcmp(method, "leapfrog") == 0) { leapfrog_step(n, x, y, z, vx, vy, vz, m, G, B, bonds, fixed, dt); } else { fprintf(stderr, "[C-engine] 未知算法: %s\n", method); exit(1); } /* 边界条件（与 Python Limit_in_box 一致） */ for (int i = 0; i < n; i++) { if (fixed[i*3+0] && fixed[i*3+1] && fixed[i*3+2]) continue; limit_in_box(&x[i], &vx[i], -box_a, box_a); limit_in_box(&y[i], &vy[i], -box_a, box_a); limit_in_box(&z[i], &vz[i], -box_a, box_a); } /* 逐自由度固定约束（与 Python apply_fixed_constraints 一致） */ for (int i = 0; i < n; i++) { if (fixed[i*3+0]) { x[i] = pos_0[i*3+0]; vx[i] = 0.0; } if (fixed[i*3+1]) { y[i] = pos_0[i*3+1]; vy[i] = 0.0; } if (fixed[i*3+2]) { z[i] = pos_0[i*3+2]; vz[i] = 0.0; } } } // ======================================================================== // JSON 输出 // ======================================================================== static void write_trajectory_json(const char *path, const Trajectory *traj, const SimParams *params, const AtomData *atoms, const BondData *bonds) { FILE *f = fopen(path, "w"); if (!f) die("无法写入 trajectory.txt"); fprintf(f, "{\n"); const char *names[] = {"traj_x","traj_y","traj_z","traj_vx","traj_vy","traj_vz"}; double *arrs[] = {traj->x, traj->y, traj->z, traj->vx, traj->vy, traj->vz}; for (int a = 0; a < 6; a++) { fprintf(f, " \"%s\": [\n", names[a]); for (int t = 0; t < traj->n_steps; t++) { fprintf(f, " ["); for (int i = 0; i < traj->n_atoms; i++) { fprintf(f, "%.15g", arrs[a][t * traj->n_atoms + i]); if (i < traj->n_atoms - 1) fputc(',', f); } fprintf(f, "]"); if (t < traj->n_steps - 1) fputc(',', f); fputc('\n', f); } fprintf(f, " ]"); fputc(',', f); fputc('\n', f); } /* 标量参数 */ fprintf(f, " \"NT\": %d,\n", params->NT); fprintf(f, " \"DT\": %.15g,\n", params->DT); fprintf(f, " \"NSTEP\": %d,\n", params->NSTEP); fprintf(f, " \"method\": \"%s\",\n", params->method); fprintf(f, " \"warmup_steps\": %d,\n", params->warmup_steps); fprintf(f, " \"G\": [%.15g, %.15g, %.15g],\n", params->G[0], params->G[1], params->G[2]); fprintf(f, " \"B\": [%.15g, %.15g, %.15g],\n", params->B[0], params->B[1], params->B[2]); fprintf(f, " \"atom_ids\": ["); for (int i = 0; i < atoms->n_atoms; i++) { if (i > 0) fputc(',', f); fprintf(f, "%d", atoms->atom_ids[i]); } fprintf(f, "],\n"); fprintf(f, " \"atom_masses\": ["); for (int i = 0; i < atoms->n_atoms; i++) { if (i > 0) fputc(',', f); fprintf(f, "%.15g", atoms->masses[i]); } fprintf(f, "],\n"); fprintf(f, " \"bond_pairs\": ["); for (int b = 0; b < bonds->n_bonds; b++) { if (b > 0) fputc(',', f); fprintf(f, "[%d, %d]", bonds->pairs[b*2], bonds->pairs[b*2+1]); } fprintf(f, "],\n"); fprintf(f, " \"bond_stiffness\": ["); for (int b = 0; b < bonds->n_bonds; b++) { if (b > 0) fputc(',', f); fprintf(f, "%.15g", bonds->stiffness[b]); } fprintf(f, "],\n"); fprintf(f, " \"bond_rest_lengths\": ["); for (int b = 0; b < bonds->n_bonds; b++) { if (b > 0) fputc(',', f); fprintf(f, "%.15g", bonds->rest_lengths[b]); } fprintf(f, "]\n"); fprintf(f, "}\n"); fclose(f); } // ======================================================================== // 主函数 // ======================================================================== int main(int argc, char **argv) { if (argc < 4) { fprintf(stderr, "用法: %s \n", argv[0]); return 1; } const char *input_dir = argv[1]; const char *output_dir = argv[2]; const char *param_path = argv[3]; clock_t t0 = clock(); SimParams params = read_params(param_path); AtomData atoms = read_coord(input_dir); BondData bonds = read_bonds(input_dir, &atoms); printf("[C-engine] 原子数=%d, 键数=%d, NT=%d, DT=%.6g, method=%s\n", atoms.n_atoms, bonds.n_bonds, params.NT, params.DT, params.method); int n = atoms.n_atoms; double *x = (double*)xmalloc(n * sizeof(double)); double *y = (double*)xmalloc(n * sizeof(double)); double *z = (double*)xmalloc(n * sizeof(double)); double *vx = (double*)xmalloc(n * sizeof(double)); double *vy = (double*)xmalloc(n * sizeof(double)); double *vz = (double*)xmalloc(n * sizeof(double)); for (int i = 0; i < n; i++) { x[i] = atoms.pos_0[i*3+0]; y[i] = atoms.pos_0[i*3+1]; z[i] = atoms.pos_0[i*3+2]; vx[i] = atoms.vel_0[i*3+0]; vy[i] = atoms.vel_0[i*3+1]; vz[i] = atoms.vel_0[i*3+2]; } /* 分配轨迹缓冲区：改用 record_steps */ int record_steps = params.NT - params.warmup_steps; Trajectory traj; traj.n_steps = record_steps; traj.n_atoms = n; traj.x = (double*)xmalloc(record_steps * n * sizeof(double) * 6); traj.y = traj.x + record_steps * n; traj.z = traj.y + record_steps * n; traj.vx = traj.z + record_steps * n; traj.vy = traj.vx + record_steps * n; traj.vz = traj.vy + record_steps * n; /* 预热 */ for (int s = 0; s < params.warmup_steps; s++) { apply_step(params.method, n, x, y, z, vx, vy, vz, atoms.masses, params.G, params.B, &bonds, atoms.fixed, atoms.pos_0, params.box_a, params.DT); } /* 记录 */ for (int s = 0; s < record_steps; s++) { for (int i = 0; i < n; i++) { traj.x[ s * n + i] = x[i]; traj.y[ s * n + i] = y[i]; traj.z[ s * n + i] = z[i]; traj.vx[s * n + i] = vx[i]; traj.vy[s * n + i] = vy[i]; traj.vz[s * n + i] = vz[i]; } apply_step(params.method, n, x, y, z, vx, vy, vz, atoms.masses, params.G, params.B, &bonds, atoms.fixed, atoms.pos_0, params.box_a, params.DT); } char out_path[512]; snprintf(out_path, sizeof(out_path), "%s/trajectory.txt", output_dir); write_trajectory_json(out_path, &traj, ¶ms, &atoms, &bonds); clock_t t1 = clock(); double elapsed = (double)(t1 - t0) / CLOCKS_PER_SEC; printf("[C-engine] 计算完成: %d 步, %.3f s\n", record_steps, elapsed); /* 清理 */ free(traj.x); free(x); free(y); free(z); free(vx); free(vy); free(vz); free(atoms.atom_ids); free(atoms.masses); free(atoms.radii); free(atoms.pos_0); free(atoms.vel_0); free(atoms.fixed); if (bonds.pairs) { free(bonds.pairs); free(bonds.stiffness); free(bonds.rest_lengths); } return 0; }