There is a bug in the routine to calculate the eigenvectors for the non mass-weighted Hessian Matrix (non mwHm) (IBRION=5, NWRITE=3). This routine is particularly important to get the input file of the Improved Dimer Method (IDM). The resulting input file would have the dimer direction rotated from the correct value, thus reducing its stability. I have found this behaviour in the 3.2.X, 5.3.X and 5.4.1 versions. I describe herein how to duplicate the bug with a simple model system.
# # # # # # # # # # B U G D E S C R I P T I O N # # # # # # # # # #
I put CO in a 10³ ų box, oriented in z-direction. I fully relaxed the molecule and then calculated the frequencies with NWRITE=3, to get the non mwHm. INCAR, POSCAR, and part of the OUTCAR, are attached at the end.
The eigenvectors for the non mwHm matrix reported in the outcar file are: (1) non-orthogonal, (2) non-normalized, and (3) displaced from the (+0.707,+0.707) direction which would be the right solution as described below.
Doing some reverse engineering, I found that VASP does the following (IBRION=5, NWRITE=3):
(1) Get the non mwHm, let us call it H.
(2) Symmetrize it.
H=0.5Ĥ + 0.5Ĥ^t (eq. 2)
(3) Multiply by M^(-1/2) at the left and right sides of "H" to get the mass-weighted Hessian matrix. Let us call it "W".
W = M^{-1/2} H M^{-1/2} (eq. 1)
(4) Diagonalize "W".
(5) Get the eigenvalues/eigenvectors of "H" from the eigenvalues/eigenvectors of "W" by eq. (3). "P" are the eigenvectors of "H" and "Q" the eigenvectors of "W".
P = M^{-1/2} Q (eq. 3)
Yet, the 5th step is not correct. The eigenvalues/eigenvectors of "H" and "W" are not proportional in any way, as it will be shown below. There is no relation between the eigenvalues/eigenvectors of "H" and "W". Below I prove that eq. 3, which is used by VASP, is wrong. The correct procedure would be to get the eigenvalues and eigenvectors from "H" and "W" separately. So the bux fix will be to compute the eigenvalues of "H" separately and not to use the spurious eq. 3.
# # # # # # # # # # W H A T I T S H O U L D B E # # # # # # # # # #
The non mass-weigthed Potential Energy Surface (PES) of this problem has two dimensions, corresponding to the "z" coordinates of C and O. If the harmonic approximation is used, the PES resembles a parabolic cylinder. One of the eigenvectors for the non-mwHm is zero, corresponding to the translation of the CO molecule. The other eigenvalue is related with the spring constant of C and O coming closer at the same pace, which describes the maximum increase in potential energy. This is, the eigenvectors should be exactly:
(+r,-r) for the atoms getting closer to each other, and
(+r,+r) for the translation (eigenvalue: 0)
Where r=sqrt(0.5) (either positive or negative)
As a side note, the eigenvectors of the non mwHm correctly describe the PES that is used on the IDM. The eigenvalues for the non mwHm should have the same units of "H", this is, in eV/Ų.
Now we take the (non mass-weighted) Hessian matrix, just after the flag "SECOND DERIVATIVES (NOT SYMMETRIZED)" in the OUTCAR file. We symmetrize it, and diagonalize it by hand, we would get the eigenvectors/eigenvalues that were expected theoretically:
(+0.707, -0.707), eigenvalue = 228 eV/Ų (C and O getting closer. Maximum variation in potential energy.)
(+0.707, +0.707), eigenvalue = zero (Translation: the PES is flat along this direction.)
Instead, if we use the VASP results, we get this (see OUTCAR below):
(-0.218161,+0.163620),
(-0.188845,-0.189020),
VASP only reports the eigenvalues of the mwHm. The eigenvectors reported by VASP are not normalized nor orthogonal. The first eigenvector is displaced 10º from the direction of maximum variation in the PES. However, the direction of the second eigenvector is correct.
* * *
On the other hand, the eigenvectors for the mwHm define frequencies. These vectors are not proportional to (+r,±r), because they belong to a different vectorial space, and the atomic masses of C and O are different. To get the mwHm (let us call it "W), we multiply "H" at the right and left side by "M^{-1/2}", equation (1). "M^(-1/2)" is a diagonal matrix that contains the mass of the atoms, to the power -1/2, once for each degree of freedom. The eigenvalues/eigenvectors of "W" are:
(-0.7559,+0.6546), eigenvalue: -16.625 eV/Ų·uma (Vibrational mode: stretching)
(+0.6546,+0.7559), eigenvalue: zero (Kernel. This "vibration" has no physical meaning.)
The 1st and 2nd components are related to the displacements of C and O respectively. During the vibration, the C atom will be displaced more than the O atom. Both numbers are therefore related with equation 3:
0.7559=0.6546*sqrt(16/12)
The results provided by VASP for the mass-weighted Hessian matrix match those expected theoretically:
(-0.756080,+0.654480)
(-0.654480,-0.756080)
As the PES defined by matrix "H" does not belong to the same vectorial space than matrix "W", their eigenvectors are not proportional. The only exception of this rule would be that the masses for all atoms where equal.
# # # # # # # # # # P R O V I N G W R O N G E Q . 3 # # # # # # # # # #
Let "D" and "P" be the eigenvalues/eigenvectors matrices of the non mwHm, "H". Let "E" and "Q" be the eigenvalues/eigenvectors matrices of the mass-weighted Hessian matrix, "W". Because "H" and "W" are symmetric matrices, then:
D = P^{-1} H P (eq. 4a)
E = Q^{-1} W Q (eq. 5a)
The superscript "-1" indicates inverse. Because "H" and "W" are symmetric matrices, "P" and "Q" are orthogonal, then:
P^t P = I (eq. 6)
Q^t Q = I (eq. 7)
The superscript "t" indicates transposition. And eq. 4a and 5a can be rewritten as:
D = P^t H P (eq. 4b)
E = Q^t W Q (eq. 5b)
By definition, the mass-weighted Hessian matrix, "W", can be written as a function of "H":
W = M^{-1/2} H M^{-1/2} (eq. 1)
Where "M^{-1/2}" is a diagonal matrix, whose components are the inverse of the square roots of the masses. Eq. (5) can be rewriten as:
E = Q^t M^{-1/2} H M^{-1/2} Q (eq.
Here it comes the bug. VASP gets "P" from "Q" by doing this:
P = M^{-1/2} Q (eq. 3) # wrong #
Which comes by comparing the terms at the right of "H" in eq. (4b) and (8). However, this is not correct because "E" and "D" contains different eigenvalues and are not equal to each other. Besides, Eq. (4a) requires that "M^{-1/2} Q" were symmetric, to fulfil eq. (6). This condition is also not fulfilled, which can be proved by "reductio ad absurdum". By substituting equation (3) on (6) we get:
Q^t M^{-1/2} = [M^{-1/2} Q]^{-1}
Q^t M^{-1/2} = Q^{-1} [M^{-1/2}]^{-1} [Property: (AB)^{-1} = B^{-1} A^{-1}]
Q^t M^{-1/2} = Q^t [M^{-1/2}]^{-1} [Eq. 7]
Q^t M^{-1/2} ≠ Q^t M^{+1/2} [Applying the definition of M, we reach an "absurdum" ]
This means that the right side of eq. (3) is not orthogonal, so it is not correct to get the eigenvectors of "H" (i.e., "P") from the eigenvectors of "W" (i.e. "Q"), and both matrices need to be diagonalized separately.
