Bug: eigenvectors of non massweighted Hessian matrix
Moderators: Global Moderator, Moderator

 Newbie
 Posts: 6
 Joined: Wed Nov 23, 2016 6:55 pm
 License Nr.: 458
Bug: eigenvectors of non massweighted Hessian matrix
Dear all,
There is a bug in the routine to calculate the eigenvectors for the non massweighted 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 zdirection. 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) nonorthogonal, (2) nonnormalized, 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 massweighted 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 massweigthed 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 nonmwHm 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 massweighted) 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 massweighted 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 massweighted 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 massweighted 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.
There is a bug in the routine to calculate the eigenvectors for the non massweighted 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 zdirection. 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) nonorthogonal, (2) nonnormalized, 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 massweighted 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 massweigthed 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 nonmwHm 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 massweighted) 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 massweighted 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 massweighted 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 massweighted 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.

 Administrator
 Posts: 2921
 Joined: Tue Aug 03, 2004 8:18 am
 License Nr.: 458
Re: Bug: eigenvectors of non massweighted Hessian matrix
The Hessian matrix is calculated, but printed is only dynamical matrix.
Cf. "Eigenvectors and eigenvalues of the dynamical matrix".
Standard units are mass weighted cartesian coordinates,
using NWRITE=3 also in cartesian coordinates.
Cf. "Eigenvectors and eigenvalues of the dynamical matrix".
Standard units are mass weighted cartesian coordinates,
using NWRITE=3 also in cartesian coordinates.

 Newbie
 Posts: 6
 Joined: Wed Nov 23, 2016 6:55 pm
 License Nr.: 458
Re: Bug: eigenvectors of non massweighted Hessian matrix
Hi there!
Are you planning to fix this bug anytime soon?
Kind regards,
Rodrigo García Muelas.
Are you planning to fix this bug anytime soon?
Kind regards,
Rodrigo García Muelas.

 Administrator
 Posts: 2921
 Joined: Tue Aug 03, 2004 8:18 am
 License Nr.: 458
Re: Bug: eigenvectors of non massweighted Hessian matrix
The starting point of your reasoning about the calculation
of vibrational eigenvalues/eigenvectors is:
"The eigenvectors for the non mwHm matrix reported in the outcar file are:
(1) nonorthogonal, (2) nonnormalized, ..."
There are no eigenvectors of the non mwHm reported in the outcar file.
Just "Eigenvectors and eigenvalues of the dynamical matrix"
Note that the manual speaks about the calculation of the Hessian matrix
https://cms.mpi.univie.ac.at/vasp/vasp/ ... ION_6.html
but it is in outcar not written. What is written is the dynamical matrix.
of vibrational eigenvalues/eigenvectors is:
"The eigenvectors for the non mwHm matrix reported in the outcar file are:
(1) nonorthogonal, (2) nonnormalized, ..."
There are no eigenvectors of the non mwHm reported in the outcar file.
Just "Eigenvectors and eigenvalues of the dynamical matrix"
Note that the manual speaks about the calculation of the Hessian matrix
https://cms.mpi.univie.ac.at/vasp/vasp/ ... ION_6.html
but it is in outcar not written. What is written is the dynamical matrix.

 Newbie
 Posts: 6
 Joined: Wed Nov 23, 2016 6:55 pm
 License Nr.: 458
Re: Bug: eigenvectors of non massweighted Hessian matrix
Thank you for your reply.
I checked the Manual. Besides the eigenvectors/eigenvalues of the dynamical matrix, there is another set of eigenvectors printed under the flag "Eigenvectors after division by SQRT(mass)" (IBRION=5; NWRITE=3). Those vectors have no physical meaning at all; in fact, they are not even orthonormal. The purpose of these vectors was to represent the eigenvectors of the Hessian matrix (not the dynamical matrix), which are needed for the Improved Dimer Method (IDM). (IBRION=44). This is very clear in the VASP manual:
https://cms.mpi.univie.ac.at/vasp/vasp/ ... ample.html
For the IDM you actually need the eigenvectors of the Hessian Matrix (nwwHm), which are not reported anywhere. The eigenvectors of the Hessian matrix represent the curvature of the Potential Energy Surface, so its physical meaning is intensively used in the papers of Henkelman & Jonsson and Heyden et al.
You may say that, as VASP normalizes the input vector for IBRION=44, you can still use the unstable mode reported under "Eigenvectors after division by SQRT(mass)". However, this is not right, as the direction of the unstable mode can be rotated from the "true" direction. Fortunately, the Hessian Matrix is reported in the OUTCAR file under the flag "SECOND DERIVATIVES (NOT SYMMETRIZED)". After realize this, now I just take the Hessian matrix, symmeterize it, diagonalize it, and use it as the input of all my IDM. The IDM algorithm are far more stable with that input.
Below these lines, I attach a fragment of the OUTCAR file. Please notice that the direction of the first eigenvector under "Eigenvectors after division by SQRT(mass)" (0.218;+0.164) is rotated 8.13º with respect to the real value of the PES: (0.707,+0.707). This was a test example with a CO molecule in gas phase, but this problem also happens for much bigger systems.
As a summary, it will be desirable that VASP calculates the eigenvectors of the Hessian matrix for IBRION=5 and NWRITE=3, and print these results instead of the "Eigenvectors after division by SQRT(mass)". So you shall make two diagonalizations: For the dynamic matrix (used in frequencies, etc) and for the Hessian matrix (used for IDM and any method requiring the Potential Energy Surface).
Best,
Rodrigo.
SECOND DERIVATIVES (NOT SYMMETRIZED)

1Z 2Z
1Z 114.005552 114.005552
2Z 113.795271 113.795271
Eigenvectors and eigenvalues of the dynamical matrix

1 f = 63.702530 THz 400.254804 2PiTHz 2124.887630 cm1 263.452584 meV
X Y Z dx dy dz
5.000000 5.000000 5.571739 0 0 0.756080
5.000000 5.000000 4.428261 0 0 0.654480
2 f/i= 0.029098 THz 0.182829 2PiTHz 0.970612 cm1 0.120341 meV
X Y Z dx dy dz
5.000000 5.000000 5.571739 0 0 0.654480
5.000000 5.000000 4.428261 0 0 0.756080
Eigenvectors after division by SQRT(mass)
Eigenvectors and eigenvalues of the dynamical matrix

1 f = 63.702530 THz 400.254804 2PiTHz 2124.887630 cm1 263.452584 meV
X Y Z dx dy dz
5.000000 5.000000 5.571739 0 0 0.218161
5.000000 5.000000 4.428261 0 0 0.163620
2 f/i= 0.029098 THz 0.182829 2PiTHz 0.970612 cm1 0.120341 meV
X Y Z dx dy dz
5.000000 5.000000 5.571739 0 0 0.188845
5.000000 5.000000 4.428261 0 0 0.189020
I checked the Manual. Besides the eigenvectors/eigenvalues of the dynamical matrix, there is another set of eigenvectors printed under the flag "Eigenvectors after division by SQRT(mass)" (IBRION=5; NWRITE=3). Those vectors have no physical meaning at all; in fact, they are not even orthonormal. The purpose of these vectors was to represent the eigenvectors of the Hessian matrix (not the dynamical matrix), which are needed for the Improved Dimer Method (IDM). (IBRION=44). This is very clear in the VASP manual:
https://cms.mpi.univie.ac.at/vasp/vasp/ ... ample.html
For the IDM you actually need the eigenvectors of the Hessian Matrix (nwwHm), which are not reported anywhere. The eigenvectors of the Hessian matrix represent the curvature of the Potential Energy Surface, so its physical meaning is intensively used in the papers of Henkelman & Jonsson and Heyden et al.
You may say that, as VASP normalizes the input vector for IBRION=44, you can still use the unstable mode reported under "Eigenvectors after division by SQRT(mass)". However, this is not right, as the direction of the unstable mode can be rotated from the "true" direction. Fortunately, the Hessian Matrix is reported in the OUTCAR file under the flag "SECOND DERIVATIVES (NOT SYMMETRIZED)". After realize this, now I just take the Hessian matrix, symmeterize it, diagonalize it, and use it as the input of all my IDM. The IDM algorithm are far more stable with that input.
Below these lines, I attach a fragment of the OUTCAR file. Please notice that the direction of the first eigenvector under "Eigenvectors after division by SQRT(mass)" (0.218;+0.164) is rotated 8.13º with respect to the real value of the PES: (0.707,+0.707). This was a test example with a CO molecule in gas phase, but this problem also happens for much bigger systems.
As a summary, it will be desirable that VASP calculates the eigenvectors of the Hessian matrix for IBRION=5 and NWRITE=3, and print these results instead of the "Eigenvectors after division by SQRT(mass)". So you shall make two diagonalizations: For the dynamic matrix (used in frequencies, etc) and for the Hessian matrix (used for IDM and any method requiring the Potential Energy Surface).
Best,
Rodrigo.
SECOND DERIVATIVES (NOT SYMMETRIZED)

1Z 2Z
1Z 114.005552 114.005552
2Z 113.795271 113.795271
Eigenvectors and eigenvalues of the dynamical matrix

1 f = 63.702530 THz 400.254804 2PiTHz 2124.887630 cm1 263.452584 meV
X Y Z dx dy dz
5.000000 5.000000 5.571739 0 0 0.756080
5.000000 5.000000 4.428261 0 0 0.654480
2 f/i= 0.029098 THz 0.182829 2PiTHz 0.970612 cm1 0.120341 meV
X Y Z dx dy dz
5.000000 5.000000 5.571739 0 0 0.654480
5.000000 5.000000 4.428261 0 0 0.756080
Eigenvectors after division by SQRT(mass)
Eigenvectors and eigenvalues of the dynamical matrix

1 f = 63.702530 THz 400.254804 2PiTHz 2124.887630 cm1 263.452584 meV
X Y Z dx dy dz
5.000000 5.000000 5.571739 0 0 0.218161
5.000000 5.000000 4.428261 0 0 0.163620
2 f/i= 0.029098 THz 0.182829 2PiTHz 0.970612 cm1 0.120341 meV
X Y Z dx dy dz
5.000000 5.000000 5.571739 0 0 0.188845
5.000000 5.000000 4.428261 0 0 0.189020

 Administrator
 Posts: 2921
 Joined: Tue Aug 03, 2004 8:18 am
 License Nr.: 458
Re: Bug: eigenvectors of non massweighted Hessian matrix
The developers promised to modify the code in the way that IBRION=5 and NWRITE=3
will print eigenvalues and eigenvectors of the Hessian matrix.
will print eigenvalues and eigenvectors of the Hessian matrix.

 Newbie
 Posts: 6
 Joined: Wed Nov 23, 2016 6:55 pm
 License Nr.: 458
Re: Bug: eigenvectors of non massweighted Hessian matrix
Thank you very much for your kind support!
Rodrigo.
Rodrigo.