The paper https://arxiv.org/pdf/2004.01469.pdf seems to present an interesting way to setup a general Lindblad equation.
There may still be surprises, e.g., that you do not end up in thermal equilibrium, but at least you can reduce the problem of setting up a Lindblad operator to the problem of setting up the spectral density.
Implemented; only implemented a Lorentzian as functor; converted MolVibration/OH/3 to use this scheme with spectral density 1/E.
I am not quite sure yet about the usage of functors for spectral densities, which is why only a single Lorentzian was implemented. Needs some more playtesting to get a good feeling for what you might want to do.
Along the way, I am also not quite sure what to do with the terms that are not easy to calculate (time-dependent Lindblad, Lamb shift); altogether, using a Lindblad Liouvillian feels again awkward like what you do not want to use. Maybe I should try out the exact non-Markovian approach of Tannor etc.
Diff:
Diff: