Thank you for that explanation ... that is the way I understood the math to work in the case where the density matrix is diagonal. The issue with the example that I gave is that the off-diagonal elements of the matrix are time dependent ... i.e. the magnitude of the anharmonic cross-couplings between the vibrational modes are dependent on the population of those modes (by population I mean the number of vibrational quanta in a given mode). I believe the ramifications of this are that you cannot choose a unique basis that diagonalizes the density matrix at all times, but like I said, I am not that well-versed in the details of density matrices, so I may not have that completely correct.
However, I have an even more fundamental problem because I do not understand from a physical point of view how the entropy of the isolated vibrational system in the example I gave can possibly be time-invariant. It seems intuitively clear to me from the phenomenology of intramolecular vibrational redistribution that the process is entropically driven. In other words, the energy starts out localized as a single quantum of excitation in a single vibrational mode, and then becomes "randomized" as one or more quanta of excitation in multiple vibrational modes with lower energy. The total internal energy of the system remains constant, but the probability of the energy finding its way back into the mode that was initially excited is (I think) vanishingly small. That seems consistent with the evolution of the state from low entropy (all the energy in a single mode) to higher entropy (the energy redistributed among many modes).
A possible counter-argument to the description I gave above might be that, even though the energy is "randomized", at any instant in time it is described by a unique "pure state" of the system, which would have the same von Neumann entropy (zero) as the initial state with a single quantum of excitation in a single mode. This argument is probably valid for small molecules where the density of states is low, however molecules belonging to symmetric point groups with degenerate irreducible representations have formal degeneracies that would give rise to non-zero entropies for particular combinations of vibrational quanta.
I would appreciate any insights you have on this ...