com evaluated by com_chunks and by fix_rigid are not the same

I am trying to calculate a rotation angle of a molecule.

To do that I need the com of the molecule:

compute chunks all chunk/atom molecule
compute molcom all com/chunk chunks
compute xu all property/atom xu
compute yu all property/atom yu
variable dX atom c_xu-v_molcom[1][1]
variable dY atom c_yu-v_molcom[1][2]
variable theta atom atan2(v_dY,v_dX)*(mol==1)
fix theta0 all store/state 0 v_theta
run 0
variable dtheta atom (atan2(v_dY,v_dX)-f_theta0)*(mol==1)
fix dtheta all ave/chunk 1 1 1 chunks f_dtheta

In this way in f_dtheta I can access the molecule rotation angle.

However, I have noticed that the com of the molecule results different when evaluated by fix_rigid:

fix printout all ave/time 1 1 1 c_molcom[1] c_molcom[2] mode vector file molcom.dat
run 0
unfix printout

fix nve all rigid/nve molecule
fix printout all ave/time 1 1 1 f_nve[1] f_nve[2] mode vector file molcom_rigid.dat
run 0
unfix printout
unfix nve

In molcom.dat the com coordinates are (1.292, -0.215905), while in molcom_rigid.dat they are (1.292, 1.9841), the latter clearly being wrapped.

Such ambiguity is a problem since dX,dY should not depend on the molecule position.

In the fix_rigid manual I read:
The center of mass (COM) for each body is similar to unwrapped coordinates written to a dump file. It will always be inside (or slightly outside) the simulation box.

I find this statement unclear: unwrapped coordinates are not always inside the simulation box. Wrapped coords are.

In the compute_com_chunk manual page I read instead:
This compute calculates the x,y,z coordinates of the center-of-mass for each chunk, which includes all effects due to atoms passing through periodic boundaries.

What is the way to guarantee that dX,dY are calculated correctly, accounting for PBC?

Thank you in advance,
Roberto

The behavior of chunks with respect to coordinate unwrapping is known to be inconsistent.
https://github.com/lammps/lammps/issues/1710

This is because it combines different “consumers” of chunks with different expectations.

We’ve had discussions about this, but as you can see from the age of the original pull request. Nobody has yet found the time to seriously look into it and implement a sufficiently general solution.

Axel.

Dear Axel,

thank you for the reply and for the info.

From what I see, the fix_rigid do wraps the com coordinates (contrary to what documentation says). If I define the distance from com by the following

fix nve all rigid/nve molecule
variable dX atom x-f_nve[1]
variable dY atom y-f_nve[2]

would it be consistent with PBC if the ghost atoms are enough to cover the molecule size?

Thank you in advance,
Roberto

Dear Axel,

thank you for the reply and for the info.

From what I see, the fix_rigid do wraps the com coordinates (contrary to what documentation says). If I define the distance from com by the following

fix nve all rigid/nve molecule
variable dX atom x-f_nve[1]
variable dY atom y-f_nve[2]

would it be consistent with PBC if the ghost atoms are enough to cover the molecule size?

Thank you in advance,
Roberto

you are only looking at the wrapped COM coordinates of the rigid bodies. you also need to check out the image flags, which are computed per rigid body, too, and applied when needed.

axel.

p.s.: can you please correct the e-mail address for your mailing list subscription? for every response i get an annoying auto-response that your e-mail address is no longer current.

Quoting from So I guess the correct way is to provide unwrapped coordinates in input, and use done thank you, Roberto

Why guess?

Set up some small test input with atoms/molecules that straddle PBC boundaries or have as a group traveled multiple times through the simulation box (i.e. image flags are > 1) and compare the output LAMMPS produces against what it should produce.
To give you an authoritative answer I would have to do the same thing. LAMMPS is by far too large a code with too many features to be able to memorize the exact behavior of every detail all the time.

Axel.