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

fix: C/C++/Fortran 引擎补齐 wrap_position 和 t=0 驱动力

Browse Source 
- 三引擎均新增 wrap_position 周期边界回绕（调用在 limit_in_box 后）
- 三引擎均新增 t=0 驱动力初始调用（在预热循环前）
- 至此三引擎算法与 Python 完全一致
This commit is contained in:
Ying-Li Niu
2026-06-11 18:57:41 +08:00
parent 2ce8ded482
commit b9ec622808
3 changed files with 59 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
engines/c/main.c
+16
View File
@@ -507,6 +507,12 @@ static void limit_in_box(double *pos, double *vel, double lo, double hi) {
if (*pos < lo) { *pos = lo; *vel = -*vel; }
}
/* 周期边界回绕：超出边界的位置绕到对侧（与 Python wrap_position 一致）*/
static void wrap_position(double *pos, double lo, double hi) {
if (*pos > hi) *pos = lo;
if (*pos < lo) *pos = hi;
}
/* ── 显式欧拉法 ──────────── */
static void explicit_euler_step(
int n, double *x, double *y, double *z,
@@ -752,6 +758,13 @@ static void apply_step(
limit_in_box(&z[i], &vz[i], -box_a, box_a);
}
/* 周期边界回绕（与 Python wrap_position 一致）*/
for (int i = 0; i < n; i++) {
wrap_position(&x[i], -box_a, box_a);
wrap_position(&y[i], -box_a, box_a);
wrap_position(&z[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; }
@@ -899,6 +912,9 @@ int main(int argc, char **argv) {
traj.vz = traj.vy + record_steps * n;
/* 预热 */
/* 初始时刻 t=0 驱动力（与 Python run_simulation 一致）*/
if (params.driving_force) apply_driving_force(n, x, y, z, vx, vy, vz, 0.0, 0, params.DT, &drivers);
for (int s = 0; s < params.warmup_steps; s++) {
double tw = (s + 1) * params.DT;
if (params.driving_force) apply_driving_force(n, x, y, z, vx, vy, vz, tw, s, params.DT, &drivers);
engines/cpp/main.cpp
+17
View File
@@ -386,6 +386,12 @@ static void limit_in_box(double &pos, double &vel, double lo, double hi) {
if (pos < lo) { pos = lo; vel = -vel; }
}
/* 周期边界回绕（与 Python wrap_position 一致）*/
static void wrap_position(double &pos, double lo, double hi) {
if (pos > hi) pos = lo;
if (pos < lo) pos = hi;
}
// ========================================================================
// 四种积分方法（只做位置/速度更新，不含边界条件）
// 与 Python: Explicit_Euler_Method / Implicit_Euler_Method /
@@ -598,6 +604,13 @@ static void apply_step(
limit_in_box(z[i], vz[i], -box_a, box_a);
}
// 周期边界回绕（与 Python wrap_position 一致）
for (int i = 0; i < n; i++) {
wrap_position(x[i], -box_a, box_a);
wrap_position(y[i], -box_a, box_a);
wrap_position(z[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]; vx[i] = 0.0; }
@@ -811,6 +824,10 @@ int main(int argc, char **argv) {
std::vector<double> traj_vz(record_steps * n);
// 预热
// 初始时刻 t=0 驱动力（与 Python run_simulation 一致）
if (params.driving_force)
apply_driving_force(n, x.data(), y.data(), z.data(), vx.data(), vy.data(), vz.data(), 0.0, 0, params.DT, drivers);
for (int s = 0; s < params.warmup_steps; s++) {
double tw = (s + 1) * params.DT;
if (params.driving_force)
engines/fortran/main.f90
+26
View File
@@ -108,6 +108,17 @@ program dynamics_f90
allocate(traj_vx(record_steps, n), traj_vy(record_steps, n), traj_vz(record_steps, n))
! 预热
! 初始时刻 t=0 驱动力（与 Python run_simulation 一致）
if (driving_force /= 0 .and. n_drivers > 0) then
call apply_driving(n, x, y, z, vx, vy, vz, 0.0d0, 0, DT, &
n_drivers, drv_atom_idx, &
drv_amp_x, drv_amp_y, drv_amp_z, &
drv_freq_x, drv_freq_y, drv_freq_z, &
drv_phi_x, drv_phi_y, drv_phi_z, &
drv_has_period, drv_period_cycles, &
drv_freeze_x, drv_freeze_y, drv_freeze_z)
end if
do s = 1, warmup_steps
if (driving_force /= 0 .and. n_drivers > 0) then
tw = (s * 1.0d0) * DT
@@ -512,6 +523,14 @@ subroutine limit_in_box(pos, vel, lo, hi)
end if
end subroutine limit_in_box
! 周期边界回绕（与 Python wrap_position 一致）
subroutine wrap_position(pos, lo, hi)
double precision, intent(inout) :: pos
double precision, intent(in) :: lo, hi
if (pos > hi) pos = lo
if (pos < lo) pos = hi
end subroutine wrap_position
! ── 显式欧拉法 ────────────
subroutine explicit_euler_step(n, x, y, z, vx, vy, vz, m, G, B, &
nb, bp, bk, br, fixed, dt, &
@@ -730,6 +749,13 @@ subroutine apply_step(method, n, x, y, z, vx, vy, vz, m, G, B, &
call limit_in_box(z(i), vz(i), -box_a, box_a)
end do
! 周期边界回绕（与 Python wrap_position 一致）
do i = 1, n
call wrap_position(x(i), -box_a, box_a)
call wrap_position(y(i), -box_a, box_a)
call wrap_position(z(i), -box_a, box_a)
end do
! 逐自由度固定约束（与 Python apply_fixed_constraints 一致）
do i = 1, n
if (fixed(i,1) /= 0) then; x(i) = pos_0(i, 1); vx(i) = 0.0d0; end if