freud cluster center of mass vs normally calculating com

[jude vishnu]

Hi everyone,

I have been using the real position of the particles from box.unwrap(pos,images) to calculate the center of mass . I then wrap it so that it lies inside the box (for binning purposes)

However when I do the same thing with the help of freud.cluster function, I get a different result . Do you know what may be the issue? . I think this is causing the error later in gyration tensor calculations. So can someone please help me out in understanding what am I doing wrong. These are the codes for the 2 approaches:

segment_length=102

overlapping_segments = np.lib.stride_tricks.sliding_window_view(cluster_keys, (segment_length,),axis=1)
 overlapping_segments = overlapping_segments.reshape(-1,segment_length)

1st approach :

realposition = box.unwrap(position,image)

realpos = realposition[overlapping_segments]

realpos = realpos.reshape(-1,segment_length,3)

center_of_mass = np.nanmean(realpos, axis=1,keepdims=True)

comwrap = box.wrap(center_of_mass)

differences = realpos -center_of_mass

gyration_tensor = np.nanmean(np.einsum('ijk,ijl->ijkl', differences, differences), axis=1)

2nd approach:

snap =traj[frame]
position = snap.particles.position
image = snap.particles.image
sel_pos = position[overlapping_segments]
sel_pos = sel_pos.reshape(-1,segment_length,3)
gyration_tensor,comwrap,rg2=Rgfunc(sel_pos,box)

def Rgfunc(points,box):

     cluster_id =                     np.arange(0,len(points)).repeat(segment_length,axis=0).reshape(-1,segment_length)


    system = (box,points.reshape(-1,3))
    clp = freud.cluster.ClusterProperties()
    clp.compute(system,cluster_id.reshape(-1))

    
    return clp.gyrations, clp.centers,(clp.radii_of_gyration)**2

best,

Jude

[Tommy Waltmann]

Hi Jude,

The reason the two approaches give different results is because in approach 1, the points are unwrapped before the calculations are done while in approach 2 they are not. The `ClusterProperties` class will not unwrap the input points. Try unwrapping the points before giving them to the freud `ClusterProperties` class and see if the results from the two different approaches match.

Best,

Tommy

[Joshua Anderson]

Freud's ClusterProperties computes the center of mass accounting for periodic boundary conditions:

https://github.com/glotzerlab/freud/blob/39ea38d1e23874f4e49337c42b7d23a7c9ca6aa1/cpp/box/Box.h#L384-L396

and then locally wraps vectors as needed to compute properties about that center of mass:

https://github.com/glotzerlab/freud/blob/39ea38d1e23874f4e49337c42b7d23a7c9ca6aa1/cpp/cluster/ClusterProperties.cc#L73-L81

This is far more accurate than starting from unwrapped coordinates. Floating point values have only so many digits of accuracy. If all of those digits are used to represent the unwrapped position, then are none left over for accuracy in computed delta r values. See https://en.wikipedia.org/wiki/Catastrophic_cancellation.

------

Joshua A. Anderson, Ph.D.

Research Area Specialist, Chemical Engineering, University of Michigan

-- 
You received this message because you are subscribed to the Google Groups "freud-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to freud-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/freud-users/d91bc33a-684c-47a5-a6e3-610c19be50e8n%40googlegroups.com.

[Bradley Dice]

I tried to reply to this earlier and lost my draft. In addition to what Tommy and Josh said, using `np.mean` (or `np.nanmean` as in your snippet) will just give you the wrong answer. Consider a 1-d system with a box going from (-1, 1] -- inclusive on the right boundary. If you have unwrapped points at -1.1 and 1.1, your wrapped points would be at -0.9 and +0.9. Their mean will be zero. But in the periodic system, the proper center of mass should actually be on the boundary at 1 (not zero). Even without considering the numerical issues of large floats, you have to use the right formula to get this result, as mentioned in the freud docs. There are lots of cases where this has a real influence on the center of mass calculation, not just in "constructed" cases. See links below:

https://freud.readthedocs.io/en/latest/modules/box.html#freud.box.Box.center_of_mass

https://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions