How to calculate moment of inertia for coarse-grain model

[Vinnarasi Saravanan]

Hi,

I want to calculate principal moment of inertia for coarse-grain (Martini FF in NAMD) protein system. Is it possible to calculate by MDAnalysis?

Thanks

[jona...@barnoud.net]

Hi,

The "AtomGroup" class has a "moment_of_inertia" method <https://docs.mdanalysis.org/stable/documentation_pages/core/groups.html?highlight=atomgroup#MDAnalysis.core.groups.AtomGroup.moment_of_inertia>. Make sure the masses are correct, you can check them with the "masses" attribute of the "AtomGroup", pay special attention to the virtual sites. The masses are the most likely thing to get wrong when reading a coarse-grained structure in MDAnalysis. If the masses are read correctly, or if you fixed them in your universe, it should not matter that the universe is about a coarse-grained simulation.

Cheers,

Jonathan

[Vinnarasi Saravanan]

Hi,

I verified the "masses" of my coarse-grained protein system, it seems correct. According to your suggestion, I tried to use MDAnalysis API functions principal_axes and moment_of_inertia to calculate these matrices for a group of selected atoms for my coarse-grained system using the following source code

import MDAnalysis
import numpy as np
u = MDAnalysis.Universe('CG.psf', 'CG.dcd')
CA = u.select_atoms("protein and name CA")
I = CA.moment_of_inertia()
UT = CA.principal_axes()
# transpose the row-vector layout UT = [p1, p2, p3]
U = UT.T
# test that U diagonalizes I
Lambda = U.T.dot(I.dot(U))
print(Lambda)
# check that it is diagonal (to machine precision)
print(np.allclose(Lambda - np.diag(np.diagonal(Lambda)), 0))

My coarse-grained (Martini force field using NAMD) system is an icosahedral capsid protein and it has 181234 atoms. Generally, protein is represented as a bead type according to the Martini force field. Hence, I wrote "protein" in the selection line instead of ("protein and name CA"). Is that correct? Maybe I'm wrong. Then, I got the "FALSE" statement after verifying "print(np.allclose(Lambda - np.diag(np.diagonal(Lambda)), 0))". Could you help me through this problem? 

Thanks

 

With regards

Vinnarasi Saravanan,

Postdoctoral fellow,

Chemical Engineering department,

Indian Institute of Technology Bombay,

Mumbai, India.

--
You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/ir33ScV3hvc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/c0cd654b-3297-465d-9060-b6bba93c8887n%40googlegroups.com.

[Oliver Beckstein]

Hi Vinnarasi Saravanan,

Can you show us the actual code that you ran (with the correct selection string)?

How many atoms are in group CA?

Did you look at Lambda, I, U? Print them all and show them in such a way that it is obvious to us what is what.

Oliver

You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAA81hZ94yk7kgQX%3DaRXF3Qd3W44qAJ2HP6KoHz3LbtvMhbFxEQ%40mail.gmail.com.

--

Oliver Beckstein (he/his/him)

orbe...@gmail.com