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

fix: driven atoms now oscillate around their initial equilibrium position

Browse Source 
Previously, apply_driving_force set absolute position to A*cos(2π f t + φ),
ignoring the atom's initial coordinates. For atoms not at the origin (e.g.,
atom 120 at x=119), this incorrectly forced them back toward the origin each
step, causing severe distortion and numerical explosion.

Fix: store each driven atom's initial position as eq_pos/eq_x/eq_y/eq_z at
load time; position is now eq + A*cos(2π f t + φ) in all four engines
(Python, C, C++, Fortran).

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Ying-Li Niu
2026-06-13 07:23:25 +08:00
parent e62e536cee
commit d489222eaf
4 changed files with 98 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files
engines/c/main.c
+12 -4
View File
@@ -84,6 +84,7 @@ typedef struct {
double *phi_x, *phi_y, *phi_z; /* radians */
int *has_period; /* 0=all, 1=limited cycles */
double *period_cycles; /* number of cycles */
double *eq_x, *eq_y, *eq_z; /* 平衡位置（初始坐标） */
double *freeze_x, *freeze_y, *freeze_z;
} DriverData;
@@ -129,12 +130,16 @@ static DriverData read_driver(const char *input_dir, const AtomData *atoms) {
d.phi_z = (double*)xmalloc(n_lines * sizeof(double));
d.has_period = (int*)xmalloc(n_lines * sizeof(int));
d.period_cycles = (double*)xmalloc(n_lines * sizeof(double));
d.eq_x = (double*)xmalloc(n_lines * sizeof(double));
d.eq_y = (double*)xmalloc(n_lines * sizeof(double));
d.eq_z = (double*)xmalloc(n_lines * sizeof(double));
d.freeze_x = (double*)xmalloc(n_lines * sizeof(double));
d.freeze_y = (double*)xmalloc(n_lines * sizeof(double));
d.freeze_z = (double*)xmalloc(n_lines * sizeof(double));
/* 初始化 freeze 数组 */
/* 初始化 freeze/eq 数组 */
for (int i = 0; i < n_lines; i++) {
d.eq_x[i] = d.eq_y[i] = d.eq_z[i] = 0.0;
d.freeze_x[i] = d.freeze_y[i] = d.freeze_z[i] = 0.0;
}
@@ -177,6 +182,9 @@ static DriverData read_driver(const char *input_dir, const AtomData *atoms) {
if (ii < 0) continue;
d.atom_idx[idx] = ii;
d.eq_x[idx] = atoms->pos_0[ii*3+0];
d.eq_y[idx] = atoms->pos_0[ii*3+1];
d.eq_z[idx] = atoms->pos_0[ii*3+2];
d.amp_x[idx] = amp_x;
d.amp_y[idx] = amp_y;
d.amp_z[idx] = amp_z;
@@ -738,9 +746,9 @@ static void apply_driving_force(
}
}
double px = drivers->amp_x[d] * cos(2*M_PI*drivers->freq_x[d]*t + drivers->phi_x[d]);
double py = drivers->amp_y[d] * cos(2*M_PI*drivers->freq_y[d]*t + drivers->phi_y[d]);
double pz = drivers->amp_z[d] * cos(2*M_PI*drivers->freq_z[d]*t + drivers->phi_z[d]);
double px = drivers->eq_x[d] + drivers->amp_x[d] * cos(2*M_PI*drivers->freq_x[d]*t + drivers->phi_x[d]);
double py = drivers->eq_y[d] + drivers->amp_y[d] * cos(2*M_PI*drivers->freq_y[d]*t + drivers->phi_y[d]);
double pz = drivers->eq_z[d] + drivers->amp_z[d] * cos(2*M_PI*drivers->freq_z[d]*t + drivers->phi_z[d]);
double vpx = -drivers->amp_x[d]*2*M_PI*drivers->freq_x[d]*sin(2*M_PI*drivers->freq_x[d]*t + drivers->phi_x[d]);
double vpy = -drivers->amp_y[d]*2*M_PI*drivers->freq_y[d]*sin(2*M_PI*drivers->freq_y[d]*t + drivers->phi_y[d]);
double vpz = -drivers->amp_z[d]*2*M_PI*drivers->freq_z[d]*sin(2*M_PI*drivers->freq_z[d]*t + drivers->phi_z[d]);
engines/cpp/main.cpp
+7 -3
View File
@@ -86,6 +86,7 @@ struct DriverData {
std::vector<double> phi_x, phi_y, phi_z; // radians
std::vector<int> has_period; // 0=all, 1=limited
std::vector<double> period_cycles;
std::vector<double> eq_x, eq_y, eq_z; // 平衡位置（初始坐标）
std::vector<double> freeze_x, freeze_y, freeze_z;
};
@@ -314,6 +315,9 @@ static DriverData read_driver(const std::string &input_dir, const AtomData &atom
continue;
}
d.atom_idx.push_back(idx);
d.eq_x.push_back(atoms.pos_0[idx*3+0]);
d.eq_y.push_back(atoms.pos_0[idx*3+1]);
d.eq_z.push_back(atoms.pos_0[idx*3+2]);
d.amp_x.push_back(ax); d.amp_y.push_back(ay); d.amp_z.push_back(az);
d.freq_x.push_back(fx); d.freq_y.push_back(fy); d.freq_z.push_back(fz);
// Convert degrees to radians
@@ -700,9 +704,9 @@ static void apply_driving_force(
}
const double TWO_PI = 2.0 * M_PI;
double px = drivers.amp_x[d] * std::cos(TWO_PI * drivers.freq_x[d] * t + drivers.phi_x[d]);
double py = drivers.amp_y[d] * std::cos(TWO_PI * drivers.freq_y[d] * t + drivers.phi_y[d]);
double pz = drivers.amp_z[d] * std::cos(TWO_PI * drivers.freq_z[d] * t + drivers.phi_z[d]);
double px = drivers.eq_x[d] + drivers.amp_x[d] * std::cos(TWO_PI * drivers.freq_x[d] * t + drivers.phi_x[d]);
double py = drivers.eq_y[d] + drivers.amp_y[d] * std::cos(TWO_PI * drivers.freq_y[d] * t + drivers.phi_y[d]);
double pz = drivers.eq_z[d] + drivers.amp_z[d] * std::cos(TWO_PI * drivers.freq_z[d] * t + drivers.phi_z[d]);
double vpx = -drivers.amp_x[d] * TWO_PI * drivers.freq_x[d] * std::sin(TWO_PI * drivers.freq_x[d] * t + drivers.phi_x[d]);
double vpy = -drivers.amp_y[d] * TWO_PI * drivers.freq_y[d] * std::sin(TWO_PI * drivers.freq_y[d] * t + drivers.phi_y[d]);
double vpz = -drivers.amp_z[d] * TWO_PI * drivers.freq_z[d] * std::sin(TWO_PI * drivers.freq_z[d] * t + drivers.phi_z[d]);
engines/fortran/main.f90
+22 -5
View File
@@ -41,6 +41,7 @@ program dynamics_f90
double precision, allocatable :: drv_phi_x(:), drv_phi_y(:), drv_phi_z(:)
integer, allocatable :: drv_has_period(:)
double precision, allocatable :: drv_period_cycles(:)
double precision, allocatable :: drv_eq_x(:), drv_eq_y(:), drv_eq_z(:)
double precision, allocatable :: drv_freeze_x(:), drv_freeze_y(:), drv_freeze_z(:)
! 运行时位置/速度
@@ -82,11 +83,12 @@ program dynamics_f90
! 读取驱动力
n_drivers = 0
if (driving_force /= 0) then
call read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
call read_driver(input_dir, n_atoms, atom_ids, pos_0, 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_eq_x, drv_eq_y, drv_eq_z, &
drv_freeze_x, drv_freeze_y, drv_freeze_z)
end if
@@ -134,6 +136,7 @@ program dynamics_f90
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_eq_x, drv_eq_y, drv_eq_z, &
drv_freeze_x, drv_freeze_y, drv_freeze_z)
end if
@@ -146,6 +149,7 @@ program dynamics_f90
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_eq_x, drv_eq_y, drv_eq_z, &
drv_freeze_x, drv_freeze_y, drv_freeze_z)
end if
call apply_step(method, n, x, y, z, vx, vy, vz, masses, G, B, &
@@ -171,6 +175,7 @@ program dynamics_f90
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_eq_x, drv_eq_y, drv_eq_z, &
drv_freeze_x, drv_freeze_y, drv_freeze_z)
end if
traj_x(s, :) = x; traj_y(s, :) = y; traj_z(s, :) = z
@@ -205,6 +210,7 @@ program dynamics_f90
deallocate(drv_freq_x, drv_freq_y, drv_freq_z)
deallocate(drv_phi_x, drv_phi_y, drv_phi_z)
deallocate(drv_has_period, drv_period_cycles)
deallocate(drv_eq_x, drv_eq_y, drv_eq_z)
deallocate(drv_freeze_x, drv_freeze_y, drv_freeze_z)
end if
@@ -804,14 +810,16 @@ end subroutine apply_step
! ========================================================================
! 读取驱动力 driver.txt
! ========================================================================
subroutine read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
subroutine read_driver(input_dir, n_atoms, atom_ids, pos_0, 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_eq_x, drv_eq_y, drv_eq_z, &
drv_freeze_x, drv_freeze_y, drv_freeze_z)
character(len=*), intent(in) :: input_dir
integer, intent(in) :: n_atoms, atom_ids(n_atoms)
double precision, intent(in) :: pos_0(n_atoms, 3)
integer, intent(out) :: n_drivers
integer, allocatable, intent(out) :: drv_atom_idx(:)
double precision, allocatable, intent(out) :: drv_amp_x(:), drv_amp_y(:), drv_amp_z(:)
@@ -819,6 +827,7 @@ subroutine read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
double precision, allocatable, intent(out) :: drv_phi_x(:), drv_phi_y(:), drv_phi_z(:)
integer, allocatable, intent(out) :: drv_has_period(:)
double precision, allocatable, intent(out) :: drv_period_cycles(:)
double precision, allocatable, intent(out) :: drv_eq_x(:), drv_eq_y(:), drv_eq_z(:)
double precision, allocatable, intent(out) :: drv_freeze_x(:), drv_freeze_y(:), drv_freeze_z(:)
character(len=512) :: path, line, period_str
@@ -831,6 +840,7 @@ subroutine read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
double precision :: pxx_tmp(MX), pyy_tmp(MX), pzz_tmp(MX)
double precision :: pc_tmp(MX)
double precision :: fzx_tmp(MX), fzy_tmp(MX), fzz_tmp2(MX)
double precision :: eqx_tmp(MX), eqy_tmp(MX), eqz_tmp(MX)
n_drivers = 0
path = trim(input_dir) // '/driver.txt'
@@ -859,6 +869,9 @@ subroutine read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
if (idx < 0) cycle
idx_tmp(n_drivers) = idx - 1 ! 0-based index
eqx_tmp(n_drivers) = pos_0(idx, 1)
eqy_tmp(n_drivers) = pos_0(idx, 2)
eqz_tmp(n_drivers) = pos_0(idx, 3)
ax_tmp(n_drivers) = ax; ay_tmp(n_drivers) = ay; az_tmp(n_drivers) = az
fxx_tmp(n_drivers) = fx; fyy_tmp(n_drivers) = fy; fzz_tmp(n_drivers) = fz
! degrees to radians
@@ -886,6 +899,7 @@ subroutine read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
allocate(drv_freq_x(n_drivers), drv_freq_y(n_drivers), drv_freq_z(n_drivers))
allocate(drv_phi_x(n_drivers), drv_phi_y(n_drivers), drv_phi_z(n_drivers))
allocate(drv_has_period(n_drivers), drv_period_cycles(n_drivers))
allocate(drv_eq_x(n_drivers), drv_eq_y(n_drivers), drv_eq_z(n_drivers))
allocate(drv_freeze_x(n_drivers), drv_freeze_y(n_drivers), drv_freeze_z(n_drivers))
drv_atom_idx = idx_tmp(1:n_drivers)
@@ -894,6 +908,7 @@ subroutine read_driver(input_dir, n_atoms, atom_ids, n_drivers, &
drv_phi_x = pxx_tmp(1:n_drivers); drv_phi_y = pyy_tmp(1:n_drivers); drv_phi_z = pzz_tmp(1:n_drivers)
drv_has_period = per_tmp(1:n_drivers)
drv_period_cycles = pc_tmp(1:n_drivers)
drv_eq_x = eqx_tmp(1:n_drivers); drv_eq_y = eqy_tmp(1:n_drivers); drv_eq_z = eqz_tmp(1:n_drivers)
drv_freeze_x = fzx_tmp(1:n_drivers); drv_freeze_y = fzy_tmp(1:n_drivers); drv_freeze_z = fzz_tmp2(1:n_drivers)
write(*, '("[Fortran-engine] 已加载驱动力: ", i0, " 条定义")') n_drivers
@@ -906,6 +921,7 @@ subroutine apply_driving(n, x, y, z, vx, vy, vz, t, step, dt, &
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_eq_x, drv_eq_y, drv_eq_z, &
drv_freeze_x, drv_freeze_y, drv_freeze_z)
integer, intent(in) :: n, step, n_drivers
integer, intent(in) :: drv_atom_idx(n_drivers), drv_has_period(n_drivers)
@@ -915,6 +931,7 @@ subroutine apply_driving(n, x, y, z, vx, vy, vz, t, step, dt, &
double precision, intent(in) :: drv_freq_x(n_drivers), drv_freq_y(n_drivers), drv_freq_z(n_drivers)
double precision, intent(in) :: drv_phi_x(n_drivers), drv_phi_y(n_drivers), drv_phi_z(n_drivers)
double precision, intent(in) :: drv_period_cycles(n_drivers)
double precision, intent(in) :: drv_eq_x(n_drivers), drv_eq_y(n_drivers), drv_eq_z(n_drivers)
double precision, intent(inout) :: drv_freeze_x(n_drivers), drv_freeze_y(n_drivers), drv_freeze_z(n_drivers)
integer :: d, idx, period_steps
@@ -941,9 +958,9 @@ subroutine apply_driving(n, x, y, z, vx, vy, vz, t, step, dt, &
end if
end if
px = drv_amp_x(d) * cos(TWO_PI * drv_freq_x(d) * t + drv_phi_x(d))
py = drv_amp_y(d) * cos(TWO_PI * drv_freq_y(d) * t + drv_phi_y(d))
pz = drv_amp_z(d) * cos(TWO_PI * drv_freq_z(d) * t + drv_phi_z(d))
px = drv_eq_x(d) + drv_amp_x(d) * cos(TWO_PI * drv_freq_x(d) * t + drv_phi_x(d))
py = drv_eq_y(d) + drv_amp_y(d) * cos(TWO_PI * drv_freq_y(d) * t + drv_phi_y(d))
pz = drv_eq_z(d) + drv_amp_z(d) * cos(TWO_PI * drv_freq_z(d) * t + drv_phi_z(d))
vpx = -drv_amp_x(d) * TWO_PI * drv_freq_x(d) * sin(TWO_PI * drv_freq_x(d) * t + drv_phi_x(d))
vpy = -drv_amp_y(d) * TWO_PI * drv_freq_y(d) * sin(TWO_PI * drv_freq_y(d) * t + drv_phi_y(d))
vpz = -drv_amp_z(d) * TWO_PI * drv_freq_z(d) * sin(TWO_PI * drv_freq_z(d) * t + drv_phi_z(d))