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

General magnetostatic shape–shape interactions


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.


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)
            if dist < radius:
                    D_r[xind,yind,zind] = 1

"Define the frequency vectors"
unitsPerPixel = 1/radius #Rescale into units of particle radii

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"
R= 1 #Sphere radius
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:

                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) \

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


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?

  • $\begingroup$ 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. $\endgroup$
    – ConvexHull
    Commented Jun 14 at 8:22


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.