[lammps-users] Fixing overlaps when populating region with rigid molecules

Hello to whomever it may concern,
I am trying to run a simulation where I read in a molecule file, then shrink the box while keeping the particles the same size to simulate particle growth. I start with highly dispersed initial conditions with random placement and orientation of around 500 molecules. I keep noticing that some molecules (around 1%) are overlapping with each other initially and they stay that way throughout the simulation. I am using fix rigid to make the molecules have a rigid structure, so I don’t think I can use a minimization step. I have considered simply running a separate simulation that randomly places down the same number of spheres and then runs a minimization to handle overlaps, and then reading that data in as position data that would be sure not to have molecule overlaps. I am concerned about ensuring random orientations for the molecules after this data is then read in with this approach. My question is: do you have an idea how to ensure my molecules do not overlap when I place them or do you know how to ensure they have random orientations from some position data to be read in? I have attached a recent version of my code below. Thank you in advance for your time and help in this matter.


dimension 2
atom_style sphere
fix prop all property/atom mol ghost yes
atom_modify map array
boundary p p p
newton off
comm_modify vel yes

#–Define Region–#

region packing block -335.4101966249685 335.4101966249685 -335.4101966249685 335.4101966249685 -0.5 0.5 units box
create_box 2 packing extra/special/per/atom 2
neighbor 1 bin
neigh_modify delay 0

#–Define Molecule–#

molecule 1 pyparts/molecule_template_input.txt
create_atoms 1 random 500 56 packing mol 1 489

#–Define Particle Interactions–#

variable kn equal 200000.0
variable kt equal 0.0
variable gamma_n equal 50.0
variable gamma_t equal 50.0
variable xmu equal 500.0
variable dampflag equal 0
pair_style gran/hertz/history {kn} {kt} {gamma_n} & {gamma_t} {xmu} {dampflag}
pair_coeff * *
timestep 0.01


fix rigfix all rigid molecule

neigh_modify exclude molecule/intra all

#–Visc Damp–#

variable Tempstart equal 0.0
variable Tempstop equal 0.0
variable damp equal 0.01 #high val equal low visc, 50= 50 time steps to damp completely
fix viscdamp all langevin {Tempstart} {Tempstop} ${damp} 478098 #zero yes


fix delbox all deform 1 x trate -3.859342747599233e-05 y trate -3.859342747599233e-05 remap x
fix thedimension all enforce2d
variable a1 equal atoms
variable v1 equal vol
compute t1 all temp/sphere
compute p1 all pressure t1

thermo 10000
thermo_style custom step time vol density atoms press temp
thermo_modify lost ignore norm no temp t1 press p1
fix globalout all ave/time 10 5 10000 v_v1 v_a1 c_p1 c_t1 mode scalar ave one &
file file_place/file_name_global.txt

#–Dump values–#

dump pos all custom 10000 file_place/file_name_dump.txt id x y fx fy mol
run 10000000

i see three possible options (in order of ease of use)

  1. don’t use random positions but a simple lattice to place your molecules but make the spacing large enough so that no overlaps from random orientations happen. initialize with random velocities and then use an MD run with fix deform to slowly compress the system to the desired final density.
  2. start with using pair style soft with a “flat” potential (so that particles may overlap as much as they wish without a problem) and then run an MD with fix adapt to “grow” the repulsive soft potential until you reach a value similar to the diameter of your granular model and then switch pair style
  3. populate the region not with create atoms but with fix deposit. for that you can indicate a minimum distance via the near keyword.