Dear all,
I am working in a project to run simulations in LAMMPS for particles each of them consists of 36 beads in the ring with hexagonal shape and 24 beads into 12 chains (two beads for each chain). The ring needs to be rigid and the chains are flexible. In the input file I use the commands:
group ring type 1
group chain type 2
fix 1 ring rigid/npt molecule temp 2.00 2.00 0.50 iso 0.05 0.05 1.0
fix 2 chain npt temp 2.00 2.00 0.50 iso 0.05 0.05 1.0
I used also the time step of 0.001 in LJ units. The problem is that in the output file (log file) I received NaN and I tried to fix the problem playing with the value of time step. But the problem remains. Is it something wrong in my input file, please?
I appreciate any help! Thank you very much!
Kind Regards,
Nikoletta Pakalidou
Getting NaN is typically a sign of bad dynamics. Is there a lot of overlap in your initial setup? If it is too large, then a small time step can’t save you I guess.
What interaction potential do you have between your particles? You can try replacing it with a soft potential and equilibrating your system first to reduce overlap.
Also, it helps to send all of your data files and input scripts.
If you want to perform NPT dynamics (barostatting), you cannot use both fix npt and fix rigid/npt !
Follow this link :
http://lammps.sandia.gov/doc/fix_rigid.html
Good catch SaeeD!
There also seem to be some problems with your data file I think. You made all particles part of molecule 0, which, I think, is not what you want, since your fix treats all atoms in the same molecule as one rigid body. (Can you use molecule ID == 0 at all?).
Furthermore, you have a lot of bonds between atoms of type 1, which make up the ring. If you use fix rigid, there is no need to define bonds between them if the rigid fix uses molecule ids to determine the clusters.
Maybe you should start with a single ring, run a simulation (without both fix rigid/npt and fix npt!) to see if everything behaves as expected, and then replicate that system?
Hi Stefan,
I’m sorry about my stupid question, but what do you mean that all particles are made part of molecule 0?
Also, you are right. I have a lot of bonds between atoms of type 1 (atoms into the ring) so perhaps I should not use the fix rigid command for something extra. But if I want to use it, it will be wrong?
Check this page to see how LAMMPS reads in data files: http://lammps.sandia.gov/doc/read_data.html
The second column for atom_style angle is the molecule id, i.e., all atoms that have the same molecule id belong to the same molecule. If you tell fix rigid to work on molecules, then it will treat all atoms that belong to the same molecule as one rigid body. I don’t think that is what you want.
Using fix rigid on atoms that are also bonded together is not “wrong” if you account for it everywhere where it matters (interpreting energie and the like), but it is computationally not very efficient and it also makes debugging your whole system a lot more complicated. Keep in mind that the bond will not do anything as the particles that are rigid can not move with respect to each other.
Again, create a simpler data file with a single molecule, check that it behaves as it should and then replicate that system, that is much more easy.
Hi all,
Thank you so much for your help. You are right. The problem was at the molecule id. Thank you again!
Best,
Nikoletta
From the fix rigid doc page:
For bodystyle molecule, each set of atoms in the fix group with a different molecule ID is treated as a rigid body. This option is allowed for fix rigid and fix rigid/small, and their sub-styles. Note that atoms with a molecule ID = 0 will be treated as a single rigid body. For a system with atomic solvent (typically this is atoms with molecule ID = 0) surrounding rigid bodies, this may not be what you want. Thus you should be careful to use a fix group that only includes atoms you want to be part of rigid bodies.
Steve