Let us have a force F as a long distance attraction force: \[f^i(r) = -\frac{r^i}{|r|}\phi(|r|)\] Where \(r\) is the vector between the atracted molecules, \(i\) is x or y, and \(\phi\) is the dependence on the distance.

Now, let us have a lattice L with a density field \(\rho\), we have a summed force: \[F([x\ y]) = \sum_{k,l} f([x\ y]-[k\ l])\rho([k\ l])\] Let us now say that: \[f_{a,b} = f([a\ b])\] \[F_{a,b} = F([a\ b])\] \[\rho_{a,b} = \rho([a\ b])\] We have: \[F_{x,y} = \sum_{k,l} f_{x-k,y-l}\rho_{k,l}\] This is a convolution. Let \(\hat f\) be the FFT of \(f\): \[\hat f_{s,t} = \sum_{a,b} e^{i\frac{sa+tb}{n}2\pi}f_{a,b}\] We have: \[\hat F_{s,t} = \sum_{a,b} e^{i\frac{sa+tb}{n}2\pi}F_{a,b} = \] \[ = \sum_{a,b} e^{i\frac{sa+tb}{n}2\pi}\sum_{k,l} f_{a-k,b-l}\rho_{k,l} = \] \[ = \sum_{a,b,k,l} e^{i\frac{sa+tb}{n}2\pi} f_{a-k,b-l}\rho_{k,l} = \] We change the addition variables with \(q=a-k\) and \(p=b-l\): \[ = \sum_{p,q,k,l} e^{i\frac{s(q+k)+t(p+l)}{n}2\pi} f_{q,p}\rho_{k,l} = \] \[ = \sum_{p,q} e^{i\frac{sq+tp}{n}2\pi}f_{q,p}\sum_{k,l} e^{i\frac{sk+tl}{n}2\pi}\rho_{k,l} = \] \[ = \hat f_{s,t}\hat\rho_{s,t}\] This means that the FFT of the total force field \(F\) is the point-wise multiplication of FFT of the density and the FFT of the distance dependent f.

Example:

Let us have a lattice 64x64 and rho field:

Lx = 64
Ly = 64
rho = matrix(0.0,Lx,Ly)
rho[1:64,9:56] = 0.5
image(rho,asp=1,main="Density field")

Let us now prepare the calculation of f:

phi = function(r) {
  ifelse(r > 6, 0, -24*0.23*(2*(5/r)^13-(5/r)^7))
}
x0 = c(0:(Lx/2),(-Lx/2+1):(-1)) # "stencil" in 1D
p = expand.grid(x=x0,y=x0) # "stencil" in 2D
r = sqrt(rowSums(p^2)) # calculating distance r
f = -p/r * phi(r)
f[r < 1e-10,] = 0
plot(p,type="n",asp=1,xlim=c(-5,5),ylim=c(-5,5),main="Kernel force f")
arrows(p[,1],p[,2],(p+f)[,1],(p+f)[,2],angle = 15,length = 0.1)

fx = matrix(f[,1],Lx,Ly) #transforming f into lattice fields
fy = matrix(f[,2],Lx,Ly)

Now we can calculate the total force F with FFT:

hat_rho = fft(rho)
hat_fx  = fft(fx)
hat_fy = fft(fy)
hat_Fx  = hat_rho * hat_fx
hat_Fy  = hat_rho * hat_fy
Fx = fft(hat_Fx, inverse = TRUE)/length(hat_Fx)
Fy = fft(hat_Fy, inverse = TRUE)/length(hat_Fx)

Let’s see F:

image(rho,asp=1,main="Density")

image(Re(Fx),asp=1,main="Force X")

image(Re(Fy),asp=1,main="Force Y")

Fx0 = matrix(0,Lx,Ly)
Fy0 = matrix(0,Lx,Ly)
for (x1 in 1:64) {
  for (y1 in 1:64) {
    for (x2 in 1:64) {
      for (y2 in 1:64) {
        v = c(x1-x2,y1-y2)
        r = sqrt(sum(v^2))
        f0 = - v/r * phi(r)
        if (r == 0) f0 = c(0,0)
        Fx0[x1,y1] = Fx0[x1,y1] + f0[1]*rho[x2,y2]
        Fy0[x1,y1] = Fy0[x1,y1] + f0[2]*rho[x2,y2]
      }
    }
  }
}
image(Fx0,asp=1)

image(Fy0,asp=1)

Let us compare the error:

Fx = Re(Fx)
Fy = Re(Fy)
range(Fx)
## [1] 0 0
range(Fx0)
## [1] -6844683934  6844683934
range(Fx-Fx0)
## [1] -6844683934  6844683934