Anomalous values of the first six normal modes

GROMACS version: 2021.1
GROMACS modification: No
Hi. I am simulating a lysozyme molecule in explicit flexible water molecules with Gromacs. I filter the trajectory to eliminate the roto-translational movement of the protein; then I compute the Hessian matrix from the PDB structure minimized in vacuo, and the normal modes via the nmeig module. The first six eigenvalues are very small (but not 0). When I project the trajectory on the normal modes, I expect the projections on the first six nm, corresponding to the roto-translational degrees of freedom, to be negligible. The opposite happens: the values of the first six projections are two orders of magnitude larger than those of the projections on the remaining eigenvectors. Moreover, the anomalous values are not stable. If I slightly change the normal modes, by changing the energy minimization path of the PDB structure, the anomalous projections change a lot, but remain very large. Does anybody have an explication for these results?

Do you use GROMACS compiled with double precision?

Thank you for your attention. Yes, I use double precision, absolutely.

Hmm, I never calculated normal modes for molecules this large (and never with nmeig), but maybe I can help anyway.

First, what confuses me is that you

filter the trajectory to eliminate the roto-translational movement of the protein

Is that why you expect there to be no rot-trans normal modes?

Second, I think this could all be just a numerical issue (not perfectly minimized for example), that’s why i asked about double precision.
It is not unexpected to me, that the eigenvector of one of the first six normal modes should have large values. Think of a harmonic oscillator with a very small frequency but constant energy (1/2 k T). It will have large amplitude. Are the eigenvectors of nmeig normalized? Do your first six normal modes look at least a bit like translation/rotation (many atoms moving in the same direction)?

(i) I filter the trajectory to eliminate the roto-translational movement of the protein because I am interested in analysing the energy exchange among internal vibrational degrees of freedom of the protein. As a matter of fact, the unfiltered and filtered 10 ns trajectories I am using are almost equal, the roto-translational motion is negligible (the translational motion had already been compensated for during the computation of the trajectory).

(ii) The eigenvectors produced by nmeig are (almost) orthogonal, but NOT normalized: their norm is in the range 0.37 - 0.40 u^0.5 nm, because nmeig computes them mass-weighted in default. Also the first six eigenvectors have norms in that range. I tried to have the eigenvectors computed non mass-weighted, with the option -nom, but apparently nmeig assumes always the option -m . From the Gromacs manual, nmeig module: “-m Divide elements of Hessian by product of sqrt(mass) of involved atoms prior to diagonalization. This should be used for ‘Normal Modes’ analysis”.

(iii) The eigenvalues of the first six eigenvectors are small. The first ten eigenvalues are as follows:

 1      -0.0392398
 2      -0.0207044
 3      -0.0118039
 4     -0.00603985
 5     -0.00274185
 6        0.019326
 7         1.20237
 8         1.67662
 9         1.92187
10         3.05489

(iv) The projections of the trajectory on the first six eigenvectors are abnormally large, but the amplitude of their oscillation is smaller than that of the projection on the following “well-behaved” normal modes. Unfortunately, I am told that as a “new user” I cannot upload attachments (?). If you give me your email I will send you two figures with a sample of the projections on eigenvectors 1 -10. Thank you.

Ok, I did not really get the first time that you are confused about the projections of the trajectory on the first six normal modes. I get now that you expect those to be zero, when the trajectory was filtered before.

Seeing the eigenvalues not being precisely zero, my best guess is that your structure is not fully minimized. I guess it is pretty hard to properly minimize a big protein. Maybe you could try minimizing with more strict settings and see if things improve.

You are welcome to send me figures.

(email removed again)

We had some more conversation over mail. To summarize:

  1. The eigenvalues of the normal modes not being zero comes likely from the molecule (protein) being large and numerical inaccuracies arising.
  2. To visualize what a normal mode represent it might be useful to make an animation from the eigenvector (see e.g. on this page, using e.g. this VMD tool)

Yes, the energy minimization is not perfect; I made various attempts to improve it, but the final result is more or less the same. As a consequence the eigenvalues of the first six normal modes are small, but not zero. Moreover, a still visualization of these six modes at T = 0 yields partial and fragmented structures of the protein. This explains the anomalous high average value of the projection of the trajectory on each of these six modes, as a very strong distortion of the fragmented structures is needed to approach the complete running structure of the protein. Visualizing these six projections with VMD produces still frames. My conclusion: the first six modes are not computed properly; as I am not intersted in the rotation and translation of the protein I will discard them.