Implementing a soft repulsion force for a polydisperse system

[gkerdelian]

Hello,

I'm trying to implement a soft (linear) repulsive force in a hoomd simulation of the following form:
F_ij = k(a_i + a_j -r_ij ) if r_ij < a_i + a_j
= 0 otherwise
where r_ij is the distance between the centers of two particles with radii a_i and a_j , and the force is along the relative position vector of the two particles. In other words, particles repel when overlapping and don't interact otherwise.

For a monodisperse system (all particles have the same radius a), I used hoomd.md.pair.Table with the code below.

# Defining cut-off radius
k = 100
a_i = a_j = a
r_cut = a_i + a_j 
n_points = 1000 

# Computing force and potential
r = np.linspace(1e-4, r_cut, n_points)
F = k * (r_cut - r)
U = 0.5 * k * (r_cut - r)**2 

# Appending 0s for continuity
r = np.append(r, r_cut + 1e-4)
F = np.append(F, 0)
U = np.append(U, 0)

# Defining force and appending to integrator
nl = hoomd.md.nlist.Cell(buffer = a)
soft_repulsive_force = hoomd.md.pair.Table(nlist = nl)
soft_repulsive_force.params[('A','A')] = dict(r_min = 0, U = U, F = F)
soft_repulsive_force.r_cut[('A', 'A')] = a_i+a_j
integrator.forces.append(soft_repulsive_force)

For a system of particles with radii sampled from a uniform distribution [0.9a, 1.1a], how do I make sure that the correct cut-off radius is used every time? I have seen suggestions to define a typeid for each particle size and loop over all possible combinations of typeids and particle sizes, but I imagine that will be very inefficient for systems of thousands of particles.

Thanks in advance.

[joaander]

Unfortunately, the typing system in HOOMD-blue is the only way to achieve polydisperse r_cut values. HOOMD will issue a warning about performance, but on CPUs and modern GPUs the performance loss is manageable (see #2037 (comment)). This does result in O(N^2) memory usage for the r_cut and potential parameter matrix.

FYI, you can obtain better performance by using the compiled potential DPDConservative(https://hoomd-blue.readthedocs.io/en/v5.2.0/hoomd/md/pair/dpdconservative.html) which has the same functional form as yours except for a shift in the energy. DPDConservative will store far less memory per type pair than Table.

I am working on a new code that will allow users to define custom per particle attributes and compute potentials that depend on these attributes. That code is not ready to share open source yet, but will be soon. Watch this discussion board for an announcement (https://github.com/glotzerlab/hoomd-blue/discussions).

[gkerdelian]

Alright, thank you for the response and the tip on DPDConservative!