I'm calculating normal vibrational modes in a large molecular system. My goal is to obtain, for each normal mode, the vibrational frequency, the list of displacement vectors and the reduced mass.
I'm using the nmodes() routine from the Amber Molecular Dynamics suite, which allows me to get directly:
- a file (in vecs format) that contains the frequencies of the normal modes of my system and information from which I can easily get the displacement vectors,
- another file which contains all the masses for all the atoms in my system.
I would like to calculate the reduced mass of each normal mode using this data. I already asked in the Amber MD mailing list about this, and obtained this answer, pasted below, which suggests that this should be a straightforward calculation:
Amber doesn't compute reduced masses for normal modes. I guess you
would have to write a script yourself to do this: the masses are in the
prmtop file, and the eigenvectors are in the vecs file.
For my problem, I need precisely the concept of reduced mass (for polyatomic systems) as described in the Gaussian manual. There I do find a section Calculate reduced mass, force constants and cartesian displacements but I don't really follow it. It seems to work with the inertia tensor which I think I would be able to calculate, but I'm not sure in the details of how to use this. I do read that as a consequence of employing elements of the mass-weighted inertia matrix:
the sum of the squares of the cartesian displacements is 1
I also looked the paper Van Vlijmen, H. W. T.; Karplus, M., “Analysis of Calculated Normal Modes of a Set of Native and Partially Unfolded Proteins,” J. Phys. Chem. B, 1999, 103(15), 3009-3021, where the equations 1 and 2 seemed to be appliable, but I'm still not sure how:
$|\nabla^2E-\lambda M| = 0$ (1)
$|M^{-1/2}(\nabla^2E)M^{-1/2}-\lambda I|=0$ (2)
I come from a biochemistry background, so I'm not familiarised with normal mode calculations. Can anyone tell me wich will be the equation to use in this case? All help is appreciated, thanks.
Best Answer
It's been a long time since I dealt with normal modes (and it was on a completely different system), so I hope I won't get this completely wrong :-)
First of all, the two equations you write at the end of your question is essentially the same equation written differently. By solving Eq.(2) you get eigenvalues (which are linked to the frequencies) and eigenvectors (which are linked to the atom displacements, weighted by their mass). However, I had a look at the Amber docs and it looks like they output non-mass-weighted quantities. In this case you can use these displacements to compute the reduced mass by generalising the well-known relation
$$ \frac{1}{\mu} = \frac{1}{m_1} + \frac{1}{m_2} + \cdots \qquad (1) $$
to a case where the mass of each atom is weighted by its displacements. Let's take the $k$-th normal mode, which is associated to the eigenvector $\mathbf{E}_k$. The latter should be normalised with the sum of the squares, that is, by imposing $\sum_i E_{i,k}^2 = A$, where $A$ is a constant and $i$ runs over all the degrees of freedom ($x$, $y$ and $z$ coordinates for each atom). The constant $A$ sort of sets the coordinate system you want to use (cfr. Gaussian's docs). For example, Gaussian uses $A = 1$. Regardless of the value of $A$, we can always define our reduced mass by generalising Eq. (1) and write:
$$ \frac{1}{\mu_k} = \sum_i \frac{E_{i,k}^2}{m_i} \qquad (2) $$
where $m_i$ is the mass associated to the $i$-th degree of freedom.
Let's try to understand whether Eq.(2) makes sense with two simple examples: