Hi all,
I am trying to simulate a bunch of patchy particles in Lammps. I have defined the patchy particle with two attractive patches (using small particles bonded to the central larger particle) placed on the opposite poles.
Now, I use a fix rigid/small/langevin for the thermostat and use a fix defrom to do shear simulations.
I find that using the rigid style blows up the simulation after some steps and some particles are lost. Given below is an excerpt from the input script:
I have used 1000 particles (3000 total counting the patch particles also) in a simulation box. The input data file is such that each particle has a different molecule ID. I have tried decreasing the timestep also. Its not working.
compute t_def all temp/deform
fix 1 all rigid/nvt/small molecule temp 1.0 1.0 (100.0*dt)
fix_modify 1 temp t_def
fix 2 all deform 1 xy erate {xyrate} remap v
I get the following error:
ERROR: Lost atoms: original 3000 current 2082 (…/thermo.cpp:438)
Last command: run ${eq_step}
Can someone help?
Regards.
There is not enough information here to and no (simple, small) example to reproduce the issue, so it is not possible to provide any specific recommendation.
What makes you conclude that the cause of the crash lies in fix rigid/small and not - for instance - the choice of model and parameters or the use of fix deform or something else?
Thanks for providing an input and data file.
They clearly confirm, that it is not fix rigid/small that is causing the problem directly, but it is your data file, that is not conforming to what fix rigid/small needs.
The output and error message clearly shows what is going wrong. You have:
fix 1 all rigid/nvt/small molecule temp 1.0 1.0 $(100.0*dt)
fix 1 all rigid/nvt/small molecule temp 1.0 1.0 0.5
create bodies CPU = 0.001 seconds
1 rigid bodies with 3000 atoms
41.134664 = max distance from body owner to body atom
and
WARNING: Bonds are defined but no bond style is set (src/force.cpp:191)
WARNING: Likewise 1-2 special neighbor interactions != 1.0 (src/force.cpp:193)
WARNING: Using compute temp/deform with no fix deform defined (src/compute_temp_deform.cpp:78)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.62
ghost atom cutoff = 3.62
binsize = 1.3375, bins = 35 35 35
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/multi/newton/tri
stencil: half/multi/3d/tri
bin: multi
WARNING: Communication cutoff adjusted to 3.62 (src/comm.cpp:732)
ERROR: Rigid body extent > ghost cutoff - use comm_modify cutoff (src/RIGID/fix_rigid_small.cpp:645)
That means, that your data file does not indicate individual rigid objects by molecule IDs as required by fix rigid/small and thus the entire system is treated as a single rigid object.
Apart from the fact, that this is not what you want, it would work with fix rgid, but not with fix rigid small, where all constituent atoms of a rigid object have to be within the subdomain plug ghost atom cutoff.
Please also note the warning about having bonds defined without a bond style. If you want to have bonds present, but do not want to compute bonded interactions yet want to have the exclusions applied as if bonds were present, then you can use bond_style zero.
TL;DR,
you need to correct your data file or correct the procedure how you created it.
axel.
Thanks a lot Dr. Kohlmayer.
I modified the data file by changing the molecule IDs of each of the particles and entirely removing all the bonds. I still get an error :
I have attached the new data file.
Lost atoms: original 3000 current 2082 (…/thermo.cpp:438)
Last command: run ${eq_step}
data_new.packed (133 KB)
That data file works for me. It ran for at least 500000 steps with the provided input file (until I stopped).
What LAMMPS version do you use?
BTW: your neighbor list setting is extremely wasteful. commenting it out and thus sticking with the default setting will make your input run over 30% faster for me.
I am using the lammps 24 Jan2020 version. The crash occurs after around 11,00,000 steps for dt = 0.0002. I did some calculations and it seems like this time step corresponds to the one when the box boundaries are remapped after a single shear cycle i.e when the box vector reaches L0/2 and then is rescaled back to its initial position.
Is there a workaround for this?
None that I know. I’ve never done such a simulation and thus never have run across this kind of issue.
I suspect you would need to change the source code and trigger that fix rigid/small rebuilds the rigid objects at that step.
I am using the lammps 24 Jan2020 version. The crash occurs after around 11,00,000 steps for dt = 0.0002. I did some calculations and it seems like this time step corresponds to the one when the box boundaries are remapped after a single shear cycle i.e when the box vector reaches L0/2 and then is rescaled back to its initial position.
Is there a workaround for this?
Thanks a lot!!
Thanks a lot Dr. Kohlmayer. I will try to work on that.