scipy.linalg.eig_banded appears to give unreliable results inside solve_schrodinger_equation

[amarkpayne]

Bug Description

This was discovered in the course of debugging #1682. I am opening this up as a separate issue to refocus our efforts, and because it is possible that this is not the only issue affecting the failing pdep test.

Inside HinderedRotor.solve_schrodinger_equation (rmgpy.statmech.torsion.pyx) we construct a complex-valued Hamiltonian matrix banded in lower triangular form (i.e. the full Hamiltonian is a 401x401 matrix, but here we only store the diagonals that are non-zero that appear on and below the main diagonal. Therefore, H is only 6x401, but represents a full 401x401 Hamiltonian). This banded representation then gets passed to scipy.linalg.eig_banded, which takes as input a complex-valued Hamiltonian represented in this way and returns the eigenvalues.

Here is where the problem occurs. When running two different versions of RMG-Py (see below for more details), the exact same Hamiltonian input into scipy.linalg.eig_banded yields completely different eigenvalues, where some of them don't even match to 1 significant figure. As far as we can tell the Hamiltonian in question (which I have attached by outputting as a list to a YAML file) is not poorly conditioned. It is worth noting that this can be hard to reproduce, because it appears to be system dependent. I am seeing this issue based on current master, but others aren't necessarily.

What is reproducible though is that the eigenvalues are different if you request scipy.linalg.eig_banded to return the eigenvectors as well, which as noted in the documentation for the underlying FORTRAN LAPACK function chbevd.f at http://www.netlib.org/lapack/explore-html/d5/db8/chbevd_8f_source.html will cause a different algorithm to be used. It is not surprising that the eigenvalues will be different if a different algorithm is used, but it is surprising if they are so different that they don't even agree to 1 significant digit. The eigenvalues were different when the eigenvectors were requested in both RMG-Py versions, and worse they were still different across the different RMG-Py versions.

This is potentially a major bug, because we use these eigenvalues (energies) to calculate the partition function of the hindered rotor, which of course is the most important quantity in all of statmech which we derive everything else from. The resulting partition functions can differ in as much as the second significant digit, so it is important that we figure out what is going on here.

How To Reproduce

Finding different versions of RMG-Py where the output of scipy.linalg.eig_banded differs can be difficult to find. It appears to be system dependent, and we really don't have a good understanding of why at this point.

Instead, I recommend changing the eigvals_only keyword argument from False to True and comparing the output.

Expected Behavior

The returned eigenvalues should have more precision than they appear to have

Installation Information

• OS (include version if known): Ubuntu 18.04

• Installation method: source

• RMG-Py version "py3 (nearly current master)" information:

  • RMG-Py: b691b45
  • RMG-database: e409c4

• RMG-Py version "py2 (v2.5.0)" information:

  • RMG-Py: 068ec11
  • RMG-database: e409c4

Additional Context

On an interesting note, I tried making the following modification to HinderedRotor.solve_schrodinger_equation to convert the banded representation to a dense numpy array, and then use np.linalg.eig to find the eigenvalues, which I think calls a different LAPACK function. Using the code below, I get different eigenvalues again, but both versions of RMG-Py now agree with each other. The eigenvalues found appear to have repeat entries (they come in pairs of two I think for the most part), which is interesting.

I am not sure at all if this is correct (I could have messed up converting to a dense matrix, and it is also possible that the eigenvalues found from the numpy call are also unreliable), but if this is actually the correct result then both RMG-Py codes yield the same rate coefficient in the pdep test, it is just that they are now a factor of 10 smaller than what the test number is.

# OUTDATED. Please see the comment below

# Populate Hamiltonian matrix (banded in lower triangular form)
        H = self.get_hamiltonian(n_basis)
        H = diags(H, np.arange(H.shape[0])*-1)
        H = H.todense()
        # The overlap matrix is the identity matrix, i.e. this is a standard
        # eigenvalue problem

        # Find the eigenvalues and eigenvectors of the Hamiltonian matrix
        #E, _ = scipy.linalg.eig_banded(H, lower=True, eigvals_only=False, overwrite_a_band=False)
        E = np.real(np.sort(np.linalg.eig(H)[0]))
        # Don't consider zero-point energy here
        self.energies = E - np.min(E)
        print('Here are the energies:\n{0}'.format(self.energies))

        # Return the eigenvalues
        return self.energies

I would be very interested to see what other codes (what does Julia say?) get for the eigenvalues for this problem. What even is the correct answer?

Final note: I had to change the file extension for the hamiltonian, but it is a YAML file
hamiltonian.txt

[amarkpayne]

Update: We realized that my code was incorrect for reconstructing the dense Hamiltonian. We forgot to include the complex conjugate on the upper diags. We could have also gotten arround that by using np.linalg.eigh though. Here is the updated code:

# Populate Hamiltonian matrix (banded in lower triangular form)
        H = self.get_hamiltonian(n_basis)
        m = H.shape[0]
        indices = [0]+list(range(1, m))+[-i for i in range(1,m)]
        H = np.concatenate((H, np.conj(H[1:, :])))
        H = diags(H, np.array(indices))
        H = H.todense()
        # The overlap matrix is the identity matrix, i.e. this is a standard
        # eigenvalue problem

        # Find the eigenvalues and eigenvectors of the Hamiltonian matrix
        #E, _ = scipy.linalg.eig_banded(H, lower=True, eigvals_only=False, overwrite_a_band=False)
        E = np.real(np.sort(np.linalg.eig(H)[0]))
        # Don't consider zero-point energy here
        self.energies = E - np.min(E)

        # Return the eigenvalues
        return self.energies

[github-actions]

This issue is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant issue, otherwise it will automatically be closed in 30 days.