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