如何实现扩散系数的计算使用周期性边界条件（PBC）？
我使用一组分子动力学有位置的粒子在一个长度的盒子里. 中国人民银行的实施是
do k = 1, 3
pos(k,i) = modulo(pos(k,i),length)
end do
现在我正在使用以下代码
do it=1,nstep
diff=0
do i = 1,nbody
pos2(:)=pos2(:)+pos(:,i)
diff=diff+dot_product(pos(:,i),pos(:,i))
end do
diff=diff/nbody-dot_product(pos2(:),pos2(:))/nbody**2
end do
diff=diff/nstep/6
我认为它对应于
但我不太确定 PBC 是否以正确的方式考虑在内。
有人能帮我吗？
谢谢马特奥
模拟完成后，您需要展开位置。我使用以下子程序：
subroutine unfold_positions(L, X, Xu)
! Unwinds periodic positions. It is using the positions of particles from the
! first time step as a reference and then tracks them as they evolve possibly
! outside of this box. The X and Xu arrays are of the type X(:, j, i), which
! are the (x,y,z) coordinates of the j-th particle in the i-th time step.
real(dp), intent(in) :: L ! Box length
! Positions in [0, L]^3 with possible jumps:
real(dp), intent(in) :: X(:, :, :)
! Unwinded positions in (-oo, oo)^3 with no jumps (continuous):
real(dp), intent(out) :: Xu(:, :, :)
real(dp) :: d(3), Xj(3)
integer :: i, j
Xu(:, :, 1) = X(:, :, 1)
do i = 2, size(X, 3)
do j = 1, size(X, 2)
Xj = X(:, j, i-1) - X(:, j, i) + [L/2, L/2, L/2]
Xj = Xj - L*floor(Xj/L)
d = [L/2, L/2, L/2] - Xj
Xu(:, j, i) = Xu(:, j, i-1) + d
end do
end do
end subroutine
Xj = Xj - L*floor(Xj/L)我认为你可以用, 好点来替换行Xj = modulo(Xj, L)，我没有意识到你可以使用modulo浮点数。
我解决了我自己的问题。正如我在问题中所说，我正在使用 PBC：
do k = 1, 3
pos(k,i) = modulo(pos(k,i),length)
end do
并且粒子存在动力学
pos(k,i)=pos(k,i)+vel(k,i)*dt+0.5*force(k,i)*dt**2
为了计算扩散系数，我使用了第二个位置变量，比如说。这个新变量像原来的变量一样演变：
pos0(k,i)=pos0(k,i)+vel(k,i)*dt+0.5*force(k,i)*dt**2
但对他们我不适用PBC。因此，扩散系数的计算类似于：
do it=1,nstep
diff=0
do i = 1,nbody
pos2(:)=pos2(:)+pos0(:,i)
diff=diff+dot_product(pos0(:,i),pos0(:,i))
end do
diff=diff/nbody-dot_product(pos2(:),pos2(:))/nbody**2
end do
diff=diff/nstep/6
相反，其他数量是使用原始变量计算的，PBC 对此持有：
do k = 1, 3
pos(k,i) = modulo(pos(k,i),length)
end do
我希望这个答案对正在寻找类似东西的人有所帮助。