# Numerical integration in Fourier space over 3D grid

I am attempting to implement a model outlined in this paper:

Background

This model allows the calculation of magnetostatic interaction energies between objects of arbitrary shape. In the model, an object is defined as a three-dimensional array $$D(\pmb{r})$$ which is equal to zero where there is not a particle and one where there is, this is called the "shape function". The calculation occurs in Fourier space, the Fourier transform of the shape function is the shape amplitude, $$D(\pmb{k})$$, a three-dimensional grid with associated three other three-dimensional grids $$k_x$$, $$k_y$$, and $$k_z$$ which assign frequencies in the $$x$$, $$y$$, and $$z$$ directions at every grid point in $$D(\pmb{k})$$.

With this background, I want to evaluate (4) from the paper:

$$E_m=\dfrac{\overline{K}_d}{4\pi^3}Re\left[\int d^3\pmb{k} D_1(\pmb{k})D_1^*(\pmb{k})\times (\pmb{\hat{m}_{1}}\cdot \pmb{\hat{k}}) (\pmb{\hat{m}_{2}}\cdot \pmb{\hat{k}})e^{i\pmb{k}\cdot\pmb{\rho}}\right]$$

where the variable $$\pmb{\rho}$$ is just the displacement of the two objects being considered, $$\pmb{\rho}=[\Delta x, \Delta y, \Delta z]$$, $$\pmb{\hat{m}_{i}}$$ is the magnetization unit vector of the object, and $$\pmb{\hat{k}}$$ is the unit frequency vector. $$Re$$ refers to the real part being taken only. This equation has a nondimensionalized, discrete representation:

$$\tilde{E}_m=\dfrac{1}{L^3} \sum_{\pmb{k_z}} \sum_{\pmb{k_y}} \sum_{\pmb{k_x}} \left[ \dfrac{D_i(\pmb{k})D_{j}^*(\pmb{k})}{V_i V_j} \right] \left[ m_i^\alpha \dfrac{k^\alpha k^\beta}{k^2} m_j^\beta \right] e^{i (\pmb{k}\cdot\pmb{\rho})}$$

Where $$m_{i/j}^{\alpha/\beta}$$ refers to the components of the unit magnetization vectors (denoted by $$\alpha$$ or $$\beta$$) of the particle (denoted by $$i$$ and $$j$$). $$L$$ is the prefactor of integration across each direction and is equal to $$N\delta$$ where $$N$$ is the size of one dimension of the grid (assumed to be a cube) and $$\delta$$ is the size of each grid point. Summation across all $$\alpha/\beta$$ is implied.

If particles are uniformly magnetized along just one coordinate axis, say the z-axis and the particles are identical, the sum simplifies:

$$\tilde{E}_m=\dfrac{1}{L^3} \sum_{\pmb{k_z}} \sum_{\pmb{k_y}} \sum_{\pmb{k_x}} \dfrac{|D(\pmb{k})|^2}{V^2} \dfrac{k_z^2}{k^2} e^{i (\pmb{k}\cdot\pmb{\rho})}$$

This is the sum I am evaluating. I compare my result to a simple nondimensionalized analytical result for spheres magnetized in the z-axis and separated in the x-axis given in the paper in (27):

$$\tilde{E}(\rho_x)=\dfrac{1}{2\pi \rho_x^3}$$

Using length-units of particle radii, I wrote a Python script to execute the integration.

Code


import numpy as np

#Create the shape function of a circle

gridSize = 128 #Speed up FFT algo with a power of 2 size
D_r = np.zeros([gridSize,gridSize,gridSize]) #Preallocate grid

radius = 64 #Define sphere size

for xind in range(gridSize):
x = xind - (gridSize-1)/2

for yind in range(gridSize):
y = yind - (gridSize-1)/2

for zind in range(gridSize):
z = zind - (gridSize-1)/2
dist =  np.sqrt((x)**2+(y)**2+(z)**2)

D_r[xind,yind,zind] = 1

"Define the frequency vectors"

freqs = np.fft.fftshift(np.fft.fftfreq(gridSize,d=unitsPerPixel))
#Stretch these out into grids to define the x, y, and z frequencies at every point of the grid D(k)
kx, ky, kz = np.meshgrid(freqs, freqs, freqs, indexing='ij')

"Define the shape amplitudes"
D_k = np.fft.fftshift(np.fft.fftn(D_r)) #Particle 1
D_k2 = D_k #Particle 2, identical in this case

"Define the magnetization and separation vectors"
m1 = np.array([0,0,1]) #particle 1 magnetization
m2 = np.array([0,0,1]) #particle 2 magnetization
rho = np.array([11,0,0]) #paricle separaion, must be >= 2

"Calculate the analytical solution for refeence"
V = 4/3*np.pi*R**3 #Sphere volume
Ea = 2/(4*np.pi*rho[0]**3) #Dimensionless interaction energy

"Computer the numerical solution"
Etest = 0

for xind in range(gridSize):
for yind in range(gridSize):
for zind in range(gridSize):
#Build the frequency vector
k = np.array([kx[xind,yind,zind],ky[xind,yind,zind],kz[xind,yind,zind]])

#Check if this is the zero point
if k[0] == 0 and k[1]== 0 and k[2]== 0:
continue
else:

Etest += 1/V**2*D_k[xind,yind,zind]*np.conjugate(D_k2[xind,yind,zind]) \
*k[2]*k[2]/(k[0]**2+k[1]**2+k[2]**2) \
*np.exp(1j*np.dot(rho,k))

L = gridSize*unitsPerPixel
print('The numerical solution is:')
print(np.real(Etest * 1 / L**3))
print('The analytical solution is:')
print(Ea)


Questions

Something is clearly going seriously wrong. The two problems are that my output is far too large and my result isn't actually a strong function of the object separation, $$\rho$$.

1. For the order of magnitude issue, I think this may stem from my definition of $$L$$. But I am directly using the author's definition. Is my definition of $$L$$ correct?

2. For the $$\rho$$ dependence issue, I worry this stems from my definition of the frequency vectors using NumPy's fftfreq function as $$\rho$$ only shows up in the exponential Fourier shift factor. Does this sound right? Is my implementation of the frequency grids correct?

• I might be wrong, but is the integral not simply an inverse Fourier transform? You never call the inverse FFT. That said, you don't need the last three loops + the sum. Commented Jun 14 at 8:22