Matrix diagonalization algorithm

[Spathi]

TL;DR Summary

I need to compute the vibrational frequencies of a molecule when the matrix of force constants (second derivative of the energy by the Cartesian coordinates) is provided. For such computation, this matrix must be diagonalized.

I need to compute the vibrational frequencies of a molecule when the matrix of force constants (second derivative of the energy by the Cartesian coordinates) is provided. For such computation, this matrix must be diagonalized. Here is an example of a matrix which must be diagonalized:

-0.0001 0.0000 0.0000 0.0001 0.0000 0.0000 0.0001 0.0000 0.0000

0.0000 0.3533 0.0000 0.0000 -0.1766 0.1068 0.0000 -0.1766 -0.1068

0.0000 0.0000 0.3671 0.0000 0.1181 -0.1836 0.0000 -0.1181 -0.1836

0.0001 0.0000 0.0000 -0.0001 0.0000 0.0000 0.0000 0.0000 0.0000

0.0000 -0.1766 0.1181 0.0000 0.2811 -0.1124 0.0000 -0.1044 -0.0056

0.0000 0.1068 -0.1836 0.0000 -0.1124 0.1238 0.0000 0.0056 0.0597

0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0000

0.0000 -0.1766 -0.1181 0.0000 -0.1044 0.0056 0.0000 0.2811 0.1124

0.0000 -0.1068 -0.1836 0.0000 -0.0056 0.0597 0.0000 0.1124 0.1238

Here is the same matrix again, with better precision:

-0.000126656264 3.64590948E-12 -3.87537004E-12 6.33281292E-5 -5.91919046E-12 3.22027804E-12 6.33281284E-5 -5.0934172E-15 -3.82247162E-13

3.64590948E-12 0.353287782 -1.25257096E-12 3.70977256E-12 -0.176643891 0.10680789 6.01837625E-12 -0.176643891 -0.10680789

-3.87537004E-12 -1.25257096E-12 0.367103546 3.02013197E-12 0.118068294 -0.183551773 -1.98117073E-12 -0.118068294 -0.183551773

6.33281292E-5 3.70977256E-12 3.02013197E-12 -7.01014427E-5 1.88779069E-13 -3.06838337E-12 6.77331697E-6 -3.91275366E-12 -2.89242863E-13

-5.91919046E-12 -0.176643891 0.118068294 1.88779069E-13 0.281060184 -0.112438092 -3.91210419E-12 -0.104416293 -0.0056302017

3.22027804E-12 0.10680789 -0.183551773 -3.06838337E-12 -0.112438092 0.123835163 2.88954654E-13 0.0056302017 0.0597166103

6.33281284E-5 6.01837625E-12 -1.98117073E-12 6.77331697E-6 -3.91210419E-12 2.88954654E-13 -7.01014423E-5 1.86304372E-13 3.0670495E-12

-5.0934172E-15 -0.176643891 -0.118068294 -3.91275366E-12 -0.104416293 0.0056302017 1.86304372E-13 0.281060184 0.112438092

-3.82247162E-13 -0.10680789 -0.183551773 -2.89242863E-13 -0.0056302017 0.0597166103 3.0670495E-12 0.112438092 0.123835163

I thought I understood how to diagonalize a matrix. Firstly my program examines the first column, and finds the biggest (by absolute value) item on it. This item is -0.000126656264 (first row). Then my program adds the first row to all other rows so that all items in other rows in column 1 become zero (for example, the first row is added to second row with 3.64590948E-12/0.000126656264 coefficient). After this iteration, we get the following matrix:

-0.000126656264 3.64590948E-12 -3.87537004E-12 6.33281292E-5 -5.91919046E-12 3.22027804E-12 6.33281284E-5 -5.0934172E-15 -3.82247162E-13

0 0.353287782 -1.25257107155586E-12 5.53272721939959E-12 -0.176643891 0.10680789 7.8413308863709E-12 -0.176643891 -0.10680789

0 -1.25257107155586E-12 0.367103546 1.08244703567311E-12 0.118068294 -0.183551773 -3.91885563984886E-12 -0.118068294 -0.183551773

0 5.53272721939959E-12 1.08244703567311E-12 -3.84373794999999E-5 -2.77081603014399E-12 -1.45824442119094E-12 3.84373797700001E-5 -3.9153003684874E-12 -4.80366435549632E-13

0 -0.176643891 0.118068294 -2.77081603014399E-12 0.281060184 -0.112438092 -6.87169925175656E-12 -0.104416293 -0.0056302017

0 0.10680789 -0.183551773 -1.45824442119094E-12 -0.112438092 0.123835163 1.89909358246879E-12 0.0056302017 0.0597166103

0 7.8413308863709E-12 -3.91885563984886E-12 3.84373797700001E-5 -6.87169925175656E-12 1.89909358246879E-12 -3.84373798999999E-5 1.83757663544772E-13 2.87592592986476E-12

0 -0.176643891 -0.118068294 -3.9153003684874E-12 -0.104416293 0.0056302017 1.83757663544772E-13 0.281060184 0.112438092

0 -0.10680789 -0.183551773 -4.80366435549632E-13 -0.0056302017 0.0597166103 2.87592592986476E-12 0.112438092 0.123835163

Then we check the second column. We find the biggest item in rows 2-9: it is 0.353287782 (row 2). Then we add row 2 to rows 1 and 3-9 so that the whole column 2 becomes filled with zeros (for example, we add row 2 to row 1 with multiplication coefficient -3.64590948E-12/0.353287782). After this iteration, the matrix looks like follows:

-0.000126656264 0 -3.87537003998707E-12 6.33281292E-5 -4.09623572E-12 2.11802679178164E-12 6.33281284E-5 1.8178613228E-12 7.20004086218364E-13

0 0.353287782 -1.25257107155586E-12 5.53272721939959E-12 -0.176643891 0.10680789 7.8413308863709E-12 -0.176643891 -0.10680789

0 0 0.367103546 1.08244703569273E-12 0.118068293999374 -0.183551772999621 -3.91885563982106E-12 -0.118068294000626 -0.183551773000379

0 0 1.08244703569273E-12 -3.84373794999999E-5 -4.4524204441942E-15 -3.13092870397103E-12 3.84373797700001E-5 -1.1489367587876E-12 1.19231784723045E-12

0 -1.35525271560688E-20 0.118068293999374 -4.4524204441942E-15 0.1927382385 -0.059034147 -2.95103380857111E-12 -0.1927382385 -0.0590341467

0 0 -0.183551772999621 -3.13092870397103E-12 -0.059034147 0.0915444188885829 -4.71540358008398E-13 0.0590341467 0.0920073544114171

0 0 -3.91885563982106E-12 3.84373797700001E-5 -2.95103380857111E-12 -4.71540358008398E-13 -3.84373798999999E-5 4.10442310673022E-12 5.24655987034195E-12

0 0 -0.118068294000626 -1.1489367587876E-12 -0.1927382385 0.0590341467 4.10442310673022E-12 0.1927382385 0.059034147

0 0 -0.183551773000379 1.19231784723045E-12 -0.0590341467 0.0920073544114171 5.24655987034195E-12 0.059034147 0.0915444188885829

After repeating the steps until the last column, we get this final matrix:

-0.000126656264 -5.63831667531326E-33 -2.81915833765663E-33 0 0 0 0.000126656257814183 -4.17866507795235E-12 -2.44793127062406E-12

0 0.353287782 -7.73422120937375E-21 0 0 0 -2.30568319218387E-9 -0.35328792014169 -1.37915956891113E-7

0 1.03390773611606E-20 0.367103546 0 6.7762635780344E-21 0 3.98006690721581E-9 2.37398846246612E-7 -0.36710330863164

0 -1.48047016107681E-31 -7.40235080538405E-32 3.84373797700001E-5 0 0 -3.84373798999998E-5 1.15339244229211E-12 8.56167022525175E-13

0 -1.35525271560688E-20 -6.7762635780344E-21 0 0.154764968560726 0 -2.04323699675084E-12 -0.154764968560324 2.99373715655519E-10

0 2.62537764105279E-29 1.3126888205264E-29 0 0 0.000231467911417079 5.02067342532649E-12 2.99373715018316E-10 -0.000231467611795766

0 -1.78922510440251E-31 -8.94612552201253E-32 0 0 0 1.40000420028219E-13 6.60441071561968E-18 6.54836315815767E-18

0 -1.35525271560335E-20 -6.77626357801675E-21 0 0 0 6.63578079646145E-18 3.88915968412662E-16 -1.2521817723127E-12

0 2.62156556997463E-29 1.31078278498732E-29 -1.97215226305253E-31 0 0 -6.46286665078964E-18 -1.25295941923224E-12 5.99999604757664E-10

So, the diagonal values are:

-0.000126656264 0.353287782 0.367103546 3.84373797700001E-5 0.154764968560726 0.000231467911417079 1.40000420028219E-13 3.88915968412662E-16 5.99999604757664E-10

Is this correct? I will write further why I suppose that this is wrong.

 

[Arjan82]

What you are describing is the Gaussian elimination algorithm.

Then we check the second column. We find the biggest item in rows 2-9: it is 0.353287782 (row 2). Then we add row 2 to rows 1 and 3-9

Three remarks:
- You should also swap the rows such that the row that you are going to use the first element of is at the top most position of the rows you are working on. So, if you are working on the second column, the first row should remain untouched and you should swap the row with the one at the second row (in your case the largest element is already at the second row, but that is just luck). Sorry if this is unclear, look at the Wikipedia article.
- The largest element should be the larges elements of row 2-9, so excluding row 1!
- You should not anything to row 1 if you are working on the second column. Also you should not add anything to row 1 or 2 if you are working on the third column etc.


Usually doing this by itself does not make sense. Usually you are solving something like Mx = b where your matrix is M, or you are inverting a matrix (or a bunch of other stuff). In both cases you should remember all operations you do because that will give you either x or M^-1, depending on what you are looking for.

 

[pasmith]

TL;DR Summary: I need to compute the vibrational frequencies of a molecule when the matrix of force constants (second derivative of the energy by the Cartesian coordinates) is provided. For such computation, this matrix must be diagonalized.

Are these frequencies related to the eigenvalues of the matrix? In this case, you don't need to "diagonalize" the matrix. Whatever programming language you are using probably has a linear algebra library which will include a function to calculate the eigenvalues (and eigenvectors, if needed) of a given matrix.

 

[Arjan82]

That is actually a good point by @pasmith:

Unless you are trying to learn about the algorithm and/or computer science (programming) itself, you should indeed find some library or tool alike to do this for you. Think about Python + numpy, Matlab/Octave, LAPACK, etc. etc. There is a whole list here:

https://en.wikipedia.org/wiki/Comparison_of_linear_algebra_libraries

 

[Spathi]

- You should also swap the rows such that the row that you are going to use the first element of is at the top most position of the rows you are working on. So, if you are working on the second column, the first row should remain untouched and you should swap the row with the one at the second row (in your case the largest element is already at the second row, but that is just luck). Sorry if this is unclear, look at the Wikipedia article.

Sorry, your English is not perfect, and I think I found something strange in the Wikipedia article cited. Some more questions:
1) You want to say that I should not add row 2 to row1 for eliminating the item in Column2 Row1 (after the first iteration which eliminates the items in column 1, rows 2-9)?
2) You want to say that I need (in some cases) to switch the first row with some other row at this iteration?

 

[Spathi]

Are these frequencies related to the eigenvalues of the matrix? In this case, you don't need to "diagonalize" the matrix. Whatever programming language you are using probably has a linear algebra library which will include a function to calculate the eigenvalues (and eigenvectors, if needed) of a given matrix.

Possibly I will need to create one more thread in other section; my task it to solve the nuclear Schrödinger equation for the "rigid rotor - harmonic oscillator" approximation. Personally I find strange that in the books it is written that the matrix should be diagonalized, and it is unclear for me, where are the values at the right column related to linear equations (sorry I don't remember how they are called in English).