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
# MSD
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.
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.
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?
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.