Radius of gyration from gyration tensor

[bemde42]

Hello fellow amepers,

Question regarding the gyration tensor function and radius of gyration:
I would expect that the sum of the eigenvalues of the gyration tensor equals the radius of gyration squared.
However, using the functions amep.cluster.radius_of_gyration(...) and amep.cluster.gyration_tensor(...) do lead to different results in that manner.

Does anyone know what is going on or how this can be explained?

I'm happy to supply some example data with a mwe if necessary.

[kay-ro]

Hi Anna,
thanks for your message.
@aritra-mukhopadhyay implemented the gyration_tensor functionality. we're looking into that and will get back to you as soon as possible. :)
Best,
Kay-Ro

[kay-ro]

Hi @bemde42 ,
can you supply a mwe?
Best,
Kay-Ro

[bemde42]

gyrationstuffex.zip

Code:

#get biggest clusters
def biggest_clusters(clusters,ids, cut):
    sizes = amep.cluster.sizes(clusters)
    bigclustersizes = [x for x in sizes if x > cut]
    bigclusters = clusters[:len(bigclustersizes)]
    bigclusterids = ids[:len(bigclustersizes)]
    return bigclusters, bigclusterids, bigclustersizes

sc, scids, scs = biggest_clusters(clusters, ids, cut)

def cluster_geometry(clusters):
    #gyration tensor, geometric center, gyration radius
    gc = amep.cluster.geometric_center(coords, box, pbc=True, clusters=clusters)
    gt = amep.cluster.gyration_tensor(coords, box, masses, pbc=True, clusters=clusters)
    rg = amep.cluster.radius_of_gyration(coords, box, masses, pbc=True, clusters=clusters)
    return gc, gt, rg

def cluster_roundness(clusters):
    # sanity check radius of gyration first
    gcs,gts,rgs = cluster_geometry(clusters)
    print("how many clusters",len(clusters))
    for i, _ in enumerate(clusters):
        print("gyration tensor\n",gts[i])
        ev = np.linalg.eigvalsh(gts[i])
        print("ev", ev)
        print("sum ev",np.sum(ev))
        print("sqrt(sum(ev))", np.sqrt(np.sum(ev)))
        print("rgs", rgs[i])
        print("sqrt(sum(ev))-rgs",np.sqrt(np.sum(ev))-rgs[i])
    return("done")   

Example-output for last frame in trajectory (first 8 clusters, cutoff set to zero):```

gyration tensor
[[19205.38266726 10240.5301499 0. ]
[10240.5301499 20255.21533317 0. ]
[ 0. 0. 0. ]]
ev [ 0. 9476.32440947 29984.27359095]
sum ev 39460.59800042739
sqrt(sum(ev)) 198.64691792330277
rgs 26.24259484545399
sqrt(sum(ev))-rgs 172.40432307784877
gyration tensor
[[26284.94515836 11510.56577301 0. ]
[11510.56577301 21563.7030148 0. ]
[ 0. 0. 0. ]]
ev [ 0. 12174.18977843 35674.45839472]
sum ev 47848.64817315752
sqrt(sum(ev)) 218.74333858007543
rgs 24.56543289451019
sqrt(sum(ev))-rgs 194.17790568556524
gyration tensor
[[24253.79650475 14325.45322042 0. ]
[14325.45322042 24684.37515868 0. ]
[ 0. 0. 0. ]]
ev [ 0. 10142.01497041 38796.15669302]
sum ev 48938.171663428824
sqrt(sum(ev)) 221.2197361526065
rgs 22.991205092706306
sqrt(sum(ev))-rgs 198.22853105990018
gyration tensor
[[21666.24422751 11962.13894285 0. ]
[11962.13894285 22927.42272653 0. ]
[ 0. 0. 0. ]]
ev [ 0. 10318.08517473 34275.58177931]
sum ev 44593.666954045
sqrt(sum(ev)) 211.17212636625365
rgs 17.661072914566905
sqrt(sum(ev))-rgs 193.51105345168673
gyration tensor
[[21816.68124992 9240.00351879 0. ]
[ 9240.00351879 22648.75597964 0. ]
[ 0. 0. 0. ]]
ev [ 0. 12983.35365733 31482.08357223]
sum ev 44465.437229560855
sqrt(sum(ev)) 210.86829356155195
rgs 13.608387328688766
sqrt(sum(ev))-rgs 197.2599062328632
gyration tensor
[[21445.7562812 12950.53001716 0. ]
[12950.53001716 23807.69038417 0. ]
[ 0. 0. 0. ]]
ev [ 0. 9622.45822842 35630.98843696]
sum ev 45253.44666537577
sqrt(sum(ev)) 212.7285751030542
rgs 14.175698356027354
sqrt(sum(ev))-rgs 198.55287674702686

[kay-ro]

for future reference: please make a working miminal example out of the code.

#need masses for gyration tensor
mass = 5e-2
masses = np.empty(traj[0].n(), dtype=np.double)
masses.fill(mass)

coords=traj[-1].coords()
box_boundary=traj[-1].box
clusters, ids=amep.cluster.identify(coords, box_boundary)

def cluster_geometry(coords, clusters, box_boundary, masses):
    #gyration tensor, geometric center, gyration radius
    gc = amep.cluster.geometric_center(coords, box_boundary, pbc=True, clusters=clusters)
    gt = amep.cluster.gyration_tensor(coords, box_boundary, masses, pbc=True, clusters=clusters)
    rg = amep.cluster.radius_of_gyration(coords, box_boundary, masses, pbc=True, clusters=clusters)
    return gc, gt, rg

def cluster_roundness(coords, clusters, box_boundary, masses):
    # sanity check radius of gyration first
    gcs,gts,rgs = cluster_geometry(coords,clusters,box_boundary, masses)
    print("how many clusters",len(clusters))
    for i, _ in enumerate(clusters):
        print(i)
#         print("gyration tensor\n",gts[i])
        ev = np.linalg.eigvalsh(gts[i])
#         print("ev", ev)
#         print("sum ev",np.sum(ev))
        print("sqrt(sum(ev))", np.sqrt(np.sum(ev))*0.05)
        print("rgs", rgs[i])
#         print("sqrt(sum(ev))-rgs",np.sqrt(np.sum(ev))-rgs[i])
cluster_roundness(coords, clusters[:3],box_boundary, masses)

[kay-ro]

there seems to be a problem in the distance calculation. this will be elevated to an issue.

[kay-ro]

issue #78 created

[kay-ro]

fixed with issues #78 #79. will be published with amep version 1.1.0