PArticle jumping out of the box while initializing and initial warm up. #2006
Replies: 4 comments 2 replies
-
|
Just to have every script in the same place === Initialize HOOMD (CPU) === === Define Parameters === === Compute Safe Box Size === === Adjust Droplet Separation to Fit Inside the Box === print(f"Using box size: {box_size:.2f} x {box_size:.2f} x {box_z:.2f}") === Generate Random Positions for Two Non-Overlapping Circular Droplets === def generate_particles_in_circle(center_x, center_y, radius, num_particles, type_A_limit): ==Generate first droplet (left) == Generate second droplet (right) == Combine lists === Apply Periodic Boundary Conditions === === Create HOOMD Snapshot === === Apply snapshot to the simulation === Neighbor List === === Step 1: WCA Warm-up === wca = hoomd.md.pair.LJ(nlist=nlist) integrator.forces.append(wca) print("Running WCA Warm-up (5,000 steps)...") === Step 2: Lennard-Jones Equilibration === lj = hoomd.md.pair.LJ(nlist=nlist) integrator.forces.append(lj) print("Running Lennard-Jones Equilibration (5,000 steps)...") ~ |
Beta Was this translation helpful? Give feedback.
-
|
Uniform random placement of particles is almost certain to produce overlapping particles with initial forces so large that time integration can never be stable. You should either a) place particles so that the initial forces are small or b) start with a potential that has bounded forces (DPDConservative or Gaussian) and run until the particles separate before activating switching to the Lennard Jones potential. For future reference, you can use markdown formatting to make code blocks easier to read: https://docs.github.com/en/get-started/writing-on-github/getting-started-with-writing-and-formatting-on-github/basic-writing-and-formatting-syntax?#quoting-code |
Beta Was this translation helpful? Give feedback.
-
|
Thank you for the information but I am facing the overlapping issue now. Do you have any suggestion how to resolve overlapping of particle? here is the updated code :import hoomd === Initialize HOOMD (CPU) ===device = hoomd.device.CPU() === Define Parameters ===N = 5000 print(f"✅ Using box size: {box_size:.2f} x {box_size:.2f} x {box_z:.2f}") === Generate Random Placement with Monte Carlo Sampling for Liquid Droplets ===def generate_random_droplet(center_x, center_y, radius, num_A, num_B, sigma): Generate two droplets with non-overlapping placementpos1, type1 = generate_random_droplet(-separation / 2, 0, radius, int(0.65 * N // 2), int(0.35 * N // 2), sigma) positions = np.array(pos1 + pos2) Apply Periodic Boundary Conditionspositions[:, 0] = (positions[:, 0] + box_size / 2) % box_size - box_size / 2 === Create Initial HOOMD Snapshot ===snapshot = hoomd.snapshot.Snapshot() === Step 1: Gaussian Warm-up Phase ===nlist = hoomd.md.nlist.Cell(buffer=buffer) print("Running Gaussian warm-up (100,000 steps)...") === Step 2: Lennard-Jones Equilibration ===lj = hoomd.md.pair.LJ(nlist=nlist, default_r_cut=r_cut) print("Running Lennard-Jones equilibration (200,000 steps)...") === Step 3: DPD Equilibration ===dpd = hoomd.md.pair.DPD(nlist=nlist, kT=kT_LJ, default_r_cut=1.5) print("Running DPD equilibration (200,000 steps)...") |
Beta Was this translation helpful? Give feedback.
-
|
@Bharat-github6 you can use Packmol to randomly fill a sphere while avoiding overlapping particles. mBuild is a python package that does lots of things, but one of these is a useful python wrapper around Packmol, and can save structures to the .gsd file format. Here is a quick example of creating a sphere with an A and B mixture. You can check out the parameters in
|
Beta Was this translation helpful? Give feedback.
-
|
Hi Chris,
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 == "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?
*Bharat Poudel*
*Graduate Student *
*Department of Material Science*
*University of Vermont*
*Burlington, VT, USA*
*601-466-6774*
…On Tue, Feb 25, 2025 at 3:50 PM Chris Jones ***@***.***> wrote:
@Bharat-github6 <https://github.com/Bharat-github6> you can use Packmol
<https://m3g.github.io/packmol/> to randomly fill a sphere while avoiding
overlapping particles.
mBuild <https://github.com/mosdef-hub/mbuild> is a python package that
does lots of things, but one of these is a useful python wrapper around
Packmol, and can save structures to the .gsd file format.
Here is a quick example of creating a sphere with an A and B mixture. You
can check out the parameters in mbuild.packing.fill_sphere to change
things like sphere size, particle spacing (i.e. radius), etc..
import mbuild as mb
import random
bead = mb.Compound(name="X")
n_compounds = 1000
A_ratio = 0.65
num_A = int(A_ratio * n_compounds)
sphere = mb.fill_sphere(compound=bead, n_compounds=n_compounds, sphere=[5, 5, 5, 2])
particles = [p for p in sphere.particles()]
random.shuffle(particles)
for particle in particles[0:num_A + 1]:
particle.name = "A"
for particle in particles[num_A+1:-1]:
particle.name = "B"
sphere.save("sphere.gsd")
sphere.visualize().show()
—
Reply to this email directly, view it on GitHub
<#2006 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AN64NKDTI3UMUIBNTRCOCXT2RTXUVAVCNFSM6AAAAABXW6Q7PKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTEMZRHA3TQMA>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
Beta Was this translation helpful? Give feedback.
-
|
[image: initial_gaussian.png]
*for your reference, I have attached the snapshot after gaussian and the
trajectory video after lennard jones , I want the final structure to be
similar to the png I have attached *
*Bharat Poudel*
*Graduate Student *
*Department of Material Science*
*University of Vermont*
*Burlington, VT, USA*
*601-466-6774*
On Tue, Mar 18, 2025 at 10:49 AM bharat poudel ***@***.***>
wrote:
… Hi Chris,
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 == "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?
*Bharat Poudel*
*Graduate Student *
*Department of Material Science*
*University of Vermont*
*Burlington, VT, USA*
*601-466-6774*
On Tue, Feb 25, 2025 at 3:50 PM Chris Jones ***@***.***>
wrote:
> @Bharat-github6 <https://github.com/Bharat-github6> you can use Packmol
> <https://m3g.github.io/packmol/> to randomly fill a sphere while
> avoiding overlapping particles.
>
> mBuild <https://github.com/mosdef-hub/mbuild> is a python package that
> does lots of things, but one of these is a useful python wrapper around
> Packmol, and can save structures to the .gsd file format.
>
> Here is a quick example of creating a sphere with an A and B mixture. You
> can check out the parameters in mbuild.packing.fill_sphere to change
> things like sphere size, particle spacing (i.e. radius), etc..
>
> import mbuild as mb
> import random
>
> bead = mb.Compound(name="X")
> n_compounds = 1000
> A_ratio = 0.65
> num_A = int(A_ratio * n_compounds)
>
> sphere = mb.fill_sphere(compound=bead, n_compounds=n_compounds, sphere=[5, 5, 5, 2])
> particles = [p for p in sphere.particles()]
> random.shuffle(particles)
>
> for particle in particles[0:num_A + 1]:
> particle.name = "A"
> for particle in particles[num_A+1:-1]:
> particle.name = "B"
>
> sphere.save("sphere.gsd")
> sphere.visualize().show()
>
> —
> Reply to this email directly, view it on GitHub
> <#2006 (comment)>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AN64NKDTI3UMUIBNTRCOCXT2RTXUVAVCNFSM6AAAAABXW6Q7PKVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTEMZRHA3TQMA>
> .
> You are receiving this because you were mentioned.Message ID:
> ***@***.***
> com>
>
|
Beta Was this translation helpful? Give feedback.
Uh oh!
There was an error while loading. Please reload this page.
-
Hi, I am new to hoomd-blue and started with initializing the two coalescence droplets with certain distance apart with the certain density of particles.
Few particles moved out of thebox
Cartesian coordinates:
x: -4.54029e+07 y: 1.95933e+07 z: 0
Fractional coordinates:
f.x: -433368 f.y: 187018 f.z: 0.5
Local box lo: (-52.3837, -52.3837, -3.1)
hi: (52.3837, 52.3837, 3.1)
Can you please suggest how to fix this PBC conditions.
Script is given below:
import hoomd
import hoomd.md
import numpy as np
=== Initialize HOOMD (CPU) ===
device = hoomd.device.CPU()
sim = hoomd.Simulation(device=device, seed=42)
=== Define Parameters ===
N = 5000 # Total number of particles
fraction_A = 0.65
fraction_B = 0.35
N_A = int(N * fraction_A)
N_B = N - N_A
density = 1.16 # Target density
sigma = 1.0 # Particle diameter
dt_initial = 1e-6 # Very small step for stability
dt_intermediate = 0.0001 # Slightly larger after minimization
dt_final = 0.001 # Final time step for LJ
kT = 0.35 # Lennard-Jones temperature
r_cut = 2.5 * sigma # Lennard-Jones cutoff
buffer = 0.1 # Neighbor list buffer
rmax = r_cut + buffer # Neighbor search radius
=== Compute Safe Box Size ===
radius = np.sqrt(N / (2 * np.pi * density)) # Radius of each droplet
box_size = 4 * radius # Box size adjusted to fit both droplets
box_z = 2 * rmax + 1.0 # Small z-dimension for 2D simulation
=== Adjust Droplet Separation to Fit Inside the Box ===
max_separation = box_size / 2 - radius - rmax # Prevent droplets from touching box edges
separation = min(2 * radius + 3 * rmax, max_separation)
print(f"Using box size: {box_size:.2f} x {box_size:.2f} x {box_z:.2f}")
print(f"Droplet separation: {separation:.2f}")
=== Generate Random Positions for Two Non-Overlapping Circular Droplets ===
positions = []
types = []
def generate_particles_in_circle(center_x, center_y, radius, num_particles, type_A_limit):
local_positions = []
local_types = []
Generate first droplet (left)
pos1, type1 = generate_particles_in_circle(-separation / 2, 0, radius, N // 2, N_A // 2)
Generate second droplet (right)
pos2, type2 = generate_particles_in_circle(separation / 2, 0, radius, N - len(pos1), N_A - len(type1))
Combine lists
positions = np.array(pos1 + pos2)
types = type1 + type2
=== Apply Periodic Boundary Conditions ===
positions[:, 0] = (positions[:, 0] + box_size / 2) % box_size - box_size / 2
positions[:, 1] = (positions[:, 1] + box_size / 2) % box_size - box_size / 2
positions[:, 2] = 0
=== Create HOOMD Snapshot ===
snapshot = hoomd.Snapshot()
snapshot.configuration.box = [box_size, box_size, box_z, 0, 0, 0]
snapshot.particles.N = N
snapshot.particles.types = ['A', 'B']
snapshot.particles.typeid[:] = [0 if t == 'A' else 1 for t in types]
snapshot.particles.position[:] = positions
snapshot.particles.velocity[:] = np.zeros((N, 3)) # Set all initial velocities to zero
Apply snapshot to the simulation
sim.create_state_from_snapshot(snapshot)
print("Simulation state successfully initialized!")
=== Neighbor List ===
nlist = hoomd.md.nlist.Cell(buffer=buffer, exclusions=[])
=== Step 1: WCA Warm-up ===
integrator = hoomd.md.Integrator(dt=dt_initial)
cv = hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All())
integrator.methods.append(cv)
wca = hoomd.md.pair.LJ(nlist=nlist)
wca.params[('A', 'A')] = dict(epsilon=0.1, sigma=1.0)
wca.params[('B', 'B')] = dict(epsilon=0.1, sigma=1.0)
wca.params[('A', 'B')] = dict(epsilon=0.1, sigma=1.0)
wca.r_cut[('A', 'A')] = 2**(1/6) * sigma
wca.r_cut[('B', 'B')] = 2**(1/6) * sigma
wca.r_cut[('A', 'B')] = 2**(1/6) * sigma
wca.shift_mode = 'xplor'
integrator.forces.append(wca)
sim.operations.integrator = integrator
print("Running WCA Warm-up (5,000 steps)...")
sim.run(5000)
print("WCA Warm-up Completed!")
=== Step 2: Lennard-Jones Equilibration ===
integrator.dt = dt_final
langevin = hoomd.md.methods.Langevin(filter=hoomd.filter.All(), kT=kT)
integrator.methods[0] = langevin
lj = hoomd.md.pair.LJ(nlist=nlist)
lj.params[('A', 'A')] = dict(epsilon=1.0, sigma=1.0)
lj.params[('B', 'B')] = dict(epsilon=1.0, sigma=1.0)
lj.params[('A', 'B')] = dict(epsilon=1.0, sigma=1.0)
lj.r_cut[('A', 'A')] = r_cut
lj.r_cut[('B', 'B')] = r_cut
lj.r_cut[('A', 'B')] = r_cut
lj.shift_mode = 'xplor'
integrator.forces.append(lj)
sim.run(5000)
~
Beta Was this translation helpful? Give feedback.
All reactions