如何在Fortran每个时间步定义新函数以降低相互作用计算量

问题解答 #

首先明确两个核心结论：

1. 你设想的封装函数的方式无法将单时间步计算量从N²降到2N
   你现在的routine_interaction子程序内部本身就包含了遍历所有N个粒子的循环，不管你是直接调用子程序，还是把它封装为你说的new_function，在3号do循环的N次迭代里，每调用一次这个函数/子程序都要执行N次内部计算，总计算量依然是N*N=N²，封装本身不会减少实际执行的运算次数。

2. Fortran语法上支持你把现有的子程序改写为返回受力值的函数，但仅改变写法不改变计算量
   你可以把routine_interaction改写为如下函数形式：

function interaction_force(ii, N, positions0) result(force)
    integer, intent(in) :: ii, N
    real, intent(in) :: positions0(N)
    real :: force, r, s
    integer :: i4
    force = 0.0
    do i4 = 1, N
        if (ii /= i4) then
            r = abs(positions0(ii) - positions0(i4))
            s = sign(1.0, positions0(i4) - positions0(ii))
            force = force - (1.0/r**4 - 1.0/r**2)*s
        end if
    end do
end function interaction_force

之后在3号循环里可以直接用force = interaction_force(i3, N, positions0)调用，写法更简洁，但计算量和原来完全一致。

可选优化方案 #

如果要降低计算量，可以利用相互作用力的牛顿第三定律对称性优化：两个粒子之间的作用力大小相等、方向相反，不需要每个粒子都单独遍历所有其他粒子算一遍，只需要遍历所有i<j的粒子对，一次计算同时给两个粒子赋值受力即可，总计算量可以降到N*(N-1)/2，接近原来的一半，优化后的时间步逻辑参考：

! 先定义存储所有粒子受力的数组
real :: force_arr(N)
! 时间步内先算所有粒子对的作用力
force_arr(:) = 0.0
do i = 1, N-1
    do j = i+1, N
        r = abs(positions0(i) - positions0(j))
        s = sign(1.0, positions0(j) - positions0(i))
        f = -(1.0/r**4 - 1.0/r**2)*s
        force_arr(i) = force_arr(i) + f
        force_arr(j) = force_arr(j) - f
    end do
end do
! 再循环更新每个粒子的位置
do 3 i3=1,N
    param=param+100
    noise=(ran2(param)-0.5)*L/10
    positions(i2,i3)=positions0(i3)+dt*(noise+force_arr(i3))
3 continue

注意这个优化还是属于O(N²)复杂度的范畴，你设想的O(N)复杂度的全对全相互作用计算目前不存在可行方案。

内容的提问来源于stack exchange，提问作者J.A