Project

General

Profile

Bug #743

Angular momentum removal incorrect

Added by David van der Spoel over 8 years ago. Updated over 7 years ago.

Status:
Closed
Priority:
Normal
Category:
mdrun
Target version:
Affected version - extra info:
all
Affected version:
Difficulty:
uncategorized
Close

Description

Hi, gromacs developers,
It seems that there is a bug regarding the angular moment removal. This may not be very important for many users since most of them are doing PBC where it is irrelevant, but if you are doing cluster or if you really care the angular momentum like I
do(I want the angular momentum to be exactly zero since this corresponds to some quantum resolved state), then it is important.
The angular momentum removal part in mdlib/vcm.c where it evaluate moment of inertia is wrong, subroutine update_tensor.
*Specifically, the diagonal part is not calculated right. According to either */Classical Mechancs, by Goldstein or Theoretical Mechanics of Particles and Continua, by A.L.Fetter and J.D.Walecka.
I_{xx}=\sum_{i} m_{i}
(r_{i}2-x_{i}2), where x can be x,y,z.
But the /**update_tensor gives
*/I_{xx}=\sum_{i} m_{i}(x_{i}2).
/*/I initially thought there could be some other routine to correct update_tensor, but apparently, it is calculated according to that, and I manually calculated them and turns out my theory is right, the update_tensor is wrong as it stands.
/Again, this would not affect, say, 99.999% of the user, but it will always be good to have a clean code.
Thanks
/

/Xiaohu/
/

History

#1 Updated by Rossen Apostolov almost 8 years ago

  • Target version changed from 4.5.5 to 4.5.6

#2 Updated by Berk Hess over 7 years ago

  • Status changed from New to Closed

The code below is in check_cm_grp, so the term you are missing is added later.
It has to be done this way to avoid multiple reductions over all MPI ranks (expensive).

/* Subtract the center of mass contribution from the inertia 
       * tensor (this is not as trivial as it seems, but due to
       * some cancellation we can still do it, even in parallel).
*/
clear_mat(Icm);
update_tensor(vcm->group_x[g],tm,Icm);
m_sub(vcm->group_i[g],Icm,vcm->group_i[g]);

Also available in: Atom PDF