How can we simulate the interactions among polydisperse particles using ESPResSo?

[ywem]

When dealing with a system composed of polydisperse particles, and assuming that the particle diameters adhere to a distribution such as a normal (Gaussian) distribution, how should one input these particles into the system?

[RudolfWeeber]

t present, this can only be achieved approximately by using particle types..
I.e., using some 20 "size classes", one particle type for each size class and a Lennard-Jones interaction for all pairs of size classes.
This is also, where you would decide on the mixing rule, i.e., the sigma and epsilon-parameter of the Lennard-Jones interaction between particles of size class i and size class j.

One could use np.histogram to sample the distributio into the size classes.

If this approximation is not viable for the probelm in quesiton, one would have to

• add a particle property which sotres a "radius"
• make a modified Lennard-Jones (or other) non-bonded interaction which makes use of the "radii" of the invovled particles and applies a mixing rule.

If you would like to implement that and have some C++ knowledge, I can point out the places to start form in the code.

[ywem]

Dear RudolfWeeber,
Thanks for your reply.
For the first approach, I have some additional queries. How do we ascertain the diameter of particles in different classes? And when it comes to other types of interactions, such as electrostatic interaction, how do we define and characterize these interactions between particles in different classes?

Best regards

[RudolfWeeber]

Hi,

here's a sketch (not tested, just a rough idea)

n_classes = 11 # odd number for symmetry
assert n_classes%2==1
r_min =1.0 # smallerst particle size
r_max = 3.0 # largest particle size
rs = np.linspace(r_min, r_max, n_classes) 

# Mixing rule (this is a physics decisin)
# for example
sigma = lambda i,j : rs[i] + rs[j]
epsilon = lambda i,j: 1.0 # but can depend on size class, too
# and similarly for the cutoff.


for i in range(n_classes):
  for j in range(i,n_classes):
    system.lennard_jones.non_bonded_inter[i,j].set_params(
      epsilon=epsilon(i,j), sigma=sigma(i,j), cutoff=....) 

For electrostatics, you can set the charge for each particle individually.

Hope that helps.
Rudolf

[ywem]

Dear Rudolf,
Thank you very much for your detailed explanation. I'll give it a try.

Best wishes,
Ywem