我是一个有用的助手,可以翻译文本。
这里是最终答案和一些时间记录。
0) 普通版本(非常慢)
In [16]: %timeit total_energy(points)
1 loops, best of 3: 14.9 s per loop
1)SciPy版本
In [9]: %timeit scipy_LJ(points)
10 loops, best of 3: 44 ms per loop
1-2) Numpy 版本
%timeit sum( potential(np.sqrt((x[i]-x[:i])**2 + (y[i]-y[:i])**2 + (z[i] - z[:i])**2)).sum() for i in range(N-1))
10 loops, best of 3: 126 ms per loop
2)极快的Fortran版本(!表示注释)
subroutine EnergyForces(Pos, PEnergy, Dim, NAtom)
implicit none
integer, intent(in) :: Dim, NAtom
real(8), intent(in), dimension(0:NAtom-1, 0:Dim-1) :: Pos
real(8), intent(out) :: PEnergy
real(8), dimension(Dim) :: rij, Posi
real(8) :: d2, id2, id6, id12
real(8) :: rc2, Shift
integer :: i, j
PEnergy = 0.
do i = 0, NAtom - 1
Posi = Pos(i,:)
do j = i + 1, NAtom - 1
rij = Pos(j,:) - Posi
d2 = sum(rij * rij)
id2 = 1. / d2
id6 = id2 * id2 * id2
id12 = id6 * id6
PEnergy = PEnergy + 4. * (id12 - id6)
enddo
enddo
end subroutine
调用之后
In [14]: %timeit ljlib.energyforces(points.transpose(), 3, N)
10000 loops, best of 3: 61 us per loop
3) 结论 Fortran比scipy快1000倍,比numpy快3000倍,比纯Python快数百万倍。这是因为Scipy版本创建了一个差异矩阵然后分析它,而Fortran版本则一切都在运行时完成。
i
只到N-1
?在这种情况下,您将错过最后一个条目(N-1
被排除在外)。x
是否总是一维的,还是需要进行泛化处理? - eickenberg