droplet particle placed in groups with gaps

[Bharat-github6]

I am running the initialization for 10000 particle with 5000 on each droplet using the code below:
import mbuild as mb
import numpy as np
import random
import hoomd
import hoomd.md
import gsd.hoomd
import os

Define Parameters for 10,000 particles

n_compounds = 5000 # Particles per droplet (5000 in each droplet)
total_particles = 2 * n_compounds # Total particles (10,000)
A_ratio = 0.65 # Fraction of type A
density = 1.16 # Target density for liquid-like behavior

Larger particle sizes for better interaction

radius_A = 0.5 # Type A size
radius_B = 0.5 # Type B size

Compute droplet radius and reduce it slightly to avoid overlap

radius = np.sqrt(n_compounds / (np.pi * density)) * 0.9 # Smaller droplet radius (10% reduction)

Increase gap further (encourage coalescence, avoid overlap)

gap = radius * 0.4 # Increased gap to avoid overlap but still close

=== Generate 2D Circular Droplets ===

def generate_circular_droplet(center, droplet_radius, num_particles, radius_A, radius_B):
"""Generates a proper circular droplet with particles randomly placed inside."""
droplet = mb.Compound()
count = 0

while count < num_particles:
    r = droplet_radius * np.sqrt(random.uniform(0, 1))  # Uniformly distributed inside circle
    theta = random.uniform(0, 2 * np.pi)  # Random angle
    x = center[0] + r * np.cos(theta)
    y = center[1] + r * np.sin(theta)
    pos = np.array([x, y, 0])  # 2D placement

    particle = mb.Compound(name="A" if count < int(A_ratio * num_particles) else "B", pos=pos)
    particle.radii = radius_A if [particle.name](http://particle.name/) == "A" else radius_B  # Assign larger radii
    droplet.add(particle)
    count += 1

return droplet

Move droplets even closer together

droplet1_center = np.array([-radius - gap, 0, 0])
droplet2_center = np.array([radius + gap, 0, 0])

Generate two 2D circular liquid droplets

droplet1 = generate_circular_droplet(droplet1_center, radius, n_compounds, radius_A, radius_B)
droplet2 = generate_circular_droplet(droplet2_center, radius, n_compounds, radius_A, radius_B)

Combine both droplets into one system

system = mb.Compound()
system.add(droplet1)
system.add(droplet2)

Define HOOMD-recognized particle types

system.particle_types = ["A", "B"]

Expand box size based on bounding box (reduced box size for tighter packing)

bounding_box = system.get_boundingbox()
box_size = max(bounding_box.lengths[0], bounding_box.lengths[1]) * 2.5 # Reduced from 3 to 2.5 for tighter packing
r_cut = 2.5 # Increased Lennard-Jones cutoff
box_z = 1 # make z dimension to 0 to ensure a 2D system

Set the simulation box

system.box = mb.Box(lengths=[box_size, box_size, box_z])

print(f"Final Box Size: ({box_size}, {box_size}, {box_z})")

Remove existing file before saving

if os.path.exists("droplets_2D_10000.gsd"):
os.remove("droplets_2D_10000.gsd")

Save to file

system.save("droplets_2D_10000.gsd")# Initialize HOOMD Simulation
device = hoomd.device.CPU()
sim = hoomd.Simulation(device=device, seed=42)
sim.create_state_from_gsd(filename="droplets_2D_10000.gsd")
box = hoomd.Box(Lx=box_size, Ly=box_size, Lz=0)
sim.state.set_box(box)

Step 1: Gaussian Warm-up (Increased to 500,000 steps for better equilibration)

nlist = hoomd.md.nlist.Cell(buffer=0.2)
gaussian = hoomd.md.pair.Gaussian(nlist=nlist, default_r_cut=1.5)
gaussian.params[('A', 'A')] = dict(epsilon=5.0, sigma=0.5)
gaussian.params[('B', 'B')] = dict(epsilon=5.0, sigma=0.5)
gaussian.params[('A', 'B')] = dict(epsilon=5.0, sigma=0.5)

integrator_gs = hoomd.md.Integrator(dt=0.0005) # Reduced time step
brownian = hoomd.md.methods.Brownian(filter=hoomd.filter.All(), kT=0.0) # Correct class name here
integrator_gs.methods.append(brownian)
integrator_gs.forces.append(gaussian)
sim.operations.integrator = integrator_gs

print("Running Gaussian warm-up (500,000 steps)...") # Increased warm-up steps
sim.run(500000)
print("Gaussian warm-up completed!")
hoomd.write.GSD.write(state=sim.state, filename="aftergaussian_10000.gsd", mode='wb')

Step 3: Lennard-Jones Equilibration (Reverted LJ parameters)

lj = hoomd.md.pair.LJ(nlist=nlist, default_r_cut=r_cut)
brownian_lj = hoomd.md.methods.Brownian(filter=hoomd.filter.All(), kT=0.35) # Correct class name here

Revert the Lennard-Jones parameters to those that worked before

lj.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0) # Epsilon for A-A interaction
lj.params[('B', 'B')] = dict(epsilon=0.5, sigma=0.88) # Epsilon for B-B interaction
lj.params[('A', 'B')] = dict(epsilon=1.5, sigma=0.8) # Epsilon for A-B interaction

Reduce time step to allow for better stability in the Lennard-Jones equilibration

integrator_lj = hoomd.md.Integrator(dt=0.0001) # Smaller time step for better stability

integrator_lj.forces = [lj]
integrator_lj.methods = [brownian_lj]
sim.operations.integrator = integrator_lj

gsd_writer = hoomd.write.GSD(filename="LJEquilibration_10000.gsd", trigger=hoomd.trigger.Periodic(20000), filter=hoomd.filter.All(), mode='wb')
sim.operations.writers.append(gsd_writer)

print("Running Lennard-Jones equilibration (500,000 steps)...") # Extended equilibration time
sim.run(5000000)
print("Lennard-Jones equilibration completed!")

Define Positive and Negative Filters

class PositiveFilter(hoomd.filter.CustomFilter):
def hash(self):
return hash("PositiveFilter")

def __eq__(self, other):
    return isinstance(other, PositiveFilter)

def __call__(self, state):
    with state.cpu_local_snapshot as snap:
        xposition = snap.particles.position[:, 0]
        indices = (xposition >= 0)
        return snap.particles.tag[indices]

class NegativeFilter(hoomd.filter.CustomFilter):
def hash(self):
return hash("NegativeFilter")

def __eq__(self, other):
    return isinstance(other, NegativeFilter)

def __call__(self, state):
    with state.cpu_local_snapshot as snap:
        xposition = snap.particles.position[:, 0]

class NegativeFilter(hoomd.filter.CustomFilter):
def hash(self):
return hash("NegativeFilter")

def __eq__(self, other):
    return isinstance(other, NegativeFilter)

def __call__(self, state):
    with state.cpu_local_snapshot as snap:
        xposition = snap.particles.position[:, 0]

class NegativeFilter(hoomd.filter.CustomFilter):
def hash(self):
return hash("NegativeFilter")

def __eq__(self, other):
    return isinstance(other, NegativeFilter)

def __call__(self, state):
    with state.cpu_local_snapshot as snap:
        xposition = snap.particles.position[:, 0]
        indices = (xposition < 0)
        return snap.particles.tag[indices]

Slowly move the two droplets close to each other (same logic as before)

positive_filter = PositiveFilter()
negative_filter = NegativeFilter()

positiveforce = hoomd.md.force.Constant(filter=negative_filter)
negativeforce = hoomd.md.force.Constant(filter=positive_filter)
positiveforce.constant_force["A", "B"] = (0.05, 0, 0)
negativeforce.constant_force["A", "B"] = (-0.05, 0, 0)

integrator_lj.forces.append(positiveforce)
integrator_lj.forces.append(negativeforce)
sim.run(400000)

print("Approaching completed!")

Save Final State

From the above code, the particles are arranged well after gaussian warm up but after lennard jones the output file opened in ovito shows some groupism between particles but not overlapping and huge gaps in between each each groupism .How to reduce it. I USD the same gaussian and lennard jones parameter for 5000 particle but they were placed well so I think the parameters are good enough and I increased box sizes for 10000 particle than 5000 particle. i generated initial structure for 5000 parcile but having problem with this 10000 particle. i want to avoid those gaps in the 10000 particle how to do it?

[joaander]

Your code is very long and not well formatted. This makes it very hard to read. See https://docs.github.com/en/get-started/writing-on-github/getting-started-with-writing-and-formatting-on-github to learn how to format your posts.

If I understand it correctly, you place particles near each other randomly, push them apart with a repulsive potential, and then continue the simulation with an attractive potential. If this is the case, I would expect small aggregates to form initially. Depending on the size of your simulation box and the diffusion rates of the particles, it could require 10's of millions to many billions (or more) steps to form a single aggregate. Note that Brownian integration will have the slowest diffusion rate of all the integration methods.

This discussion board is to answer questions regarding the technical details of using HOOMD-blue. You should work with your advisor on questions related to your simulation outcomes and how those relate to your research goals.