ylniu/dynamics
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

refactor: 引擎源码移至 src/，编译产物统一至 release/

Browse Source 
- 源码: engines/c/ → engines/src/c/
- 源码: engines/cpp/ → engines/src/cpp/
- 源码: engines/fortran/ → engines/src/fortran/
- 编译产物: engines/release/ (静态编译)
  - dynamics_c.exe   (516KB, C 引擎)
  - dynamics_cpp.exe (3.3MB, C++ 引擎)
  - dynamics_f90.exe (988KB, Fortran 引擎)
- compute.py: engine_map 指向 release/
- compute.py: param.json 和校准缓存移至 release/
- engine_dll.py: DLL 搜索路径指向 release/
- release/ 目录纳入 git 版本管理，用户克隆后直接可用
This commit is contained in:
Ying-Li Niu
2026-06-17 15:48:30 +08:00
parent 0e636e275d
commit a16d2239a1
13 changed files with 4937 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files
engines/src/c/Makefile
+82
View File
@@ -0,0 +1,82 @@
# engines/c/Makefile
# 跨平台编译：make → 本地系统编译
# make linux → Linux 交叉编译（需 x86_64-linux-gnu-gcc）
# make windows → Windows 交叉编译（需 x86_64-w64-mingw32-gcc）
# make macos → macOS 交叉编译（需 osxcross 工具链）
CC = gcc
CFLAGS = -O3 -march=native -Wall -Wextra
LDFLAGS = -lm
SRCS = main.c
LIB_SRC = dynamics_lib.c
# 自动检测系统
UNAME_S := $(shell uname -s 2>/dev/null || echo Windows)
# 目标文件名：统一使用 .exe 后缀（方便 Python 跨平台调用）
TARGET = build/dynamics_c.exe
# DLL 目标（平台自动选择后缀）
ifeq ($(UNAME_S),Linux)
DLL_TARGET = build/dynamics_c.so
DLL_FLAGS = -shared -fPIC
else ifeq ($(UNAME_S),Darwin)
DLL_TARGET = build/dynamics_c.dylib
DLL_FLAGS = -dynamiclib
else
DLL_TARGET = build/dynamics_c.dll
DLL_FLAGS = -shared
endif
# ── 本地编译 ─────────────────────────────────
.PHONY: all dll clean linux windows macos
all: $(TARGET)
dll: $(DLL_TARGET)
$(TARGET): $(SRCS) | build
$(CC) $(CFLAGS) -o $@ $(SRCS) $(LDFLAGS)
@echo " === C engine built: $@ ==="
$(DLL_TARGET): $(LIB_SRC) | build
$(CC) $(CFLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC) $(LDFLAGS)
@echo " === C DLL built: $@ ==="
build:
mkdir -p build
# ── 交叉编译 ─────────────────────────────────
# Linux → Linux (x86_64)
linux: CROSS_PREFIX = x86_64-linux-gnu-
linux: CC = $(CROSS_PREFIX)gcc
linux: CFLAGS = -O3 -march=x86-64 -Wall -Wextra
linux: $(SRCS) | build
$(CC) $(CFLAGS) -o build/dynamics_c_linux.exe $(SRCS) $(LDFLAGS)
@echo " === Linux binary: build/dynamics_c_linux.exe ==="
# 任意平台 → Windows (x86_64)
# 需要安装 MinGW 交叉编译器：
# apt install mingw-w64 (Debian/Ubuntu)
# brew install mingw-w64 (macOS)
windows: CROSS_PREFIX = x86_64-w64-mingw32-
windows: CC = $(CROSS_PREFIX)gcc
windows: CFLAGS = -O3 -march=x86-64 -Wall -Wextra
windows: $(SRCS) | build
$(CC) $(CFLAGS) -o build/dynamics_c_win.exe $(SRCS) $(LDFLAGS)
@echo " === Windows binary: build/dynamics_c_win.exe ==="
# 任意平台 → macOS (x86_64)
# 需要安装 osxcross 工具链
macos: CROSS_PREFIX = x86_64-apple-darwin-
macos: CC = $(CROSS_PREFIX)gcc
macos: CFLAGS = -O3 -march=x86-64 -Wall -Wextra
macos: $(SRCS) | build
$(CC) $(CFLAGS) -o build/dynamics_c_mac.exe $(SRCS) $(LDFLAGS)
@echo " === macOS binary: build/dynamics_c_mac.exe ==="
# ── 编译所有平台 ──────────────────────────────
all-platforms: linux windows macos
clean:
rm -rf build *.o
engines/src/c/dynamics_lib.c
+554
View File
@@ -0,0 +1,554 @@
/**
* engines/c/dynamics_lib.c
* -------------------------
* 纯计算 DLL：无文件 I/O，所有数据由 Python 以 NumPy 数组传入，
* 结果直接写入 Python 预分配的输出数组。
* 算法与 main.c 和 compute.py 保持完全一致。
*
* 编译（Windows DLL）:
* gcc -O3 -march=native -shared -o build/dynamics_c.dll dynamics_lib.c -lm
* 编译（Linux .so）:
* gcc -O3 -march=native -shared -fPIC -o build/dynamics_c.so dynamics_lib.c -lm
* 编译（macOS .dylib）:
* gcc -O3 -march=native -dynamiclib -o build/dynamics_c.dylib dynamics_lib.c -lm
*/
#ifdef _WIN32
# define EXPORT __declspec(dllexport)
#else
# define EXPORT __attribute__((visibility("default")))
#endif
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
/* ── 驱动力结构体 ─────────────────────────────────────────── */
typedef struct {
int n_drivers;
const int *idx; /* [n_drivers] 0-based local atom index */
const double *amp; /* [n_drivers*3] (ax,ay,az) interleaved */
const double *freq; /* [n_drivers*3] */
const double *phi; /* [n_drivers*3] radians */
const double *eq; /* [n_drivers*3] equilibrium positions */
const double *ncycles; /* [n_drivers] 0=unlimited */
const int *has_period; /* [n_drivers] */
/* mutable freeze positions (allocated internally) */
double *freeze; /* [n_drivers*3] */
} Drivers;
/* ── 加速度：保守力（弹簧键 + 均匀重力场）────────────────── */
static void accel_conservative(
int n, const double *x, const double *y, const double *z,
const double *m,
double Gx, double Gy, double Gz,
int gravity_field, int elastic_force,
int n_bonds, const int *bond_pairs,
const double *bond_k, const double *bond_r0,
double *ax, double *ay, double *az)
{
for (int i = 0; i < n; i++) {
ax[i] = gravity_field ? Gx : 0.0;
ay[i] = gravity_field ? Gy : 0.0;
az[i] = gravity_field ? Gz : 0.0;
}
if (!elastic_force || n_bonds == 0) return;
for (int b = 0; b < n_bonds; b++) {
int ii = bond_pairs[b*2];
int jj = bond_pairs[b*2+1];
double dx = x[jj] - x[ii];
double dy = y[jj] - y[ii];
double dz = z[jj] - z[ii];
double dist = sqrt(dx*dx + dy*dy + dz*dz);
if (dist < 1e-12) continue;
double k = bond_k[b];
double r0 = bond_r0[b];
double fac = k * (dist - r0) / dist;
double fx = fac * dx, fy = fac * dy, fz_b = fac * dz;
ax[ii] += fx / m[ii]; ay[ii] += fy / m[ii]; az[ii] += fz_b / m[ii];
ax[jj] -= fx / m[jj]; ay[jj] -= fy / m[jj]; az[jj] -= fz_b / m[jj];
}
}
/* ── 完整加速度（含阻尼）────────────────────────────────── */
static void accel_full(
int n, const double *x, const double *y, const double *z,
const double *vx, const double *vy, const double *vz,
const double *m,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bond_pairs,
const double *bond_k, const double *bond_r0,
double *ax, double *ay, double *az)
{
accel_conservative(n, x, y, z, m, Gx, Gy, Gz,
gravity_field, elastic_force,
n_bonds, bond_pairs, bond_k, bond_r0,
ax, ay, az);
if (damping_force) {
for (int i = 0; i < n; i++) {
ax[i] -= Bx * vx[i] / m[i];
ay[i] -= By * vy[i] / m[i];
az[i] -= Bz * vz[i] / m[i];
}
}
}
/* ── 边界：反弹（与 main.c limit_in_box 一致）────────────── */
static inline void _limit1(double *p, double *v, double lo, double hi) {
if (*p > hi) { *p = hi; *v = -fabs(*v); }
if (*p < lo) { *p = lo; *v = fabs(*v); }
}
/* ── 边界：回绕（与 main.c wrap_position 一致）──────────── */
static inline void _wrap1(double *p, double lo, double hi) {
if (*p > hi) *p = lo;
if (*p < lo) *p = hi;
}
/* ── 边界 + 固定约束（与 main.c apply_step 末尾一致）──────── */
static void apply_boundary_and_constraints(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const int *fixed, const double *pos_init,
double box_a)
{
double lo = -box_a, hi = box_a;
/* 反弹 */
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue;
_limit1(&x[i], &vx[i], lo, hi);
_limit1(&y[i], &vy[i], lo, hi);
_limit1(&z[i], &vz[i], lo, hi);
}
/* 回绕 */
for (int i = 0; i < n; i++) {
_wrap1(&x[i], lo, hi);
_wrap1(&y[i], lo, hi);
_wrap1(&z[i], lo, hi);
}
/* 逐自由度固定约束：与 main.c 和 Python apply_fixed_constraints 一致 */
for (int i = 0; i < n; i++) {
if (fixed[i*3+0]) { x[i] = pos_init[i*3+0]; vx[i] = 0.0; }
if (fixed[i*3+1]) { y[i] = pos_init[i*3+1]; vy[i] = 0.0; }
if (fixed[i*3+2]) { z[i] = pos_init[i*3+2]; vz[i] = 0.0; }
}
}
/* ══════════════════════════════════════════════════════════
* 蛙跳法（与 main.c leapfrog_step 完全一致）
* x(t), v(t-dt/2) → x(t+dt), v(t+dt/2)
* 无阻尼：纯辛积分。有阻尼：半隐式处理 α = B·dt/(2m)
* ══════════════════════════════════════════════════════════ */
static void leapfrog_step(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const double *m, const int *fixed,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bp, const double *bk, const double *br0,
double dt)
{
double *ax = (double*)alloca(n*sizeof(double)*3);
double *ay = ax+n; double *az = ay+n;
accel_conservative(n, x, y, z, m, Gx, Gy, Gz,
gravity_field, elastic_force,
n_bonds, bp, bk, br0, ax, ay, az);
int has_damp = damping_force && (Bx != 0.0 || By != 0.0 || Bz != 0.0);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue;
if (has_damp) {
double ax_ = Bx*dt/(2.0*m[i]);
double ay_ = By*dt/(2.0*m[i]);
double az_ = Bz*dt/(2.0*m[i]);
vx[i] = (vx[i]*(1.0-ax_) + ax[i]*dt) / (1.0+ax_);
vy[i] = (vy[i]*(1.0-ay_) + ay[i]*dt) / (1.0+ay_);
vz[i] = (vz[i]*(1.0-az_) + az[i]*dt) / (1.0+az_);
} else {
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;
}
}
/* ══════════════════════════════════════════════════════════
* 显式欧拉法（与 main.c explicit_euler_step 一致）
* ══════════════════════════════════════════════════════════ */
static void euler_step(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const double *m, const int *fixed,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bp, const double *bk, const double *br0,
double dt)
{
double *ax = (double*)alloca(n*sizeof(double)*3);
double *ay = ax+n; double *az = ay+n;
accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz,
gravity_field, elastic_force, damping_force,
n_bonds, bp, bk, br0, ax, ay, az);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && 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;
}
}
/* ══════════════════════════════════════════════════════════
* 隐式欧拉法（与 main.c implicit_euler_step 完全一致）
*
* main.c 逻辑：
* 1. 用 v_next ≈ (v + G·dt)/(1 + γ·dt) 预测（只含重力+阻尼，不含弹簧）
* 2. 用 (x, v_next) 计算完整加速度 a_next
* 3. v += a_next·dt; x += v·dt
* ══════════════════════════════════════════════════════════ */
static void implicit_euler_step(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const double *m, const int *fixed,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bp, const double *bk, const double *br0,
double dt)
{
double *vxn = (double*)alloca(n*sizeof(double)*3);
double *vyn = vxn+n; double *vzn = vyn+n;
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) {
vxn[i] = vyn[i] = vzn[i] = 0.0; continue;
}
double gx = Bx / m[i], gy = By / m[i], gz = Bz / m[i];
vxn[i] = (vx[i] + Gx*dt) / (1.0 + gx*dt);
vyn[i] = (vy[i] + Gy*dt) / (1.0 + gy*dt);
vzn[i] = (vz[i] + Gz*dt) / (1.0 + gz*dt);
}
double *ax = (double*)alloca(n*sizeof(double)*3);
double *ay = ax+n; double *az = ay+n;
accel_full(n, x, y, z, vxn, vyn, vzn, m, Gx, Gy, Gz, Bx, By, Bz,
gravity_field, elastic_force, damping_force,
n_bonds, bp, bk, br0, ax, ay, az);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && 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;
}
}
/* ══════════════════════════════════════════════════════════
* 中点法（与 main.c midpoint_step 完全一致）
*
* main.c 逻辑：
* 1. a = accel(x, v)
* 2. xm = x + 0.5·v·dt; vm = v + 0.5·a·dt
* 3. x = x + vm·dt (位置更新用 vm，即中点速度)
* 4. am = accel(xm, vm)
* 5. v = v + am·dt
* ══════════════════════════════════════════════════════════ */
static void midpoint_step(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const double *m, const int *fixed,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bp, const double *bk, const double *br0,
double dt)
{
/* Allocate in one block for cache locality */
double *buf = (double*)alloca(n*sizeof(double)*9);
double *ax = buf;
double *ay = ax+n; double *az = ay+n;
double *xm = az+n; double *ym = xm+n; double *zm = ym+n;
double *vxm = zm+n; double *vym = vxm+n; double *vzm = vym+n;
accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz,
gravity_field, elastic_force, damping_force,
n_bonds, bp, bk, br0, ax, ay, az);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) {
xm[i]=x[i]; ym[i]=y[i]; zm[i]=z[i];
vxm[i]=vym[i]=vzm[i]=0.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;
/* position updated with midpoint velocity (same as main.c) */
x[i] = x[i] + vxm[i]*dt;
y[i] = y[i] + vym[i]*dt;
z[i] = z[i] + vzm[i]*dt;
}
double *axm = (double*)alloca(n*sizeof(double)*3);
double *aym = axm+n; double *azm = aym+n;
accel_full(n, xm, ym, zm, vxm, vym, vzm, m, Gx, Gy, Gz, Bx, By, Bz,
gravity_field, elastic_force, damping_force,
n_bonds, bp, bk, br0, axm, aym, azm);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue;
vx[i] += axm[i]*dt;
vy[i] += aym[i]*dt;
vz[i] += azm[i]*dt;
}
}
/* ── 驱动力（与 main.c apply_driving_force 一致）────────── */
static void apply_driving(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
double t, int step, double dt, Drivers *drv)
{
(void)n;
if (!drv || drv->n_drivers == 0) return;
const double TWO_PI = 2.0 * 3.14159265358979323846;
for (int d = 0; d < drv->n_drivers; d++) {
int idx = drv->idx[d];
double fx = drv->freq[d*3+0];
double fy = drv->freq[d*3+1];
double fz = drv->freq[d*3+2];
if (drv->has_period[d]) {
double mf = fabs(fx) > fabs(fy) ? fabs(fx) : fabs(fy);
if (fabs(fz) > mf) mf = fabs(fz);
int period_steps = 0;
if (mf > 1e-12)
period_steps = (int)(drv->ncycles[d] / mf / dt);
if (step > period_steps) {
x[idx] = drv->freeze[d*3+0];
y[idx] = drv->freeze[d*3+1];
z[idx] = drv->freeze[d*3+2];
vx[idx] = vy[idx] = vz[idx] = 0.0;
continue;
}
double px = drv->eq[d*3+0] + drv->amp[d*3+0]*cos(TWO_PI*fx*t + drv->phi[d*3+0]);
double py = drv->eq[d*3+1] + drv->amp[d*3+1]*cos(TWO_PI*fy*t + drv->phi[d*3+1]);
double pz = drv->eq[d*3+2] + drv->amp[d*3+2]*cos(TWO_PI*fz*t + drv->phi[d*3+2]);
if (step == period_steps) {
drv->freeze[d*3+0] = px;
drv->freeze[d*3+1] = py;
drv->freeze[d*3+2] = pz;
}
}
x[idx] = drv->eq[d*3+0] + drv->amp[d*3+0]*cos(TWO_PI*fx*t + drv->phi[d*3+0]);
y[idx] = drv->eq[d*3+1] + drv->amp[d*3+1]*cos(TWO_PI*fy*t + drv->phi[d*3+1]);
z[idx] = drv->eq[d*3+2] + drv->amp[d*3+2]*cos(TWO_PI*fz*t + drv->phi[d*3+2]);
vx[idx] = -drv->amp[d*3+0]*TWO_PI*fx*sin(TWO_PI*fx*t + drv->phi[d*3+0]);
vy[idx] = -drv->amp[d*3+1]*TWO_PI*fy*sin(TWO_PI*fy*t + drv->phi[d*3+1]);
vz[idx] = -drv->amp[d*3+2]*TWO_PI*fz*sin(TWO_PI*fz*t + drv->phi[d*3+2]);
}
}
/* ══════════════════════════════════════════════════════════
* 导出函数：run_dynamics
*
* 与 main.c 的计算顺序完全一致：
* 1. leapfrog 初始化 v(-dt/2)
* 2. 初始驱动 t=0
* 3. 预热循环（不记录）
* 4. 记录循环：drive → record → step → boundary → constraints
*
* 参数说明（所有数组均为 C-contiguous 行优先 float64/int32）：
* n_atoms 原子数
* pos_init 初始位置 [n_atoms*3] x0,y0,z0, x1,y1,z1, ...
* vel_init 初始速度 [n_atoms*3]
* masses 质量 [n_atoms]
* fixed 自由度约束 [n_atoms*3] int32, 1=固定
* n_bonds 键数
* bond_pairs 键对 [n_bonds*2] int32, 0-based local index
* bond_k 刚度 [n_bonds]
* bond_r0 平衡键长 [n_bonds]
* box_a 盒子半边长
* dt 时间步长
* NT 总步数（含预热）
* NSTEP 抽帧间隔
* warmup_steps 预热步数
* method_id 0=euler 1=implicit 2=midpoint 3=leapfrog
* Gx/Gy/Gz 均匀重力场加速度分量
* Bx/By/Bz 阻尼系数分量
* gravity_field / elastic_force / damping_force 力开关
* gravity_strength 原子间引力强度（暂未实现，留接口）
* n_drivers 驱动原子数
* drv_idx 驱动原子局部索引 [n_drivers] int32
* drv_amp 振幅 [n_drivers*3]
* drv_freq 频率 [n_drivers*3]
* drv_phi 初相（弧度）[n_drivers*3]
* drv_eq 平衡位置 [n_drivers*3]
* drv_ncycles 周期数 [n_drivers] 0=不限
* drv_has_period [n_drivers] int32
* n_frames 输出帧数（Python 预计算：(NT-warmup)/NSTEP 向上取整）
* out_x/y/z/vx/vy/vz 输出数组 [n_frames*n_atoms] 由 Python 预分配
* progress_cb 进度回调（可为 NULL）
*
* 返回：0=成功，负数=错误
* ══════════════════════════════════════════════════════════ */
EXPORT int run_dynamics(
int n_atoms,
const double *pos_init,
const double *vel_init,
const double *masses,
const int *fixed,
int n_bonds,
const int *bond_pairs,
const double *bond_k,
const double *bond_r0,
double box_a, double dt,
int NT, int NSTEP, int warmup_steps, int method_id,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
double gravity_strength,
int n_drivers,
const int *drv_idx,
const double *drv_amp,
const double *drv_freq,
const double *drv_phi,
const double *drv_eq,
const double *drv_ncycles,
const int *drv_has_period,
int n_frames,
double *out_x, double *out_y, double *out_z,
double *out_vx, double *out_vy, double *out_vz,
void (*progress_cb)(int step, int total))
{
(void)gravity_strength; /* 原子间引力暂未实现 */
int n = n_atoms;
/* ── 工作数组 ── */
double *x = (double*)malloc(n*sizeof(double));
double *y = (double*)malloc(n*sizeof(double));
double *z = (double*)malloc(n*sizeof(double));
double *vx = (double*)malloc(n*sizeof(double));
double *vy = (double*)malloc(n*sizeof(double));
double *vz = (double*)malloc(n*sizeof(double));
if (!x||!y||!z||!vx||!vy||!vz) return -1;
for (int i = 0; i < n; i++) {
x[i]=pos_init[i*3+0]; y[i]=pos_init[i*3+1]; z[i]=pos_init[i*3+2];
vx[i]=vel_init[i*3+0]; vy[i]=vel_init[i*3+1]; vz[i]=vel_init[i*3+2];
}
/* ── 驱动结构 ── */
Drivers drv;
drv.n_drivers = n_drivers;
drv.idx = drv_idx;
drv.amp = drv_amp;
drv.freq = drv_freq;
drv.phi = drv_phi;
drv.eq = drv_eq;
drv.ncycles = drv_ncycles;
drv.has_period = drv_has_period;
drv.freeze = NULL;
if (n_drivers > 0) {
drv.freeze = (double*)calloc(n_drivers*3, sizeof(double));
if (!drv.freeze) { free(x);free(y);free(z);free(vx);free(vy);free(vz); return -2; }
}
/* ── 内联步进宏 ── */
#define DO_STEP() do { \
switch (method_id) { \
case 0: euler_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \
gravity_field,elastic_force,damping_force, \
n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \
case 1: implicit_euler_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \
gravity_field,elastic_force,damping_force, \
n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \
case 2: midpoint_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \
gravity_field,elastic_force,damping_force, \
n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \
default: leapfrog_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \
gravity_field,elastic_force,damping_force, \
n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \
} \
apply_boundary_and_constraints(n,x,y,z,vx,vy,vz,fixed,pos_init,box_a); \
} while(0)
/* ── 蛙跳法：初始化 v(-dt/2) = v(0) - 0.5·a_c(0)·dt ── */
if (method_id == 3) {
double *ax0 = (double*)alloca(n*sizeof(double)*3);
double *ay0 = ax0+n; double *az0 = ay0+n;
accel_conservative(n, x, y, z, masses, Gx, Gy, Gz,
gravity_field, elastic_force,
n_bonds, bond_pairs, bond_k, bond_r0,
ax0, ay0, az0);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue;
vx[i] -= 0.5*ax0[i]*dt;
vy[i] -= 0.5*ay0[i]*dt;
vz[i] -= 0.5*az0[i]*dt;
}
}
/* ── 初始驱动 t=0（与 main.c 一致：leapfrog init 之后施加）── */
if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, 0.0, 0, dt, &drv);
/* ── 预热（不记录）── */
for (int s = 0; s < warmup_steps; s++) {
double tw = (s + 1) * dt;
if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, tw, s, dt, &drv);
DO_STEP();
}
/* ── 记录循环 ── */
int record_steps = NT - warmup_steps;
int prog_interval = record_steps / 100;
if (prog_interval < 1) prog_interval = 1;
int frame_idx = 0;
for (int s = 0; s < record_steps; s++) {
if (progress_cb && s % prog_interval == 0 && s > 0)
progress_cb(s, record_steps);
double t = (s + warmup_steps) * dt;
if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, t, s, dt, &drv);
/* 抽帧记录（drive 之后，step 之前，与 main.c 一致）*/
if (s % NSTEP == 0 && frame_idx < n_frames) {
int base = frame_idx * n;
for (int i = 0; i < n; i++) {
out_x [base+i] = x[i]; out_y [base+i] = y[i]; out_z [base+i] = z[i];
out_vx[base+i] = vx[i]; out_vy[base+i] = vy[i]; out_vz[base+i] = vz[i];
}
frame_idx++;
}
DO_STEP();
}
#undef DO_STEP
free(x); free(y); free(z);
free(vx); free(vy); free(vz);
if (drv.freeze) free(drv.freeze);
return 0;
}
engines/src/c/main.c
+1114
View File
File diff suppressed because it is too large Load Diff
engines/src/cpp/Makefile
+49
View File
@@ -0,0 +1,49 @@
# engines/cpp/Makefile
CXX = g++
SRCS = main.cpp
LIB_SRC = dynamics_lib.cpp
UNAME_S := $(shell uname -s 2>/dev/null || echo Windows)
CXXFLAGS = -O3 -march=native -std=c++17 -Wall -Wextra -D_USE_MATH_DEFINES
# Windows 下静态链接运行时，避免 libstdc++-6.dll / libgcc_s_seh-1.dll 版本冲突
ifeq ($(UNAME_S),Windows)
STATIC_FLAGS = -static-libgcc -static-libstdc++
else
STATIC_FLAGS =
endif
TARGET = build/dynamics_cpp.exe
ifeq ($(UNAME_S),Linux)
DLL_TARGET = build/dynamics_cpp.so
DLL_FLAGS = -shared -fPIC
else ifeq ($(UNAME_S),Darwin)
DLL_TARGET = build/dynamics_cpp.dylib
DLL_FLAGS = -dynamiclib
else
DLL_TARGET = build/dynamics_cpp.dll
DLL_FLAGS = -shared
endif
.PHONY: all dll clean
all: $(TARGET)
dll: $(DLL_TARGET)
$(TARGET): $(SRCS) | build
$(CXX) $(CXXFLAGS) $(STATIC_FLAGS) -o $@ $(SRCS)
@echo " === C++ engine built: $@ ==="
$(DLL_TARGET): $(LIB_SRC) | build
$(CXX) $(CXXFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC)
@echo " === C++ DLL built: $@ ==="
build:
mkdir -p build
clean:
rm -rf build *.o
engines/src/cpp/dynamics_lib.cpp
+450
View File
@@ -0,0 +1,450 @@
/**
* engines/cpp/dynamics_lib.cpp
* -----------------------------
* 纯计算 DLL（C++ 版）：无文件 I/O，所有数据由 Python 以 NumPy 数组传入。
* 算法与 main.cpp / compute.py 保持完全一致。
*
* 编译（Windows）:
* g++ -O3 -march=native -std=c++17 -shared -o build/dynamics_cpp.dll dynamics_lib.cpp
* 编译（Linux）:
* g++ -O3 -march=native -std=c++17 -shared -fPIC -o build/dynamics_cpp.so dynamics_lib.cpp
* 编译（macOS）:
* g++ -O3 -march=native -std=c++17 -dynamiclib -o build/dynamics_cpp.dylib dynamics_lib.cpp
*/
#ifdef _WIN32
# define EXPORT extern "C" __declspec(dllexport)
#else
# define EXPORT extern "C" __attribute__((visibility("default")))
#endif
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <vector>
/* ── 驱动力结构体 ─────────────────────────────────────────── */
struct Drivers {
int n_drivers = 0;
const int *idx = nullptr;
const double *amp = nullptr;
const double *freq = nullptr;
const double *phi = nullptr;
const double *eq = nullptr;
const double *ncycles = nullptr;
const int *has_period = nullptr;
std::vector<double> freeze; /* [n_drivers*3] 冻结位置（period 结束时锁定）*/
};
/* ── 加速度：保守力（弹簧键 + 均匀重力场）────────────────── */
static void accel_conservative(
int n, const double *x, const double *y, const double *z,
const double *m,
double Gx, double Gy, double Gz,
int gravity_field, int elastic_force,
int n_bonds, const int *bond_pairs,
const double *bond_k, const double *bond_r0,
double *ax, double *ay, double *az)
{
for (int i = 0; i < n; i++) {
ax[i] = gravity_field ? Gx : 0.0;
ay[i] = gravity_field ? Gy : 0.0;
az[i] = gravity_field ? Gz : 0.0;
}
if (!elastic_force || n_bonds == 0) return;
for (int b = 0; b < n_bonds; b++) {
int ii = bond_pairs[b*2];
int jj = bond_pairs[b*2+1];
double dx = x[jj]-x[ii], dy = y[jj]-y[ii], dz = z[jj]-z[ii];
double dist = std::sqrt(dx*dx + dy*dy + dz*dz);
if (dist < 1e-12) continue;
double fac = bond_k[b] * (dist - bond_r0[b]) / dist;
double fx = fac*dx, fy = fac*dy, fz_b = fac*dz;
ax[ii] += fx/m[ii]; ay[ii] += fy/m[ii]; az[ii] += fz_b/m[ii];
ax[jj] -= fx/m[jj]; ay[jj] -= fy/m[jj]; az[jj] -= fz_b/m[jj];
}
}
/* ── 完整加速度（含阻尼）────────────────────────────────── */
static void accel_full(
int n, const double *x, const double *y, const double *z,
const double *vx, const double *vy, const double *vz,
const double *m,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bond_pairs,
const double *bond_k, const double *bond_r0,
double *ax, double *ay, double *az)
{
accel_conservative(n, x, y, z, m, Gx, Gy, Gz,
gravity_field, elastic_force,
n_bonds, bond_pairs, bond_k, bond_r0,
ax, ay, az);
if (damping_force) {
for (int i = 0; i < n; i++) {
ax[i] -= Bx * vx[i] / m[i];
ay[i] -= By * vy[i] / m[i];
az[i] -= Bz * vz[i] / m[i];
}
}
}
/* ── 边界：反弹 ──────────────────────────────────────────── */
static inline void _limit1(double &p, double &v, double lo, double hi) {
if (p > hi) { p = hi; v = -std::fabs(v); }
if (p < lo) { p = lo; v = std::fabs(v); }
}
/* ── 边界：回绕 ──────────────────────────────────────────── */
static inline void _wrap1(double &p, double lo, double hi) {
if (p > hi) p = lo;
if (p < lo) p = hi;
}
/* ── 边界 + 固定约束 ────────────────────────────────────── */
static void apply_boundary_and_constraints(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const int *fixed, const double *pos_init, double box_a)
{
double lo = -box_a, hi = box_a;
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue;
_limit1(x[i], vx[i], lo, hi);
_limit1(y[i], vy[i], lo, hi);
_limit1(z[i], vz[i], lo, hi);
}
for (int i = 0; i < n; i++) {
_wrap1(x[i], lo, hi);
_wrap1(y[i], lo, hi);
_wrap1(z[i], lo, hi);
}
for (int i = 0; i < n; i++) {
if (fixed[i*3+0]) { x[i] = pos_init[i*3+0]; vx[i] = 0.0; }
if (fixed[i*3+1]) { y[i] = pos_init[i*3+1]; vy[i] = 0.0; }
if (fixed[i*3+2]) { z[i] = pos_init[i*3+2]; vz[i] = 0.0; }
}
}
/* ══════════════════════════════════════════════════════════
* 蛙跳法（与 main.cpp leapfrog_step 完全一致）
* ══════════════════════════════════════════════════════════ */
static void leapfrog_step(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const double *m, const int *fixed,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bp, const double *bk, const double *br0, double dt)
{
std::vector<double> buf(n * 3);
double *ax = buf.data(), *ay = ax+n, *az = ay+n;
accel_conservative(n, x, y, z, m, Gx, Gy, Gz,
gravity_field, elastic_force,
n_bonds, bp, bk, br0, ax, ay, az);
bool has_damp = damping_force && (Bx != 0.0 || By != 0.0 || Bz != 0.0);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue;
if (has_damp) {
double ax_ = Bx*dt/(2.0*m[i]);
double ay_ = By*dt/(2.0*m[i]);
double az_ = Bz*dt/(2.0*m[i]);
vx[i] = (vx[i]*(1.0-ax_) + ax[i]*dt) / (1.0+ax_);
vy[i] = (vy[i]*(1.0-ay_) + ay[i]*dt) / (1.0+ay_);
vz[i] = (vz[i]*(1.0-az_) + az[i]*dt) / (1.0+az_);
} else {
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 euler_step(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const double *m, const int *fixed,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bp, const double *bk, const double *br0, double dt)
{
std::vector<double> buf(n * 3);
double *ax = buf.data(), *ay = ax+n, *az = ay+n;
accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz,
gravity_field, elastic_force, damping_force,
n_bonds, bp, bk, br0, ax, ay, az);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && 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;
}
}
/* ══════════════════════════════════════════════════════════
* 隐式欧拉法（与 main.cpp implicit_euler_step 完全一致）
* ══════════════════════════════════════════════════════════ */
static void implicit_euler_step(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const double *m, const int *fixed,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bp, const double *bk, const double *br0, double dt)
{
std::vector<double> vbuf(n * 3), abuf(n * 3);
double *vxn = vbuf.data(), *vyn = vxn+n, *vzn = vyn+n;
double *ax = abuf.data(), *ay = ax+n, *az = ay+n;
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) {
vxn[i] = vyn[i] = vzn[i] = 0.0; continue;
}
double gx = Bx/m[i], gy = By/m[i], gz = Bz/m[i];
vxn[i] = (vx[i] + Gx*dt) / (1.0 + gx*dt);
vyn[i] = (vy[i] + Gy*dt) / (1.0 + gy*dt);
vzn[i] = (vz[i] + Gz*dt) / (1.0 + gz*dt);
}
accel_full(n, x, y, z, vxn, vyn, vzn, m, Gx, Gy, Gz, Bx, By, Bz,
gravity_field, elastic_force, damping_force,
n_bonds, bp, bk, br0, ax, ay, az);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && 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;
}
}
/* ══════════════════════════════════════════════════════════
* 中点法（与 main.cpp midpoint_step 完全一致）
* ══════════════════════════════════════════════════════════ */
static void midpoint_step(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
const double *m, const int *fixed,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
int n_bonds, const int *bp, const double *bk, const double *br0, double dt)
{
std::vector<double> buf(n * 9);
double *ax = buf.data();
double *ay = ax+n; double *az = ay+n;
double *xm = az+n; double *ym = xm+n; double *zm = ym+n;
double *vxm = zm+n; double *vym = vxm+n; double *vzm = vym+n;
accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz,
gravity_field, elastic_force, damping_force,
n_bonds, bp, bk, br0, ax, ay, az);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) {
xm[i]=x[i]; ym[i]=y[i]; zm[i]=z[i];
vxm[i]=vym[i]=vzm[i]=0.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;
}
std::vector<double> abuf(n * 3);
double *axm = abuf.data(), *aym = axm+n, *azm = aym+n;
accel_full(n, xm, ym, zm, vxm, vym, vzm, m, Gx, Gy, Gz, Bx, By, Bz,
gravity_field, elastic_force, damping_force,
n_bonds, bp, bk, br0, axm, aym, azm);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue;
vx[i] += axm[i]*dt;
vy[i] += aym[i]*dt;
vz[i] += azm[i]*dt;
}
}
/* ── 驱动力 ─────────────────────────────────────────────── */
static void apply_driving(
int n, double *x, double *y, double *z,
double *vx, double *vy, double *vz,
double t, int step, double dt, Drivers &drv)
{
(void)n;
if (drv.n_drivers == 0) return;
constexpr double TWO_PI = 2.0 * 3.14159265358979323846;
for (int d = 0; d < drv.n_drivers; d++) {
int idx = drv.idx[d];
double fx = drv.freq[d*3+0];
double fy = drv.freq[d*3+1];
double fz = drv.freq[d*3+2];
if (drv.has_period[d]) {
double mf = std::fabs(fx) > std::fabs(fy) ? std::fabs(fx) : std::fabs(fy);
if (std::fabs(fz) > mf) mf = std::fabs(fz);
int period_steps = 0;
if (mf > 1e-12)
period_steps = (int)(drv.ncycles[d] / mf / dt);
if (step > period_steps) {
x[idx] = drv.freeze[d*3+0];
y[idx] = drv.freeze[d*3+1];
z[idx] = drv.freeze[d*3+2];
vx[idx] = vy[idx] = vz[idx] = 0.0;
continue;
}
double px = drv.eq[d*3+0] + drv.amp[d*3+0]*std::cos(TWO_PI*fx*t + drv.phi[d*3+0]);
double py = drv.eq[d*3+1] + drv.amp[d*3+1]*std::cos(TWO_PI*fy*t + drv.phi[d*3+1]);
double pz = drv.eq[d*3+2] + drv.amp[d*3+2]*std::cos(TWO_PI*fz*t + drv.phi[d*3+2]);
if (step == period_steps) {
drv.freeze[d*3+0] = px;
drv.freeze[d*3+1] = py;
drv.freeze[d*3+2] = pz;
}
}
x[idx] = drv.eq[d*3+0] + drv.amp[d*3+0]*std::cos(TWO_PI*fx*t + drv.phi[d*3+0]);
y[idx] = drv.eq[d*3+1] + drv.amp[d*3+1]*std::cos(TWO_PI*fy*t + drv.phi[d*3+1]);
z[idx] = drv.eq[d*3+2] + drv.amp[d*3+2]*std::cos(TWO_PI*fz*t + drv.phi[d*3+2]);
vx[idx] = -drv.amp[d*3+0]*TWO_PI*fx*std::sin(TWO_PI*fx*t + drv.phi[d*3+0]);
vy[idx] = -drv.amp[d*3+1]*TWO_PI*fy*std::sin(TWO_PI*fy*t + drv.phi[d*3+1]);
vz[idx] = -drv.amp[d*3+2]*TWO_PI*fz*std::sin(TWO_PI*fz*t + drv.phi[d*3+2]);
}
}
/* ══════════════════════════════════════════════════════════
* 导出函数：run_dynamics（接口与 C 版完全相同）
* ══════════════════════════════════════════════════════════ */
EXPORT int run_dynamics(
int n_atoms,
const double *pos_init,
const double *vel_init,
const double *masses,
const int *fixed,
int n_bonds,
const int *bond_pairs,
const double *bond_k,
const double *bond_r0,
double box_a, double dt,
int NT, int NSTEP, int warmup_steps, int method_id,
double Gx, double Gy, double Gz,
double Bx, double By, double Bz,
int gravity_field, int elastic_force, int damping_force,
double gravity_strength,
int n_drivers,
const int *drv_idx,
const double *drv_amp,
const double *drv_freq,
const double *drv_phi,
const double *drv_eq,
const double *drv_ncycles,
const int *drv_has_period,
int n_frames,
double *out_x, double *out_y, double *out_z,
double *out_vx, double *out_vy, double *out_vz,
void (*progress_cb)(int step, int total))
{
(void)gravity_strength;
int n = n_atoms;
std::vector<double> xv(n), yv(n), zv(n);
std::vector<double> vxv(n), vyv(n), vzv(n);
for (int i = 0; i < n; i++) {
xv[i]=pos_init[i*3+0]; yv[i]=pos_init[i*3+1]; zv[i]=pos_init[i*3+2];
vxv[i]=vel_init[i*3+0]; vyv[i]=vel_init[i*3+1]; vzv[i]=vel_init[i*3+2];
}
double *x=xv.data(), *y=yv.data(), *z=zv.data();
double *vx=vxv.data(), *vy=vyv.data(), *vz=vzv.data();
Drivers drv;
drv.n_drivers = n_drivers;
drv.idx = drv_idx;
drv.amp = drv_amp;
drv.freq = drv_freq;
drv.phi = drv_phi;
drv.eq = drv_eq;
drv.ncycles = drv_ncycles;
drv.has_period = drv_has_period;
if (n_drivers > 0)
drv.freeze.assign(n_drivers * 3, 0.0);
#define DO_STEP() do { \
switch (method_id) { \
case 0: euler_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \
gravity_field,elastic_force,damping_force, \
n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \
case 1: implicit_euler_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \
gravity_field,elastic_force,damping_force, \
n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \
case 2: midpoint_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \
gravity_field,elastic_force,damping_force, \
n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \
default: leapfrog_step(n,x,y,z,vx,vy,vz,masses,fixed,Gx,Gy,Gz,Bx,By,Bz, \
gravity_field,elastic_force,damping_force, \
n_bonds,bond_pairs,bond_k,bond_r0,dt); break; \
} \
apply_boundary_and_constraints(n,x,y,z,vx,vy,vz,fixed,pos_init,box_a); \
} while(0)
/* 蛙跳法：初始化 v(-dt/2) */
if (method_id == 3) {
std::vector<double> ibuf(n * 3);
double *ax0=ibuf.data(), *ay0=ax0+n, *az0=ay0+n;
accel_conservative(n, x, y, z, masses, Gx, Gy, Gz,
gravity_field, elastic_force,
n_bonds, bond_pairs, bond_k, bond_r0,
ax0, ay0, az0);
for (int i = 0; i < n; i++) {
if (fixed[i*3] && fixed[i*3+1] && fixed[i*3+2]) continue;
vx[i] -= 0.5*ax0[i]*dt;
vy[i] -= 0.5*ay0[i]*dt;
vz[i] -= 0.5*az0[i]*dt;
}
}
/* 初始驱动 t=0 */
if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, 0.0, 0, dt, drv);
/* 预热 */
for (int s = 0; s < warmup_steps; s++) {
double tw = (s + 1) * dt;
if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, tw, s, dt, drv);
DO_STEP();
}
/* 记录循环 */
int record_steps = NT - warmup_steps;
int prog_interval = std::max(1, record_steps / 100);
int frame_idx = 0;
for (int s = 0; s < record_steps; s++) {
if (progress_cb && s % prog_interval == 0 && s > 0)
progress_cb(s, record_steps);
double t = (s + warmup_steps) * dt;
if (n_drivers > 0) apply_driving(n, x, y, z, vx, vy, vz, t, s, dt, drv);
if (s % NSTEP == 0 && frame_idx < n_frames) {
int base = frame_idx * n;
for (int i = 0; i < n; i++) {
out_x [base+i] = x[i]; out_y [base+i] = y[i]; out_z [base+i] = z[i];
out_vx[base+i] = vx[i]; out_vy[base+i] = vy[i]; out_vz[base+i] = vz[i];
}
frame_idx++;
}
DO_STEP();
}
#undef DO_STEP
return 0;
}
engines/src/cpp/main.cpp
+1025
View File
File diff suppressed because it is too large Load Diff
engines/src/fortran/Makefile
+47
View File
@@ -0,0 +1,47 @@
# engines/fortran/Makefile
FC = gfortran
FFLAGS = -O3 -march=native -Wall -Wextra
SRCS = main.f90
LIB_SRC = dynamics_lib.f90
UNAME_S := $(shell uname -s 2>/dev/null || echo Windows)
ifeq ($(UNAME_S),Windows)
STATIC_FLAGS = -static-libgcc -static-libgfortran -static-libquadmath
else
STATIC_FLAGS =
endif
TARGET = build/dynamics_f90.exe
ifeq ($(UNAME_S),Linux)
DLL_TARGET = build/dynamics_f90.so
DLL_FLAGS = -shared -fPIC
else ifeq ($(UNAME_S),Darwin)
DLL_TARGET = build/dynamics_f90.dylib
DLL_FLAGS = -dynamiclib
else
DLL_TARGET = build/dynamics_f90.dll
DLL_FLAGS = -shared -fPIC
endif
.PHONY: all dll clean
all: $(TARGET)
dll: $(DLL_TARGET)
$(TARGET): $(SRCS) | build
$(FC) $(FFLAGS) $(STATIC_FLAGS) -o $@ $(SRCS)
@echo " === Fortran engine built: $@ ==="
$(DLL_TARGET): $(LIB_SRC) | build
$(FC) $(FFLAGS) $(STATIC_FLAGS) $(DLL_FLAGS) -o $@ $(LIB_SRC)
@echo " === Fortran DLL built: $@ ==="
build:
mkdir -p build
clean:
rm -rf build *.o *.mod
engines/src/fortran/dynamics_lib.f90
+483
View File
@@ -0,0 +1,483 @@
! engines/fortran/dynamics_lib.f90
! ---------------------------------
! 纯计算 DLL（Fortran 版）：无文件 I/O，由 Python ctypes 调用。
! 算法与 main.f90 / compute.py 完全一致。
! 使用 iso_c_binding 导出 C 兼容接口。
!
! 编译（Windows）:
! gfortran -O3 -march=native -shared -fPIC -o build/dynamics_f90.dll dynamics_lib.f90
! 编译（Linux）:
! gfortran -O3 -march=native -shared -fPIC -o build/dynamics_f90.so dynamics_lib.f90
! 编译（macOS）:
! gfortran -O3 -march=native -dynamiclib -o build/dynamics_f90.dylib dynamics_lib.f90
module dynamics_dll
use iso_c_binding, only: c_int, c_double, c_funptr, c_f_procpointer, c_associated
implicit none
private
real(c_double), parameter :: TWO_PI = 2.0d0 * 3.14159265358979323846d0
public :: run_dynamics
contains
! ── 保守加速度 ───────────────────────────────────────────────
subroutine accel_conservative(n, x, y, z, m, Gx, Gy, Gz, &
gravity_field, elastic_force, &
n_bonds, bond_pairs, bond_k, bond_r0, &
ax, ay, az)
integer, intent(in) :: n, gravity_field, elastic_force, n_bonds
real(c_double), intent(in) :: x(n), y(n), z(n), m(n)
real(c_double), intent(in) :: Gx, Gy, Gz
integer, intent(in) :: bond_pairs(2, n_bonds)
real(c_double), intent(in) :: bond_k(n_bonds), bond_r0(n_bonds)
real(c_double), intent(out) :: ax(n), ay(n), az(n)
integer :: b, ii, jj
real(c_double) :: dx, dy, dz, dist, fac, fx, fy, fz_b
if (gravity_field /= 0) then
ax = Gx; ay = Gy; az = Gz
else
ax = 0.0d0; ay = 0.0d0; az = 0.0d0
end if
if (elastic_force == 0 .or. n_bonds == 0) return
do b = 1, n_bonds
ii = bond_pairs(1, b) + 1 ! 0-based → 1-based
jj = bond_pairs(2, b) + 1
dx = x(jj)-x(ii); dy = y(jj)-y(ii); dz = z(jj)-z(ii)
dist = sqrt(dx*dx + dy*dy + dz*dz)
if (dist < 1.0d-12) cycle
fac = bond_k(b) * (dist - bond_r0(b)) / dist
fx = fac*dx; fy = fac*dy; fz_b = fac*dz
ax(ii) = ax(ii) + fx/m(ii); ay(ii) = ay(ii) + fy/m(ii); az(ii) = az(ii) + fz_b/m(ii)
ax(jj) = ax(jj) - fx/m(jj); ay(jj) = ay(jj) - fy/m(jj); az(jj) = az(jj) - fz_b/m(jj)
end do
end subroutine
! ── 完整加速度（含阻尼）──────────────────────────────────────
subroutine accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az)
integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds
real(c_double), intent(in) :: x(n), y(n), z(n), vx(n), vy(n), vz(n), m(n)
real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz
integer, intent(in) :: bond_pairs(2, n_bonds)
real(c_double), intent(in) :: bond_k(n_bonds), bond_r0(n_bonds)
real(c_double), intent(out) :: ax(n), ay(n), az(n)
integer :: i
call accel_conservative(n, x, y, z, m, Gx, Gy, Gz, &
gravity_field, elastic_force, &
n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az)
if (damping_force /= 0) then
do i = 1, n
ax(i) = ax(i) - Bx*vx(i)/m(i)
ay(i) = ay(i) - By*vy(i)/m(i)
az(i) = az(i) - Bz*vz(i)/m(i)
end do
end if
end subroutine
! ── 边界 + 固定约束 ──────────────────────────────────────────
subroutine apply_bc(n, x, y, z, vx, vy, vz, fixed, pos_init, box_a)
integer, intent(in) :: n
real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
integer, intent(in) :: fixed(3, n)
real(c_double), intent(in) :: pos_init(3, n), box_a
integer :: i
real(c_double) :: lo, hi
lo = -box_a; hi = box_a
! 反弹
do i = 1, n
if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) cycle
if (x(i)>hi) then; x(i)=hi; vx(i)=-abs(vx(i)); end if
if (x(i)<lo) then; x(i)=lo; vx(i)= abs(vx(i)); end if
if (y(i)>hi) then; y(i)=hi; vy(i)=-abs(vy(i)); end if
if (y(i)<lo) then; y(i)=lo; vy(i)= abs(vy(i)); end if
if (z(i)>hi) then; z(i)=hi; vz(i)=-abs(vz(i)); end if
if (z(i)<lo) then; z(i)=lo; vz(i)= abs(vz(i)); end if
end do
! 回绕
do i = 1, n
if (x(i)>hi) x(i)=lo; if (x(i)<lo) x(i)=hi
if (y(i)>hi) y(i)=lo; if (y(i)<lo) y(i)=hi
if (z(i)>hi) z(i)=lo; if (z(i)<lo) z(i)=hi
end do
! 逐自由度固定约束
do i = 1, n
if (fixed(1,i)/=0) then; x(i)=pos_init(1,i); vx(i)=0.0d0; end if
if (fixed(2,i)/=0) then; y(i)=pos_init(2,i); vy(i)=0.0d0; end if
if (fixed(3,i)/=0) then; z(i)=pos_init(3,i); vz(i)=0.0d0; end if
end do
end subroutine
! ── 蛙跳法 ───────────────────────────────────────────────────
subroutine leapfrog_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt)
integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds
real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
real(c_double), intent(in) :: m(n), bond_k(n_bonds), bond_r0(n_bonds)
integer, intent(in) :: fixed(3,n), bond_pairs(2,n_bonds)
real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz, dt
real(c_double) :: ax(n), ay(n), az(n), ax_, ay_, az_
logical :: has_damp
integer :: i
call accel_conservative(n, x, y, z, m, Gx, Gy, Gz, &
gravity_field, elastic_force, &
n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az)
has_damp = (damping_force/=0) .and. (abs(Bx)+abs(By)+abs(Bz) > 0.0d0)
do i = 1, n
if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) cycle
if (has_damp) then
ax_ = Bx*dt/(2.0d0*m(i)); ay_ = By*dt/(2.0d0*m(i)); az_ = Bz*dt/(2.0d0*m(i))
vx(i) = (vx(i)*(1.0d0-ax_) + ax(i)*dt)/(1.0d0+ax_)
vy(i) = (vy(i)*(1.0d0-ay_) + ay(i)*dt)/(1.0d0+ay_)
vz(i) = (vz(i)*(1.0d0-az_) + az(i)*dt)/(1.0d0+az_)
else
vx(i) = vx(i)+ax(i)*dt; vy(i) = vy(i)+ay(i)*dt; vz(i) = vz(i)+az(i)*dt
end if
x(i) = x(i)+vx(i)*dt; y(i) = y(i)+vy(i)*dt; z(i) = z(i)+vz(i)*dt
end do
end subroutine
! ── 显式欧拉法 ───────────────────────────────────────────────
subroutine euler_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt)
integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds
real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
real(c_double), intent(in) :: m(n), bond_k(n_bonds), bond_r0(n_bonds)
integer, intent(in) :: fixed(3,n), bond_pairs(2,n_bonds)
real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz, dt
real(c_double) :: ax(n), ay(n), az(n)
integer :: i
call accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az)
do i = 1, n
if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) cycle
x(i) = x(i)+vx(i)*dt; y(i) = y(i)+vy(i)*dt; z(i) = z(i)+vz(i)*dt
vx(i)= vx(i)+ax(i)*dt; vy(i)= vy(i)+ay(i)*dt; vz(i)= vz(i)+az(i)*dt
end do
end subroutine
! ── 隐式欧拉法 ───────────────────────────────────────────────
subroutine implicit_euler_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt)
integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds
real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
real(c_double), intent(in) :: m(n), bond_k(n_bonds), bond_r0(n_bonds)
integer, intent(in) :: fixed(3,n), bond_pairs(2,n_bonds)
real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz, dt
real(c_double) :: vxn(n), vyn(n), vzn(n), ax(n), ay(n), az(n)
real(c_double) :: gamma_x, gamma_y, gamma_z
integer :: i
do i = 1, n
if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) then
vxn(i)=0.0d0; vyn(i)=0.0d0; vzn(i)=0.0d0; cycle
end if
gamma_x = Bx/m(i); gamma_y = By/m(i); gamma_z = Bz/m(i)
vxn(i) = (vx(i)+Gx*dt)/(1.0d0+gamma_x*dt)
vyn(i) = (vy(i)+Gy*dt)/(1.0d0+gamma_y*dt)
vzn(i) = (vz(i)+Gz*dt)/(1.0d0+gamma_z*dt)
end do
call accel_full(n, x, y, z, vxn, vyn, vzn, m, Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az)
do i = 1, n
if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) cycle
vx(i)=vx(i)+ax(i)*dt; vy(i)=vy(i)+ay(i)*dt; vz(i)=vz(i)+az(i)*dt
x(i) =x(i) +vx(i)*dt; y(i) =y(i) +vy(i)*dt; z(i) =z(i) +vz(i)*dt
end do
end subroutine
! ── 中点法 ───────────────────────────────────────────────────
subroutine midpoint_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt)
integer, intent(in) :: n, gravity_field, elastic_force, damping_force, n_bonds
real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
real(c_double), intent(in) :: m(n), bond_k(n_bonds), bond_r0(n_bonds)
integer, intent(in) :: fixed(3,n), bond_pairs(2,n_bonds)
real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz, dt
real(c_double) :: ax(n), ay(n), az(n)
real(c_double) :: xm(n), ym(n), zm(n), vxm(n), vym(n), vzm(n)
real(c_double) :: axm(n), aym(n), azm(n)
integer :: i
call accel_full(n, x, y, z, vx, vy, vz, m, Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, ax, ay, az)
do i = 1, n
if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) then
xm(i)=x(i); ym(i)=y(i); zm(i)=z(i)
vxm(i)=0.0d0; vym(i)=0.0d0; vzm(i)=0.0d0; cycle
end if
xm(i) = x(i) +0.5d0*vx(i)*dt; ym(i) = y(i) +0.5d0*vy(i)*dt; zm(i) = z(i) +0.5d0*vz(i)*dt
vxm(i) = vx(i)+0.5d0*ax(i)*dt; vym(i) = vy(i)+0.5d0*ay(i)*dt; vzm(i) = vz(i)+0.5d0*az(i)*dt
x(i) = x(i) +vxm(i)*dt; y(i) = y(i) +vym(i)*dt; z(i) = z(i) +vzm(i)*dt
end do
call accel_full(n, xm, ym, zm, vxm, vym, vzm, m, Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, axm, aym, azm)
do i = 1, n
if (fixed(1,i)/=0 .and. fixed(2,i)/=0 .and. fixed(3,i)/=0) cycle
vx(i)=vx(i)+axm(i)*dt; vy(i)=vy(i)+aym(i)*dt; vz(i)=vz(i)+azm(i)*dt
end do
end subroutine
! ── 驱动力施加 ───────────────────────────────────────────────
subroutine apply_drive(n, x, y, z, vx, vy, vz, t, step, dt, &
nd, drv_idx, drv_amp, drv_freq, drv_phi, &
drv_eq, drv_ncycles, drv_has_period, freeze)
integer, intent(in) :: n, nd, step
real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
real(c_double), intent(in) :: t, dt
integer, intent(in) :: drv_idx(nd), drv_has_period(nd)
real(c_double), intent(in) :: drv_amp(3,nd), drv_freq(3,nd)
real(c_double), intent(in) :: drv_phi(3,nd), drv_eq(3,nd)
real(c_double), intent(in) :: drv_ncycles(nd)
real(c_double), intent(inout) :: freeze(3,nd)
integer :: d, idx, ps
real(c_double) :: fx, fy, fz, mf, px, py, pz
do d = 1, nd
idx = drv_idx(d) + 1 ! 0-based → 1-based
fx = drv_freq(1,d); fy = drv_freq(2,d); fz = drv_freq(3,d)
if (drv_has_period(d) /= 0) then
mf = max(abs(fx), max(abs(fy), abs(fz)))
ps = 0
if (mf > 1.0d-12) ps = int(drv_ncycles(d)/mf/dt)
if (step > ps) then
x(idx)=freeze(1,d); y(idx)=freeze(2,d); z(idx)=freeze(3,d)
vx(idx)=0.0d0; vy(idx)=0.0d0; vz(idx)=0.0d0
cycle
end if
px = drv_eq(1,d)+drv_amp(1,d)*cos(TWO_PI*fx*t+drv_phi(1,d))
py = drv_eq(2,d)+drv_amp(2,d)*cos(TWO_PI*fy*t+drv_phi(2,d))
pz = drv_eq(3,d)+drv_amp(3,d)*cos(TWO_PI*fz*t+drv_phi(3,d))
if (step == ps) then
freeze(1,d)=px; freeze(2,d)=py; freeze(3,d)=pz
end if
end if
x(idx) = drv_eq(1,d)+drv_amp(1,d)*cos(TWO_PI*fx*t+drv_phi(1,d))
y(idx) = drv_eq(2,d)+drv_amp(2,d)*cos(TWO_PI*fy*t+drv_phi(2,d))
z(idx) = drv_eq(3,d)+drv_amp(3,d)*cos(TWO_PI*fz*t+drv_phi(3,d))
vx(idx) = -drv_amp(1,d)*TWO_PI*fx*sin(TWO_PI*fx*t+drv_phi(1,d))
vy(idx) = -drv_amp(2,d)*TWO_PI*fy*sin(TWO_PI*fy*t+drv_phi(2,d))
vz(idx) = -drv_amp(3,d)*TWO_PI*fz*sin(TWO_PI*fz*t+drv_phi(3,d))
end do
end subroutine
! ══════════════════════════════════════════════════════════════
! 导出函数：run_dynamics（C 兼容接口，bind(C)）
! 接口与 C/C++ DLL 完全相同（扁平 C-contiguous 数组）。
! ══════════════════════════════════════════════════════════════
integer(c_int) function run_dynamics( &
n_atoms, pos_init, vel_init, masses, fixed, &
n_bonds, bond_pairs, bond_k, bond_r0, &
box_a, dt, NT, NSTEP, warmup_steps, method_id, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, gravity_strength, &
n_drivers, drv_idx, drv_amp, drv_freq, drv_phi, drv_eq, &
drv_ncycles, drv_has_period, &
n_frames, out_x, out_y, out_z, out_vx, out_vy, out_vz, &
progress_cb) &
bind(C, name="run_dynamics")
integer(c_int), value, intent(in) :: n_atoms, n_bonds, NT, NSTEP
integer(c_int), value, intent(in) :: warmup_steps, method_id
integer(c_int), value, intent(in) :: gravity_field, elastic_force, damping_force
integer(c_int), value, intent(in) :: n_drivers, n_frames
real(c_double), value, intent(in) :: box_a, dt
real(c_double), value, intent(in) :: Gx, Gy, Gz, Bx, By, Bz
real(c_double), value, intent(in) :: gravity_strength
! 扁平数组：Python 传入 C-contiguous int32/float64
! Fortran 以列优先解释，维度反转：(3,n) 对应 C 的 n×3
real(c_double), intent(in) :: pos_init(3, n_atoms)
real(c_double), intent(in) :: vel_init(3, n_atoms)
real(c_double), intent(in) :: masses(n_atoms)
integer(c_int), intent(in) :: fixed(3, n_atoms)
integer(c_int), intent(in) :: bond_pairs(2, n_bonds)
real(c_double), intent(in) :: bond_k(n_bonds), bond_r0(n_bonds)
integer(c_int), intent(in) :: drv_idx(n_drivers)
real(c_double), intent(in) :: drv_amp(3, n_drivers)
real(c_double), intent(in) :: drv_freq(3, n_drivers)
real(c_double), intent(in) :: drv_phi(3, n_drivers)
real(c_double), intent(in) :: drv_eq(3, n_drivers)
real(c_double), intent(in) :: drv_ncycles(n_drivers)
integer(c_int), intent(in) :: drv_has_period(n_drivers)
real(c_double), intent(out) :: out_x(n_atoms, n_frames)
real(c_double), intent(out) :: out_y(n_atoms, n_frames)
real(c_double), intent(out) :: out_z(n_atoms, n_frames)
real(c_double), intent(out) :: out_vx(n_atoms, n_frames)
real(c_double), intent(out) :: out_vy(n_atoms, n_frames)
real(c_double), intent(out) :: out_vz(n_atoms, n_frames)
type(c_funptr), value, intent(in) :: progress_cb
! 进度回调接口
abstract interface
subroutine cb_iface(step, total) bind(C)
use iso_c_binding
integer(c_int), value :: step, total
end subroutine
end interface
procedure(cb_iface), pointer :: cb_ptr
integer :: n, s, frame_idx, record_steps, prog_interval, nd
real(c_double) :: t, tw
real(c_double), allocatable :: x(:), y(:), z(:), vx(:), vy(:), vz(:)
real(c_double), allocatable :: ax0(:), ay0(:), az0(:)
real(c_double), allocatable :: freeze(:,:)
logical :: has_cb
n = n_atoms
nd = n_drivers
allocate(x(n), y(n), z(n), vx(n), vy(n), vz(n))
do s = 1, n
x(s) = pos_init(1,s); y(s) = pos_init(2,s); z(s) = pos_init(3,s)
vx(s) = vel_init(1,s); vy(s) = vel_init(2,s); vz(s) = vel_init(3,s)
end do
allocate(freeze(3, max(nd,1)))
freeze = 0.0d0
has_cb = c_associated(progress_cb)
if (has_cb) call c_f_procpointer(progress_cb, cb_ptr)
! ── 蛙跳法：初始化 v(-dt/2) ─────────────────────────────
if (method_id == 3) then
allocate(ax0(n), ay0(n), az0(n))
call accel_conservative(n, x, y, z, masses, Gx, Gy, Gz, &
gravity_field, elastic_force, &
n_bonds, bond_pairs, bond_k, bond_r0, ax0, ay0, az0)
do s = 1, n
if (fixed(1,s)/=0 .and. fixed(2,s)/=0 .and. fixed(3,s)/=0) cycle
vx(s)=vx(s)-0.5d0*ax0(s)*dt
vy(s)=vy(s)-0.5d0*ay0(s)*dt
vz(s)=vz(s)-0.5d0*az0(s)*dt
end do
deallocate(ax0, ay0, az0)
end if
! ── 初始驱动 t=0 ─────────────────────────────────────────
if (nd > 0) call apply_drive(n, x, y, z, vx, vy, vz, 0.0d0, 0, dt, &
nd, drv_idx, drv_amp, drv_freq, drv_phi, &
drv_eq, drv_ncycles, drv_has_period, freeze)
! ── 预热 ─────────────────────────────────────────────────
do s = 0, warmup_steps-1
tw = (s+1)*dt
if (nd>0) call apply_drive(n, x, y, z, vx, vy, vz, tw, s, dt, &
nd, drv_idx, drv_amp, drv_freq, drv_phi, &
drv_eq, drv_ncycles, drv_has_period, freeze)
call do_step(n, x, y, z, vx, vy, vz, masses, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt, method_id, &
pos_init, box_a)
end do
! ── 记录循环 ─────────────────────────────────────────────
record_steps = NT - warmup_steps
prog_interval = max(1, record_steps/100)
frame_idx = 0
do s = 0, record_steps-1
if (has_cb .and. mod(s, prog_interval)==0 .and. s>0) call cb_ptr(s, record_steps)
t = (s+warmup_steps)*dt
if (nd>0) call apply_drive(n, x, y, z, vx, vy, vz, t, s, dt, &
nd, drv_idx, drv_amp, drv_freq, drv_phi, &
drv_eq, drv_ncycles, drv_has_period, freeze)
if (mod(s, NSTEP)==0 .and. frame_idx<n_frames) then
frame_idx = frame_idx+1
out_x(:, frame_idx) = x
out_y(:, frame_idx) = y
out_z(:, frame_idx) = z
out_vx(:, frame_idx) = vx
out_vy(:, frame_idx) = vy
out_vz(:, frame_idx) = vz
end if
call do_step(n, x, y, z, vx, vy, vz, masses, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt, method_id, &
pos_init, box_a)
end do
deallocate(x, y, z, vx, vy, vz, freeze)
run_dynamics = 0
contains
subroutine do_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt, method_id, &
pos_init, box_a)
integer, intent(in) :: n, gravity_field, elastic_force, damping_force
integer, intent(in) :: n_bonds, method_id
real(c_double), intent(inout) :: x(n), y(n), z(n), vx(n), vy(n), vz(n)
real(c_double), intent(in) :: m(n), bond_k(n_bonds), bond_r0(n_bonds)
integer, intent(in) :: fixed(3,n), bond_pairs(2,n_bonds)
real(c_double), intent(in) :: Gx, Gy, Gz, Bx, By, Bz, dt, box_a
real(c_double), intent(in) :: pos_init(3, n)
select case (method_id)
case (0)
call euler_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt)
case (1)
call implicit_euler_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt)
case (2)
call midpoint_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt)
case default
call leapfrog_step(n, x, y, z, vx, vy, vz, m, fixed, &
Gx, Gy, Gz, Bx, By, Bz, &
gravity_field, elastic_force, damping_force, &
n_bonds, bond_pairs, bond_k, bond_r0, dt)
end select
call apply_bc(n, x, y, z, vx, vy, vz, fixed, pos_init, box_a)
end subroutine do_step
end function run_dynamics
end module dynamics_dll
engines/src/fortran/main.f90
+1071
View File
File diff suppressed because it is too large Load Diff