! 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)hi) then; y(i)=hi; vy(i)=-abs(vy(i)); end if if (y(i)hi) then; z(i)=hi; vz(i)=-abs(vz(i)); end if if (z(i)hi) x(i)=lo; if (x(i)hi) y(i)=lo; if (y(i)hi) z(i)=lo; if (z(i) 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