Hard sphere with brownian dynamics

Hi, I am new to lammps and I was trying to perform a brownian dynamics simulation of an hard sphere, but I always get the ‘Lost atom’ error. I am sure this is something that I am overlooking, but I don’t know how to solve it. The input file I am using is slightly different from the one provided in the example folder of the brownian dynamics pakcage:

variable        rng string gaussian
variable        seed string 198098
variable        temp string 1.0
variable        gamma_t string 5.0 
variable        gamma_r string 0.7
variable        params string ${rng}_${temp}_${gamma_r}_${gamma_t}

units           lj
atom_style      hybrid dipole sphere
dimension       3
boundary 	p p p 
newton off

lattice         sc 0.4
region          box block 0 10 0 10 0 10
create_box      1 box
create_atoms    1 random 1000 12345 box

pair_style lj/cut 1.12

pair_coeff 1 1 1.0 1.0

mass            * 1.0
set             type  * dipole/random ${seed} 1.0
velocity        all create 1.0 1 loop geom

neighbor        1.0 bin
neigh_modify    every 1 delay 1 check yes

fix 1 all brownian/sphere ${temp} ${seed} rng ${rng} gamma_r ${gamma_r} gamma_t ${gamma_t}

#initialisation for the main run

compute         msd  all msd
fix 		MSD_output all ave/correlate 25 100 5000 c_msd[*] file msd.out ave running  
thermo_style    custom step ke pe c_msd[*]

dump            1 all custom 1000 dump_3d.lammpstrj id type &
                x y z xu yu zu mux muy muz fx fy fz
dump_modify     1 first yes sort id

timestep        1.0
thermo          10000

# main run
run             300000000

Am I missing some major physics element or I am making some huge mistake somewhere?
Thank you all in advance!

When reporting errors, please always also report which version of LAMMPS you are using.

With this choice of pair style you do not have a “hard sphere” interaction, but rather a repulsive-only Lennard-Jones interaction. To simulate real hard spheres is rather tricky for an MD code like LAMMPS with discretized time and requires implementing a different integrator.

This timestep is far too large for your choice of units, mass, and interactions to result in stable numerics during time integration.

Sorry I forgot to mention it, I am on LAMMPS 8 Apr 2021.
Yes I understood that it is not so easy to have a real ‘hard sphere’, but should not be the source of the error, am I correct? Or should I change it drastically?

Unfortunately even if I make the time steps way smaller, I still get the error. I am able to simulate only if I drastically reduce the number of particles in the system

timestep=1 is too large, I ran your system with step 1e-6 and add this line
“minimize 1.0e-4 1.0e-6 1000 10000”
and its running fine, didnt finish it though.

In comparison to a real hard sphere interaction, the Lennard-Jones potential is rather soft and squishy.

There have been many discussions about the different causes of lost atoms errors in the past. You can search and look them up in the forum since we have archived all discussions since 2005.

How much smaller?

LAMMPS will do what you ask it to do. It will compute the forces resulting from your model and use the forces to integrate Newton’s equation of motion numerically using the Velocity Verlet algorithm. In that way, your input is syntactically correct. However, that does not say anything about the scientific relevance of your model or whether the choice of parameters and settings is compatible with each other to lead to a numerically stable calculation. There are several inconsistencies in your input:

  • the default time step for reduced units is 0.005 tau, which is quite different from 1.0 tau
  • you are using a per-atom sphere diameter property, but your potential does not use it
  • you are using a per-atom dipole property, but you potential does not use it and thus setting the dipoles or changing them does not affect the trajectory in any way. It would be purely random.
  • your use of fix ave/correlate in combination with compute msd makes no sense: 8.3.7. Calculate diffusion coefficients — LAMMPS documentation

None of these do concern the LAMMPS code - it will just compute what you ask - but they affect the relevance and meaning of your simulations. Discussing that is something that you should be doing with your adviser and others that have a vested interest in your research outcome.

Is it possible to simulate hard spheres with fix_gcmc, which can perform MC moves?

Apart from the fact that there is currently no way to enter a hard sphere potential, how should that work with how fix gcmc works? The energy would be either infinity when the spheres overlap or zero if they don’t. Fix gcmc requires a finite and comparable energy to accept/reject moves.

The issue is not so much about MD versus MC but the fact that you need to propagate hard spheres differently. With discretized time you will have to reconstruct where on the timeline two spheres had contact and then apply reflection condition for that point in time and determine the velocity vector and remaining displacement afterwards. This is quite similar to what fix wall/reflect is doing (check out its source code). With “soft” interactions you don’t have a hard reflection but just a finite force pushing or pulling and you ignore that for the discretized time the force changes a little bit over that time for as long as that is a small change and you have sufficient error cancellation. This is a different way of looking at why you cannot use arbitrary large timesteps.

Thanks, the part about lack of “hard spheres” in MD makes sense and something that I had realized quite recently.

The energy would be either infinity when the spheres overlap or zero if they don’t. Fix gcmc requires a finite and comparable energy to accept/reject moves

From what I can understand, with normal hard particle MC simulations, when spheres overlap, the energy is calculated as infinity. The move is then automatically rejected because the Boltzmann factor will always be larger than a random number between 0 and 1 (not 100% sure if I got the order of that right). So essentially, any move with overlapping spheres will just be rejected 100% of the time, and a different trial move will need to be made.

You forgot about the second part of my comment. Since the energy will be zero for all other cases how can you do anything else but just random displacements?

In my personal case, I am hoping to put in long-range electrostatic interactions into my MC simulation

What has that to do with the discussion topic of hard spheres?
Please start a new discussion for a new topic.
“Hijacking” somebody else’s discussion is generally considered a violation of forum etiquette.