Overflow errors in simulating systems with Lennard-Jones and magnetostatics for polydispersed particles

[ywem]

Dear users,
I am trying to simulating polydispersed particles following the tutorial of "Ferrofluid". When I modified the program to incorporate polydisperse particle sizes, the following error occurred during the execution of tune.skin.

RuntimeError: Overflow in the image box count while folding a particle coordinate into the primary simulation box. Maybe a particle experienced a huge force.

In this simulation, I segregate the particles into 8 groups, with each group giving a dipolar moment using the program below.

N_classes = 8
dp = [1.0, 1.2, 1.6, 2.0, 2.5, 3.1, 3.9, 4.9]
Nd_classes = [1205, 2307, 1037, 327, 90, 24, 6, 2]   #Number of each type of diameter
m = math.sqrt(KT)*np.sqrt(Lambda)/np.sqrt(LJ_EPSILON)
np.random.seed(seed=1)
dip_phi = 2. * np.pi * np.random.random((N_PART, 1))
dip_cos_theta = 2 * np.random.random((N_PART, 1)) - 1
dip_sin_theta = np.sin(np.arccos(dip_cos_theta))
dip = np.hstack((
    dip_sin_theta * np.sin(dip_phi),
    dip_sin_theta * np.cos(dip_phi),
    dip_cos_theta))
Nd = 0
for i in range(N_classes):
    for j in range(Nd_classes[i]):
        dip[j+Nd,0:2] = m[i]*dp[i]**(3/2)*dip[j+Nd,0:2]
    Nd = Nd + Nd_classes[i]

And defined their interactions by employing the following approach

for i in range(N_classes):
    for j in range(N_classes):
        LJ_SIGMA[i,j] = dp[i] + dp[j]
        LJ_CUT[i,j] = 2**(-5./6.) * (dp[i]+dp[j])
for i in range(N_classes):
    for j in range(i,N_classes):
        system.non_bonded_inter[i, j].lennard_jones.set_params(epsilon=LJ_EPSILON, sigma=LJ_SIGMA[i,j], cutoff=LJ_CUT[i,j], shift="auto")

I have also attempted to segregate the particles into two groups, with each group sharing identical parameters. However, the aforementioned problem persists.
I wonder if any of the users have encountered a similar issue.

[jngrad]

Hi,

The error message is relatively transparent. Two particles overlapped with each other and went ballistic, and the image box counter has become undefined. I can only speculate on the origin of the large forces, since your MWE isn't entirely self-contained. Maybe the time step was not adjusted according to the new epsilon values; with $\sigma = 9.8$ the force will be $\sigma^{12} \approx 7\cdot 10^{11}$ times larger than in the original tutorial, thus the time step must get 12 orders of magnitude smaller. It could also be that you adjusted the time step but didn't adjust the steepest descent algorithm's f_max parameter, whose value is a function of the time step (the charged matter tutorial shows its expression), or max_displacement parameter, whose value is a function of sigma (often 1% of the largest sigma).

[ywem]

Dear jngrad,

I have revisited the program and confirmed that fmax was indeed not set correctly. Now, the program is functioning properly. I sincerely appreciate your response.