Symmetrize the dynmat before diagonalize it

[litman90]

The 'np.linalg.eigh' algorithm is designed for symmetric (hermitian) matrix. Due to numerical inaccuracy the dynmat is not symmetric.

Example for a saddle point in a water dimer, frequencies in cm^-1: Previous: -175.2 -20.3 -11.9 -6.0 -5.5 9.8 21.2 120.3 159.9 201.4
With the correction
-176.1 - 2.5E-05 -1.4E-05 1.0E-05 1.7E-05 4.0E-05 5.2E-05 122.8 159.8 201.8

(Note that results with the correction are the same if we use np.linalg.eig, which is for a general matrix. Although, in this case, a real output is not guaranteed and tiny imaginary part could appear. Of course np.linalg.eigh is faster )

[mahrossi]

Just to give the heads up: I had discussed this with Yair and indeed we were getting "much crappier"(they were not so bad) frequencies with i-pi than with basically any other code. So this symmetrization (before using the algorithm for symmetric matrices) helps a lot. Could also be relevant in other contexts where instabilities were previously seen.

[venkatkapil24]

Interesting, since the differences are not small at all. So far I had benchmarked phonons calculations on forcefields for which there was not much numerical noise. So I guess I did not use it everywhere.

However, instead of doing it in the print all function I think is better to do it in each of the phonons methods (fd,nmfd and ennmfd) so that even the restart files also contain the noise-free symmetric dynamtrix.

[mahrossi]

@venkatkapil24 I agree

[litman90]

@venkatkapil24 What do you mean with "to do it in each of the phonons methods". I can imagine two things. 1_ To symmetrize the dynmatrix and the refdynmatrix just before calling "self.applyasr(self.refdynmatrix)" 2 Instead of updating the hessian row by row, we can also update at the same time the corresponding column.

[venkatkapil24]

@litman90 : I meant inside the DummyPhononCalculator. It contains a function called transform which is executed once the dynmatrix (or the refdynmatrix) is evaluated.

I have added a line that symmetrises the dynmatrix (or the refdynmatrix) . I also added a line just after apply_asr().

[litman90]

I tested again. Everything is ok.