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}$ refer to the components of the unit magnetization vectors (denoted by $\alpha$ or $\beta$) of 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?

Thank you for reading, any help is greatly appreciated.



Your Answer

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