Simulation of polydisperse system

[edwinb-ai]

Hello everyone.

I am attempting to perform simulations of polydisperse systems. The most straightforward way to implement this, that I could think of, was to assign each particle a different diameter drawn from a continuous probability distribution, and as such each particle has assigned a unique type.

Since I am using hard-like interactions (particles interact via the WCA potential), I see that sometimes particles are touching each other, something that should not be the case. I was wondering if there is a better way of doing this, since most likely I am incorrectly assigning the interactions between particles. I provide here the code that I am using to do the simulations.

import itertools
from pathlib import Path

import gsd.hoomd
import hoomd
import numpy as np

PROJECT_DIR = Path().resolve()


def from_snapshot_to_frame(snapshot):
    """A function that takes a Snapshot from the current state of the simulation,
    and converts it to a Frame, in order to preserve data like diameters of particles."""
    frame = gsd.hoomd.Frame()
    frame.particles.N = snapshot.particles.N
    frame.particles.position = snapshot.particles.position
    frame.particles.orientation = snapshot.particles.orientation
    frame.particles.typeid = snapshot.particles.typeid
    frame.particles.types = snapshot.particles.types
    frame.particles.diameter = snapshot.particles.diameter
    frame.configuration.box = snapshot.configuration.box

    return frame


def main():
    m = 32
    total_particles = 2 * m**2
    n_particles = 1024
    rng = np.random.default_rng()
    print(f"Number of particles: {n_particles}")

    # Particles' positioning
    spacing = 1.5
    K = np.ceil(np.sqrt(total_particles)).astype(np.int64)
    L = K * spacing
    x = np.linspace(-L / 2, L / 2, K, endpoint=False)
    # Repeat twice for a 2D system
    position = np.array(list(itertools.product(x, repeat=2)))
    position = position[:n_particles]
    position = np.column_stack((position, np.zeros((n_particles,))))
    # "4" is the size of the quaternion for setting the orientation
    orientation = [(1, 0, 0, 0)] * n_particles

    frame = gsd.hoomd.Frame()
    frame.particles.N = n_particles
    frame.particles.position = position
    frame.particles.orientation = orientation

    # Diameters and types are now going to be continuous using polydispersity
    polydispersity = 0.09
    diameters = rng.normal(loc=1.0, scale=polydispersity, size=n_particles)
    frame.particles.diameter = diameters
    frame.particles.mass = diameters**2
    frame.particles.typeid = np.zeros((n_particles,))
    frame.particles.types = [f"A{i}" for i in range(1, n_particles + 1)]

    # Set the box dimensions to be 2D
    frame.configuration.box = [L, L, 0, 0, 0, 0]

    file_path = PROJECT_DIR.joinpath("lattice.gsd")
    with gsd.hoomd.open(name=file_path, mode="w") as f:
        f.append(frame)

    # Create the simulation object and set a random number
    device = hoomd.device.CPU()
    simulation = hoomd.Simulation(device=device, seed=rng.integers(2**16))
    simulation.create_state_from_snapshot(snapshot=frame)

    # Integration and pair interactions
    integrator = hoomd.md.Integrator(dt=0.00001)
    cell = hoomd.md.nlist.Cell(buffer=0.5)
    lj = hoomd.md.pair.LJ(nlist=cell)

    # Interactions are now split into identical and cross-term
    types = frame.particles.types

    # First, we do identical terms
    for t, d in zip(types, diameters):
        lj.params[(t, t)] = dict(epsilon=1, sigma=d)
        lj.r_cut[(t, t)] = 2 ** (1.0 / 6.0) * d

    # Now we do cross terms, we double loop
    for i in range(n_particles - 1):
        for j in range(i + 1, n_particles):
            cross_sigma = (diameters[i] + diameters[j]) / 2.0
            lj.params[(types[i], types[j])] = dict(epsilon=1, sigma=cross_sigma)
            lj.r_cut[(types[i], types[j])] = 2 ** (1.0 / 6.0) * cross_sigma

    current_packing = np.pi * np.sum(diameters) / (4.0 * L**2)
    device.notice(current_packing)

    integrator.forces.append(lj)

    # Set the temperature and the integration method
    temperature = 1.5
    brownian = hoomd.md.methods.Brownian(filter=hoomd.filter.All(), kT=temperature)
    integrator.methods.append(brownian)
    simulation.operations.integrator = integrator
    simulation.state.thermalize_particle_momenta(
        filter=hoomd.filter.All(), kT=temperature
    )

    simulation.run(1e6)

    snapshot = simulation.state.get_snapshot()
    frame = from_snapshot_to_frame(snapshot)

    with gsd.hoomd.open(name=PROJECT_DIR.joinpath("random.gsd"), mode="w") as f:
        f.append(frame)

    # Compression step
    initial_volume_fraction = simulation.state.N_particles / simulation.state.box.volume
    device.notice(f"Initial density: {initial_volume_fraction}")
    final_packing_fraction = 0.74
    final_rho = final_packing_fraction * 4.0 / np.pi
    final_volume = simulation.state.N_particles / final_rho
    inverse_volume_ramp = hoomd.variant.box.InverseVolumeRamp(
        initial_box=simulation.state.box,
        final_volume=final_volume,
        t_start=simulation.timestep,
        t_ramp=100_000,
    )
    box_resize = hoomd.update.BoxResize(
        trigger=hoomd.trigger.Periodic(10),
        box=inverse_volume_ramp,
    )
    simulation.operations.updaters.append(box_resize)
    simulation.run(5e5)
    device.notice(
        f"Final density: {simulation.state.N_particles / simulation.state.box.volume}"
    )
    simulation.operations.updaters.remove(box_resize)

    snapshot = simulation.state.get_snapshot()
    frame = from_snapshot_to_frame(snapshot)

    with gsd.hoomd.open(name=PROJECT_DIR.joinpath("compressed.gsd"), mode="w") as f:
        f.append(frame)

if __name__ == "__main__":
   main()

I also provide here a snapshot of my system, which shows ordering, something that should also not happen at this state point (temperature and density).

[joaander]

You assign the same type id to all particles: frame.particles.typeid = np.zeros((n_particles,)). Therefore, they will all interact with the LJ sigma you assigned for the pair (types[0], types[0]).

You should assign different type ids to the particles that you want to have a different size.

[edwinb-ai]

Thank you! Seems to be working now.