Topic Options
#22685 - 10/31/09 11:31 PM Center of mass bug in misc/quicka.src
slaw Offline
Forum Member

Registered: 09/11/06
Posts: 146
Hi,

I found a bug in the way that quicka.src is calculating the center of mass for a given atom selection (around line 400). Even though the "MASS" keyword is used, it is not actually mass weighting the coordinates. It should calculate the following:

SUM (coordinate x mass) / SUM (mass)

Instead, it is doing:

SUM (coordinate) / SUM (mass)

almost as if doing a center-of-geometry and center-of-mass hybrid.

This is fixed by replacing the content around line 400 with:

IF (QMASS) THEN
XT=XT+X(I)*AM
YT=YT+Y(I)*AM
ZT=ZT+Z(I)*AM
ELSE
XT=XT+X(I)
YT=YT+Y(I)
ZT=ZT+Z(I)
ENDIF


Sean

Top
#26123 - 01/03/11 11:24 AM Re: Center of mass bug in misc/quicka.src [Re: slaw]
patsalov Offline
Forum Member

Registered: 02/12/07
Posts: 12
This is still an issue with c35b3:

Here's a trivial example, calculating the center of mass and geometric center of a single atom (which should be identical)

CHARMM> quick sele segid A .and. resid 1 .and. type CA end
SELRPN> 1 atoms have been selected out of 3410
QUICKA: Center of GEOM of 1 selected atoms is at: 4.17100 13.68000 13.99100


CHARMM> quick sele segid A .and. resid 1 .and. type CA end mass
SELRPN> 1 atoms have been selected out of 3410
QUICKA: Center of MASS of 1 selected atoms is at: 0.34727 1.13896 1.16485

The coordinates have been divided by 12.01, the mass of CA. : (

Please, developers, incorporate this trivial bugfix into the codebase!

Top
#26124 - 01/03/11 12:05 PM Re: Center of mass bug in misc/quicka.src [Re: patsalov]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4711
Loc: ~ 59N, 15E
It is in there, in the next release (c36). The fix is however easy to implement yourself in routine GETNEXTP in source file source/misc/quicka.src. In two places you have to multiply x,y,z coordinates by AM in a sum, like this:
ENDIF
XT=XT+XCOMP(I)*AM
YT=YT+YCOMP(I)*AM
ZT=ZT+ZCOMP(I)*AM
ELSE
IF(.NOT.INITIA(I,X,Y,Z)) THEN
IF(PRNLEV > 2) THEN
CALL ATOMID(I,SIDDN,RIDDN,RESDN,ACDN)
WRITE(OUTU,35) I, &
SIDDN(1:idleng),RIDDN(1:idleng), &
RESDN(1:idleng),ACDN(1:idleng)
ENDIF
BAD=.TRUE.
RETURN
ENDIF
XT=XT+X(I)*AM
YT=YT+Y(I)*AM
ZT=ZT+Z(I)*AM

You may also be able to use a q&d workaround:
coor copy comp ! save
scalar x prod mass
scalar y prod mass
scalar z prod mass
quick .... !NO MASS KEYWORD HERE
coor copy !restore
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#27709 - 06/21/11 11:04 AM Re: Center of mass bug in misc/quicka.src [Re: lennart]
Randy Bin Lin Offline
Forum Member

Registered: 03/03/10
Posts: 2
Loc: Houston, TX | Baltimore, MD
As of c36a6, the issue above is fixed:

Input:
! verify that mass option is working
! for one atom gom and com should be the same
quick select bynu 1 end
quick select bynu 1 end mass

Output:

CHARMM> quick select bynu 1 end
SELRPN> 1 atoms have been selected out of 7953
QUICKA: Center of GEOM of 1 selected atoms is at: 0.67892 0.87544 1.51325

CHARMM> quick select bynu 1 end mass
SELRPN> 1 atoms have been selected out of 7953
QUICKA: Center of MASS of 1 selected atoms is at: 0.67892 0.87544 1.51325
_________________________
I think therefore I am

Top
#27710 - 06/21/11 11:26 AM Re: Center of mass bug in misc/quicka.src [Re: Randy Bin Lin]
lennart Online   content

Forum Member

Registered: 09/25/03
Posts: 4711
Loc: ~ 59N, 15E
?
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top

Moderator:  BRBrooks, bucknerj, chmgr, John Legato