Diverging HEOM Dynamics with Fitted Correlation Function in QuTiP

[AMIT UPADHYAY]

Dear QuTiP community,

I am currently simulating the HEOM dynamics of a three-site Frenkel-exciton model, where each site interacts with an independent bosonic bath characterized by an Ohmic spectral density. To model the environment, I used the QuTiP function env.approximate("prony", tlist_corr, separate=True, Nr=8, Ni=8) to fit the exact bath correlation function for use with HEOMSolver.

However, I’m encountering divergence in the population dynamics at steady state. I have also tested other fitting methods such as "cf" and "espira-I", and increased both the number of exponential terms and the hierarchy depth (up to 4). Unfortunately, going beyond this depth is computationally infeasible due to time constraints.

I would greatly appreciate it if someone could help verify whether there is an error in my simulation setup. For reference, I’ve attached the Python script I am using.

Thank you for your time and support.

Best regards,
Amit Kumar Upadhyay

---

Let me know if you'd like a shorter or more technical version.

Ask ChatGPT

[Neill Lambert]

hi amit,

For the divergence of the heom solver, yeah this can occur when there is lack of convergence in the cut-off.  

A couple of tips that might help: 

- probably its better to avoid using ''real units'', so setting constant=1, and appropriately renormalizing the temperature, and using hbar=kb=1 are good steps (I think the cf fit method in particular might have trouble with real units).

 - second, when you define a bath from a function, you set this wMax = wc parameter. probably this needs to be a bit larger, say 10*wc  (it is, as far as i recall, the range of some internal fft's used to get the correlation function). Alternatively, you could use our predefined ohmic bath (https://qutip.readthedocs.io/en/stable/guide/environments/bosonic.html#pre-defined-environments), which has analytical expressions for the correlation functions, so something like:

env = qutip.OhmicEnvironment(T=Tempheom, alpha=lam/wc, wc=wc, s=1)

Working out which fitting method works best for a particular problem can be a bit of an art, I suggest trying 'aaa',  and 'cf', and trying to see which gives the best fit for minimal number of exponents.  alternatively, start with a poor fit with just a few exponents, and then see how the system dynamics converge as both the fit is gradually improved and the convergence parameter is increased

all the best

neill

--
You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/qutip/a1e4bc9b-c69b-45e8-abc6-c9f59baef0c6n%40googlegroups.com.

--

-----------------------------------
Senior Research Scientist
RIKEN, Japan
Research Homepage

[Neill Lambert]

Hi amit,

I had a quick play with your notebook, I made those changes I suggested, and seem to get a more reasonable result (with some caveats, not sure I rescaled everything correctly).

Here 'cf', which is just fitting the correlation functions with exponentials with scipys curve_fit, seems to do a better job than prony, not sure why.  anyways, maybe a useful basis for you to experiment further.

all the best

neill