Discrepancy in particle coordinates using rigid/small style in a linear molecule (trimer) model

LAMMPS LAMMPS Development
1

Hello everyone!

I’m currently implementing a model of a linear molecule composed of three particles, where the intramolecular interaction is defined using bond and angle harmonic potentials. The intermolecular interaction is modeled using a Lennard-Jones potential with a cutoff of 2.2. However, I’m facing an issue with the dynamics of the system after a certain simulation time when trying to fix the angle using the rigid/small style.

Upon analyzing the results using VMD, I observed a significant discrepancy in the particle coordinates when utilizing the xu yu zu and x y z styles. Specifically, the “unwrapped” coordinates show an abrupt drift of the entire simulation box in a certain direction. I have tried different ensembles, but it is intriguing that in NPT (constant number of particles, pressure, and temperature) simulations, the x y z coordinates remain centered, and the system does not exhibit any drift. Additionally, all thermodynamic parameters seem to be stable without any temporal drift.

I have attempted several variations, such as adjusting the time_step, using fix_recenter, fix_momentum, and neigh_modify exclude, but none of these approaches have resolved the issue. In all cases, using the xu yu zu styles, the system experiences unexpected behavior. However, when using rattle and fixing the angle, I have not encountered any problems, although it is not recommended for linear molecules.

I would greatly appreciate any guidance or suggestions on why this discrepancy in the particle coordinates is occurring and how I can effectively resolve it. Are there any special considerations I should be aware of when using the rigid/small style in such simple linear molecule systems?

Thank you all for your time and expertise! Any help would be highly appreciated.

Best regards,
Cristian

2

The is all too vague and general to make some assessment and try to come up with an explanation (or figure out where your expectations are not aligned with the documentation or the documentation not with the implementation). You need to provide a minimal input example to two the two different cases and explain where exactly you are looking.

I also suggest you have a close look at the documentation of fix rigid/small since there should be several paragraphs discussing image flags. If I am not mistaken, those may be changed when the rigid bodies are defined to have a consistent behavior since with the rigid/small fixes one of the constituent atoms is pegged to discriminate whether image flags need to be incremented and decremented and also determines the migration between subdomains.

3

Hello Axel! Of course, here’s the input that I’m implementing:

units lj
atom_style angle

neighbor 0.8 bin
#neighbor 0.4 bin
neigh_modify every 1 delay 1
comm_modify vel yes

pair_style lj/cut 2.2
pair_modify shift yes
bond_style harmonic
angle_style harmonic

read_data DATA

pair_modify shift yes
pair_coeff * * 1.0 1.0
bond_coeff 1 100 0.96
angle_coeff 1 50 180.0
special_bonds lj 0 0 1

timestep 0.01
#timestep 0.005

#fix MOM all momentum 100 linear 1 1 1 angular
#neigh_modify exclude molecule/intra all

fix 1 all rigid/npt/small molecule temp $x $x 1 iso 0.0 0.0 10.0
#fix 1 all rigid/nvt/small molecule temp $x $x 1
#fix REC all recenter INIT INIT INIT
#fix 1999 all rattle 0.0001 20 10000 a 1
run 1000000

From reading the forum, what I’m experiencing is a “Flying ice cube.” I have also tried resetting the velocities frequently, but to no avail – it keeps flying away!

As an additional piece of information that might help in identifying my error, it’s worth noting that the DATA file I’m reading represents an equilibrated configuration that has gone through several equilibration processes/times. Therefore, the image values are relatively high (20 <|i| < 100).

Another factor to consider is that using the same input with an equilibrated data file containing non-linear molecules (e.g., 90-degree angles), there are no issues with system drift in the unwrapped coordinates. In that case, the system’s dynamics closely resemble the system without the constraint imposed by the rigid command.

Regarding the documentation of fix rigid/small and its image_flags parameter, I understand that the issue might lie there. However, I find it challenging to grasp the precise modifications needed to avoid the system drift as described.

I keep reading! Thank you in advance for your help!!

4

That should not be a problem, since the limit for a normal LAMMPS executable (i.e. not compiled with -DLAMMPS_BIGBIG) is 512. But you can also use reset_atoms image to reset those to be as close to zero as possible.

One more item to experiment with would be the timestep.
And another is to use fix rigid/npt instead of fix rigid/npt/small.

I cannot make any promises to have a closer look until we have released the next stable version of LAMMPS, but if you cannot sort it out until then, please also make the data file available for some experimentation.

5

Dear Axel, thank you so much for your help. I have reinstalled LAMMPS to its current version. I was previously using one from 2020 (my mistake, I thought I was using this year’s version!), and I no longer get the ‘Flying ice cube’ effect. It must have been a little coding issue back then.

Best regards,
Cristian

2 Likes
6

The LAMMPS developers are always happy to hear when a problem was resolved by using a more recent version. This is what we strive for.

7

DATA_4_0.6 (1.0 MB)

Hello again, well, new trials regarding the ‘Flying ice cube’ effect; the new version of LAMMPS still generates this effect, but at much longer times compared to previous versions. On the other hand, I see that the temperature oscillates to a lesser extent in the new version.

I’m attaching the DATA file in case you want to run some tests on it.

Also, to avoid this effect in the runs, it can be achieved using the command ‘velocity all zero linear,’ as Axel and other users have suggested in various posts. Here’s an example:

variable ntotal index 1000000
variable nchunk index 200000
variable nloop index (round(v_ntotal/v_nchunk)) variable i loop {nloop}
label loop

velocity all zero linear
fix rigid all rigid/npt/small molecule temp $x x 1.0 iso 0.0 0.0 1.0 run {nchunk} post no

next i
jump SELF loop

8

Please reformat using triple backquotes (```) for quoting from input files. As such this is useless since the forum software recognizes characters like $ or * or _ as indicators for applying text style formatting.

Please see Please Read This First: Guidelines and Suggestions for posting LAMMPS questions for details on this (and other things :wink:)