diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index fbc7cfc641..ffee236f08 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -141,7 +141,7 @@ jobs: # Install the extra requirement python -m pip install pytest>=5.2 pytest-rerunfailures # tests python -m pip install ipython # ipython - python -m pip install loky tqdm # extras + python -m pip install loky tqdm mpmath # extras python -m pip install "coverage${{ matrix.coverage-requirement }}" chardet python -m pip install pytest-cov coveralls pytest-fail-slow diff --git a/doc/apidoc/apidoc.rst b/doc/apidoc/apidoc.rst index ff2d1f3390..b5d20b8902 100644 --- a/doc/apidoc/apidoc.rst +++ b/doc/apidoc/apidoc.rst @@ -13,6 +13,7 @@ complete list of QuTiP's public classes and functions. quantumobject.rst time_dep.rst solver.rst + environments.rst heom.rst piqs.rst visualization.rst diff --git a/doc/apidoc/environments.rst b/doc/apidoc/environments.rst new file mode 100644 index 0000000000..00edf32d0b --- /dev/null +++ b/doc/apidoc/environments.rst @@ -0,0 +1,52 @@ +************ +Environments +************ + +Bosonic Environments +-------------------- + +.. autoclass:: qutip.core.BosonicEnvironment + :members: + +.. autoclass:: qutip.core.DrudeLorentzEnvironment + :members: + :inherited-members: + :show-inheritance: + :exclude-members: from_correlation_function, from_power_spectrum, from_spectral_density + +.. autoclass:: qutip.core.UnderDampedEnvironment + :members: + :inherited-members: + :show-inheritance: + :exclude-members: from_correlation_function, from_power_spectrum, from_spectral_density + +.. autoclass:: qutip.core.OhmicEnvironment + :members: + :inherited-members: + :show-inheritance: + :exclude-members: from_correlation_function, from_power_spectrum, from_spectral_density + +.. autoclass:: qutip.core.CFExponent + :members: + +.. autoclass:: qutip.core.ExponentialBosonicEnvironment + :members: + :show-inheritance: + +.. autofunction:: qutip.core.environment.system_terminator + + +Fermionic Environments +---------------------- + +.. autoclass:: qutip.core.FermionicEnvironment + :members: + :exclude-members: from_correlation_function, from_power_spectrum, from_spectral_density + +.. autoclass:: qutip.core.LorentzianEnvironment + :members: + :show-inheritance: + +.. autoclass:: qutip.core.ExponentialFermionicEnvironment + :members: + :show-inheritance: diff --git a/doc/apidoc/utilities.rst b/doc/apidoc/utilities.rst index 72c69b694e..98def612c4 100644 --- a/doc/apidoc/utilities.rst +++ b/doc/apidoc/utilities.rst @@ -7,7 +7,7 @@ Utility Functions ----------------- .. automodule:: qutip.utilities - :members: n_thermal, clebsch, convert_unit + :members: n_thermal, clebsch, convert_unit, iterated_fit .. _fileio: diff --git a/doc/biblio.rst b/doc/biblio.rst index 9ae69b60df..24ab4dcf77 100644 --- a/doc/biblio.rst +++ b/doc/biblio.rst @@ -3,78 +3,86 @@ Bibliography ============ +.. Note: first letter of entries must be escaped to avoid rst parsing as enumerated list + https://docutils.sourceforge.io/docs/ref/rst/restructuredtext.html#enumerated-lists + .. [BCSZ08] - W. Bruzda, V. Cappellini, H.-J. Sommers, K. Życzkowski, *Random Quantum Operations*, Phys. Lett. A **373**, 320-324 (2009). :doi:`10.1016/j.physleta.2008.11.043`. + \W. Bruzda, V. Cappellini, H.-J. Sommers, K. Życzkowski, *Random Quantum Operations*, Phys. Lett. A **373**, 320-324 (2009). :doi:`10.1016/j.physleta.2008.11.043`. .. [Hav03] - Havel, T. *Robust procedures for converting among Lindblad, Kraus and matrix representations of quantum dynamical semigroups*. Journal of Mathematical Physics **44** 2, 534 (2003). :doi:`10.1063/1.1518555`. + \T. Havel, *Robust procedures for converting among Lindblad, Kraus and matrix representations of quantum dynamical semigroups*, J. Math. Phys. **44** 2, 534 (2003). :doi:`10.1063/1.1518555`. .. [Wat13] - Watrous, J. |theory-qi|_, lecture notes. - -.. The trick with |text|_ is to get an italic link, and is described in the - Docutils FAQ at https://docutils.sourceforge.net/FAQ.html#is-nested-inline-markup-possible. - -.. |theory-qi| replace:: *Theory of Quantum Information* -.. _theory-qi: https://cs.uwaterloo.ca/~watrous/TQI-notes/ + \J. Watrous, |theory-qi|_, lecture notes. .. [Mez07] - F. Mezzadri, *How to generate random matrices from the classical compact groups*, Notices of the AMS **54** 592-604 (2007). :arxiv:`math-ph/0609050`. + \F. Mezzadri, *How to generate random matrices from the classical compact groups*, Notices of the AMS **54**, 592-604 (2007). :arxiv:`math-ph/0609050`. .. [Moh08] - M. Mohseni, A. T. Rezakhani, D. A. Lidar, *Quantum-process tomography: Resource analysis of different strategies*, Phys. Rev. A **77**, 032322 (2008). :doi:`10.1103/PhysRevA.77.032322`. + \M. Mohseni, A. T. Rezakhani, D. A. Lidar, *Quantum-process tomography: Resource analysis of different strategies*, Phys. Rev. A **77**, 032322 (2008). :doi:`10.1103/PhysRevA.77.032322`. .. [Gri98] - M. Grifoni, P. Hänggi, *Driven quantum tunneling*, Physics Reports **304**, 299 (1998). :doi:`10.1016/S0370-1573(98)00022-2`. + \M. Grifoni, P. Hänggi, *Driven quantum tunneling*, Phys. Rep. **304**, 299 (1998). :doi:`10.1016/S0370-1573(98)00022-2`. .. [Gar03] - Gardineer and Zoller, *Quantum Noise* (Springer, 2004). + Gardiner and Zoller, *Quantum Noise* (Springer, 2004). .. [Bre02] H.-P. Breuer and F. Petruccione, *The Theory of Open Quantum Systems* (Oxford, 2002). .. [Coh92] - C. Cohen-Tannoudji, J. Dupont-Roc, G. Grynberg, *Atom-Photon Interactions: Basic Processes and Applications*, (Wiley, 1992). + \C. Cohen-Tannoudji, J. Dupont-Roc, G. Grynberg, *Atom-Photon Interactions: Basic Processes and Applications*, (Wiley, 1992). .. [WBC11] - C. Wood, J. Biamonte, D. G. Cory, *Tensor networks and graphical calculus for - open quantum systems*. :arxiv:`1111.6950` + \C. Wood, J. Biamonte, D. G. Cory, *Tensor networks and graphical calculus for open quantum systems*. :arxiv:`1111.6950` .. [AKN98] - D. Aharonov, A. Kitaev, and N. Nisan, *Quantum circuits with mixed states*, in Proceedings of the - thirtieth annual ACM symposium on Theory of computing, 20-30 (1998). :arxiv:`quant-ph/9806029` + \D. Aharonov, A. Kitaev, N. Nisan, *Quantum circuits with mixed states*, in Proceedings of the Thirtieth Annual ACM STOC, 20-30 (1998). :arxiv:`quant-ph/9806029` .. [dAless08] - D. d’Alessandro, *Introduction to Quantum Control and Dynamics*, (Chapman & Hall/CRC, 2008). + \D. d’Alessandro, *Introduction to Quantum Control and Dynamics*, (Chapman & Hall/CRC, 2008). .. [Byrd95] - R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, *A Limited Memory Algorithm for Bound Constrained Optimization*, SIAM J. Sci. Comput. **16**, 1190 (1995). :doi:`10.1137/0916069` + \R. H. Byrd, P. Lu, J. Nocedal, C. Zhu, *A Limited Memory Algorithm for Bound Constrained Optimization*, SIAM J. Sci. Comput. **16**, 1190 (1995). :doi:`10.1137/0916069` .. [Flo12] - F. F. Floether, P. de Fouquieres, and S. G. Schirmer, *Robust quantum gates for open systems via optimal control: Markovian versus non-Markovian dynamics*, New J. Phys. **14**, 073023 (2012). :doi:`10.1088/1367-2630/14/7/073023` + \F. F. Floether, P. de Fouquieres, S. G. Schirmer, *Robust quantum gates for open systems via optimal control: Markovian versus non-Markovian dynamics*, New J. Phys. **14**, 073023 (2012). :doi:`10.1088/1367-2630/14/7/073023` .. [Lloyd14] - S. Lloyd and S. Montangero, *Information theoretical analysis of quantum optimal control*, Phys. Rev. Lett. **113**, 010502 (2014). :doi:`10.1103/PhysRevLett.113.010502` + \S. Lloyd, S. Montangero, *Information theoretical analysis of quantum optimal control*, Phys. Rev. Lett. **113**, 010502 (2014). :doi:`10.1103/PhysRevLett.113.010502` .. [Doria11] - P. Doria, T. Calarco & S. Montangero, *Optimal Control Technique for Many-Body Quantum Dynamics*, Phys. Rev. Lett. **106**, 190501 (2011). :doi:`10.1103/PhysRevLett.106.190501` + \P. Doria, T. Calarco, S. Montangero, *Optimal Control Technique for Many-Body Quantum Dynamics*, Phys. Rev. Lett. **106**, 190501 (2011). :doi:`10.1103/PhysRevLett.106.190501` .. [Caneva11] - T. Caneva, T. Calarco, & S. Montangero, *Chopped random-basis quantum optimization*, Phys. Rev. A **84**, 022326 (2011). :doi:`10.1103/PhysRevA.84.022326` + \T. Caneva, T. Calarco, S. Montangero, *Chopped random-basis quantum optimization*, Phys. Rev. A **84**, 022326 (2011). :doi:`10.1103/PhysRevA.84.022326` .. [Rach15] - N. Rach, M. M. Müller, T. Calarco, and S. Montangero, *Dressing the chopped-random-basis optimization: A bandwidth-limited access to the trap-free landscape*, Phys. Rev. A. **92**, 062343 (2015). :doi:`10.1103/PhysRevA.92.062343` + \N. Rach, M. M. Müller, T. Calarco, S. Montangero, *Dressing the chopped-random-basis optimization: A bandwidth-limited access to the trap-free landscape*, Phys. Rev. A. **92**, 062343 (2015). :doi:`10.1103/PhysRevA.92.062343` .. [Wis09] - - Wiseman, H. M. & Milburn, G. J. *Quantum Measurement and Control*, (Cambridge University Press, 2009). + \H. M. Wiseman, G. J. Milburn, *Quantum Measurement and Control*, (Cambridge University Press, 2009). .. [NKanej] - - N Khaneja et. al. *Optimal control of coupled spin dynamics: Design of NMR pulse sequences by gradient ascent algorithms.* J. Magn. Reson. **172**, 296–305 (2005). :doi:`10.1016/j.jmr.2004.11.004` + \N. Khaneja *et al.*, *Optimal control of coupled spin dynamics: Design of NMR pulse sequences by gradient ascent algorithms.* J. Magn. Reson. **172**, 296–305 (2005). :doi:`10.1016/j.jmr.2004.11.004` .. [Donvil22] - B. Donvil, P. Muratore-Ginanneschi, *Quantum trajectory framework for general time-local master equations*, Nat Commun **13**, 4140 (2022). :doi:`10.1038/s41467-022-31533-8`. + \B. Donvil, P. Muratore-Ginanneschi, *Quantum trajectory framework for general time-local master equations*, Nat Commun **13**, 4140 (2022). :doi:`10.1038/s41467-022-31533-8`. .. [Abd19] - M. Abdelhafez, D. I. Schuster, J. Koch, *Gradient-based optimal control of open quantum systems using quantumtrajectories and automatic differentiation*, Phys. Rev. A **99**, 052327 (2019). :doi:`10.1103/PhysRevA.99.052327`. + \M. Abdelhafez, D. I. Schuster, J. Koch, *Gradient-based optimal control of open quantum systems using quantumtrajectories and automatic differentiation*, Phys. Rev. A **99**, 052327 (2019). :doi:`10.1103/PhysRevA.99.052327`. + +.. [BoFiN23] + \N. Lambert, T. Raheja, S. Cross, P. Menczel, S. Ahmed, A. Pitchford, D. Burgarth, F. Nori, *QuTiP-BoFiN: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics*, Phys. Rev. Research **5**, 013181 (2023). :doi:`10.1103/PhysRevResearch.5.013181`. + +.. [Lambert19] + \N. Lambert, S. Ahmed, M. Cirio, F. Nori, *Modelling the ultra-strongly coupled spin-boson model with unphysical modes*, Nat Commun **10**, 3721 (2019). :doi:`10.1038/s41467-019-11656-1`. + + + +.. The trick with |text|_ is to get an italic link, and is described in the + Docutils FAQ at https://docutils.sourceforge.net/FAQ.html#is-nested-inline-markup-possible. + This is at the bottom of the source file to avoid extra whitespace. + +.. |theory-qi| replace:: *Theory of Quantum Information* +.. _theory-qi: https://cs.uwaterloo.ca/~watrous/TQI-notes/ \ No newline at end of file diff --git a/doc/changes/2534.feature b/doc/changes/2534.feature new file mode 100644 index 0000000000..acad251d7c --- /dev/null +++ b/doc/changes/2534.feature @@ -0,0 +1 @@ +This change introduces the qutip.core.environment module, which contains classes that characterize bosonic and fermionic thermal environments. \ No newline at end of file diff --git a/doc/guide/guide-environments.rst b/doc/guide/guide-environments.rst new file mode 100644 index 0000000000..292baf682c --- /dev/null +++ b/doc/guide/guide-environments.rst @@ -0,0 +1,491 @@ +.. _environments guide: + +************************************ +Environments of Open Quantum Systems +************************************ + +*written by* |pm|_ *and* |gs|_ + +.. _pm: https://www.menczel.net/ +.. |pm| replace:: *Paul Menczel* +.. _gs: https://gsuarezr.github.io/ +.. |gs| replace:: *Gerardo Suarez* +.. (this is a workaround for italic links in rst) + +QuTiP can describe environments of open quantum systems. +They can be passed to various solvers, where their influence is taken into account exactly or approximately. +In the following, we will discuss bosonic and fermionic thermal environments. +In our definitions, we follow [BoFiN23]_. + +Note that currently, we only support a single coupling term per environment. +If a more generalized coupling would be useful to you, please let us know on GitHub. + + +.. _bosonic environments guide: + +Bosonic Environments +-------------------- + +A bosonic environment is described by a continuum of harmonic oscillators. +The open quantum system and its environment are thus described by the Hamiltonian + +.. math:: + + H = H_{\rm s} + Q \otimes X + \sum_k \hbar\omega_k\, a_k^\dagger a_k , \qquad X = \sum_k g_k (a_k + a_k^\dagger) + +(in the case of a single bosonic environment). +Here, :math:`H_{\rm s}` is the Hamiltonian of the open system and :math:`Q` a system coupling operator. +The sums enumerate the environment modes with frequencies :math:`\omega_k` and couplings :math:`g_k`. +The couplings are described by the spectral density + +.. math:: + + J(\omega) = \pi \sum_k g_k^2 \delta(\omega - \omega_k) . + +Equivalently, a bosonic environment can be described by its auto-correlation function + +.. math:: + + C(t) = \langle X(t) X \rangle = \sum_{k} g_{k}^{2} \Big( \cos(\omega_{k} t) + \underbrace{( 2 n_{k}+1)}_{\coth(\frac{\beta \omega_{k}}{2})} + - i \sin(\omega_{k} t) \Big) + +(where :math:`X(t)` is the interaction picture operator) or its power spectrum + +.. math:: + + S(\omega) = \int_{-\infty}^\infty \mathrm dt\, C(t)\, \mathrm e^{\mathrm i\omega t} , + +which is the inverse Fourier transform of the correlation function. +The correlation function satisfies the symmetry relation :math:`C(-t) = C(t)^\ast`. + +Assuming that the initial state is thermal with inverse temperature +:math:`\beta`, in the continuum limit the correlation function and the power +spectrum can be calculated from the spectral density via + +.. math:: + :label: cfandps + + \begin{aligned} + C(t) &= \int_0^\infty \frac{\mathrm d\omega}{\pi}\, J(\omega) \Bigl[ \coth\Bigl( \frac{\beta\omega}{2} \Bigr) \cos\bigl( \omega t \bigr) - \mathrm i \sin\bigl( \omega t \bigr) \Bigr] , \\ + S(\omega) &= \operatorname{sign}(\omega)\, J(|\omega|) \Bigl[ \coth\Bigl( \frac{\beta\omega}{2} \Bigr) + 1 \Bigr] . + \end{aligned} + +Here, :math:`\operatorname{sign}(\omega) = \pm 1` depending on the sign of :math:`\omega`. +At zero temperature, these equations become :math:`C(t) = \int_0^\infty \frac{\mathrm d\omega}{\pi} J(\omega) \mathrm e^{-\mathrm i\omega t}` and :math:`S(\omega) = 2 \Theta(\omega) J(|\omega|)`, where :math:`\Theta` is the Heaviside function. + +If the environment is coupled weakly to the open system, the environment induces quantum jumps with transition rates :math:`\gamma \propto S(-\Delta\omega)`, where :math:`\Delta\omega` is the energy change in the system corresponding to the quantum jump. +In the strong coupling case, QuTiP provides exact integration methods based on multi-exponential decompositions of the correlation function, see below. + +.. note:: + We generally assume that the frequencies :math:`\omega_k` are all positive and hence :math:`J(\omega) = 0` for :math:`\omega \leq 0`. + To handle a spectral density :math:`J(\omega)` with support at negative frequencies, one can use an effective spectral density :math:`J'(\omega) = J(\omega) - J(-\omega)`. + This process might result in the desired correlation function, because + + .. math:: + + \int_0^\infty \frac{\mathrm d\omega}{\pi}\, J'(\omega) \bigl[ \cdots \bigr] = \int_{-\infty}^\infty \frac{\mathrm d\omega}{\pi}\, J(\omega) \bigl[ \cdots \bigr] , + + where :math:`[\cdots]` stands for the square brackets in Eq. :eq:`cfandps`. + Note however that the derivation of Eq. :eq:`cfandps` is not valid in this situation, since the environment does not have a thermal state. + +.. note:: + Note that the expressions provided above for :math:`S(\omega)` are ill-defined at :math:`\omega=0`. + For zero temperature, we simply have :math:`S(0) = 0`. + + For non-zero temperature, one obtains + + .. math:: + + S(0) = 2\beta^{-1}\, J'(0) . + + Hence, :math:`S(0)` diverges if the spectral density is sub-ohmic. + + +Pre-defined Environments +------------------------ + +Ohmic Environment +^^^^^^^^^^^^^^^^^ + +Ohmic environments can be constructed in QuTiP using the class :class:`.OhmicEnvironment`. +They are characterized by spectral densities of the form + +.. math:: + :label: ohmicf + + J(\omega) = \alpha \frac{\omega^s}{\omega_c^{s-1}} e^{-\omega / \omega_c} , + +where :math:`\alpha` is a dimensionless parameter that indicates the coupling strength, +:math:`\omega_{c}` is the cutoff frequency, and :math:`s` is a parameter that determines the low-frequency behaviour. +Ohmic environments are usually classified according to this parameter as + +* Sub-Ohmic (:math:`s<1`) +* Ohmic (:math:`s=1`) +* Super-Ohmic (:math:`s>1`). + +.. note:: + In the literature, the Ohmic spectral density can often be found as :math:`J(\omega) = \alpha \frac{\omega^s}{\omega_c^{s-1}} f(\omega)`, + where :math:`f(\omega)` with :math:`\lim\limits_{\omega \to \infty} f(\omega) = 0` is known as the cutoff function. + The cutoff function ensures that the spectral density and its integrals (for example :eq:`cfandps`) do not diverge. + Sometimes, with sub-Ohmic spectral densities, an infrared cutoff is used as well so that :math:`\lim\limits_{\omega \to 0} J(\omega) = 0`. + This pre-defined Ohmic environment class is restricted to an exponential cutoff function, which is one of the most commonly used in the literature. + Other cutoff functions can be used in QuTiP with user-defined environments as explained below. + +Substituting the Ohmic spectral density :eq:`ohmicf` into :eq:`cfandps`, the correlation function can be computed analytically: + +.. math:: + C(t)= \frac{\alpha}{\pi} w_{c}^{1-s} \beta^{-(s+1)} \Gamma(s+1) + \left[ \zeta\left(s+1,\frac{1+\beta w_{c} -i w_{c} t}{\beta w_{c}} + \right) +\zeta\left(s+1,\frac{1+ i w_{c} t}{\beta w_{c}}\right) + \right] , + +where :math:`\beta` is the inverse temperature, :math:`\Gamma` the Gamma function, and :math:`\zeta` the Hurwitz zeta function. +The zero temperature case can be obtained by taking the limit :math:`\beta \to \infty`, which results in + +.. math:: + C(t) = \frac{\alpha}{\pi} \omega_c^2\, \Gamma(s+1) (1+ i \omega_{c} t)^{-(s+1)} . + +The evaluation of the zeta function for complex arguments requires `mpmath`, so certain features of the Ohmic enviroment are +only available if `mpmath` is installed. + +Multi-exponential approximations to Ohmic environments can be obtained through +the fitting procedures :meth:`approx_by_cf_fit<.BosonicEnvironment.approx_by_cf_fit>` +and :meth:`approx_by_sd_fit<.BosonicEnvironment.approx_by_sd_fit>`. +The following example shows how to create a sub-Ohmic environment, and how to use +:meth:`approx_by_cf_fit<.BosonicEnvironment.approx_by_cf_fit>` to fit the real and imaginary parts +of the correlation function with two exponential terms each. + +.. plot:: + :context: reset + :nofigs: + + import numpy as np + import qutip as qt + import matplotlib.pyplot as plt + + # Define a sub-Ohmic environment with the given temperature, coupling strength and cutoff + env = qt.OhmicEnvironment(T=0.1, alpha=1, wc=3, s=0.7) + + # Fit the correlation function with three exponential terms + tlist = np.linspace(0, 3, 250) + approx_env, info = env.approx_by_cf_fit(tlist, target_rsme=None, Nr_max=3, Ni_max=3, maxfev=1e8) + +The environment `approx_env` created here could be used, for example, with the :ref:`HEOM solver`. +The variable `info` contains info about the convergence of the fit; here, we will just plot the fit together with +the analytical correlation function. Note that a larger number of exponential terms would have yielded a better result. + +.. plot:: + :context: + + plt.plot(tlist, np.real(env.correlation_function(tlist)), label='Real part (analytic)') + plt.plot(tlist, np.real(approx_env.correlation_function(tlist)), '--', label='Real part (fit)') + + plt.plot(tlist, np.imag(env.correlation_function(tlist)), label='Imag part (analytic)') + plt.plot(tlist, np.imag(approx_env.correlation_function(tlist)), '--', label='Imag part (fit)') + + plt.xlabel('Time') + plt.ylabel('Correlation function') + plt.tight_layout() + plt.legend() + + +.. _dl env guide: + +Drude-Lorentz Environment +^^^^^^^^^^^^^^^^^^^^^^^^^ + +Drude-Lorentz environments, also known as overdamped environments, can be constructed in QuTiP +using the class :class:`.DrudeLorentzEnvironment`. They are characterized by spectral densities of the form + +.. math:: + J(\omega) = \frac{2 \lambda \gamma \omega}{\gamma^{2}+\omega^{2}} , + +where :math:`\lambda` is a coupling strength (with the dimension of energy) and :math:`\gamma` the cutoff frequency. + +To compute the corresponding correlation function, one can apply the Matsubara expansion: + +.. math:: + C(t) = \sum_{k=0}^{\infty} c_k e^{- \nu_k t} + +The coefficients of this expansion are + +.. math:: + + \nu_{k} = \begin{cases} + \gamma & k = 0\\ + {2 \pi k} / {\beta} & k \geq 1\\ + \end{cases} \;, \qquad + c_k = \begin{cases} + \lambda \gamma [\cot(\beta \gamma / 2) - i] & k = 0\\ + \frac{4 \lambda \gamma \nu_k }{ (\nu_k^2 - \gamma^2)\beta} & k \geq 1\\ + \end{cases} \;. + +The function :meth:`approx_by_matsubara<.DrudeLorentzEnvironment.approx_by_matsubara>` creates a multi-exponential +approximation to the Drude-Lorentz environment by truncating this series at a finite index :math:`N_k`. +This approximation can then be used with the HEOM solver, for example. +The :ref:`HEOM section` of this guide contains further examples using the Drude-Lorentz enviroment. + +Similarly, the function :meth:`approx_by_pade<.DrudeLorentzEnvironment.approx_by_pade>` can be used to apply +and truncate the numerically more efficient Padé expansion. + +Underdamped Environment +^^^^^^^^^^^^^^^^^^^^^^^ + +Underdamped environments can be constructed in QuTiP +using the class :class:`.UnderDampedEnvironment`. They are characterized by spectral densities of the form + +.. math:: + J(\omega) = \frac{\lambda^{2} \Gamma \omega}{(\omega_0^{2}- + \omega^{2})^{2}+ \Gamma^{2} \omega^{2}} , + +where :math:`\lambda`, :math:`\Gamma` and :math:`\omega_0` are the coupling strength +(with dimension :math:`(\text{energy})^{3/2}`), the cutoff frequency and the resonance frequency. + +Similar to the Drude-Lorentz environment, the correlation function can be approximated by a +Matsubara expansion. This functionality is available with the +:meth:`approx_by_matsubara<.UnderDampedEnvironment.approx_by_matsubara>` function. + +For small temperatures, the Matsubara expansion converges slowly. It is recommended to instead use a fitting procedure +for the Matsubara contribution as described in [Lambert19]_. + + +User-Defined Environments +------------------------- + +As stated in the introduction, a bosonic environment is fully characterized +by its temperature and spectral density (SD), or alternatively by its correlation function (CF) +or its power spectrum (PS). QuTiP allows for the creation of an user-defined environment by +specifying either the spectral density, the correlation function, or the power spectrum. + +QuTiP then computes the other two functions based on the provided one. To do so, it converts between +the SD and the PS using the formula +:math:`S(\omega) = \operatorname{sign}(\omega)\, J(|\omega|) \bigl[ \coth( \beta\omega / 2 ) + 1 \bigr]` +introduced earlier, and between the PS and the CF using the fast Fourier transform. +The former conversion requires the bath temperature to be specified; the latter requires a cutoff frequency (or cutoff time) +to be provided together with the specified function (SD, CF or PS). +In this way, all characteristic functions can be computed from the specified one. + +The following example manually creates an environment with an underdamped spectral density. +It then compares the correlation function obtained via fast Fourier transformation with the Matsubara expansion. +The slow convergence of the Matsubara expansion is visible around :math:`t=0`. + +.. plot:: + :context: close-figs + + # Define underdamped environment parameters + T = 0.1 + lam = 1 + gamma = 2 + w0 = 5 + + # User-defined environment based on SD + def underdamped_sd(w): + return lam**2 * gamma * w / ((w**2 - w0**2)**2 + (gamma*w)**2) + env = qt.BosonicEnvironment.from_spectral_density(underdamped_sd, wMax=50, T=T) + + tlist = np.linspace(-2, 2, 250) + plt.plot(tlist, np.real(env.correlation_function(tlist)), label='FFT') + + # Pre-defined environment and Matsubara approximations + env2 = qt.UnderDampedEnvironment(T, lam, gamma, w0) + for Nk in range(0, 11, 2): + approx_env = env2.approx_by_matsubara(Nk) + plt.plot(tlist, np.real(approx_env.correlation_function(tlist)), label=f'Nk={Nk}') + + plt.xlabel('Time') + plt.ylabel('Correlation function (real part)') + plt.tight_layout() + plt.legend() + + +Multi-Exponential Approximations +-------------------------------- + +Many approaches to simulating the dynamics of an open quantum system strongly coupled to an environment +assume that the environment correlation function can be approximated by a multi-exponential expansion like + +.. math:: + C(t) = C_R(t) + \mathrm i C_I(t) , \qquad + C_{R,I}(t) = \sum_{k=1}^{N_{R,I}} c^{R,I}_k \exp[-\nu^{R,I}_k t] + +with small numbers :math:`N_{R,I}` of exponents. +Note that :math:`C_R(t)` and :math:`C_I(t)` are the real and imaginary parts of the correlation function, +but the coefficients :math:`c^{R,I}_k` and exponents :math:`\nu^{R,I}_k` are not required to be real in general. + +In the previous sections, various methods of obtaining multi-exponential approximations were introduced. +The output of these approximation functions are :class:`.ExponentialBosonicEnvironment` objects. +An :class:`.ExponentialBosonicEnvironment` is basically a collection of :class:`.CFExponent` s, which store (in the bosonic case) +the coefficient, the exponent, and whether the exponent contributes to the real part, the imaginary part, or both. +As we have already seen above, one can then compute the spectral density, correlation function and power spectrum corresponding +to the exponents, in order to compare them to the original, exact environment. + +Let :math:`c_k \mathrm e^{-\nu_k t}` be a term in the correlation function (i.e., :math:`c_k = c^R_k` or :math:`c_k = \mathrm i c^I_k`). +The corresponding term in the power spectrum is + +.. math:: + S_k(\omega) = 2\Re\Bigr[ \frac{c_k}{\nu_k - \mathrm i\omega} \Bigr] + +and, if a temperature has been specified, the corresponding term in the spectral density can be computed as described above. + +The following example shows how to manually create an :class:`.ExponentialBosonicEnvironment` for the simple example +:math:`C(t) = c \mathrm e^{-\nu t}` with real :math:`c`, :math:`\nu`. The power spectrum then is a Lorentzian, +:math:`S(\omega) = 2c\nu / (\nu^2 + \omega^2)`. + +.. plot:: + :context: close-figs + + c = 1 + nu = 2 + wlist = np.linspace(-3, 3, 250) + + env = qt.ExponentialBosonicEnvironment([c], [nu], [], []) + + plt.figure(figsize=(4, 3)) + plt.plot(wlist, env.power_spectrum(wlist)) + plt.plot(wlist, 2 * c * nu / (nu**2 + wlist**2), '--') + plt.xlabel('Frequency') + plt.ylabel('Power spectrum') + plt.tight_layout() + + +.. _fermionic environments guide: + +Fermionic environments +---------------------- + +The implementation of fermionic environments in QuTiP is not yet as advanced as the bosonic environments. +Currently, user-defined fermionic environments and fitting are not implemented. + +However, the overall structure of fermionic environments in QuTiP is analogous to the bosonic environments. +There is one pre-defined fermionic environment, the Lorentzian environment, and multi-exponential fermionic environments. +Lorentzian environments can be approximated by multi-exponential Matsubara or Padé expansions. + +In the fermionic case, we consider the number-conserving Hamiltonian + +.. math:: + H = H_s + (B^\dagger c + c^\dagger B) + \sum_k \hbar\omega_k\, b^\dagger_k b_k , \qquad + B = \sum_k f_k b_k , + +where the bath operators :math:`b_k` and the system operator :math:`c` obey fermionic anti-commutation relations. +In analogy to the bosonic case, we define the spectral density + +.. math:: + J(\omega) = 2\pi \sum_k f_k^2\, \delta[\omega - \omega_k] , + +which may however now be defined for all (including negative) frequencies, since the spectrum of each mode is bounded. + +The fermionic environment is characterized either by its spectral density, inverse temperature :math:`\beta` and chemical potential :math:`\mu`, +or equivalently by two correlation functions or by two power spectra. The correlation functions are + +.. math:: + C^\sigma(t) = \langle B^\sigma(t) B^{-\sigma} \rangle + = \int_{-\infty}^\infty \frac{\mathrm d\omega}{2\pi}\, J(\omega)\, + \mathrm e^{\sigma \mathrm i\omega t}\, f_F(\sigma \beta[\omega - \mu]) , + +where :math:`\sigma = \pm 1`, :math:`B^+ = B^\dagger` and :math:`B^- = B`. +Further, :math:`f_F(x) = (\mathrm e^x + 1)^{-1}` is the Fermi-Dirac function. +Note that we still have :math:`C^\sigma(-t) = C^\sigma(t)^\ast`. +The power spectra are the Fourier transformed correlation functions, + +.. math:: + S^\sigma(\omega) = \int_{-\infty}^\infty \mathrm dt\, C^\sigma(t)\, \mathrm e^{-\sigma \mathrm i\omega t} + = J(\omega) f_F(\sigma\beta[\omega - \mu]) . + +Since :math:`f_F(x) + f_F(-x) = 1`, we have :math:`S^+(\omega) + S^-(\omega) = J(\omega)`. + +.. note:: + The relationship between the spectral density and the two power spectra (or the two correlation functions) is not one-to-one. + A pair of functions :math:`S^\pm(\omega)` is physical if they satisfy the condition + + .. math:: + S^-(\omega) = \mathrm e^{\beta(\omega - \mu)}\, S^+(\omega) . + + For the correlation functions, the condition becomes :math:`C^-(t) = \mathrm e^{-\beta\mu}\, C^+(t - \mathrm i\beta)^\ast`. + For flexibility, we do not enforce the power spectra / correlation functions to be physical in this sense. + +.. _lorentzian env guide: + +Lorentzian Environment +^^^^^^^^^^^^^^^^^^^^^^ + +Fermionic Lorentzian environments are represented by the class :class:`.LorentzianEnvironment`. +They are characterized by spectral densities of the form + +.. math:: + J(\omega) = \frac{\gamma W^2}{(\omega - \omega_0)^2 + W^2} , + +where :math:`\gamma` is the coupling strength, :math:`W` the spectral width and :math:`\omega_0` the resonance frequency. +Often, the resonance frequency is taken to be equal to the chemical potential of the environment. + +As with the bosonic Drude-Lorentz environments, multi-exponential approximations of the correlation functions, + +.. math:: + C^\sigma(t) \approx \sum_{k=0}^{N_k} c^\sigma_k e^{- \nu^\sigma_k t} , + +can be obtained using the Matsubara or Padé expansions. +The functions :meth:`approx_by_matsubara<.LorentzianEnvironment.approx_by_matsubara>` and +:meth:`approx_by_pade<.LorentzianEnvironment.approx_by_pade>` implement these approximations in QuTiP, +yielding approximated environments that can be used, for example, with the HEOM solver. +Note that for this type of environment, the Matsubara expansion is very inefficient, converging much more slowly than the Padé expansion. +Typically, at least :math:`N_k \geq 20` is required for good convergence. + +For reference, we tabulate the values of the coefficients and exponents in the following. +For the Matsubara expansion, they are + +.. math:: + + \nu^\sigma_{k} = \begin{cases} + W - \mathrm i \sigma\, \omega_0 & k = 0\\ + \frac{(2k - 1) \pi}{\beta} - \mathrm i \sigma\, \mu & k \geq 1\\ + \end{cases} \;, \qquad + c^\sigma_k = \begin{cases} + \frac{\gamma W}{2} f_F[\sigma\beta(\omega_0 - \mu) + \mathrm i\, \beta W] & k = 0\\ + \frac{\mathrm i \gamma W^2}{\beta} \frac{1}{[\mathrm i \sigma (\omega_0 - \mu) + (2k - 1) \pi / \beta]^2 - W^2} & k \geq 1\\ + \end{cases} \;. + +The Padé decomposition approximates the Fermi distribution as: + +.. math:: + + f_F(x) \approx f_F^{\mathrm{approx}}(x) = \frac{1}{2} - \sum_{k=0}^{N_k} \frac{2\kappa_k x}{x^2 + \epsilon_k^2} + +where :math:`\kappa_k` and :math:`\epsilon_k` are coefficients that depend on :math:`N_k` and are defined in +`J. Chem Phys 133, "Efficient on the fly calculation of time correlation functions in computer simulations" `_. +This approach results in the exponents + +.. math:: + + \nu^\sigma_{k} = \begin{cases} + W - \mathrm i \sigma\, \omega_0 & k = 0\\ + \frac{\epsilon_k}{\beta} - \mathrm i \sigma\, \mu & k \geq 1\\ + \end{cases} \;, \qquad + c^\sigma_k = \begin{cases} + \frac{\gamma W}{2} f_F^{\mathrm{approx}}[\sigma\beta(\omega_0 - \mu) + \mathrm i\, \beta W] & k = 0\\ + \frac{\mathrm i\, \kappa_k \gamma W^2}{\beta} \frac{1}{[\mathrm i\sigma(\omega_0 - \mu) + \epsilon_k / \beta]^2 - W^2} & k \geq 1\\ + \end{cases} \;. + + +Multi-Exponential Fermionic Environment +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Analogous to the :class:`.ExponentialBosonicEnvironment` in the bosonic case, the :class:`.ExponentialFermionicEnvironment` describes fermionic environments +where the correlation functions are given by multi-exponential decompositions, + +.. math:: + C^\sigma(t) \approx \sum_{k=0}^{N_k^\sigma} c^\sigma_k e^{-\nu^\sigma_k t} . + +Like in the bosonic case, the class allows us to automatically compute the spectral density and power spectra that correspond to the +multi-exponential correlation functions. +In this case, they are + +.. math:: + S^\sigma(\omega) = \sum_{k=0}^{N_k^\sigma} 2\Re\Bigr[ \frac{c_k^\sigma}{\nu_k^\sigma + \mathrm i \sigma\, \omega} \Bigr] + +and :math:`J(\omega) = S^+(\omega) + S^-(\omega)`. + + +.. plot:: + :context: reset + :include-source: false + :nofigs: \ No newline at end of file diff --git a/doc/guide/guide.rst b/doc/guide/guide.rst index e450ca737e..dace78b10b 100644 --- a/doc/guide/guide.rst +++ b/doc/guide/guide.rst @@ -13,6 +13,7 @@ Users Guide guide-tensor.rst guide-super.rst guide-dynamics.rst + guide-environments.rst guide-heom.rst guide-steady.rst guide-piqs.rst diff --git a/doc/guide/heom/bosonic.rst b/doc/guide/heom/bosonic.rst index fa56357320..a7af0efa2f 100644 --- a/doc/guide/heom/bosonic.rst +++ b/doc/guide/heom/bosonic.rst @@ -10,7 +10,7 @@ spectral density, :math:`J_D`, are H_{sys} &= \frac{\epsilon \sigma_z}{2} + \frac{\Delta \sigma_x}{2} - J_D &= \frac{2\lambda \gamma \omega}{(\gamma^2 + \omega^2)}, + J_D(\omega) &= \frac{2\lambda \gamma \omega}{\gamma^2 + \omega^2}, We will demonstrate how to describe the bath using two different expansions of the spectral density correlation function (Matsubara's expansion and @@ -23,23 +23,28 @@ the two bath expansions, :class:`~qutip.solver.heom.DrudeLorentzBath` and truncated expansion and show how to include an approximation to all of the remaining terms in the bath expansion. -Afterwards, we will show how to calculate the bath expansion coefficients and to +.. admonition:: Environment API + + We will also explain how to achieve the same results using the :class:`.DrudeLorentzEnvironment` + that was introduced in the :ref:`section on environments `. + The "bath" classes are part of an older API that is less powerful than the "environment" API, + but often more convenient to use when one only uses the HEOM solver and does not need any of the new features. + +Afterwards, we will show how to calculate the correlation function expansion coefficients and to use those coefficients to construct your own bath description so that you can -implement your own bosonic baths. +implement your own bosonic baths / environments. Finally, we will demonstrate how to simulate a system coupled to multiple independent baths, as occurs, for example, in certain photosynthesis processes. -A notebook containing a complete example similar to this one implemented in -BoFiN can be found in -`example notebook 1a `__. +A tutorial notebook containing a complete example similar to this one is the +`HEOM example notebook 1a `_. Describing the system and bath ------------------------------ -First, let us construct the system Hamiltonian, :math:`H_{sys}`, and the initial -system state, ``rho0``: +First, let us construct the system Hamiltonian ``H_sys`` and the initial system state ``rho0``: .. plot:: :context: reset @@ -69,9 +74,8 @@ Now let us describe the bath properties: # System-bath coupling operator: Q = sigmaz() -where :math:`\gamma` (``gamma``), :math:`\lambda` (``lam``) and :math:`T` are -the parameters of a Drude-Lorentz bath, and ``Q`` is the coupling operator -between the system and the bath. +where :math:`\gamma` (``gamma``), :math:`\lambda` (``lam``) and the temperature :math:`T` are +the parameters of a Drude-Lorentz bath, and ``Q`` is the coupling operator between the system and the bath. We may the pass these parameters to either :class:`~qutip.solver.heom.DrudeLorentzBath` or @@ -94,8 +98,28 @@ the bath correlations: # Padé expansion: bath = DrudeLorentzPadeBath(Q, lam, gamma, T, Nk) -Where ``Nk`` is the number of terms to retain within the expansion of the -bath. +Here, ``Nk`` is the number of terms to retain within the expansion of the bath. + +.. admonition:: Environment API + + Using the environment API, we first create an abstract :class:`.DrudeLorentzEnvironment` describing the bath, + and then use its functions to create exponential expansions such as the Matsubara and Pade ones: + + .. plot:: + :context: + :nofigs: + + from qutip.core.environment import DrudeLorentzEnvironment + + env = DrudeLorentzEnvironment(T, lam, gamma) + + # Matsubara expansion: + approx = env.approx_by_matsubara(Nk) + + # Padé expansion: + approx = env.approx_by_pade(Nk) + + Note that the coupling operator ``Q`` is not part of the environment objects. .. _heom-bosonic-system-and-bath-dynamics: @@ -124,26 +148,23 @@ and to calculate the system evolution as a function of time: result = solver.run(rho0, tlist) The ``max_depth`` parameter determines how many levels of the hierarchy to -retain. As a first approximation hierarchy depth may be thought of as similar +retain. As a first approximation, hierarchy depth may be thought of as similar to the order of Feynman Diagrams (both classify terms by increasing number of interactions). The ``result`` is a standard QuTiP results object with the attributes: -- ``times``: the times at which the state was evaluated (i.e. ``tlist``) -- ``states``: the system states at each time -- ``expect``: a list with the values of each ``e_ops`` at each time -- ``e_data``: a dictionary with the values of each ``e_op`` at each time -- ``ado_states``: see below (an instance of - :class:`~qutip.solver.heom.HierarchyADOsState`) +- ``times``: The times at which the state was evaluated (i.e. ``tlist``). +- ``states``: The system states at each time. +- ``expect``: A list with the values of each expectation operator at each time. +- ``e_data``: A dictionary with the values of each expectation operator at each time. +- ``ado_states``: See below (a list of instances of :class:`~qutip.solver.heom.HierarchyADOsState`). If ``ado_return=True`` is passed to ``.run(...)`` the full set of auxilliary density operators (ADOs) that make up the hierarchy at each time will be -returned as ``.ado_states``. We will describe how to use these to determine -other properties, such as system-bath currents, later in the fermionic guide -(see :ref:`heom-determining-currents`). - -If one has a full set of ADOs from a previous call of ``.run(...)`` you may +returned as ``result.ado_states``. We will describe how to use these to determine +other properties, such as system-bath currents, later in the :ref:`fermionic guide `. +If one has a full set of ADOs from a previous call of ``.run(...)``, one may supply it as the initial state of the solver by calling ``.run(result.ado_states[-1], tlist, ado_init=True)``. @@ -151,6 +172,18 @@ As with other QuTiP solvers, if expectation operators or functions are supplied using ``.run(..., e_ops=[...])`` the expectation values are available in ``result.expect`` and ``result.e_data``. +.. admonition:: Environment API + + When using the environment API, one needs to pass the coupling operator + to the HEOM solver together with the approximated environment: + + .. plot:: + :context: + :nofigs: + + solver = HEOMSolver(H_sys, (approx, Q), max_depth=max_depth, options=options) + + Below we run the solver again, but use ``e_ops`` to store the expectation values of the population of the system states and the coherence: @@ -171,14 +204,14 @@ values of the population of the system states and the coherence: result = solver.run(rho0, tlist, e_ops={"11": P11p, "22": P22p, "12": P12p}) # Plot the results: - fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8)) - axes.plot(result.times, result.e_data["11"], 'b', linewidth=2, label="P11") - axes.plot(result.times, result.e_data["12"], 'r', linewidth=2, label="P12") - axes.set_xlabel(r't', fontsize=28) - axes.legend(loc=0, fontsize=12) + fig, axes = plt.subplots(1, 1, sharex=True, figsize=(6, 6)) + axes.plot(result.times, np.real(result.e_data["11"]), 'b', linewidth=2, label="P11") + axes.plot(result.times, np.real(result.e_data["12"]), 'r', linewidth=2, label="P12") + axes.set_xlabel(r't', fontsize=16) + axes.legend(loc=0, fontsize=16) -Steady-state +Steady state ------------ Using the same solver, we can also determine the steady state of the @@ -191,9 +224,9 @@ combined system and bath using: steady_state, steady_ados = solver.steady_state() where ``steady_state`` is the steady state of the system and ``steady_ados`` -if the steady state of the full hierarchy. The ADO states are -described more fully in :ref:`heom-determining-currents` and -:class:`~qutip.solver.heom.HierarchyADOsState`. +is the steady state of the full hierarchy. The ADO states are +described more fully in the section on :ref:`determining currents ` +and in the API documentation for :class:`~qutip.solver.heom.HierarchyADOsState`. Matsubara Terminator @@ -234,62 +267,48 @@ calculating the terminator for a given expansion: This captures the Markovian effect of the remaining terms in the expansion without having to fully model many more terms. -The value ``delta`` is an approximation to the strength of the effect of +The terminator amplitude ``delta`` is an approximation to the strength of the effect of the remaining terms in the expansion (i.e. how strongly the terminator is coupled to the rest of the system). +.. admonition:: Environment API -Matsubara expansion coefficients --------------------------------- + Here, the terminator amplitude can be returned directly by the ``approx_by_matsubara`` and ``approx_by_pade`` methods used earlier. + Based on it, the special function ``environment.system_terminator`` can then be used to construct the terminator Liouvillian: -So far we have relied on the built-in -:class:`~qutip.solver.heom.DrudeLorentzBath` to construct the Drude-Lorentz -bath expansion for us. Now we will calculate the coefficients ourselves and -construct a :class:`~qutip.solver.heom.BosonicBath` directly. A similar -procedure can be used to apply :class:`~qutip.solver.heom.HEOMSolver` to any -bosonic bath for which we can calculate the expansion coefficients. + .. plot:: + :context: + :nofigs: -The real and imaginary parts of the correlation function, :math:`C(t)`, for the -bosonic bath is expanded in an expontential series: + from qutip.core.environment import system_terminator -.. math:: + # Matsubara expansion: + approx, delta = env.approx_by_matsubara(Nk, compute_delta=True) - C(t) &= C_{real}(t) + i C_{imag}(t) + # Padé expansion: + approx, delta = env.approx_by_pade(Nk, compute_delta=True) - C_{real}(t) &= \sum_{k=0}^{\infty} c_{k,real} e^{- \nu_{k,real} t} + # Add terminator to the system Liouvillian: + terminator = system_terminator(Q, delta) + HL = liouvillian(H_sys) + terminator - C_{imag}(t) &= \sum_{k=0}^{\infty} c_{k,imag} e^{- \nu_{k,imag} t} + # Construct solver + solver = HEOMSolver(HL, (approx, Q), max_depth=max_depth, options=options) -In the specific case of Matsubara expansion for the Drude-Lorentz bath, the -coefficients of this expansion are, for the real part, :math:`C_{real}(t)`: -.. math:: - - \nu_{k,real} &= \begin{cases} - \gamma & k = 0\\ - {2 \pi k} / {\beta } & k \geq 1\\ - \end{cases} - - c_{k,real} &= \begin{cases} - \lambda \gamma [\cot(\beta \gamma / 2) ] & k = 0\\ - \frac{4 \lambda \gamma \nu_k }{ (\nu_k^2 - \gamma^2)\beta} & k \geq 1\\ - \end{cases} - -and the imaginary part, :math:`C_{imag}(t)`: - -.. math:: - - \nu_{k,imag} &= \begin{cases} - \gamma & k = 0\\ - 0 & k \geq 1\\ - \end{cases} +Matsubara expansion coefficients +-------------------------------- - c_{k,imag} &= \begin{cases} - - \lambda \gamma & k = 0\\ - 0 & k \geq 1\\ - \end{cases} +So far we have relied on the built-in +:class:`~qutip.solver.heom.DrudeLorentzBath` to construct the Drude-Lorentz +bath expansion for us. Now we will calculate the coefficients ourselves and +construct a :class:`~qutip.solver.heom.BosonicBath` directly. A similar +procedure can be used to apply :class:`~qutip.solver.heom.HEOMSolver` to any +bosonic bath for which we can calculate the expansion coefficients. -And now the same numbers calculated in Python: +The Matsubara expansion of the Drude-Lorentz correlation function is detailed in +the section on the :ref:`Drude-Lorentz Environment
`. +Let us calculate the coefficients and exponents in Python: .. plot:: :context: @@ -331,8 +350,19 @@ After all that, constructing the bath is very straight forward: And we're done! -The :class:`~qutip.solver.heom.BosonicBath` can be used with the -:class:`~qutip.solver.heom.HEOMSolver` in exactly the same way as the baths +.. admonition:: Environment API + + The analogue of the ``BosonicBath`` is the ``ExponentialBosonicEnvironment``. + Its usage is very similar: + + .. code-block:: python + + from qutip.core.environment import ExponentialBosonicEnvironment + + env = ExponentialBosonicEnvironment(ck_real, vk_real, ck_imag, vk_imag) + +The :class:`~qutip.solver.heom.BosonicBath` and the :class:`.ExponentialBosonicEnvironment` can be used with the +:class:`~qutip.solver.heom.HEOMSolver` in exactly the same way as the baths and environments that we constructed previously using the built-in Drude-Lorentz bath expansions. @@ -340,16 +370,15 @@ Multiple baths -------------- The :class:`~qutip.solver.heom.HEOMSolver` supports having a system interact -with multiple environments. All that is needed is to supply a list of baths -instead of a single bath. +with multiple reservoirs. All that is needed is to supply a list of baths or environments +instead of a single bath or environment. In the example below we calculate the evolution of a small system where each basis state of the system interacts with a separate bath. Such an arrangement can model, for example, the Fenna–Matthews–Olson (FMO) -pigment-protein complex which plays an important role in photosynthesis ( -for a full FMO example see the notebook -https://github.com/tehruhn/bofin/blob/main/examples/example-2-FMO-example.ipynb -). +pigment-protein complex which plays an important role in photosynthesis +(for a full FMO example see the +`HEOM example notebook 2 `_). For each bath expansion, we also include the terminator in the system Liouvillian. @@ -399,10 +428,21 @@ occurs: # Plot populations: fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8)) for label, values in result.e_data.items(): - axes.plot(result.times, values, label=label) - axes.set_xlabel(r't', fontsize=28) - axes.set_ylabel(r"Population", fontsize=28) - axes.legend(loc=0, fontsize=12) + axes.plot(result.times, np.real(values), label=label) + axes.set_xlabel(r't', fontsize=16) + axes.set_ylabel(r"Population", fontsize=16) + axes.legend(loc=0, fontsize=16) + +.. admonition:: Environment API + + Instead of a list ``[bath1, bath2, ...]``, one can of course also pass multiple + environments with different coupling operators like + + .. code-block:: python + + HEOMSolver(Hsys, [(env1, Q1), (env2, Q2), ...], ...) + + or even a mixed list of baths and environments. .. plot:: diff --git a/doc/guide/heom/fermionic.rst b/doc/guide/heom/fermionic.rst index 3bef8bafe1..ceabc0f868 100644 --- a/doc/guide/heom/fermionic.rst +++ b/doc/guide/heom/fermionic.rst @@ -11,7 +11,7 @@ are H_{sys} &= c^{\dagger} c - J_D &= \frac{\Gamma W^2}{(w - \mu)^2 + W^2}, + J_D(\omega) &= \frac{\Gamma W^2}{(w - \mu)^2 + W^2}, We will demonstrate how to describe the bath using two different expansions of the spectral density correlation function (Matsubara's expansion and @@ -19,26 +19,29 @@ a Padé expansion), how to evolve the system in time, and how to calculate the steady state. Since our fermion is coupled to two reservoirs, we will construct two baths -- -one for each reservoir or lead -- and call them the left (:math:`L`) and right -(:math:`R`) baths for convenience. Each bath will have a different chemical +one for each reservoir or lead -- and call them the left (L) and right +(R) baths for convenience. Each bath will have a different chemical potential :math:`\mu` which we will label :math:`\mu_L` and :math:`\mu_R`. First we will do this using the built-in implementations of the bath expansions, :class:`~qutip.solver.heom.LorentzianBath` and :class:`~qutip.solver.heom.LorentzianPadeBath`. - Afterwards, we will show how to calculate the bath expansion coefficients and to use those coefficients to construct your own bath description so that you can implement your own fermionic baths. -Our implementation of fermionic baths primarily follows the definitions used by -Christian Schinabeck in his dissertation ( -https://opus4.kobv.de/opus4-fau/files/10984/DissertationChristianSchinabeck.pdf -) and related publications. +.. admonition:: Environment API + + As for the bosonic case, we will include brief intermissions showing how to + achieve the same results using the newer environment API. The "bath" classes + are part of an older API that is less powerful, but often more convenient to + use when one only uses the HEOM solver and does not need any of the new features. -A notebook containing a complete example similar to this one implemented in -BoFiN can be found in `example notebook 4b -`__. +Our implementation of fermionic baths primarily follows the definitions used in the +`dissertation of Christian Schinabeck `_ +and related publications. +A notebook containing a complete example similar to this one implemented in BoFiN can be found in +`HEOM example notebook 5a `_. Describing the system and bath @@ -67,8 +70,8 @@ Now let us describe the bath properties: :nofigs: # Shared bath properties: - gamma = 0.01 # coupling strength - W = 1.0 # cut-off + gamma = 0.01 # coupling strength + W = 1.0 # cut-off T = 0.025851991 # temperature beta = 1. / T @@ -79,10 +82,11 @@ Now let us describe the bath properties: # System-bath coupling operator: Q = destroy(2) -where :math:`\Gamma` (``gamma``), :math:`W` and :math:`T` are the parameters of +where :math:`\Gamma` (``gamma``), :math:`W` and the temperature :math:`T` are the parameters of an Lorentzian bath, :math:`\mu_L` (``mu_L``) and :math:`\mu_R` (``mu_R``) are the chemical potentials of the left and right baths, and ``Q`` is the coupling -operator between the system and the baths. +operator between the system and the baths. QuTiP assumes a number-conserving coupling term, as described +in the :ref:`section on fermionic environments `. We may the pass these parameters to either ``LorentzianBath`` or ``LorentzianPadeBath`` to construct an expansion of the bath correlations: @@ -105,13 +109,35 @@ We may the pass these parameters to either ``LorentzianBath`` or bath_L = LorentzianPadeBath(Q, gamma, W, mu_L, T, Nk, tag="L") bath_R = LorentzianPadeBath(Q, gamma, W, mu_R, T, Nk, tag="R") -Where ``Nk`` is the number of terms to retain within the expansion of the -bath. +Here, ``Nk`` is the number of terms to retain within the expansion of the bath. Note that we haved labelled each bath with a tag (either "L" or "R") so that we can identify the exponents from individual baths later when calculating the currents between the system and the bath. +.. admonition:: Environment API + + In analogy to the Drude-Lorentz environment in the previous section, + we first create a :class:`.LorentzianEnvironment` that describes the bath abstractly + and then invoke approximation methods: + + .. plot:: + :context: + :nofigs: + + from qutip.core.environment import LorentzianEnvironment + + env_L = LorentzianEnvironment(T, mu_L, gamma, W) + env_R = LorentzianEnvironment(T, mu_R, gamma, W) + + # Matsubara expansion: + approx_L = env_L.approx_by_matsubara(Nk, tag="L") + approx_R = env_R.approx_by_matsubara(Nk, tag="R") + + # Pade expansion: + approx_L = env_L.approx_by_pade(Nk, tag="L") + approx_R = env_R.approx_by_pade(Nk, tag="R") + System and bath dynamics ------------------------ @@ -139,12 +165,23 @@ and to calculate the system evolution as a function of time: As in the bosonic case, the ``max_depth`` parameter determines how many levels of the hierarchy to retain. - -As in the bosonic case, we can specify ``e_ops`` in order to retrieve the +Also here, we can specify ``e_ops`` in order to retrieve the expectation values of operators at each given time. See -:ref:`heom-bosonic-system-and-bath-dynamics` for a fuller description of +:ref:`the previous section ` for a fuller description of the returned ``result`` object. +.. admonition:: Environment API + + When using the environment API, we again have to pass the coupling operator + to the HEOM solver together with the approximated environment: + + .. plot:: + :context: + :nofigs: + + solver = HEOMSolver(H_sys, [(approx_L, Q), (approx_R, Q)], + max_depth=max_depth, options=options) + Below we run the solver again, but use ``e_ops`` to store the expectation values of the population of the system states: @@ -161,11 +198,11 @@ values of the population of the system states: result = solver.run(rho0, tlist, e_ops={"11": P11p, "22": P22p}) # Plot the results: - fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8,8)) - axes.plot(result.times, result.e_data["11"], 'b', linewidth=2, label="P11") - axes.plot(result.times, result.e_data["22"], 'r', linewidth=2, label="P22") - axes.set_xlabel(r't', fontsize=28) - axes.legend(loc=0, fontsize=12) + fig, axes = plt.subplots(1, 1, sharex=True, figsize=(6, 6)) + axes.plot(result.times, np.real(result.e_data["11"]), 'b', linewidth=2, label="P11") + axes.plot(result.times, np.real(result.e_data["22"]), 'r', linewidth=2, label="P22") + axes.set_xlabel(r't', fontsize=16) + axes.legend(loc=0, fontsize=16) The plot above is not very exciting. What we would really like to see in this case are the currents between the system and the two baths. We will plot @@ -187,12 +224,11 @@ bath is: .. math:: - \mathrm{Contribution from Exponent} = \pm i \mathrm{Tr}(Q^\pm \cdot A) + \text{Contribution from Exponent} = \pm i \mathrm{Tr}(Q^\pm \cdot A) -where the :math:`\pm` sign is the sign of the exponent (see the -description later in :ref:`heom-fermionic-pade-expansion-coefficients`) and -:math:`Q^\pm` is :math:`Q` for ``+`` exponents and :math:`Q^{\dagger}` for -``-`` exponents. +where the :math:`\pm` sign depends on whether the exponent corresponds to the :math:`C^+(t)` or the :math:`C^-(t)` correlation function +(see the explanation in the guide on :ref:`fermionic environments guide`) and +:math:`Q^\pm` is :math:`Q` for ``+`` exponents and :math:`Q^{\dagger}` for ``-`` exponents. The first-level exponents for the left bath are retrieved by calling ``.filter(tags=["L"])`` on ``ado_state`` which is an instance of @@ -250,10 +286,10 @@ baths, we can pass them to ``e_ops`` and plot the results: result.times, result.e_data["right_currents"], 'r', linewidth=2, label="Bath R", ) - axes.set_xlabel(r't', fontsize=28) - axes.set_ylabel(r'Current', fontsize=20) - axes.set_title(r'System to Bath Currents', fontsize=20) - axes.legend(loc=0, fontsize=12) + axes.set_xlabel(r't', fontsize=16) + axes.set_ylabel(r'Current', fontsize=16) + axes.set_title(r'System to Bath Currents', fontsize=16) + axes.legend(loc=0, fontsize=16) And now we have a more interesting plot that shows the currents to the left and right baths decaying towards their steady states! @@ -320,68 +356,15 @@ steady state is reached! Padé expansion coefficients --------------------------- -We now look at how to calculate the correlation expansion coefficients for the -Lorentzian spectral density ourselves. Once we have calculated the coefficients -we can construct a :class:`~qutip.solver.heom.FermionicBath` directly from -them. A similar procedure can be used to apply -:class:`~qutip.solver.heom.HEOMSolver` to any fermionic bath for which we can -calculate the expansion coefficients. - -In the fermionic case we must descriminate between the order in which -excitations are created within the bath, so we define two different correlation -functions, :math:`C_{+}(t)`, and :math:`C_{-}(t)`: - -.. math:: - - C^{\sigma}(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} d\omega e^{\sigma i \omega t} J(\omega) f_F[\sigma\beta(\omega - \mu)] - -where :math:`\sigma` is either ``+`` or ``-`` and, :math:`f_F` is the Fermi -distribution function, and :math:`J(\omega)` is the Lorentzian spectral density -we defined at the start. +As in the bosonic case, we can also use manually calculated correlation function +coefficients of fermionic environments with the HEOM solver. In the section on +fermionic environments, we :ref:`defined ` their +correlation functions and :ref:`calculated ` the Matsubara and Pade +expansions for Lorentzian environments. -The Fermi distribution function is: - -.. math:: - - f_F (x) = (\exp(x) + 1)^{-1} - -As in the bosonic case we can approximate this integral with a Matsubara or -Padé expansion. For the Lorentzian bath the Padé expansion converges much +For the Lorentzian bath, the Padé expansion converges much more quickly, so we will calculate the Padé expansion coefficients here. - -The Padé decomposition approximates the Fermi distribution as: - -.. math:: - - f_F(x) \approx f_F^{\mathrm{approx}}(x) = \frac{1}{2} - \sum_{l=0}^{Nk} \frac{2k_l x}{x^2 + \epsilon_l^2} - -where :math:`k_l` and :math:`\epsilon_l` are coefficients defined in -`J. Chem Phys 133, "Efficient on the fly calculation of time correlation functions in computer simulations" `_, -and :math:`Nk` specifies the cut-off in the expansion. - -Evaluating the integral for the correlation functions gives: - -.. math:: - - C^{\sigma}(t) \approx \sum_{l=0}^{Nk} \eta^{\sigma,l} e^{-\gamma_{\sigma,l}t} - -where: - -.. math:: - - \eta_{\sigma, l} &= \begin{cases} - \frac{\Gamma W}{2} f_F^{approx}(i\beta W) & l = 0\\ - -i\cdot \frac{k_l}{\beta} \cdot \frac{\Gamma W^2}{-\frac{\epsilon^2_l}{\beta^2} + W^2} & l \neq 0\\ - \end{cases} - - \gamma_{\sigma,l} &= \begin{cases} - W - \sigma i\mu & l = 0\\ - \frac{\epsilon_l}{\beta} - \sigma i \mu & l \neq 0\\ - \end{cases} - -and :math:`\beta = \frac{1}{T}`. - -And now we calculate the same numbers in Python: +In Python code, the calculation can be done as follows: .. plot:: :context: @@ -447,7 +430,7 @@ And now we calculate the same numbers in Python: # correlation function expansions: def C(sigma, mu, Nk): - """ Calculate the expansion coefficients for C_\sigma. """ + r""" Calculate the expansion coefficients for C_\sigma. """ beta = 1. / T ck = [0.5 * gamma * W * f_approx(1.0j * beta * W, Nk)] vk = [W - sigma * 1.0j * mu] @@ -484,6 +467,18 @@ The :class:`~qutip.solver.heom.FermionicBath` can be used with the :class:`~qutip.solver.heom.HEOMSolver` in exactly the same way as the baths we constructed previously using the built-in Lorentzian bath expansions. +.. admonition:: Environment API + + The analogue of the ``FermionicBath`` is the ``ExponentialFermionicEnvironment``. + Its usage is very similar: + + .. code-block:: python + + from qutip.core.environment import ExponentialFermionicEnvironment + + env_L = ExponentialFermionicEnvironment(ck_plus_L, vk_plus_L, ck_minus_L, vk_minus_L) + env_R = ExponentialFermionicEnvironment(ck_plus_R, vk_plus_R, ck_minus_R, vk_minus_R) + .. plot:: :context: reset diff --git a/doc/guide/heom/intro.rst b/doc/guide/heom/intro.rst index fd3c2ea624..fdf0c0e7f6 100644 --- a/doc/guide/heom/intro.rst +++ b/doc/guide/heom/intro.rst @@ -23,21 +23,23 @@ density: .. math:: - J_D = \frac{2\lambda \gamma \omega}{(\gamma^2 + \omega^2)}, + J_D = \frac{2\lambda \gamma \omega}{\gamma^2 + \omega^2}, or an under-damped Brownian motion spectral density: .. math:: - J_U = \frac{\alpha^2 \Gamma \omega}{[(\omega_c^2 - \omega^2)^2 + \Gamma^2 \omega^2]}. + J_U = \frac{\lambda^2 \Gamma \omega}{(\omega_c^2 - \omega^2)^2 + \Gamma^2 \omega^2}. Given the spectral density, the HEOM requires a decomposition of the bath -correlation functions in terms of exponentials. In :doc:`bosonic` we describe +correlation functions in terms of exponentials. +Generally, such decompositions can be generated using the methods available in the :ref:`environment module `. +We will go into some more detail in :doc:`bosonic`, describe how this is done with code examples, and how these expansions are passed to the solver. In addition to support for bosonic environments, QuTiP also provides support for -feriomic environments which is described in :doc:`fermionic`. +fermionic environments which is described in :doc:`fermionic`. Both bosonic and fermionic environments are supported via a single solver, :class:`.HEOMSolver`, that supports solving for both dynamics and steady-states. diff --git a/qutip/core/__init__.py b/qutip/core/__init__.py index 296f26fa40..1a7707c6a8 100644 --- a/qutip/core/__init__.py +++ b/qutip/core/__init__.py @@ -2,6 +2,7 @@ from .coefficient import * from .qobj import * from .cy.qobjevo import * +from .environment import * from .expect import * from .tensor import * from .states import * diff --git a/qutip/core/environment.py b/qutip/core/environment.py new file mode 100644 index 0000000000..9363fb3b4a --- /dev/null +++ b/qutip/core/environment.py @@ -0,0 +1,2747 @@ +""" +Classes that describe environments of open quantum systems +""" + +# Required for Sphinx to follow autodoc_type_aliases +from __future__ import annotations + +__all__ = ['BosonicEnvironment', + 'DrudeLorentzEnvironment', + 'UnderDampedEnvironment', + 'OhmicEnvironment', + 'ExponentialBosonicEnvironment', + 'FermionicEnvironment', + 'LorentzianEnvironment', + 'ExponentialFermionicEnvironment', + 'CFExponent', + 'system_terminator'] + +import abc +import enum +from time import time +from typing import Any, Callable, Literal, Sequence, overload +import warnings + +import numpy as np +from numpy.typing import ArrayLike +from scipy.linalg import eigvalsh +from scipy.interpolate import CubicSpline + +try: + from mpmath import mp + _mpmath_available = True +except ModuleNotFoundError: + _mpmath_available = False + +from ..utilities import (n_thermal, iterated_fit) +from .superoperator import spre, spost +from .qobj import Qobj + + +class BosonicEnvironment(abc.ABC): + """ + The bosonic environment of an open quantum system. It is characterized by + its spectral density and temperature or, equivalently, its power spectrum + or its two-time auto-correlation function. + + Use one of the classmethods :meth:`from_spectral_density`, + :meth:`from_power_spectrum` or :meth:`from_correlation_function` to + construct an environment manually from one of these characteristic + functions, or use a predefined sub-class such as the + :class:`DrudeLorentzEnvironment`, the :class:`UnderDampedEnvironment` or + the :class:`OhmicEnvironment`. + + Bosonic environments offer various ways to approximate the environment with + a multi-exponential correlation function, which can be used for example in + the HEOM solver. The approximated environment is represented as a + :class:`ExponentialBosonicEnvironment`. + + All bosonic environments can be approximated by directly fitting their + correlation function with a multi-exponential ansatz + (:meth:`approx_by_cf_fit`) or by fitting their spectral density with a sum + of Lorentzians (:meth:`approx_by_sd_fit`), which correspond to + underdamped environments with known multi-exponential decompositions. + Subclasses may offer additional approximation methods such as + :meth:`DrudeLorentzEnvironment.approx_by_matsubara` or + :meth:`DrudeLorentzEnvironment.approx_by_pade` in the case of a + Drude-Lorentz environment. + + Parameters + ---------- + T : optional, float + The temperature of this environment. + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + """ + + def __init__(self, T: float = None, tag: Any = None): + self.T = T + self.tag = tag + + @abc.abstractmethod + def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike): + """ + The spectral density of this environment. For negative frequencies, + a value of zero will be returned. See the Users Guide on + :ref:`bosonic environments ` for specifics + on the definitions used by QuTiP. + + If no analytical expression for the spectral density is known, it will + be derived from the power spectrum. In this case, the temperature of + this environment must be specified. + + If no analytical expression for the power spectrum is known either, it + will be derived from the correlation function via a fast fourier + transform. + + Parameters + ---------- + w : array_like or float + The frequencies at which to evaluate the spectral density. + """ + + ... + + @abc.abstractmethod + def correlation_function( + self, t: float | ArrayLike, *, eps: float = 1e-10 + ) -> (float | ArrayLike): + """ + The two-time auto-correlation function of this environment. See the + Users Guide on :ref:`bosonic environments ` + for specifics on the definitions used by QuTiP. + + If no analytical expression for the correlation function is known, it + will be derived from the power spectrum via a fast fourier transform. + + If no analytical expression for the power spectrum is known either, it + will be derived from the spectral density. In this case, the + temperature of this environment must be specified. + + Parameters + ---------- + t : array_like or float + The times at which to evaluate the correlation function. + + eps : optional, float + Used in case the power spectrum is derived from the spectral + density; see the documentation of + :meth:`BosonicEnvironment.power_spectrum`. + """ + + ... + + @abc.abstractmethod + def power_spectrum( + self, w: float | ArrayLike, *, eps: float = 1e-10 + ) -> (float | ArrayLike): + """ + The power spectrum of this environment. See the Users Guide on + :ref:`bosonic environments ` for specifics + on the definitions used by QuTiP. + + If no analytical expression for the power spectrum is known, it will + be derived from the spectral density. In this case, the temperature of + this environment must be specified. + + If no analytical expression for the spectral density is known either, + the power spectrum will instead be derived from the correlation + function via a fast fourier transform. + + Parameters + ---------- + w : array_like or float + The frequencies at which to evaluate the power spectrum. + + eps : optional, float + To derive the zero-frequency power spectrum from the spectral + density, the spectral density must be differentiated numerically. + In that case, this parameter is used as the finite difference in + the numerical differentiation. + """ + + ... + + # --- user-defined environment creation + + @classmethod + def from_correlation_function( + cls, + C: Callable[[float], complex] | ArrayLike, + tlist: ArrayLike = None, + tMax: float = None, + *, + T: float = None, + tag: Any = None, + args: dict[str, Any] = None, + ) -> BosonicEnvironment: + r""" + Constructs a bosonic environment from the provided correlation + function. The provided function will only be used for times + :math:`t \geq 0`. At times :math:`t < 0`, the symmetry relation + :math:`C(-t) = C(t)^\ast` is enforced. + + Parameters + ---------- + C : callable or array_like + The correlation function. Can be provided as a Python function or + as an array. When using a function, the signature should be + + ``C(t: array_like, **args) -> array_like`` + + where ``t`` is time and ``args`` is a dict containing the + other parameters of the function. + + tlist : optional, array_like + The times where the correlation function is sampled (if it is + provided as an array). + + tMax : optional, float + Specifies that the correlation function is essentially zero outside + the interval [-tMax, tMax]. Used for numerical integration + purposes. + + T : optional, float + Environment temperature. (The spectral density of this environment + can only be calculated from the correlation function if a + temperature is provided.) + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + + args : optional, dict + Extra arguments for the correlation function ``C``. + """ + return _BosonicEnvironment_fromCF(C, tlist, tMax, T, tag, args) + + @classmethod + def from_power_spectrum( + cls, + S: Callable[[float], float] | ArrayLike, + wlist: ArrayLike = None, + wMax: float = None, + *, + T: float = None, + tag: Any = None, + args: dict[str, Any] = None, + ) -> BosonicEnvironment: + r""" + Constructs a bosonic environment with the provided power spectrum. + + Parameters + ---------- + S : callable or array_like + The power spectrum. Can be provided as a Python function or + as an array. When using a function, the signature should be + + ``S(w: array_like, **args) -> array_like`` + + where ``w`` is the frequency and ``args`` is a dict containing the + other parameters of the function. + + wlist : optional, array_like + The frequencies where the power spectrum is sampled (if it is + provided as an array). + + wMax : optional, float + Specifies that the power spectrum is essentially zero outside the + interval [-wMax, wMax]. Used for numerical integration purposes. + + T : optional, float + Environment temperature. (The spectral density of this environment + can only be calculated from the powr spectrum if a temperature is + provided.) + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + + args : optional, dict + Extra arguments for the power spectrum ``S``. + """ + return _BosonicEnvironment_fromPS(S, wlist, wMax, T, tag, args) + + @classmethod + def from_spectral_density( + cls, + J: Callable[[float], float] | ArrayLike, + wlist: ArrayLike = None, + wMax: float = None, + *, + T: float = None, + tag: Any = None, + args: dict[str, Any] = None, + ) -> BosonicEnvironment: + r""" + Constructs a bosonic environment with the provided spectral density. + The provided function will only be used for frequencies + :math:`\omega > 0`. At frequencies :math:`\omega \leq 0`, the spectral + density is zero according to the definition used by QuTiP. See the + Users Guide on :ref:`bosonic environments ` + for a note on spectral densities with support at negative frequencies. + + Parameters + ---------- + J : callable or array_like + The spectral density. Can be provided as a Python function or + as an array. When using a function, the signature should be + + ``J(w: array_like, **args) -> array_like`` + + where ``w`` is the frequency and ``args`` is a tuple containing the + other parameters of the function. + + wlist : optional, array_like + The frequencies where the spectral density is sampled (if it is + provided as an array). + + wMax : optional, float + Specifies that the spectral density is essentially zero outside the + interval [-wMax, wMax]. Used for numerical integration purposes. + + T : optional, float + Environment temperature. (The correlation function and the power + spectrum of this environment can only be calculated from the + spectral density if a temperature is provided.) + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + + args : optional, dict + Extra arguments for the spectral density ``S``. + """ + return _BosonicEnvironment_fromSD(J, wlist, wMax, T, tag, args) + + # --- spectral density, power spectrum, correlation function conversions + + def _ps_from_sd(self, w, eps, derivative=None): + # derivative: value of J'(0) + if self.T is None: + raise ValueError( + "The temperature must be specified for this operation.") + + w = np.asarray(w, dtype=float) + if self.T == 0: + return 2 * np.heaviside(w, 0) * self.spectral_density(w) + + # at zero frequency, we do numerical differentiation + # S(0) = 2 J'(0) / beta + zero_mask = (w == 0) + nonzero_mask = np.invert(zero_mask) + + S = np.zeros_like(w) + if derivative is None: + S[zero_mask] = 2 * self.T * self.spectral_density(eps) / eps + else: + S[zero_mask] = 2 * self.T * derivative + S[nonzero_mask] = ( + 2 * np.sign(w[nonzero_mask]) + * self.spectral_density(np.abs(w[nonzero_mask])) + * (n_thermal(w[nonzero_mask], self.T) + 1) + ) + return S.item() if w.ndim == 0 else S + + def _sd_from_ps(self, w): + w = np.asarray(w, dtype=float) + J = np.zeros_like(w) + positive_mask = (w > 0) + power_spectrum = self.power_spectrum(w[positive_mask]) + + if self.T is None: + raise ValueError( + "The temperature must be specified for this operation.") + + J[positive_mask] = ( + power_spectrum / 2 / (n_thermal(w[positive_mask], self.T) + 1) + ) + return J.item() if w.ndim == 0 else J + + def _ps_from_cf(self, w, tMax): + w = np.asarray(w, dtype=float) + if w.ndim == 0: + wMax = np.abs(w) + elif len(w) == 0: + return np.array([]) + else: + wMax = max(np.abs(w[0]), np.abs(w[-1])) + + mirrored_result = _fft(self.correlation_function, wMax, tMax=tMax) + result = np.real(mirrored_result(-w)) + return result.item() if w.ndim == 0 else result + + def _cf_from_ps(self, t, wMax, **ps_kwargs): + t = np.asarray(t, dtype=float) + if t.ndim == 0: + tMax = np.abs(t) + elif len(t) == 0: + return np.array([]) + else: + tMax = max(np.abs(t[0]), np.abs(t[-1])) + + result_fct = _fft(lambda w: self.power_spectrum(w, **ps_kwargs), + tMax, tMax=wMax) + result = result_fct(t) / (2 * np.pi) + return result.item() if t.ndim == 0 else result + + # --- fitting + + def approx_by_cf_fit( + self, + tlist: ArrayLike, + target_rsme: float = 2e-5, + Nr_max: int = 10, + Ni_max: int = 10, + guess: list[float] = None, + lower: list[float] = None, + upper: list[float] = None, + sigma: float | ArrayLike = None, + maxfev: int = None, + full_ansatz: bool = False, + combine: bool = True, + tag: Any = None, + ) -> tuple[ExponentialBosonicEnvironment, dict[str, Any]]: + r""" + Generates an approximation to this environment by fitting its + correlation function with a multi-exponential ansatz. The number of + exponents is determined iteratively based on reducing the normalized + root mean squared error below a given threshold. + + Specifically, the real and imaginary parts are fit by the following + model functions: + + .. math:: + \operatorname{Re}[C(t)] = \sum_{k=1}^{N_r} \operatorname{Re}\Bigl[ + (a_k + \mathrm i d_k) \mathrm e^{(b_k + \mathrm i c_k) t}\Bigl] + , + \\ + \operatorname{Im}[C(t)] = \sum_{k=1}^{N_i} \operatorname{Im}\Bigl[ + (a'_k + \mathrm i d'_k) \mathrm e^{(b'_k + \mathrm i c'_k) t} + \Bigr]. + + If the parameter `full_ansatz` is `False`, :math:`d_k` and :math:`d'_k` + are set to zero and the model functions simplify to + + .. math:: + \operatorname{Re}[C(t)] = \sum_{k=1}^{N_r} + a_k e^{b_k t} \cos(c_{k} t) + , + \\ + \operatorname{Im}[C(t)] = \sum_{k=1}^{N_i} + a'_k e^{b'_k t} \sin(c'_{k} t) . + + The simplified version offers faster fits, however it fails for + anomalous spectral densities with + :math:`\operatorname{Im}[C(0)] \neq 0` as :math:`\sin(0) = 0`. + + Parameters + ---------- + tlist : array_like + The time range on which to perform the fit. + target_rmse : optional, float + Desired normalized root mean squared error (default `2e-5`). Can be + set to `None` to perform only one fit using the maximum number of + modes (`Nr_max`, `Ni_max`). + Nr_max : optional, int + The maximum number of modes to use for the fit of the real part + (default 10). + Ni_max : optional, int + The maximum number of modes to use for the fit of the imaginary + part (default 10). + guess : optional, list of float + Initial guesses for the parameters :math:`a_k`, :math:`b_k`, etc. + The same initial guesses are used for all values of k, and for + the real and imaginary parts. If `full_ansatz` is True, `guess` is + a list of size 4, otherwise, it is a list of size 3. + If none of `guess`, `lower` and `upper` are provided, these + parameters will be chosen automatically. + lower : optional, list of float + Lower bounds for the parameters :math:`a_k`, :math:`b_k`, etc. + The same lower bounds are used for all values of k, and for + the real and imaginary parts. If `full_ansatz` is True, `lower` is + a list of size 4, otherwise, it is a list of size 3. + If none of `guess`, `lower` and `upper` are provided, these + parameters will be chosen automatically. + upper : optional, list of float + Upper bounds for the parameters :math:`a_k`, :math:`b_k`, etc. + The same upper bounds are used for all values of k, and for + the real and imaginary parts. If `full_ansatz` is True, `upper` is + a list of size 4, otherwise, it is a list of size 3. + If none of `guess`, `lower` and `upper` are provided, these + parameters will be chosen automatically. + sigma : optional, float or list of float + Adds an uncertainty to the correlation function of the environment, + i.e., adds a leeway to the fit. This parameter is useful to adjust + if the correlation function is very small in parts of the time + range. For more details, see the documentation of + ``scipy.optimize.curve_fit``. + maxfev : optional, int + Number of times the parameters of the fit are allowed to vary + during the optimization (per fit). + full_ansatz : optional, bool (default False) + If this is set to False, the parameters :math:`d_k` are all set to + zero. The full ansatz, including :math:`d_k`, usually leads to + significantly slower fits, and some manual tuning of the `guesses`, + `lower` and `upper` is usually needed. On the other hand, the full + ansatz can lead to better fits with fewer exponents, especially + for anomalous spectral densities with + :math:`\operatorname{Im}[C(0)] \neq 0` for which the simplified + ansatz will always give :math:`\operatorname{Im}[C(0)] = 0`. + When using the full ansatz with default values for the guesses and + bounds, if the fit takes too long, we recommend choosing guesses + and bounds manually. + combine : optional, bool (default True) + Whether to combine exponents with the same frequency. See + :meth:`combine <.ExponentialBosonicEnvironment.combine>` for + details. + tag : optional, str, tuple or any other object + An identifier (name) for the approximated environment. If not + provided, a tag will be generated from the tag of this environment. + + Returns + ------- + approx_env : :class:`ExponentialBosonicEnvironment` + The approximated environment with multi-exponential correlation + function. + fit_info : dictionary + A dictionary containing the following information about the fit. + + "Nr" + The number of terms used to fit the real part of the + correlation function. + "Ni" + The number of terms used to fit the imaginary part of the + correlation function. + "fit_time_real" + The time the fit of the real part of the correlation function + took in seconds. + "fit_time_imag" + The time the fit of the imaginary part of the correlation + function took in seconds. + "rmse_real" + Normalized mean squared error obtained in the fit of the real + part of the correlation function. + "rmse_imag" + Normalized mean squared error obtained in the fit of the + imaginary part of the correlation function. + "params_real" + The fitted parameters (array of shape Nx3 or Nx4) for the real + part of the correlation function. + "params_imag" + The fitted parameters (array of shape Nx3 or Nx4) for the + imaginary part of the correlation function. + "summary" + A string that summarizes the information about the fit. + """ + + # Process arguments + if tag is None and self.tag is not None: + tag = (self.tag, "CF Fit") + + if full_ansatz: + num_params = 4 + else: + num_params = 3 + + if target_rsme is None: + target_rsme = 0 + Nr_min, Ni_min = Nr_max, Ni_max + else: + Nr_min, Ni_min = 1, 1 + + clist = self.correlation_function(tlist) + if guess is None and lower is None and upper is None: + guess_re, lower_re, upper_re = _default_guess_cfreal( + tlist, np.real(clist), full_ansatz) + guess_im, lower_im, upper_im = _default_guess_cfimag( + np.imag(clist), full_ansatz) + else: + guess_re, lower_re, upper_re = guess, lower, upper + guess_im, lower_im, upper_im = guess, lower, upper + + # Fit real part + start_real = time() + rmse_real, params_real = iterated_fit( + _cf_real_fit_model, num_params, tlist, np.real(clist), target_rsme, + Nr_min, Nr_max, guess=guess_re, lower=lower_re, upper=upper_re, + sigma=sigma, maxfev=maxfev + ) + end_real = time() + fit_time_real = end_real - start_real + + # Fit imaginary part + start_imag = time() + rmse_imag, params_imag = iterated_fit( + _cf_imag_fit_model, num_params, tlist, np.imag(clist), target_rsme, + Ni_min, Ni_max, guess=guess_im, lower=lower_im, upper=upper_im, + sigma=sigma, maxfev=maxfev + ) + end_imag = time() + fit_time_imag = end_imag - start_imag + + # Generate summary + Nr = len(params_real) + Ni = len(params_imag) + full_summary = _cf_fit_summary( + params_real, params_imag, fit_time_real, fit_time_imag, + Nr, Ni, rmse_real, rmse_imag, n=num_params + ) + + fit_info = {"Nr": Nr, "Ni": Ni, "fit_time_real": fit_time_real, + "fit_time_imag": fit_time_imag, "rmse_real": rmse_real, + "rmse_imag": rmse_imag, "params_real": params_real, + "params_imag": params_imag, "summary": full_summary} + + # Finally, generate environment and return + ckAR = [] + vkAR = [] + for term in params_real: + if full_ansatz: + a, b, c, d = term + else: + a, b, c = term + d = 0 + ckAR.extend([(a + 1j * d) / 2, (a - 1j * d) / 2]) + vkAR.extend([-b - 1j * c, -b + 1j * c]) + + ckAI = [] + vkAI = [] + for term in params_imag: + if full_ansatz: + a, b, c, d = term + else: + a, b, c = term + d = 0 + ckAI.extend([-1j * (a + 1j * d) / 2, 1j * (a - 1j * d) / 2]) + vkAI.extend([-b - 1j * c, -b + 1j * c]) + + approx_env = ExponentialBosonicEnvironment( + ckAR, vkAR, ckAI, vkAI, combine=combine, T=self.T, tag=tag) + return approx_env, fit_info + + def approx_by_sd_fit( + self, + wlist: ArrayLike, + Nk: int = 1, + target_rmse: float = 5e-6, + Nmax: int = 10, + guess: list[float] = None, + lower: list[float] = None, + upper: list[float] = None, + sigma: float | ArrayLike = None, + maxfev: int = None, + combine: bool = True, + tag: Any = None, + ) -> tuple[ExponentialBosonicEnvironment, dict[str, Any]]: + r""" + Generates an approximation to this environment by fitting its spectral + density with a sum of underdamped terms. Each underdamped term + effectively acts like an underdamped environment. We use the known + exponential decomposition of the underdamped environment, keeping `Nk` + Matsubara terms for each. The number of underdamped terms is determined + iteratively based on reducing the normalized root mean squared error + below a given threshold. + + Specifically, the spectral density is fit by the following model + function: + + .. math:: + J(\omega) = \sum_{k=1}^{N} \frac{2 a_k b_k \omega}{\left(\left( + \omega + c_k \right)^2 + b_k^2 \right) \left(\left( + \omega - c_k \right)^2 + b_k^2 \right)} + + Parameters + ---------- + wlist : array_like + The frequency range on which to perform the fit. + Nk : optional, int + The number of Matsubara terms to keep in each mode (default 1). + target_rmse : optional, float + Desired normalized root mean squared error (default `5e-6`). Can be + set to `None` to perform only one fit using the maximum number of + modes (`Nmax`). + Nmax : optional, int + The maximum number of modes to use for the fit (default 10). + guess : optional, list of float + Initial guesses for the parameters :math:`a_k`, :math:`b_k` and + :math:`c_k`. The same initial guesses are used for all values of + k. + If none of `guess`, `lower` and `upper` are provided, these + parameters will be chosen automatically. + lower : optional, list of float + Lower bounds for the parameters :math:`a_k`, :math:`b_k` and + :math:`c_k`. The same lower bounds are used for all values of + k. + If none of `guess`, `lower` and `upper` are provided, these + parameters will be chosen automatically. + upper : optional, list of float + Upper bounds for the parameters :math:`a_k`, :math:`b_k` and + :math:`c_k`. The same upper bounds are used for all values of + k. + If none of `guess`, `lower` and `upper` are provided, these + parameters will be chosen automatically. + sigma : optional, float or list of float + Adds an uncertainty to the spectral density of the environment, + i.e., adds a leeway to the fit. This parameter is useful to adjust + if the spectral density is very small in parts of the frequency + range. For more details, see the documentation of + ``scipy.optimize.curve_fit``. + maxfev : optional, int + Number of times the parameters of the fit are allowed to vary + during the optimization (per fit). + combine : optional, bool (default True) + Whether to combine exponents with the same frequency. See + :meth:`combine <.ExponentialBosonicEnvironment.combine>` for + details. + tag : optional, str, tuple or any other object + An identifier (name) for the approximated environment. If not + provided, a tag will be generated from the tag of this environment. + + Returns + ------- + approx_env : :class:`ExponentialBosonicEnvironment` + The approximated environment with multi-exponential correlation + function. + fit_info : dictionary + A dictionary containing the following information about the fit. + + "N" + The number of underdamped terms used in the fit. + "Nk" + The number of Matsubara modes included per underdamped term. + "fit_time" + The time the fit took in seconds. + "rmse" + Normalized mean squared error obtained in the fit. + "params" + The fitted parameters (array of shape Nx3). + "summary" + A string that summarizes the information about the fit. + """ + + # Process arguments + if tag is None and self.tag is not None: + tag = (self.tag, "SD Fit") + + if target_rmse is None: + target_rmse = 0 + Nmin = Nmax + else: + Nmin = 1 + + jlist = self.spectral_density(wlist) + if guess is None and lower is None and upper is None: + guess, lower, upper = _default_guess_sd(wlist, jlist) + + # Fit + start = time() + rmse, params = iterated_fit( + _sd_fit_model, 3, wlist, jlist, target_rmse, Nmin, Nmax, + guess=guess, lower=lower, upper=upper, sigma=sigma, maxfev=maxfev + ) + end = time() + fit_time = end - start + + # Generate summary + N = len(params) + summary = _fit_summary( + fit_time, rmse, N, "the spectral density", params + ) + fit_info = { + "N": N, "Nk": Nk, "fit_time": fit_time, "rmse": rmse, + "params": params, "summary": summary} + + ckAR, vkAR, ckAI, vkAI = [], [], [], [] + # Finally, generate environment and return + for a, b, c in params: + lam = np.sqrt(a + 0j) + gamma = 2 * b + w0 = np.sqrt(c**2 + b**2) + + env = UnderDampedEnvironment(self.T, lam, gamma, w0) + coeffs = env._matsubara_params(Nk) + ckAR.extend(coeffs[0]) + vkAR.extend(coeffs[1]) + ckAI.extend(coeffs[2]) + vkAI.extend(coeffs[3]) + + approx_env = ExponentialBosonicEnvironment( + ckAR, vkAR, ckAI, vkAI, combine=combine, T=self.T, tag=tag) + return approx_env, fit_info + + +class _BosonicEnvironment_fromCF(BosonicEnvironment): + def __init__(self, C, tlist, tMax, T, tag, args): + super().__init__(T, tag) + self._cf = _complex_interpolation( + C, tlist, 'correlation function', args) + if tlist is not None: + self.tMax = max(np.abs(tlist[0]), np.abs(tlist[-1])) + else: + self.tMax = tMax + + def correlation_function(self, t, **kwargs): + t = np.asarray(t, dtype=float) + result = np.zeros_like(t, dtype=complex) + positive_mask = (t >= 0) + non_positive_mask = np.invert(positive_mask) + + result[positive_mask] = self._cf(t[positive_mask]) + result[non_positive_mask] = np.conj( + self._cf(-t[non_positive_mask]) + ) + return result.item() if t.ndim == 0 else result + + def spectral_density(self, w): + return self._sd_from_ps(w) + + def power_spectrum(self, w, **kwargs): + if self.tMax is None: + raise ValueError('The support of the correlation function (tMax) ' + 'must be specified for this operation.') + return self._ps_from_cf(w, self.tMax) + + +class _BosonicEnvironment_fromPS(BosonicEnvironment): + def __init__(self, S, wlist, wMax, T, tag, args): + super().__init__(T, tag) + self._ps = _real_interpolation(S, wlist, 'power spectrum', args) + if wlist is not None: + self.wMax = max(np.abs(wlist[0]), np.abs(wlist[-1])) + else: + self.wMax = wMax + + def correlation_function(self, t, **kwargs): + if self.wMax is None: + raise ValueError('The support of the power spectrum (wMax) ' + 'must be specified for this operation.') + return self._cf_from_ps(t, self.wMax) + + def spectral_density(self, w): + return self._sd_from_ps(w) + + def power_spectrum(self, w, **kwargs): + w = np.asarray(w, dtype=float) + ps = self._ps(w) + return ps.item() if w.ndim == 0 else self._ps(w) + + +class _BosonicEnvironment_fromSD(BosonicEnvironment): + def __init__(self, J, wlist, wMax, T, tag, args): + super().__init__(T, tag) + self._sd = _real_interpolation(J, wlist, 'spectral density', args) + if wlist is not None: + self.wMax = max(np.abs(wlist[0]), np.abs(wlist[-1])) + else: + self.wMax = wMax + + def correlation_function(self, t, *, eps=1e-10): + if self.T is None: + raise ValueError('The temperature must be specified for this ' + 'operation.') + if self.wMax is None: + raise ValueError('The support of the spectral density (wMax) ' + 'must be specified for this operation.') + return self._cf_from_ps(t, self.wMax, eps=eps) + + def spectral_density(self, w): + w = np.asarray(w, dtype=float) + + result = np.zeros_like(w) + positive_mask = (w > 0) + result[positive_mask] = self._sd(w[positive_mask]) + + return result.item() if w.ndim == 0 else result + + def power_spectrum(self, w, *, eps=1e-10): + return self._ps_from_sd(w, eps) + + +class DrudeLorentzEnvironment(BosonicEnvironment): + r""" + Describes a Drude-Lorentz bosonic environment with the spectral density + + .. math:: + + J(\omega) = \frac{2 \lambda \gamma \omega}{\gamma^{2}+\omega^{2}} + + (see Eq. 15 in [BoFiN23]_). + + Parameters + ---------- + T : float + Environment temperature. + + lam : float + Coupling strength. + + gamma : float + Spectral density cutoff frequency. + + Nk : optional, int, defaults to 10 + The number of Pade exponents to be used for the calculation of the + correlation function. + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + """ + + def __init__( + self, T: float, lam: float, gamma: float, *, + Nk: int = 10, tag: Any = None + ): + super().__init__(T, tag) + + self.lam = lam + self.gamma = gamma + self.Nk = Nk + + def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike): + """ + Calculates the Drude-Lorentz spectral density. + + Parameters + ---------- + w : array_like or float + Energy of the mode. + """ + + w = np.asarray(w, dtype=float) + result = np.zeros_like(w) + + positive_mask = (w > 0) + w_mask = w[positive_mask] + result[positive_mask] = ( + 2 * self.lam * self.gamma * w_mask / (self.gamma**2 + w_mask**2) + ) + + return result.item() if w.ndim == 0 else result + + def correlation_function( + self, t: float | ArrayLike, Nk: int = None, **kwargs + ) -> (float | ArrayLike): + """ + Calculates the two-time auto-correlation function of the Drude-Lorentz + environment. The calculation is performed by summing a large number of + exponents of the Pade expansion. + + Parameters + ---------- + t : array_like or float + The time at which to evaluate the correlation function. + Nk : int, optional + The number of exponents to use. If not provided, then the value + that was provided when the class was instantiated is used. + """ + + if self.T == 0: + raise ValueError("The Drude-Lorentz correlation function diverges " + "at zero temperature.") + + t = np.asarray(t, dtype=float) + abs_t = np.abs(t) + Nk = Nk or self.Nk + ck_real, vk_real, ck_imag, vk_imag = self._pade_params(Nk) + + def C(c, v): + return np.sum([ck * np.exp(-np.asarray(vk * abs_t)) + for ck, vk in zip(c, v)], axis=0) + result = C(ck_real, vk_real) + 1j * C(ck_imag, vk_imag) + + result = np.asarray(result, dtype=complex) + result[t < 0] = np.conj(result[t < 0]) + return result.item() if t.ndim == 0 else result + + def power_spectrum( + self, w: float | ArrayLike, **kwargs + ) -> (float | ArrayLike): + """ + Calculates the power spectrum of the Drude-Lorentz environment. + + Parameters + ---------- + w : array_like or float + The frequency at which to evaluate the power spectrum. + """ + + sd_derivative = 2 * self.lam / self.gamma + return self._ps_from_sd(w, None, sd_derivative) + + @overload + def approx_by_matsubara( + self, Nk: int, combine: bool = ..., + compute_delta: Literal[False] = False, tag: Any = ... + ) -> ExponentialBosonicEnvironment: ... + + @overload + def approx_by_matsubara( + self, Nk: int, combine: bool = ..., + compute_delta: Literal[True] = True, tag: Any = ... + ) -> tuple[ExponentialBosonicEnvironment, float]: ... + + def approx_by_matsubara( + self, Nk, combine=True, compute_delta=False, tag=None + ): + """ + Generates an approximation to this environment by truncating its + Matsubara expansion. + + Parameters + ---------- + Nk : int + Number of Matsubara terms to include. In total, the real part of + the correlation function will include `Nk+1` terms and the + imaginary part `1` term. + + combine : bool, default `True` + Whether to combine exponents with the same frequency. + + compute_delta : bool, default `False` + Whether to compute and return the approximation discrepancy + (see below). + + tag : optional, str, tuple or any other object + An identifier (name) for the approximated environment. If not + provided, a tag will be generated from the tag of this environment. + + Returns + ------- + approx_env : :class:`ExponentialBosonicEnvironment` + The approximated environment with multi-exponential correlation + function. + + delta : float, optional + The approximation discrepancy. That is, the difference between the + true correlation function of the Drude-Lorentz environment and the + sum of the ``Nk`` exponential terms is approximately ``2 * delta * + dirac(t)``, where ``dirac(t)`` denotes the Dirac delta function. + It can be used to create a "terminator" term to add to the system + dynamics to take this discrepancy into account, see + :func:`.system_terminator`. + """ + + if self.T == 0: + raise ValueError("The Drude-Lorentz correlation function diverges " + "at zero temperature.") + if tag is None and self.tag is not None: + tag = (self.tag, "Matsubara Truncation") + + lists = self._matsubara_params(Nk) + approx_env = ExponentialBosonicEnvironment( + *lists, T=self.T, combine=combine, tag=tag) + + if not compute_delta: + return approx_env + + delta = 2 * self.lam * self.T / self.gamma - 1j * self.lam + for exp in approx_env.exponents: + delta -= exp.coefficient / exp.exponent + + return approx_env, delta + + @overload + def approx_by_pade( + self, Nk: int, combine: bool = ..., + compute_delta: Literal[False] = False, tag: Any = ... + ) -> ExponentialBosonicEnvironment: ... + + @overload + def approx_by_pade( + self, Nk: int, combine: bool = ..., + compute_delta: Literal[True] = True, tag: Any = ... + ) -> tuple[ExponentialBosonicEnvironment, float]: ... + + def approx_by_pade( + self, Nk, combine=True, compute_delta=False, tag=None + ): + """ + Generates an approximation to this environment by truncating its + Pade expansion. + + Parameters + ---------- + Nk : int + Number of Pade terms to include. In total, the real part of + the correlation function will include `Nk+1` terms and the + imaginary part `1` term. + + combine : bool, default `True` + Whether to combine exponents with the same frequency. + + compute_delta : bool, default `False` + Whether to compute and return the approximation discrepancy + (see below). + + tag : optional, str, tuple or any other object + An identifier (name) for the approximated environment. If not + provided, a tag will be generated from the tag of this environment. + + Returns + ------- + approx_env : :class:`ExponentialBosonicEnvironment` + The approximated environment with multi-exponential correlation + function. + + delta : float, optional + The approximation discrepancy. That is, the difference between the + true correlation function of the Drude-Lorentz environment and the + sum of the ``Nk`` exponential terms is approximately ``2 * delta * + dirac(t)``, where ``dirac(t)`` denotes the Dirac delta function. + It can be used to create a "terminator" term to add to the system + dynamics to take this discrepancy into account, see + :func:`.system_terminator`. + """ + + if self.T == 0: + raise ValueError("The Drude-Lorentz correlation function diverges " + "at zero temperature.") + if tag is None and self.tag is not None: + tag = (self.tag, "Pade Truncation") + + ck_real, vk_real, ck_imag, vk_imag = self._pade_params(Nk) + approx_env = ExponentialBosonicEnvironment( + ck_real, vk_real, ck_imag, vk_imag, + T=self.T, combine=combine, tag=tag + ) + + if not compute_delta: + return approx_env + + delta = 2 * self.lam * self.T / self.gamma - 1j * self.lam + for exp in approx_env.exponents: + delta -= exp.coefficient / exp.exponent + + return approx_env, delta + + def _pade_params(self, Nk): + eta_p, gamma_p = self._corr(Nk) + + ck_real = [np.real(eta) for eta in eta_p] + vk_real = [gam for gam in gamma_p] + # There is only one term in the expansion of the imaginary part of the + # Drude-Lorentz correlation function. + ck_imag = [np.imag(eta_p[0])] + vk_imag = [gamma_p[0]] + + return ck_real, vk_real, ck_imag, vk_imag + + def _matsubara_params(self, Nk): + """ Calculate the Matsubara coefficients and frequencies. """ + ck_real = [self.lam * self.gamma / np.tan(self.gamma / (2 * self.T))] + ck_real.extend([ + (8 * self.lam * self.gamma * self.T * np.pi * k * self.T / + ((2 * np.pi * k * self.T)**2 - self.gamma**2)) + for k in range(1, Nk + 1) + ]) + vk_real = [self.gamma] + vk_real.extend([2 * np.pi * k * self.T for k in range(1, Nk + 1)]) + + ck_imag = [-self.lam * self.gamma] + vk_imag = [self.gamma] + + return ck_real, vk_real, ck_imag, vk_imag + + # --- Pade approx calculation --- + + def _corr(self, Nk): + kappa, epsilon = self._kappa_epsilon(Nk) + + eta_p = [self.lam * self.gamma * + (self._cot(self.gamma / (2 * self.T)) - 1.0j)] + gamma_p = [self.gamma] + + for ll in range(1, Nk + 1): + eta_p.append( + (kappa[ll] * self.T) * 4 * self.lam * + self.gamma * (epsilon[ll] * self.T) + / ((epsilon[ll]**2 * self.T**2) - self.gamma**2) + ) + gamma_p.append(epsilon[ll] * self.T) + + return eta_p, gamma_p + + def _cot(self, x): + return 1. / np.tan(x) + + def _kappa_epsilon(self, Nk): + eps = self._calc_eps(Nk) + chi = self._calc_chi(Nk) + + kappa = [0] + prefactor = 0.5 * Nk * (2 * (Nk + 1) + 1) + for j in range(Nk): + term = prefactor + for k in range(Nk - 1): + term *= ( + (chi[k]**2 - eps[j]**2) / + (eps[k]**2 - eps[j]**2 + self._delta(j, k)) + ) + for k in [Nk - 1]: + term /= (eps[k]**2 - eps[j]**2 + self._delta(j, k)) + kappa.append(term) + + epsilon = [0] + eps + + return kappa, epsilon + + def _delta(self, i, j): + return 1.0 if i == j else 0.0 + + def _calc_eps(self, Nk): + alpha = np.diag([ + 1. / np.sqrt((2 * k + 5) * (2 * k + 3)) + for k in range(2 * Nk - 1) + ], k=1) + alpha += alpha.transpose() + evals = eigvalsh(alpha) + eps = [-2. / val for val in evals[0: Nk]] + return eps + + def _calc_chi(self, Nk): + alpha_p = np.diag([ + 1. / np.sqrt((2 * k + 7) * (2 * k + 5)) + for k in range(2 * Nk - 2) + ], k=1) + alpha_p += alpha_p.transpose() + evals = eigvalsh(alpha_p) + chi = [-2. / val for val in evals[0: Nk - 1]] + return chi + + +class UnderDampedEnvironment(BosonicEnvironment): + r""" + Describes an underdamped environment with the spectral density + + .. math:: + J(\omega) = \frac{\lambda^{2} \Gamma \omega}{(\omega_0^{2}- + \omega^{2})^{2}+ \Gamma^{2} \omega^{2}} + + (see Eq. 16 in [BoFiN23]_). + + Parameters + ---------- + T : float + Environment temperature. + + lam : float + Coupling strength. + + gamma : float + Spectral density cutoff frequency. + + w0 : float + Spectral density resonance frequency. + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + """ + + def __init__( + self, T: float, lam: float, gamma: float, w0: float, *, tag: Any = None + ): + super().__init__(T, tag) + + self.lam = lam + self.gamma = gamma + self.w0 = w0 + + def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike): + """ + Calculates the underdamped spectral density. + + Parameters + ---------- + w : array_like or float + Energy of the mode. + """ + + w = np.asarray(w, dtype=float) + result = np.zeros_like(w) + + positive_mask = (w > 0) + w_mask = w[positive_mask] + result[positive_mask] = ( + self.lam**2 * self.gamma * w_mask / ( + (w_mask**2 - self.w0**2)**2 + (self.gamma * w_mask)**2 + ) + ) + + return result.item() if w.ndim == 0 else result + + def power_spectrum( + self, w: float | ArrayLike, **kwargs + ) -> (float | ArrayLike): + """ + Calculates the power spectrum of the underdamped environment. + + Parameters + ---------- + w : array_like or float + The frequency at which to evaluate the power spectrum. + """ + + sd_derivative = self.lam**2 * self.gamma / self.w0**4 + return self._ps_from_sd(w, None, sd_derivative) + + def correlation_function( + self, t: float | ArrayLike, **kwargs + ) -> (float | ArrayLike): + """ + Calculates the two-time auto-correlation function of the underdamped + environment. + + Parameters + ---------- + t : array_like or float + The time at which to evaluate the correlation function. + """ + # we need an wMax so that spectral density is zero for w>wMax, guess: + wMax = self.w0 + 25 * self.gamma + return self._cf_from_ps(t, wMax) + + def approx_by_matsubara( + self, Nk: int, combine: bool = True, tag: Any = None + ) -> ExponentialBosonicEnvironment: + """ + Generates an approximation to this environment by truncating its + Matsubara expansion. + + Parameters + ---------- + Nk : int + Number of Matsubara terms to include. In total, the real part of + the correlation function will include `Nk+2` terms and the + imaginary part `2` terms. + + combine : bool, default `True` + Whether to combine exponents with the same frequency. + + tag : optional, str, tuple or any other object + An identifier (name) for the approximated environment. If not + provided, a tag will be generated from the tag of this environment. + + Returns + ------- + :class:`ExponentialBosonicEnvironment` + The approximated environment with multi-exponential correlation + function. + """ + + if tag is None and self.tag is not None: + tag = (self.tag, "Matsubara Truncation") + + lists = self._matsubara_params(Nk) + result = ExponentialBosonicEnvironment( + *lists, T=self.T, combine=combine, tag=tag) + return result + + def _matsubara_params(self, Nk): + """ Calculate the Matsubara coefficients and frequencies. """ + + if Nk > 0 and self.T == 0: + warnings.warn("The Matsubara expansion cannot be performed at " + "zero temperature. Use other approaches such as " + "fitting the correlation function.") + Nk = 0 + + Om = np.sqrt(self.w0**2 - (self.gamma / 2)**2) + Gamma = self.gamma / 2 + + z = np.inf if self.T == 0 else (Om + 1j * Gamma) / (2*self.T) + # we set the argument of the hyperbolic tangent to infinity if T=0 + ck_real = ([ + (self.lam**2 / (4 * Om)) * (1 / np.tanh(z)), + (self.lam**2 / (4 * Om)) * (1 / np.tanh(np.conjugate(z))), + ]) + + ck_real.extend([ + (-2 * self.lam**2 * self.gamma * self.T) * (2 * np.pi * k * self.T) + / ( + ((Om + 1j * Gamma)**2 + (2 * np.pi * k * self.T)**2) + * ((Om - 1j * Gamma)**2 + (2 * np.pi * k * self.T)**2) + ) + for k in range(1, Nk + 1) + ]) + + vk_real = [-1j * Om + Gamma, 1j * Om + Gamma] + vk_real.extend([ + 2 * np.pi * k * self.T + for k in range(1, Nk + 1) + ]) + + ck_imag = [ + 1j * self.lam**2 / (4 * Om), + -1j * self.lam**2 / (4 * Om), + ] + + vk_imag = [-1j * Om + Gamma, 1j * Om + Gamma] + + return ck_real, vk_real, ck_imag, vk_imag + + +class OhmicEnvironment(BosonicEnvironment): + r""" + Describes Ohmic environments as well as sub- or super-Ohmic environments + (depending on the choice of the parameter `s`). The spectral density is + + .. math:: + J(\omega) + = \alpha \frac{\omega^s}{\omega_c^{s-1}} e^{-\omega / \omega_c} . + + This class requires the `mpmath` module to be installed. + + Parameters + ---------- + T : float + Temperature of the environment. + + alpha : float + Coupling strength. + + wc : float + Cutoff parameter. + + s : float + Power of omega in the spectral density. + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + """ + + def __init__( + self, T: float, alpha: float, wc: float, s: float, *, tag: Any = None + ): + super().__init__(T, tag) + + self.alpha = alpha + self.wc = wc + self.s = s + + if _mpmath_available is False: + warnings.warn( + "The mpmath module is required for some operations on " + "Ohmic environments, but it is not installed.") + + def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike): + r""" + Calculates the spectral density of the Ohmic environment. + + Parameters + ---------- + w : array_like or float + Energy of the mode. + """ + + w = np.asarray(w, dtype=float) + result = np.zeros_like(w) + + positive_mask = (w > 0) + w_mask = w[positive_mask] + result[positive_mask] = ( + self.alpha * w_mask ** self.s + / (self.wc ** (self.s - 1)) + * np.exp(-np.abs(w_mask) / self.wc) + ) + + return result.item() if w.ndim == 0 else result + + def power_spectrum( + self, w: float | ArrayLike, **kwargs + ) -> (float | ArrayLike): + """ + Calculates the power spectrum of the Ohmic environment. + + Parameters + ---------- + w : array_like or float + The frequency at which to evaluate the power spectrum. + """ + if self.s > 1: + sd_derivative = 0 + elif self.s == 1: + sd_derivative = self.alpha + else: + sd_derivative = np.inf + return self._ps_from_sd(w, None, sd_derivative) + + def correlation_function( + self, t: float | ArrayLike, **kwargs + ) -> (float | ArrayLike): + r""" + Calculates the correlation function of an Ohmic environment using the + formula + + .. math:: + C(t)= \frac{1}{\pi} \alpha w_{c}^{1-s} \beta^{-(s+1)} \Gamma(s+1) + \left[ \zeta\left(s+1,\frac{1+\beta w_{c} -i w_{c} t}{\beta w_{c}} + \right) +\zeta\left(s+1,\frac{1+ i w_{c} t}{\beta w_{c}}\right) + \right] , + + where :math:`\Gamma` is the gamma function, and :math:`\zeta` the + Riemann zeta function. + + Parameters + ---------- + t : array_like or float + The time at which to evaluate the correlation function. + """ + t = np.asarray(t, dtype=float) + t_was_array = t.ndim > 0 + if not t_was_array: + t = np.array([t], dtype=float) + + if self.T != 0: + corr = (self.alpha * self.wc ** (1 - self.s) / np.pi + * mp.gamma(self.s + 1) * self.T ** (self.s + 1)) + z1_u = ((1 + self.wc / self.T - 1j * self.wc * t) + / (self.wc / self.T)) + z2_u = (1 + 1j * self.wc * t) / (self.wc / self.T) + result = np.asarray( + [corr * (mp.zeta(self.s + 1, u1) + mp.zeta(self.s + 1, u2)) + for u1, u2 in zip(z1_u, z2_u)], + dtype=np.cdouble + ) + else: + corr = (self.alpha * self.wc**(self.s + 1) / np.pi + * mp.gamma(self.s + 1) + * (1 + 1j * self.wc * t) ** (-self.s - 1)) + result = np.asarray(corr, dtype=np.cdouble) + + if t_was_array: + return result + return result[0] + + +class CFExponent: + """ + Represents a single exponent (naively, an excitation mode) within an + exponential decomposition of the correlation function of a environment. + + Parameters + ---------- + type : {"R", "I", "RI", "+", "-"} or one of `CFExponent.types` + The type of exponent. + + "R" and "I" are bosonic exponents that appear in the real and + imaginary parts of the correlation expansion, respectively. + + "RI" is a combined bosonic exponent that appears in both the real + and imaginary parts of the correlation expansion. The combined exponent + has a single ``vk``. The ``ck`` is the coefficient in the real + expansion and ``ck2`` is the coefficient in the imaginary expansion. + + "+" and "-" are fermionic exponents. + + ck : complex + The coefficient of the excitation term. + + vk : complex + The frequency of the exponent of the excitation term. + + ck2 : optional, complex + For exponents of type "RI" this is the coefficient of the term in the + imaginary expansion (and ``ck`` is the coefficient in the real + expansion). + + tag : optional, str, tuple or any other object + A label for the exponent (often the name of the environment). It + defaults to None. + + Attributes + ---------- + fermionic : bool + True if the type of the exponent is a Fermionic type (i.e. either + "+" or "-") and False otherwise. + + coefficient : complex + The coefficient of this excitation term in the total correlation + function (including real and imaginary part). + + exponent : complex + The frequency of the exponent of the excitation term. (Alias for `vk`.) + + All of the parameters are also available as attributes. + """ + types = enum.Enum("ExponentType", ["R", "I", "RI", "+", "-"]) + + def _check_ck2(self, type, ck2): + if type == self.types["RI"]: + if ck2 is None: + raise ValueError("RI exponents require ck2") + else: + if ck2 is not None: + raise ValueError( + "Second co-efficient (ck2) should only be specified for" + " RI exponents" + ) + + def _type_is_fermionic(self, type): + return type in (self.types["+"], self.types["-"]) + + def __init__( + self, type: str | CFExponent.ExponentType, + ck: complex, vk: complex, ck2: complex = None, tag: Any = None + ): + if not isinstance(type, self.types): + type = self.types[type] + self._check_ck2(type, ck2) + + self.type = type + self.ck = ck + self.vk = vk + self.ck2 = ck2 + + self.tag = tag + self.fermionic = self._type_is_fermionic(type) + + def __repr__(self): + return ( + f"<{self.__class__.__name__} type={self.type.name}" + f" ck={self.ck!r} vk={self.vk!r} ck2={self.ck2!r}" + f" fermionic={self.fermionic!r}" + f" tag={self.tag!r}>" + ) + + @property + def coefficient(self) -> complex: + coeff = 0 + if self.type != self.types['I']: + coeff += self.ck + else: + coeff += 1j * self.ck + if self.type == self.types['RI']: + coeff += 1j * self.ck2 + return coeff + + @property + def exponent(self) -> complex: + return self.vk + + def _can_combine(self, other, rtol, atol): + if type(self) is not type(other): + return False + if self.fermionic or other.fermionic: + return False + if not np.isclose(self.vk, other.vk, rtol=rtol, atol=atol): + return False + return True + + def _combine(self, other, **init_kwargs): + # Assumes can combine was checked + cls = type(self) + + if self.type == other.type and self.type != self.types['RI']: + # Both R or both I + return cls(type=self.type, ck=(self.ck + other.ck), + vk=self.vk, tag=self.tag, **init_kwargs) + + # Result will be RI + real_part_coefficient = 0 + imag_part_coefficient = 0 + for exp in [self, other]: + if exp.type == self.types['RI'] or exp.type == self.types['R']: + real_part_coefficient += exp.ck + if exp.type == self.types['I']: + imag_part_coefficient += exp.ck + if exp.type == self.types['RI']: + imag_part_coefficient += exp.ck2 + + return cls(type=self.types['RI'], ck=real_part_coefficient, vk=self.vk, + ck2=imag_part_coefficient, tag=self.tag, **init_kwargs) + + +class ExponentialBosonicEnvironment(BosonicEnvironment): + """ + Bosonic environment that is specified through an exponential decomposition + of its correlation function. The list of coefficients and exponents in + the decomposition may either be passed through the four lists `ck_real`, + `vk_real`, `ck_imag`, `vk_imag`, or as a list of bosonic + :class:`CFExponent` objects. + + Parameters + ---------- + ck_real : list of complex + The coefficients of the expansion terms for the real part of the + correlation function. The corresponding frequencies are passed as + vk_real. + + vk_real : list of complex + The frequencies (exponents) of the expansion terms for the real part of + the correlation function. The corresponding coefficients are passed as + ck_real. + + ck_imag : list of complex + The coefficients of the expansion terms in the imaginary part of the + correlation function. The corresponding frequencies are passed as + vk_imag. + + vk_imag : list of complex + The frequencies (exponents) of the expansion terms for the imaginary + part of the correlation function. The corresponding coefficients are + passed as ck_imag. + + exponents : list of :class:`CFExponent` + The expansion coefficients and exponents of both the real and the + imaginary parts of the correlation function as :class:`CFExponent` + objects. + + combine : bool, default True + Whether to combine exponents with the same frequency. See + :meth:`combine` for details. + + T: optional, float + The temperature of the environment. + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + """ + + _make_exponent = CFExponent + + def _check_cks_and_vks(self, ck_real, vk_real, ck_imag, vk_imag): + # all None: returns False + # all provided and lengths match: returns True + # otherwise: raises ValueError + lists = [ck_real, vk_real, ck_imag, vk_imag] + if all(x is None for x in lists): + return False + if any(x is None for x in lists): + raise ValueError( + "If any of the exponent lists ck_real, vk_real, ck_imag, " + "vk_imag is provided, all must be provided." + ) + if len(ck_real) != len(vk_real) or len(ck_imag) != len(vk_imag): + raise ValueError( + "The exponent lists ck_real and vk_real, and ck_imag and " + "vk_imag must be the same length." + ) + return True + + def __init__( + self, + ck_real: ArrayLike = None, vk_real: ArrayLike = None, + ck_imag: ArrayLike = None, vk_imag: ArrayLike = None, + *, + exponents: Sequence[CFExponent] = None, + combine: bool = True, T: float = None, tag: Any = None + ): + super().__init__(T, tag) + + lists_provided = self._check_cks_and_vks( + ck_real, vk_real, ck_imag, vk_imag) + if exponents is None and not lists_provided: + raise ValueError( + "Either the parameter `exponents` or the parameters " + "`ck_real`, `vk_real`, `ck_imag`, `vk_imag` must be provided." + ) + if exponents is not None and any(exp.fermionic for exp in exponents): + raise ValueError( + "Fermionic exponent passed to exponential bosonic environment." + ) + + exponents = exponents or [] + if lists_provided: + exponents.extend(self._make_exponent("R", ck, vk, tag=tag) + for ck, vk in zip(ck_real, vk_real)) + exponents.extend(self._make_exponent("I", ck, vk, tag=tag) + for ck, vk in zip(ck_imag, vk_imag)) + + if combine: + exponents = self.combine(exponents) + self.exponents = exponents + + @classmethod + def combine( + cls, exponents: Sequence[CFExponent], + rtol: float = 1e-5, atol: float = 1e-7 + ) -> Sequence[CFExponent]: + """ + Group bosonic exponents with the same frequency and return a + single exponent for each frequency present. + + Parameters + ---------- + exponents : list of :class:`CFExponent` + The list of exponents to combine. + + rtol : float, default 1e-5 + The relative tolerance to use to when comparing frequencies. + + atol : float, default 1e-7 + The absolute tolerance to use to when comparing frequencies. + + Returns + ------- + list of :class:`CFExponent` + The new reduced list of exponents. + """ + remaining = exponents[:] + new_exponents = [] + + while remaining: + new_exponent = remaining.pop(0) + for other_exp in remaining[:]: + if new_exponent._can_combine(other_exp, rtol, atol): + new_exponent = new_exponent._combine(other_exp) + remaining.remove(other_exp) + new_exponents.append(new_exponent) + + return new_exponents + + def correlation_function( + self, t: float | ArrayLike, **kwargs + ) -> (float | ArrayLike): + """ + Computes the correlation function represented by this exponential + decomposition. + + Parameters + ---------- + t : array_like or float + The time at which to evaluate the correlation function. + """ + + t = np.asarray(t, dtype=float) + corr = np.zeros_like(t, dtype=complex) + + for exp in self.exponents: + corr += exp.coefficient * np.exp(-exp.exponent * np.abs(t)) + corr[t < 0] = np.conj(corr[t < 0]) + + return corr.item() if t.ndim == 0 else corr + + def power_spectrum( + self, w: float | ArrayLike, **kwargs + ) -> (float | ArrayLike): + """ + Calculates the power spectrum corresponding to the multi-exponential + correlation function. + + Parameters + ---------- + w : array_like or float + The frequency at which to evaluate the power spectrum. + """ + + w = np.asarray(w, dtype=float) + S = np.zeros_like(w) + + for exp in self.exponents: + S += 2 * np.real( + exp.coefficient / (exp.exponent - 1j * w) + ) + + return S.item() if w.ndim == 0 else S + + def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike): + """ + Calculates the spectral density corresponding to the multi-exponential + correlation function. + + Parameters + ---------- + w : array_like or float + Energy of the mode. + """ + + return self._sd_from_ps(w) + + +def system_terminator(Q: Qobj, delta: float) -> Qobj: + """ + Constructs the terminator for a given approximation discrepancy. + + Parameters + ---------- + Q : :class:`Qobj` + The system coupling operator. + + delta : float + The approximation discrepancy of approximating an environment with a + finite number of exponentials, see for example + :meth:`.DrudeLorentzEnvironment.approx_by_matsubara`. + + Returns + ------- + terminator : :class:`Qobj` + A superoperator acting on the system Hilbert space. Liouvillian term + representing the contribution to the system-environment dynamics of all + neglected expansion terms. It should be used by adding it to the system + Liouvillian (i.e. ``liouvillian(H_sys)``). + """ + op = 2 * spre(Q) * spost(Q.dag()) - spre(Q.dag() * Q) - spost(Q.dag() * Q) + return delta * op + + +# --- utility functions --- + +def _real_interpolation(fun, xlist, name, args=None): + args = args or {} + if callable(fun): + return lambda w: fun(w, **args) + else: + if xlist is None or len(xlist) != len(fun): + raise ValueError("A list of x-values with the same length must be " + f"provided for the discretized function ({name})") + return CubicSpline(xlist, fun) + + +def _complex_interpolation(fun, xlist, name, args=None): + args = args or {} + if callable(fun): + return lambda t: fun(t, **args) + else: + real_interp = _real_interpolation(np.real(fun), xlist, name) + imag_interp = _real_interpolation(np.imag(fun), xlist, name) + return lambda x: real_interp(x) + 1j * imag_interp(x) + + +def _fft(f, wMax, tMax): + r""" + Calculates the Fast Fourier transform of the given function. We calculate + Fourier transformations via FFT because numerical integration is often + noisy in the scenarios we are interested in. + + Given a (mathematical) function `f(t)`, this function approximates its + Fourier transform + + .. math:: + g(\omega) = \int_{-\infty}^\infty dt\, e^{-i\omega t}\, f(t) . + + The function f is sampled on the interval `[-tMax, tMax]`. The sampling + discretization is chosen as `dt = pi / (4*wMax)` (Shannon-Nyquist + some + leeway). However, `dt` is always chosen small enough to have at least 500 + samples on the interval `[-tMax, tMax]`. + + Parameters + ---------- + wMax: float + Maximum frequency of interest + tMax: float + Support of the function f (i.e., f(t) is essentially zero for + `|t| > tMax`). + + Returns + ------- + The fourier transform of the provided function as an interpolated function. + """ + # Code adapted from https://stackoverflow.com/a/24077914 + + numSamples = int( + max(500, np.ceil(2 * tMax * 4 * wMax / np.pi + 1)) + ) + t, dt = np.linspace(-tMax, tMax, numSamples, retstep=True) + f_values = f(t) + + # Compute Fourier transform by numpy's FFT function + g = np.fft.fft(f_values) + # frequency normalization factor is 2 * np.pi / dt + w = np.fft.fftfreq(numSamples) * 2 * np.pi / dt + # In order to get a discretisation of the continuous Fourier transform + # we need to multiply g by a phase factor + g *= dt * np.exp(1j * w * tMax) + + return _complex_interpolation( + np.fft.fftshift(g), np.fft.fftshift(w), 'FFT' + ) + + +def _cf_real_fit_model(tlist, a, b, c, d=0): + return np.real((a + 1j * d) * np.exp((b + 1j * c) * np.abs(tlist))) + + +def _cf_imag_fit_model(tlist, a, b, c, d=0): + return np.sign(tlist) * np.imag( + (a + 1j * d) * np.exp((b + 1j * c) * np.abs(tlist)) + ) + + +def _default_guess_cfreal(tlist, clist, full_ansatz): + corr_abs = np.abs(clist) + corr_max = np.max(corr_abs) + tc = 2 / np.max(tlist) + + # Checks if constant array, and assigns zero + if (clist == clist[0]).all(): + if full_ansatz: + return [[0] * 4]*3 + return [[0] * 3]*3 + + if full_ansatz: + lower = [-100 * corr_max, -np.inf, -np.inf, -100 * corr_max] + guess = [corr_max, -100*corr_max, 0, 0] + upper = [100*corr_max, 0, np.inf, 100*corr_max] + else: + lower = [-20 * corr_max, -np.inf, 0] + guess = [corr_max, -tc, 0] + upper = [20 * corr_max, 0.1, np.inf] + + return guess, lower, upper + + +def _default_guess_cfimag(clist, full_ansatz): + corr_max = np.max(np.abs(clist)) + # Checks if constant array, and assigns zero + if (clist == clist[0]).all(): + if full_ansatz: + return [[0] * 4]*3 + return [[0] * 3]*3 + + if full_ansatz: + lower = [-100 * corr_max, -np.inf, -np.inf, -100 * corr_max] + guess = [0, -10 * corr_max, 0, 0] + upper = [100 * corr_max, 0, np.inf, 100 * corr_max] + else: + lower = [-20 * corr_max, -np.inf, 0] + guess = [-corr_max, -10 * corr_max, 1] + upper = [10 * corr_max, 0, np.inf] + + return guess, lower, upper + + +def _sd_fit_model(wlist, a, b, c): + return ( + 2 * a * b * wlist / ((wlist + c)**2 + b**2) / ((wlist - c)**2 + b**2) + ) + + +def _default_guess_sd(wlist, jlist): + sd_abs = np.abs(jlist) + sd_max = np.max(sd_abs) + wc = wlist[np.argmax(sd_abs)] + + if sd_max == 0: + return [0] * 3 + + lower = [-100 * sd_max, 0.1 * wc, 0.1 * wc] + guess = [sd_max, wc, wc] + upper = [100 * sd_max, 100 * wc, 100 * wc] + + return guess, lower, upper + + +def _fit_summary(time, rmse, N, label, params, + columns=['lam', 'gamma', 'w0']): + # Generates summary of fit by nonlinear least squares + if len(columns) == 3: + summary = (f"Result of fitting {label} " + f"with {N} terms: \n \n {'Parameters': <10}|" + f"{columns[0]: ^10}|{columns[1]: ^10}|{columns[2]: >5} \n ") + for k in range(N): + summary += ( + f"{k+1: <10}|{params[k][0]: ^10.2e}|{params[k][1]:^10.2e}|" + f"{params[k][2]:>5.2e}\n ") + elif len(columns) == 4: + summary = ( + f"Result of fitting {label} " + f"with {N} terms: \n \n {'Parameters': <10}|" + f"{columns[0]: ^10}|{columns[1]: ^10}|{columns[2]: ^10}" + f"|{columns[3]: >5} \n ") + for k in range(N): + summary += ( + f"{k+1: <10}|{params[k][0]: ^10.2e}|{params[k][1]:^10.2e}" + f"|{params[k][2]:^10.2e}|{params[k][3]:>5.2e}\n ") + else: + raise ValueError("Unsupported number of columns") + summary += (f"\nA normalized RMSE of {rmse: .2e}" + f" was obtained for the {label}.\n") + summary += f"The current fit took {time: 2f} seconds." + return summary + + +def _cf_fit_summary( + params_real, params_imag, fit_time_real, fit_time_imag, Nr, Ni, + rmse_real, rmse_imag, n=3 +): + # Generate nicely formatted summary with two columns for CF fit + columns = ["a", "b", "c"] + if n == 4: + columns.append("d") + summary_real = _fit_summary( + fit_time_real, rmse_real, Nr, + "the real part of\nthe correlation function", + params_real, columns=columns + ) + summary_imag = _fit_summary( + fit_time_imag, rmse_imag, Ni, + "the imaginary part\nof the correlation function", + params_imag, columns=columns + ) + + full_summary = "Correlation function fit:\n\n" + lines_real = summary_real.splitlines() + lines_imag = summary_imag.splitlines() + max_lines = max(len(lines_real), len(lines_imag)) + # Fill the shorter string with blank lines + lines_real = ( + lines_real[:-1] + + (max_lines - len(lines_real)) * [""] + [lines_real[-1]] + ) + lines_imag = ( + lines_imag[:-1] + + (max_lines - len(lines_imag)) * [""] + [lines_imag[-1]] + ) + # Find the maximum line length in each column + max_length1 = max(len(line) for line in lines_real) + max_length2 = max(len(line) for line in lines_imag) + + # Print the strings side by side with a vertical bar separator + for line1, line2 in zip(lines_real, lines_imag): + formatted_line1 = f"{line1:<{max_length1}} |" + formatted_line2 = f"{line2:<{max_length2}}" + full_summary += formatted_line1 + formatted_line2 + "\n" + return full_summary + + +# --- fermionic environments --- + +class FermionicEnvironment(abc.ABC): + r""" + The fermionic environment of an open quantum system. It is characterized by + its spectral density, temperature and chemical potential or, equivalently, + by its power spectra or its two-time auto-correlation functions. + + This class is included as a counterpart to :class:`BosonicEnvironment`, but + it currently does not support all features that the bosonic environment + does. In particular, fermionic environments cannot be constructed from + manually specified spectral densities, power spectra or correlation + functions. The only types of fermionic environment implemented at this time + are Lorentzian environments (:class:`LorentzianEnvironment`) and + environments with multi-exponential correlation functions + (:class:`ExponentialFermionicEnvironment`). + + Parameters + ---------- + T : optional, float + The temperature of this environment. + mu : optional, float + The chemical potential of this environment. + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + """ + + def __init__(self, T: float = None, mu: float = None, tag: Any = None): + self.T = T + self.mu = mu + self.tag = tag + + @abc.abstractmethod + def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike): + r""" + The spectral density of this environment. See the Users Guide on + :ref:`fermionic environments ` for + specifics on the definitions used by QuTiP. + + Parameters + ---------- + w : array_like or float + The frequencies at which to evaluate the spectral density. + """ + + ... + + @abc.abstractmethod + def correlation_function_plus( + self, t: float | ArrayLike + ) -> (float | ArrayLike): + r""" + The "+"-branch of the auto-correlation function of this environment. + See the Users Guide on + :ref:`fermionic environments ` for + specifics on the definitions used by QuTiP. + + Parameters + ---------- + t : array_like or float + The times at which to evaluate the correlation function. + """ + + ... + + @abc.abstractmethod + def correlation_function_minus( + self, t: float | ArrayLike + ) -> (float | ArrayLike): + r""" + The "-"-branch of the auto-correlation function of this environment. + See the Users Guide on + :ref:`fermionic environments ` for + specifics on the definitions used by QuTiP. + + Parameters + ---------- + t : array_like or float + The times at which to evaluate the correlation function. + """ + + ... + + @abc.abstractmethod + def power_spectrum_plus(self, w: float | ArrayLike) -> (float | ArrayLike): + r""" + The "+"-branch of the power spectrum of this environment. See the Users + Guide on :ref:`fermionic environments ` + for specifics on the definitions used by QuTiP. + + Parameters + ---------- + w : array_like or float + The frequencies at which to evaluate the power spectrum. + """ + + ... + + @abc.abstractmethod + def power_spectrum_minus( + self, w: float | ArrayLike + ) -> (float | ArrayLike): + r""" + The "-"-branch of the power spectrum of this environment. See the Users + Guide on :ref:`fermionic environments ` + for specifics on the definitions used by QuTiP. + + Parameters + ---------- + w : array_like or float + The frequencies at which to evaluate the power spectrum. + """ + + ... + + # --- user-defined environment creation + + @classmethod + def from_correlation_functions(cls, **kwargs) -> FermionicEnvironment: + r""" + User-defined fermionic environments are currently not implemented. + """ + + raise NotImplementedError("User-defined fermionic environments are " + "currently not implemented.") + + @classmethod + def from_power_spectra(cls, **kwargs) -> FermionicEnvironment: + r""" + User-defined fermionic environments are currently not implemented. + """ + + raise NotImplementedError("User-defined fermionic environments are " + "currently not implemented.") + + @classmethod + def from_spectral_density(cls, **kwargs) -> FermionicEnvironment: + r""" + User-defined fermionic environments are currently not implemented. + """ + + raise NotImplementedError("User-defined fermionic environments are " + "currently not implemented.") + + +class LorentzianEnvironment(FermionicEnvironment): + r""" + Describes a Lorentzian fermionic environment with the spectral density + + .. math:: + + J(\omega) = \frac{\gamma W^2}{(\omega - \omega_0)^2 + W^2}. + + (see Eq. 46 in [BoFiN23]_). + + Parameters + ---------- + T : float + Environment temperature. + + mu : float + Environment chemical potential. + + gamma : float + Coupling strength. + + W : float + The spectral width of the environment. + + omega0 : optional, float (default equal to ``mu``) + The resonance frequency of the environment. + + Nk : optional, int, defaults to 10 + The number of Pade exponents to be used for the calculation of the + correlation functions. + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + """ + + def __init__( + self, T: float, mu: float, gamma: float, W: float, + omega0: float = None, *, Nk: int = 10, tag: Any = None + ): + super().__init__(T, mu, tag) + + self.gamma = gamma + self.W = W + self.Nk = Nk + if omega0 is None: + self.omega0 = mu + else: + self.omega0 = omega0 + + def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike): + """ + Calculates the Lorentzian spectral density. + + Parameters + ---------- + w : array_like or float + Energy of the mode. + """ + + w = np.asarray(w, dtype=float) + return self.gamma * self.W**2 / ((w - self.omega0)**2 + self.W**2) + + def correlation_function_plus( + self, t: float | ArrayLike, Nk: int = None + ) -> (float | ArrayLike): + r""" + Calculates the "+"-branch of the two-time auto-correlation function of + the Lorentzian environment. The calculation is performed by summing a + large number of exponents of the Pade expansion. + + Parameters + ---------- + t : array_like or float + The time at which to evaluate the correlation function. + Nk : int, optional + The number of exponents to use. If not provided, then the value + that was provided when the class was instantiated is used. + """ + Nk = Nk or self.Nk + return self._correlation_function(t, Nk, 1) + + def correlation_function_minus( + self, t: float | ArrayLike, Nk: int = None + ) -> (float | ArrayLike): + r""" + Calculates the "-"-branch of the two-time auto-correlation function of + the Lorentzian environment. The calculation is performed by summing a + large number of exponents of the Pade expansion. + + Parameters + ---------- + t : array_like or float + The time at which to evaluate the correlation function. + Nk : int, optional + The number of exponents to use. If not provided, then the value + that was provided when the class was instantiated is used. + """ + Nk = Nk or self.Nk + return self._correlation_function(t, Nk, -1) + + def _correlation_function(self, t, Nk, sigma): + if self.T == 0: + raise NotImplementedError( + "Calculation of zero-temperature Lorentzian correlation " + "functions is not implemented yet.") + + t = np.asarray(t, dtype=float) + abs_t = np.abs(t) + c, v = self._corr(Nk, sigma) + + result = np.sum([ck * np.exp(-np.asarray(vk * abs_t)) + for ck, vk in zip(c, v)], axis=0) + + result = np.asarray(result, dtype=complex) + result[t < 0] = np.conj(result[t < 0]) + return result.item() if t.ndim == 0 else result + + def power_spectrum_plus(self, w: float | ArrayLike) -> (float | ArrayLike): + r""" + Calculates the "+"-branch of the power spectrum of the Lorentzian + environment. + + Parameters + ---------- + w : array_like or float + The frequency at which to evaluate the power spectrum. + """ + + return self.spectral_density(w) / (np.exp((w - self.mu) / self.T) + 1) + + def power_spectrum_minus( + self, w: float | ArrayLike + ) -> (float | ArrayLike): + r""" + Calculates the "-"-branch of the power spectrum of the Lorentzian + environment. + + Parameters + ---------- + w : array_like or float + The frequency at which to evaluate the power spectrum. + """ + + return self.spectral_density(w) / (np.exp((self.mu - w) / self.T) + 1) + + def approx_by_matsubara( + self, Nk: int, tag: Any = None + ) -> ExponentialFermionicEnvironment: + """ + Generates an approximation to this environment by truncating its + Matsubara expansion. + + Parameters + ---------- + Nk : int + Number of Matsubara terms to include. In total, the "+" and "-" + correlation function branches will include `Nk+1` terms each. + + tag : optional, str, tuple or any other object + An identifier (name) for the approximated environment. If not + provided, a tag will be generated from the tag of this environment. + + Returns + ------- + The approximated environment with multi-exponential correlation + function. + """ + if self.T == 0: + raise NotImplementedError( + "Calculation of zero-temperature Lorentzian correlation " + "functions is not implemented yet.") + + if tag is None and self.tag is not None: + tag = (self.tag, "Matsubara Truncation") + + ck_plus, vk_plus = self._matsubara_params(Nk, 1) + ck_minus, vk_minus = self._matsubara_params(Nk, -1) + + return ExponentialFermionicEnvironment( + ck_plus, vk_plus, ck_minus, vk_minus, T=self.T, mu=self.mu, + tag=tag + ) + + def approx_by_pade( + self, Nk: int, tag: Any = None + ) -> ExponentialFermionicEnvironment: + """ + Generates an approximation to this environment by truncating its + Pade expansion. + + Parameters + ---------- + Nk : int + Number of Pade terms to include. In total, the "+" and "-" + correlation function branches will include `Nk+1` terms each. + + tag : optional, str, tuple or any other object + An identifier (name) for the approximated environment. If not + provided, a tag will be generated from the tag of this environment. + + Returns + ------- + The approximated environment with multi-exponential correlation + function. + """ + if self.T == 0: + raise NotImplementedError( + "Calculation of zero-temperature Lorentzian correlation " + "functions is not implemented yet.") + + if tag is None and self.tag is not None: + tag = (self.tag, "Pade Truncation") + + ck_plus, vk_plus = self._corr(Nk, sigma=1) + ck_minus, vk_minus = self._corr(Nk, sigma=-1) + + return ExponentialFermionicEnvironment( + ck_plus, vk_plus, ck_minus, vk_minus, T=self.T, mu=self.mu, tag=tag + ) + + def _matsubara_params(self, Nk, sigma): + """ Calculate the Matsubara coefficients and frequencies. """ + def f(x): + return 1 / (np.exp(x / self.T) + 1) + + coeff_list = [( + self.W * self.gamma / 2 * + f(sigma * (self.omega0 - self.mu) + 1j * self.W) + )] + exp_list = [self.W - sigma * 1j * self.omega0] + + xk_list = [(2 * k - 1) * np.pi * self.T for k in range(1, Nk + 1)] + for xk in xk_list: + coeff_list.append( + 1j * self.gamma * self.W**2 * self.T / + ((sigma * xk - 1j * self.mu + 1j * self.omega0)**2 - self.W**2) + ) + exp_list.append( + xk - sigma * 1j * self.mu + ) + + return coeff_list, exp_list + + # --- Pade approx calculation --- + + def _corr(self, Nk, sigma): + beta = 1 / self.T + kappa, epsilon = self._kappa_epsilon(Nk) + + def f_approx(x): + f = 0.5 + for ll in range(1, Nk + 1): + f = f - 2 * kappa[ll] * x / (x**2 + epsilon[ll]**2) + return f + + eta_list = [(0.5 * self.gamma * self.W * + f_approx(beta * sigma * (self.omega0 - self.mu) + + beta * 1j * self.W))] + gamma_list = [self.W - sigma * 1.0j * self.omega0] + + for ll in range(1, Nk + 1): + eta_list.append( + -1.0j * (kappa[ll] / beta) * self.gamma * self.W**2 + / ((self.mu - self.omega0 + sigma * 1j * epsilon[ll] / beta)**2 + + self.W**2) + ) + gamma_list.append(epsilon[ll] / beta - sigma * 1.0j * self.mu) + + return eta_list, gamma_list + + def _kappa_epsilon(self, Nk): + eps = self._calc_eps(Nk) + chi = self._calc_chi(Nk) + + kappa = [0] + prefactor = 0.5 * Nk * (2 * (Nk + 1) - 1) + for j in range(Nk): + term = prefactor + for k in range(Nk - 1): + term *= ( + (chi[k]**2 - eps[j]**2) / + (eps[k]**2 - eps[j]**2 + self._delta(j, k)) + ) + for k in [Nk - 1]: + term /= (eps[k]**2 - eps[j]**2 + self._delta(j, k)) + kappa.append(term) + + epsilon = [0] + eps + + return kappa, epsilon + + def _delta(self, i, j): + return 1.0 if i == j else 0.0 + + def _calc_eps(self, Nk): + alpha = np.diag([ + 1. / np.sqrt((2 * k + 3) * (2 * k + 1)) + for k in range(2 * Nk - 1) + ], k=1) + alpha += alpha.transpose() + + evals = eigvalsh(alpha) + eps = [-2. / val for val in evals[0: Nk]] + return eps + + def _calc_chi(self, Nk): + alpha_p = np.diag([ + 1. / np.sqrt((2 * k + 5) * (2 * k + 3)) + for k in range(2 * Nk - 2) + ], k=1) + alpha_p += alpha_p.transpose() + evals = eigvalsh(alpha_p) + chi = [-2. / val for val in evals[0: Nk - 1]] + return chi + + +class ExponentialFermionicEnvironment(FermionicEnvironment): + """ + Fermionic environment that is specified through an exponential + decomposition of its correlation function. The list of coefficients and + exponents in the decomposition may either be passed through the four lists + `ck_plus`, `vk_plus`, `ck_minus`, `vk_minus`, or as a list of fermionic + :class:`CFExponent` objects. + + Alternative constructors :meth:`from_plus_exponents` and + :meth:`from_minus_exponents` are available to compute the "-" exponents + automatically from the "+" ones, or vice versa. + + Parameters + ---------- + ck_plus : list of complex + The coefficients of the expansion terms for the ``+`` part of the + correlation function. The corresponding frequencies are passed as + vk_plus. + + vk_plus : list of complex + The frequencies (exponents) of the expansion terms for the ``+`` part + of the correlation function. The corresponding coefficients are passed + as ck_plus. + + ck_minus : list of complex + The coefficients of the expansion terms for the ``-`` part of the + correlation function. The corresponding frequencies are passed as + vk_minus. + + vk_minus : list of complex + The frequencies (exponents) of the expansion terms for the ``-`` part + of the correlation function. The corresponding coefficients are passed + as ck_minus. + + exponents : list of :class:`CFExponent` + The expansion coefficients and exponents of both parts of the + correlation function as :class:`CFExponent` objects. + + T: optional, float + The temperature of the environment. + + mu: optional, float + The chemical potential of the environment. + + tag : optional, str, tuple or any other object + An identifier (name) for this environment. + """ + + def _check_cks_and_vks(self, ck_plus, vk_plus, ck_minus, vk_minus): + # all None: returns False + # all provided and lengths match: returns True + # otherwise: raises ValueError + lists = [ck_plus, vk_plus, ck_minus, vk_minus] + if all(x is None for x in lists): + return False + if any(x is None for x in lists): + raise ValueError( + "If any of the exponent lists ck_plus, vk_plus, ck_minus, " + "vk_minus is provided, all must be provided." + ) + if len(ck_plus) != len(vk_plus) or len(ck_minus) != len(vk_minus): + raise ValueError( + "The exponent lists ck_plus and vk_plus, and ck_minus and " + "vk_minus must be the same length." + ) + return True + + def __init__( + self, + ck_plus: ArrayLike = None, vk_plus: ArrayLike = None, + ck_minus: ArrayLike = None, vk_minus: ArrayLike = None, + *, + exponents: Sequence[CFExponent] = None, + T: float = None, mu: float = None, tag: Any = None + ): + super().__init__(T, mu, tag) + + lists_provided = self._check_cks_and_vks( + ck_plus, vk_plus, ck_minus, vk_minus) + if exponents is None and not lists_provided: + raise ValueError( + "Either the parameter `exponents` or the parameters " + "`ck_plus`, `vk_plus`, `ck_minus`, `vk_minus` must be " + "provided." + ) + if (exponents is not None and + not all(exp.fermionic for exp in exponents)): + raise ValueError( + "Bosonic exponent passed to exponential fermionic environment." + ) + + self.exponents = exponents or [] + if lists_provided: + self.exponents.extend(CFExponent("+", ck, vk, tag=tag) + for ck, vk in zip(ck_plus, vk_plus)) + self.exponents.extend(CFExponent("-", ck, vk, tag=tag) + for ck, vk in zip(ck_minus, vk_minus)) + + def spectral_density(self, w: float | ArrayLike) -> (float | ArrayLike): + """ + Computes the spectral density corresponding to the multi-exponential + correlation function. + + Parameters + ---------- + w : array_like or float + Energy of the mode. + """ + + return self.power_spectrum_minus(w) + self.power_spectrum_plus(w) + + def correlation_function_plus( + self, t: float | ArrayLike + ) -> (float | ArrayLike): + r""" + Computes the "+"-branch of the correlation function represented by this + exponential decomposition. + + Parameters + ---------- + t : array_like or float + The times at which to evaluate the correlation function. + """ + + return self._cf(t, CFExponent.types['+']) + + def correlation_function_minus( + self, t: float | ArrayLike + ) -> (float | ArrayLike): + r""" + Computes the "-"-branch of the correlation function represented by this + exponential decomposition. + + Parameters + ---------- + t : array_like or float + The times at which to evaluate the correlation function. + """ + + return self._cf(t, CFExponent.types['-']) + + def _cf(self, t, type): + t = np.asarray(t, dtype=float) + corr = np.zeros_like(t, dtype=complex) + + for exp in self.exponents: + if exp.type == type: + corr += exp.coefficient * np.exp(-exp.exponent * np.abs(t)) + corr[t < 0] = np.conj(corr[t < 0]) + + return corr.item() if t.ndim == 0 else corr + + def power_spectrum_plus( + self, w: float | ArrayLike + ) -> (float | ArrayLike): + r""" + Calculates the "+"-branch of the power spectrum corresponding to the + multi-exponential correlation function. + + Parameters + ---------- + w : array_like or float + The frequency at which to evaluate the power spectrum. + """ + + return self._ps(w, CFExponent.types['+'], 1) + + def power_spectrum_minus( + self, w: float | ArrayLike + ) -> (float | ArrayLike): + r""" + Calculates the "-"-branch of the power spectrum corresponding to the + multi-exponential correlation function. + + Parameters + ---------- + w : array_like or float + The frequency at which to evaluate the power spectrum. + """ + + return self._ps(w, CFExponent.types['-'], -1) + + def _ps(self, w, type, sigma): + w = np.asarray(w, dtype=float) + S = np.zeros_like(w) + + for exp in self.exponents: + if exp.type == type: + S += 2 * np.real( + exp.coefficient / (exp.exponent + sigma * 1j * w) + ) + + return S.item() if w.ndim == 0 else S diff --git a/qutip/solver/brmesolve.py b/qutip/solver/brmesolve.py index e6a11c52f0..68b5c608b6 100644 --- a/qutip/solver/brmesolve.py +++ b/qutip/solver/brmesolve.py @@ -21,6 +21,7 @@ from .options import _SolverOptions from ._feedback import _QobjFeedback, _DataFeedback from ..typing import EopsLike, QobjEvoLike, CoefficientLike +from ..core.environment import BosonicEnvironment, FermionicEnvironment def brmesolve( @@ -198,17 +199,22 @@ def brmesolve( a_op = QobjEvo(a_op, args=args, tlist=tlist) if isinstance(spectra, str): new_a_ops.append( - (a_op, coefficient(spectra, args={**args, 'w':0}))) + (a_op, coefficient(spectra, args={**args, 'w': 0}))) elif isinstance(spectra, InterCoefficient): new_a_ops.append((a_op, SpectraCoefficient(spectra))) elif isinstance(spectra, Coefficient): new_a_ops.append((a_op, spectra)) + elif isinstance(spectra, BosonicEnvironment): + spec = SpectraCoefficient( + coefficient(spectra.power_spectrum) + ) + new_a_ops.append((a_op, spec)) elif callable(spectra): sig = inspect.signature(spectra) if tuple(sig.parameters.keys()) == ("w",): spec = SpectraCoefficient(coefficient(spectra)) else: - spec = coefficient(spectra, args={**args, 'w':0}) + spec = coefficient(spectra, args={**args, 'w': 0}) new_a_ops.append((a_op, spec)) else: raise TypeError("a_ops's spectra not known") @@ -274,7 +280,7 @@ class BRSolver(Solver): name = "brmesolve" solver_options = { "progress_bar": "", - "progress_kwargs": {"chunk_size":10}, + "progress_kwargs": {"chunk_size": 10}, "store_final_state": False, "store_states": None, "normalize_output": False, @@ -331,7 +337,6 @@ def __init__( self.stats = self._initialize_stats() self.rhs._register_feedback({}, solver=self.name) - def _initialize_stats(self): stats = super()._initialize_stats() stats.update({ diff --git a/qutip/solver/heom/__init__.py b/qutip/solver/heom/__init__.py index 9a8264b073..28db052659 100644 --- a/qutip/solver/heom/__init__.py +++ b/qutip/solver/heom/__init__.py @@ -1,5 +1,5 @@ """ -This module provides solvers for system-bath evoluation using the +This module provides solvers for system-bath evolution using the HEOM (hierarchy equations of motion). See https://en.wikipedia.org/wiki/Hierarchical_equations_of_motion for a very diff --git a/qutip/solver/heom/bofin_baths.py b/qutip/solver/heom/bofin_baths.py index ff64dfc610..9cdc667c86 100644 --- a/qutip/solver/heom/bofin_baths.py +++ b/qutip/solver/heom/bofin_baths.py @@ -1,22 +1,28 @@ """ -This module provides utilities for describing baths when using the -HEOM (hierarchy equations of motion) to model system-bath interactions. - -See the ``qutip.nonmarkov.bofin_solvers`` module for the associated solver. +This module provides utilities for describing baths when using the HEOM +(hierarchy equations of motion) to model system-bath interactions. See the +``qutip.nonmarkov.bofin_solvers`` module for the associated solver. The implementation is derived from the BoFiN library (see https://github.com/tehruhn/bofin) which was itself derived from an earlier implementation in QuTiP itself. -""" -import enum +The "bath" classes in this module are closely related to the "environment" +classes in the `qutip.core.environment` module. The bath classes were +implemented first, specifically for the HEOM solver. The environment classes, +added later provide additional functionality and are designed to work also with +other QuTiP solvers. -import numpy as np -from scipy.linalg import eigvalsh +The bath classes are kept partly for backwards compatibility, partly to use as +a "shortcut" when one only wants to use the HEOM solver. + +Note that this module also contains the `BathExponent` class, which is used by +the HEOM solver internally and in the result object describing its output. +""" from qutip.core import data as _data +from qutip.core import environment from qutip.core.qobj import Qobj -from qutip.core.superoperator import spre, spost __all__ = [ "BathExponent", @@ -36,14 +42,21 @@ def _isequal(Q1, Q2, tol): return _data.iszero(_data.sub(Q1.data, Q2.data), tol=tol) -class BathExponent: +class BathExponent(environment.CFExponent): """ Represents a single exponent (naively, an excitation mode) within the decomposition of the correlation functions of a bath. + This class extends the + :class:`CFExponent ` of the + environment-module for use with the HEOM solver. In addition to the + exponent itself, the `BathExponent` keeps track of the corresponding system + coupling operator ``Q``, as well as the parameters ``dim`` and + ``sigma_bar_k_offset``. + Parameters ---------- - type : {"R", "I", "RI", "+", "-"} or BathExponent.ExponentType + type : {"R", "I", "RI", "+", "-"} or ``CFExponent.ExponentType`` The type of bath exponent. "R" and "I" are bosonic bath exponents that appear in the real and @@ -97,52 +110,32 @@ class BathExponent: All of the parameters are also available as attributes. """ - types = enum.Enum("ExponentType", ["R", "I", "RI", "+", "-"]) - - def _check_ck2(self, type, ck2): - if type == self.types["RI"]: - if ck2 is None: - raise ValueError("RI bath exponents require ck2") - else: - if ck2 is not None: - raise ValueError( - "Second co-efficient (ck2) should only be specified for RI" - " bath exponents" - ) + types = environment.CFExponent.types def _check_sigma_bar_k_offset(self, type, offset): if type in (self.types["+"], self.types["-"]): if offset is None: raise ValueError( - "+ and - bath exponents require sigma_bar_k_offset" + "+ and - type exponents require sigma_bar_k_offset" ) else: if offset is not None: raise ValueError( "Offset of sigma bar (sigma_bar_k_offset) should only be" - " specified for + and - bath exponents" + " specified for + and - type exponents" ) - def _type_is_fermionic(self, type): - return type in (self.types["+"], self.types["-"]) - def __init__( self, type, dim, Q, ck, vk, ck2=None, sigma_bar_k_offset=None, tag=None, ): - if not isinstance(type, self.types): - type = self.types[type] - self._check_ck2(type, ck2) - self._check_sigma_bar_k_offset(type, sigma_bar_k_offset) - self.type = type - self.fermionic = self._type_is_fermionic(type) + super().__init__(type, ck, vk, ck2, tag) + self.dim = dim self.Q = Q - self.ck = ck - self.vk = vk - self.ck2 = ck2 + + self._check_sigma_bar_k_offset(self.type, sigma_bar_k_offset) self.sigma_bar_k_offset = sigma_bar_k_offset - self.tag = tag def __repr__(self): dims = getattr(self.Q, "dims", None) @@ -156,6 +149,17 @@ def __repr__(self): f" tag={self.tag!r}>" ) + def _can_combine(self, other, rtol, atol): + if not super()._can_combine(other, rtol, atol): + return False + if not _isequal(self.Q, other.Q, tol=atol): + return False + return True + + def _combine(self, other): + # Assumes can combine was checked + return super()._combine(other, dim=None, Q=self.Q) + class Bath: """ @@ -163,19 +167,15 @@ class Bath: Parameters ---------- - exponents : list of BathExponent + exponents : list of :class:`.BathExponent` The exponents of the correlation function describing the bath. - - Attributes - ---------- - - All of the parameters are available as attributes. """ + def __init__(self, exponents): self.exponents = exponents -class BosonicBath(Bath): +class BosonicBath(environment.ExponentialBosonicEnvironment): """ A helper class for constructing a bosonic bath from the expansion coefficients and frequencies for the real and imaginary parts of @@ -223,119 +223,67 @@ class BosonicBath(Bath): combine : bool, default True Whether to combine exponents with the same frequency (and coupling - operator). See :meth:`combine` for details. + operator). See :meth:`.ExponentialBosonicEnvironment.combine` for + details. tag : optional, str, tuple or any other object A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from. + + Notes + ----- + This class is part of the "bath" API, which is now mirrored by the newer + "environment" API. The bath classes are kept in QuTiP for reasons of + backwards compatibility and convenience. This class is an extended version + of the :class:`.ExponentialBosonicEnvironment`, adding the parameter ``Q``, + which is not included in the newer "environment" API. """ - def _check_cks_and_vks(self, ck_real, vk_real, ck_imag, vk_imag): - if len(ck_real) != len(vk_real) or len(ck_imag) != len(vk_imag): - raise ValueError( - "The bath exponent lists ck_real and vk_real, and ck_imag and" - " vk_imag must be the same length." - ) + + def _make_exponent(self, type, ck, vk, ck2=None, tag=None): + return BathExponent(type, None, self._Q, ck, vk, ck2, tag=tag) def _check_coup_op(self, Q): if not isinstance(Q, Qobj): raise ValueError("The coupling operator Q must be a Qobj.") def __init__( - self, Q, ck_real, vk_real, ck_imag, vk_imag, combine=True, - tag=None, + self, Q, ck_real, vk_real, ck_imag, vk_imag, combine=True, tag=None ): - self._check_cks_and_vks(ck_real, vk_real, ck_imag, vk_imag) self._check_coup_op(Q) + self._Q = Q - exponents = [] - exponents.extend( - BathExponent("R", None, Q, ck, vk, tag=tag) - for ck, vk in zip(ck_real, vk_real) - ) - exponents.extend( - BathExponent("I", None, Q, ck, vk, tag=tag) - for ck, vk in zip(ck_imag, vk_imag) - ) - - if combine: - exponents = self.combine(exponents) - - super().__init__(exponents) + super().__init__(ck_real, vk_real, ck_imag, vk_imag, + combine=combine, tag=tag) @classmethod - def combine(cls, exponents, rtol=1e-5, atol=1e-7): + def from_environment(cls, env, Q, dim=None): """ - Group bosonic exponents with the same frequency and return a - single exponent for each frequency present. - - Exponents with the same frequency are only combined if they share the - same coupling operator ``.Q``. - - Note that combined exponents take their tag from the first - exponent in the group being combined (i.e. the one that occurs first - in the given exponents list). + Converts from the "environment" API to the "bath" API. A `BosonicBath` + combines the information from an `ExponentialBosonicEnvironment` and a + coupling operator. Parameters ---------- - exponents : list of BathExponent - The list of exponents to combine. - - rtol : float, default 1e-5 - The relative tolerance to use to when comparing frequencies and - coupling operators. - - atol : float, default 1e-7 - The absolute tolerance to use to when comparing frequencies and - coupling operators. - - Returns - ------- - list of BathExponent - The new reduced list of exponents. + env : :class:`.ExponentialBosonicEnvironment` + The bath. + Q : Qobj + The coupling operator for the bath. + dim : optional, int or ``None`` (default ``None``) + The maximum number of excitations for each exponent. Usually + ``None`` (i.e. unlimited). """ - groups = [] - remaining = exponents[:] - - while remaining: - e1 = remaining.pop(0) - group = [e1] - for e2 in remaining[:]: - if ( - np.isclose(e1.vk, e2.vk, rtol=rtol, atol=atol) and - _isequal(e1.Q, e2.Q, tol=atol) - ): - group.append(e2) - remaining.remove(e2) - groups.append(group) - - new_exponents = [] - for combine in groups: - exp1 = combine[0] - if (exp1.type != exp1.types.RI) and all( - exp2.type == exp1.type for exp2 in combine - ): - # the group is either type I or R - ck = sum(exp.ck for exp in combine) - new_exponents.append(BathExponent( - exp1.type, None, exp1.Q, ck, exp1.vk, tag=exp1.tag, - )) - else: - # the group includes both type I and R exponents - ck_R = ( - sum(exp.ck for exp in combine if exp.type == exp.types.R) + - sum(exp.ck for exp in combine if exp.type == exp.types.RI) - ) - ck_I = ( - sum(exp.ck for exp in combine if exp.type == exp.types.I) + - sum(exp.ck2 for exp in combine if exp.type == exp.types.RI) - ) - new_exponents.append(BathExponent( - "RI", None, exp1.Q, ck_R, exp1.vk, ck2=ck_I, - tag=exp1.tag, - )) + bath_exponents = [] + for exponent in env.exponents: + new_exponent = BathExponent( + exponent.type, dim, Q, exponent.ck, exponent.vk, + exponent.ck2, None, env.tag + ) + bath_exponents.append(new_exponent) - return new_exponents + result = cls(Q, [], [], [], [], tag=env.tag) + result.exponents = bath_exponents + return result class DrudeLorentzBath(BosonicBath): @@ -363,30 +311,42 @@ class DrudeLorentzBath(BosonicBath): combine : bool, default True Whether to combine exponents with the same frequency (and coupling - operator). See :meth:`BosonicBath.combine` for details. + operator). See :meth:`.ExponentialBosonicEnvironment.combine` for + details. tag : optional, str, tuple or any other object A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from. + + Notes + ----- + This class is part of the "bath" API, which is now mirrored by the newer + "environment" API. The bath classes are kept in QuTiP for reasons of + backwards compatibility and convenience. Creating a `DrudeLorentzBath` is + equivalent to creating a :class:`.DrudeLorentzEnvironment`, performing a + :meth:`Matsubara <.DrudeLorentzEnvironment.approx_by_matsubara>` + approximation, and finally bundling the result together with the coupling + operator ``Q`` for convenient use with the HEOM solver. """ - def __init__( - self, Q, lam, gamma, T, Nk, combine=True, tag=None, - ): - ck_real, vk_real, ck_imag, vk_imag = self._matsubara_params( - lam=lam, - gamma=gamma, - T=T, - Nk=Nk, - ) - super().__init__( - Q, ck_real, vk_real, ck_imag, vk_imag, combine=combine, tag=tag, + def __new__( + mcs, Q, lam, gamma, T, Nk, combine=True, tag=None, + ): + # Basically this makes `DrudeLorentzBath` a function + # (Q, lam, ...) -> BosonicBath + # but it is made to look like a class because it was a class in the + # initial bofin release + env = environment.DrudeLorentzEnvironment(T, lam, gamma) + matsubara_approx, delta = env.approx_by_matsubara( + Nk=Nk, combine=combine, compute_delta=True, tag=tag ) - self._dl_terminator = _DrudeLorentzTerminator( - Q=Q, lam=lam, gamma=gamma, T=T, + result = BosonicBath.from_environment(matsubara_approx, Q) + result.terminator = lambda: ( + delta, environment.system_terminator(Q, delta) ) + return result def terminator(self): """ @@ -404,29 +364,13 @@ def terminator(self): terminator : Qobj - The Matsubara terminator -- i.e. a liouvillian term representing + The Padé terminator -- i.e. a liouvillian term representing the contribution to the system-bath dynamics of all exponential expansion terms beyond ``Nk``. It should be used by adding it to the system liouvillian (i.e. ``liouvillian(H_sys)``). """ - delta, L = self._dl_terminator.terminator(self.exponents) - return delta, L - - def _matsubara_params(self, lam, gamma, T, Nk): - """ Calculate the Matsubara coefficents and frequencies. """ - ck_real = [lam * gamma / np.tan(gamma / (2 * T))] - ck_real.extend([ - (8 * lam * gamma * T * np.pi * k * T / - ((2 * np.pi * k * T)**2 - gamma**2)) - for k in range(1, Nk + 1) - ]) - vk_real = [gamma] - vk_real.extend([2 * np.pi * k * T for k in range(1, Nk + 1)]) - - ck_imag = [lam * gamma * (-1.0)] - vk_imag = [gamma] - - return ck_real, vk_real, ck_imag, vk_imag + # This is only here to keep the API doc + ... class DrudeLorentzPadeBath(BosonicBath): @@ -434,8 +378,8 @@ class DrudeLorentzPadeBath(BosonicBath): A helper class for constructing a Padé expansion for a Drude-Lorentz bosonic bath from the bath parameters (see parameters below). - A Padé approximant is a sum-over-poles expansion ( - see https://en.wikipedia.org/wiki/Pad%C3%A9_approximant). + A Padé approximant is a sum-over-poles expansion (see + https://en.wikipedia.org/wiki/Pad%C3%A9_approximant). The application of the Padé method to spectrum decompoisitions is described in "Padé spectrum decompositions of quantum distribution functions and @@ -469,32 +413,39 @@ class DrudeLorentzPadeBath(BosonicBath): combine : bool, default True Whether to combine exponents with the same frequency (and coupling - operator). See :meth:`BosonicBath.combine` for details. + operator). See :meth:`.ExponentialBosonicEnvironment.combine` for + details. tag : optional, str, tuple or any other object A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from. - """ - def __init__( - self, Q, lam, gamma, T, Nk, combine=True, tag=None - ): - eta_p, gamma_p = self._corr(lam=lam, gamma=gamma, T=T, Nk=Nk) - ck_real = [np.real(eta) for eta in eta_p] - vk_real = [gam for gam in gamma_p] - # There is only one term in the expansion of the imaginary part of the - # Drude-Lorentz correlation function. - ck_imag = [np.imag(eta_p[0])] - vk_imag = [gamma_p[0]] + Notes + ----- + This class is part of the "bath" API, which is now mirrored by the newer + "environment" API. The bath classes are kept in QuTiP for reasons of + backwards compatibility and convenience. Creating a `DrudeLorentzPadeBath` + is equivalent to creating a :class:`.DrudeLorentzEnvironment`, performing a + :meth:`Pade <.DrudeLorentzEnvironment.approx_by_pade>` approximation, and + finally bundling the result together with the coupling operator ``Q`` for + convenient use with the HEOM solver. + """ - super().__init__( - Q, ck_real, vk_real, ck_imag, vk_imag, combine=combine, tag=tag, + def __new__( + mcs, Q, lam, gamma, T, Nk, combine=True, tag=None, + ): + # See DrudeLorentzBath comment + env = environment.DrudeLorentzEnvironment(T, lam, gamma) + pade_approx, delta = env.approx_by_pade( + Nk=Nk, combine=combine, compute_delta=True, tag=tag ) - self._dl_terminator = _DrudeLorentzTerminator( - Q=Q, lam=lam, gamma=gamma, T=T, + result = BosonicBath.from_environment(pade_approx, Q) + result.terminator = lambda: ( + delta, environment.system_terminator(Q, delta) ) + return result def terminator(self): """ @@ -517,104 +468,8 @@ def terminator(self): expansion terms beyond ``Nk``. It should be used by adding it to the system liouvillian (i.e. ``liouvillian(H_sys)``). """ - delta, L = self._dl_terminator.terminator(self.exponents) - return delta, L - - def _corr(self, lam, gamma, T, Nk): - beta = 1. / T - kappa, epsilon = self._kappa_epsilon(Nk) - - eta_p = [lam * gamma * (self._cot(gamma * beta / 2.0) - 1.0j)] - gamma_p = [gamma] - - for ll in range(1, Nk + 1): - eta_p.append( - (kappa[ll] / beta) * 4 * lam * gamma * (epsilon[ll] / beta) - / ((epsilon[ll]**2 / beta**2) - gamma**2) - ) - gamma_p.append(epsilon[ll] / beta) - - return eta_p, gamma_p - - def _cot(self, x): - return 1. / np.tan(x) - - def _kappa_epsilon(self, Nk): - eps = self._calc_eps(Nk) - chi = self._calc_chi(Nk) - - kappa = [0] - prefactor = 0.5 * Nk * (2 * (Nk + 1) + 1) - for j in range(Nk): - term = prefactor - for k in range(Nk - 1): - term *= ( - (chi[k]**2 - eps[j]**2) / - (eps[k]**2 - eps[j]**2 + self._delta(j, k)) - ) - for k in [Nk - 1]: - term /= (eps[k]**2 - eps[j]**2 + self._delta(j, k)) - kappa.append(term) - - epsilon = [0] + eps - - return kappa, epsilon - - def _delta(self, i, j): - return 1.0 if i == j else 0.0 - - def _calc_eps(self, Nk): - alpha = np.diag([ - 1. / np.sqrt((2 * k + 5) * (2 * k + 3)) - for k in range(2 * Nk - 1) - ], k=1) - alpha += alpha.transpose() - evals = eigvalsh(alpha) - eps = [-2. / val for val in evals[0: Nk]] - return eps - - def _calc_chi(self, Nk): - alpha_p = np.diag([ - 1. / np.sqrt((2 * k + 7) * (2 * k + 5)) - for k in range(2 * Nk - 2) - ], k=1) - alpha_p += alpha_p.transpose() - evals = eigvalsh(alpha_p) - chi = [-2. / val for val in evals[0: Nk - 1]] - return chi - - -class _DrudeLorentzTerminator: - """ A class for calculating the terminator of a Drude-Lorentz bath - expansion. - """ - def __init__(self, Q, lam, gamma, T): - self.Q = Q - self.lam = lam - self.gamma = gamma - self.T = T - - def terminator(self, exponents): - """ Calculate the terminator for a Drude-Lorentz bath. """ - Q = self.Q - lam = self.lam - gamma = self.gamma - beta = 1 / self.T - - delta = 2 * lam / (beta * gamma) - 1j * lam - - for exp in exponents: - if exp.type == BathExponent.types["R"]: - delta -= exp.ck / exp.vk - elif exp.type == BathExponent.types["RI"]: - delta -= (exp.ck + 1j * exp.ck2) / exp.vk - else: - delta -= 1j * exp.ck / exp.vk - - op = -2*spre(Q)*spost(Q.dag()) + spre(Q.dag()*Q) + spost(Q.dag()*Q) - L_bnd = -delta * op - - return delta, L_bnd + # This is only here to keep the API doc + ... class UnderDampedBath(BosonicBath): @@ -645,67 +500,37 @@ class UnderDampedBath(BosonicBath): combine : bool, default True Whether to combine exponents with the same frequency (and coupling - operator). See :meth:`BosonicBath.combine` for details. + operator). See :meth:`.ExponentialBosonicEnvironment.combine` for + details. tag : optional, str, tuple or any other object A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from. + + Notes + ----- + This class is part of the "bath" API, which is now mirrored by the newer + "environment" API. The bath classes are kept in QuTiP for reasons of + backwards compatibility and convenience. Creating an `UnderDampedBath` is + equivalent to creating an :class:`.UnderDampedEnvironment`, performing a + :meth:`Matsubara <.UnderDampedEnvironment.approx_by_matsubara>` + approximation, and finally bundling the result together with the coupling + operator ``Q`` for convenient use with the HEOM solver. """ - def __init__( - self, Q, lam, gamma, w0, T, Nk, combine=True, tag=None, - ): - ck_real, vk_real, ck_imag, vk_imag = self._matsubara_params( - lam=lam, - gamma=gamma, - w0=w0, - T=T, - Nk=Nk, - ) - super().__init__( - Q, ck_real, vk_real, ck_imag, vk_imag, combine=combine, tag=tag, + def __new__( + mcs, Q, lam, gamma, w0, T, Nk, combine=True, tag=None, + ): + # See DrudeLorentzBath comment + env = environment.UnderDampedEnvironment(T, lam, gamma, w0) + matsubara_approx = env.approx_by_matsubara( + Nk=Nk, combine=combine, tag=tag ) + return BosonicBath.from_environment(matsubara_approx, Q) - def _matsubara_params(self, lam, gamma, w0, T, Nk): - """ Calculate the Matsubara coefficents and frequencies. """ - beta = 1/T - Om = np.sqrt(w0**2 - (gamma/2)**2) - Gamma = gamma/2. - - ck_real = ([ - (lam**2 / (4 * Om)) - * (1 / np.tanh(beta * (Om + 1.0j * Gamma) / 2)), - (lam**2 / (4*Om)) - * (1 / np.tanh(beta * (Om - 1.0j * Gamma) / 2)), - ]) - - ck_real.extend([ - (-2 * lam**2 * gamma / beta) * (2 * np.pi * k / beta) - / ( - ((Om + 1.0j * Gamma)**2 + (2 * np.pi * k/beta)**2) - * ((Om - 1.0j * Gamma)**2 + (2 * np.pi * k / beta)**2) - ) - for k in range(1, Nk + 1) - ]) - - vk_real = [-1.0j * Om + Gamma, 1.0j * Om + Gamma] - vk_real.extend([ - 2 * np.pi * k * T - for k in range(1, Nk + 1) - ]) - - ck_imag = [ - 1.0j * lam**2 / (4 * Om), - -1.0j * lam**2 / (4 * Om), - ] - vk_imag = [-1.0j * Om + Gamma, 1.0j * Om + Gamma] - - return ck_real, vk_real, ck_imag, vk_imag - - -class FermionicBath(Bath): +class FermionicBath(environment.ExponentialFermionicEnvironment): """ A helper class for constructing a fermionic bath from the expansion coefficients and frequencies for the ``+`` and ``-`` modes of @@ -755,21 +580,27 @@ class FermionicBath(Bath): A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from. + + Notes + ----- + This class is part of the "bath" API, which is now mirrored by the newer + "environment" API. The bath classes are kept in QuTiP for reasons of + backwards compatibility and convenience. This class is an extended version + of the :class:`.ExponentialFermionicEnvironment`, adding the parameter + ``Q``, which is not included in the newer "environment" API. """ def _check_cks_and_vks(self, ck_plus, vk_plus, ck_minus, vk_minus): - if len(ck_plus) != len(vk_plus) or len(ck_minus) != len(vk_minus): - raise ValueError( - "The bath exponent lists ck_plus and vk_plus, and ck_minus and" - " vk_minus must be the same length." - ) - if len(ck_plus) != len(ck_minus): + lists_provided = super()._check_cks_and_vks( + ck_plus, vk_plus, ck_minus, vk_minus) + if lists_provided and len(ck_plus) != len(ck_minus): raise ValueError( "The must be the same number of plus and minus exponents" " in the bath, and elements of plus and minus arrays" " should be arranged so that ck_plus[i] is the plus mode" " corresponding to ck_minus[i]." ) + return lists_provided def _check_coup_op(self, Q): if not isinstance(Q, Qobj): @@ -787,7 +618,55 @@ def __init__(self, Q, ck_plus, vk_plus, ck_minus, vk_minus, tag=None): exponents.append(BathExponent( "-", 2, Q, ckm, vkm, sigma_bar_k_offset=-1, tag=tag, )) - super().__init__(exponents) + super().__init__(exponents=exponents, tag=tag) + + @classmethod + def from_environment(cls, env, Q, dim=2): + """ + Converts from the "environment" API to the "bath" API. A + `FermionicBath` combines the information from an + `ExponentialFermionicEnvironment` and a coupling operator. + + Note that the "environment" API does not require fermionic exponents to + come in pairs. This method will add additional exponents with + zero coefficient in order to achieve the same amount of ``+`` and ``-`` + exponents. + + Parameters + ---------- + env : :class:`.ExponentialFermionicEnvironment` + The bath. + Q : Qobj + The coupling operator for the bath. + dim : optional, int or ``None`` (default ``2``) + The maximum number of excitations for each exponent. Usually ``2``. + """ + + # make same amount of plus and minus exponents by adding zeros + plus_exponents = [] + minus_exponents = [] + for exponent in env.exponents: + if exponent.type == BathExponent.types['+']: + plus_exponents.append(exponent) + else: + minus_exponents.append(exponent) + while len(plus_exponents) > len(minus_exponents): + minus_exponents.append(environment.CFExponent('-', 0, 0)) + while len(minus_exponents) > len(plus_exponents): + plus_exponents.append(environment.CFExponent('+', 0, 0)) + + bath_exponents = [] + for expp, expm in zip(plus_exponents, minus_exponents): + bath_exponents.append(BathExponent( + '+', dim, Q, expp.ck, expp.vk, None, 1, env.tag + )) + bath_exponents.append(BathExponent( + '-', dim, Q, expm.ck, expm.vk, None, -1, env.tag + )) + + result = cls(Q, [], [], [], [], tag=env.tag) + result.exponents = bath_exponents + return result class LorentzianBath(FermionicBath): @@ -827,36 +706,23 @@ class LorentzianBath(FermionicBath): A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from. - """ - def __init__(self, Q, gamma, w, mu, T, Nk, tag=None): - ck_plus, vk_plus = self._corr(gamma, w, mu, T, Nk, sigma=1.0) - ck_minus, vk_minus = self._corr(gamma, w, mu, T, Nk, sigma=-1.0) - - super().__init__( - Q, ck_plus, vk_plus, ck_minus, vk_minus, tag=tag, - ) - - def _corr(self, gamma, w, mu, T, Nk, sigma): - beta = 1. / T - kappa = [0.] - kappa.extend([1. for _ in range(1, Nk + 1)]) - epsilon = [0] - epsilon.extend([(2 * ll - 1) * np.pi for ll in range(1, Nk + 1)]) - - def f(x): - return 1 / (np.exp(x) + 1) - - eta_list = [0.5 * gamma * w * f(1.0j * beta * w)] - gamma_list = [w - sigma * 1.0j * mu] - for ll in range(1, Nk + 1): - eta_list.append( - -1.0j * (kappa[ll] / beta) * gamma * w**2 / - (-(epsilon[ll]**2 / beta**2) + w**2) - ) - gamma_list.append(epsilon[ll] / beta - sigma * 1.0j * mu) + Notes + ----- + This class is part of the "bath" API, which is now mirrored by the newer + "environment" API. The bath classes are kept in QuTiP for reasons of + backwards compatibility and convenience. Creating a `LorentzianBath` is + equivalent to creating a :class:`.LorentzianEnvironment`, performing a + :meth:`Matsubara <.LorentzianEnvironment.approx_by_matsubara>` + approximation, and finally bundling the result together with the coupling + operator ``Q`` for convenient use with the HEOM solver. + """ - return eta_list, gamma_list + def __new__(self, Q, gamma, w, mu, T, Nk, tag=None): + # See DrudeLorentzBath comment + env = environment.LorentzianEnvironment(T, mu, gamma, w) + mats_approx = env.approx_by_matsubara(Nk, tag) + return FermionicBath.from_environment(mats_approx, Q) class LorentzianPadeBath(FermionicBath): @@ -905,78 +771,20 @@ class LorentzianPadeBath(FermionicBath): A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from. - """ - def __init__(self, Q, gamma, w, mu, T, Nk, tag=None): - ck_plus, vk_plus = self._corr(gamma, w, mu, T, Nk, sigma=1.0) - ck_minus, vk_minus = self._corr(gamma, w, mu, T, Nk, sigma=-1.0) - super().__init__( - Q, ck_plus, vk_plus, ck_minus, vk_minus, tag=tag, - ) - - def _corr(self, gamma, w, mu, T, Nk, sigma): - beta = 1. / T - kappa, epsilon = self._kappa_epsilon(Nk) - - def f_approx(x): - f = 0.5 - for ll in range(1, Nk + 1): - f = f - 2 * kappa[ll] * x / (x**2 + epsilon[ll]**2) - return f - - eta_list = [0.5 * gamma * w * f_approx(1.0j * beta * w)] - gamma_list = [w - sigma * 1.0j * mu] + Notes + ----- + This class is part of the "bath" API, which is now mirrored by the newer + "environment" API. The bath classes are kept in QuTiP for reasons of + backwards compatibility and convenience. Creating a `LorentzianPadeBath` is + equivalent to creating a :class:`.LorentzianEnvironment`, performing a + :meth:`Pade <.LorentzianEnvironment.approx_by_pade>` approximation, and + finally bundling the result together with the coupling operator ``Q`` for + convenient use with the HEOM solver. + """ - for ll in range(1, Nk + 1): - eta_list.append( - -1.0j * (kappa[ll] / beta) * gamma * w**2 - / (-(epsilon[ll]**2 / beta**2) + w**2) - ) - gamma_list.append(epsilon[ll] / beta - sigma * 1.0j * mu) - - return eta_list, gamma_list - - def _kappa_epsilon(self, Nk): - eps = self._calc_eps(Nk) - chi = self._calc_chi(Nk) - - kappa = [0] - prefactor = 0.5 * Nk * (2 * (Nk + 1) - 1) - for j in range(Nk): - term = prefactor - for k in range(Nk - 1): - term *= ( - (chi[k]**2 - eps[j]**2) / - (eps[k]**2 - eps[j]**2 + self._delta(j, k)) - ) - for k in [Nk - 1]: - term /= (eps[k]**2 - eps[j]**2 + self._delta(j, k)) - kappa.append(term) - - epsilon = [0] + eps - - return kappa, epsilon - - def _delta(self, i, j): - return 1.0 if i == j else 0.0 - - def _calc_eps(self, Nk): - alpha = np.diag([ - 1. / np.sqrt((2 * k + 3) * (2 * k + 1)) - for k in range(2 * Nk - 1) - ], k=1) - alpha += alpha.transpose() - - evals = eigvalsh(alpha) - eps = [-2. / val for val in evals[0: Nk]] - return eps - - def _calc_chi(self, Nk): - alpha_p = np.diag([ - 1. / np.sqrt((2 * k + 5) * (2 * k + 3)) - for k in range(2 * Nk - 2) - ], k=1) - alpha_p += alpha_p.transpose() - evals = eigvalsh(alpha_p) - chi = [-2. / val for val in evals[0: Nk - 1]] - return chi + def __new__(self, Q, gamma, w, mu, T, Nk, tag=None): + # See DrudeLorentzBath comment + env = environment.LorentzianEnvironment(T, mu, gamma, w) + mats_approx = env.approx_by_pade(Nk, tag) + return FermionicBath.from_environment(mats_approx, Q) diff --git a/qutip/solver/heom/bofin_solvers.py b/qutip/solver/heom/bofin_solvers.py index 6009863228..88f544098d 100644 --- a/qutip/solver/heom/bofin_solvers.py +++ b/qutip/solver/heom/bofin_solvers.py @@ -20,10 +20,13 @@ from qutip import state_number_enumerate from qutip.core import data as _data from qutip.core.data import csr as _csr +from qutip.core.environment import ( + BosonicEnvironment, FermionicEnvironment, + ExponentialBosonicEnvironment, ExponentialFermionicEnvironment) from qutip.core import Qobj, QobjEvo from qutip.core.superoperator import liouvillian, spre, spost from .bofin_baths import ( - BathExponent, DrudeLorentzBath, + Bath, BathExponent, BosonicBath, DrudeLorentzBath, FermionicBath, ) from ..solver_base import Solver from .. import Result @@ -66,7 +69,7 @@ class HierarchyADOs: Parameters ---------- - exponents : list of BathExponent + exponents : list of :class:`.BathExponent` The exponents of the correlation function describing the bath or baths. @@ -76,7 +79,7 @@ class HierarchyADOs: Attributes ---------- - exponents : list of BathExponent + exponents : list of :class:`.BathExponent` The exponents of the correlation function describing the bath or baths. @@ -105,6 +108,7 @@ class HierarchyADOs: labels: list of tuples A list of the ADO labels within the hierarchy. """ + def __init__(self, exponents, max_depth): self.exponents = exponents self.max_depth = max_depth @@ -207,7 +211,7 @@ def exps(self, label): Returns ------- - tuple of BathExponent + tuple of :class:`.BathExponent` A tuple of BathExponents. Examples @@ -261,7 +265,7 @@ def filter(self, level=None, tags=None, dims=None, types=None): types : list of BathExponent types or list of str Filter parameter that matches the ``.type`` attribute of exponents. Types may be supplied by name (e.g. "R", "I", "+") - instead of by the actual type (e.g. ``BathExponent.types.R``). + instead of by the actual type (e.g. ``CFExponent.types.R``). Returns ------- @@ -339,6 +343,8 @@ class HierarchyADOsState: rho : Qobj The system state. + Notes + ----- In addition, all of the attributes of the hierarchy description, i.e. ``HierarchyADOs``, are provided directly on this class for convenience. E.g. one can access ``.labels``, or ``.exponents`` or @@ -371,7 +377,7 @@ def extract(self, idx_or_label): Returns ------- Qobj - A :obj:`~qutip.Qobj` representing the state of the specified ADO. + A :obj:`.Qobj` representing the state of the specified ADO. """ if isinstance(idx_or_label, int): idx = idx_or_label @@ -428,7 +434,8 @@ def heomsolve( Hierarchical Equations of Motion (HEOM) solver that supports multiple baths. - The baths must be all either bosonic or fermionic baths. + Each bath must be either bosonic or fermionic, but bosonic and fermionic + baths can be mixed. If you need to run many evolutions of the same system and bath, consider using :class:`HEOMSolver` directly to avoid having to continually @@ -438,16 +445,20 @@ def heomsolve( ---------- H : :obj:`.Qobj`, :obj:`.QobjEvo` Possibly time-dependent system Liouvillian or Hamiltonian as a Qobj or - QobjEvo. list of [:obj:`.Qobj`, :obj:`.Coefficient`] or callable that + QobjEvo. List of [:obj:`.Qobj`, :obj:`.Coefficient`], or callable that can be made into :obj:`.QobjEvo` are also accepted. - bath : Bath or list of Bath - A :obj:`Bath` containing the exponents of the expansion of the + bath : Bath specification or list of Bath specifications + A bath containing the exponents of the expansion of the bath correlation funcion and their associated coefficients and coupling operators, or a list of baths. - If multiple baths are given, they must all be either fermionic - or bosonic baths. + Each bath can be specified as *either* an object of type + :class:`.Bath`, :class:`.BosonicBath`, :class:`.FermionicBath`, or + their subtypes, *or* as a tuple ``(env, Q)``, where ``env`` is an + :class:`.ExponentialBosonicEnvironment` or an + :class:`.ExponentialFermionicEnvironment` and ``Q`` the system coupling + operator. max_depth : int The maximum depth of the heirarchy (i.e. the maximum number of bath @@ -536,9 +547,10 @@ def heomsolve( * ``ado_states``: the full ADO state at each time (only available if the results option ``ado_return`` was set to ``True``). - Each element is an instance of :class:`HierarchyADOsState`. + Each element is an instance of :class:`.HierarchyADOsState`. The state of a particular ADO may be extracted from - ``result.ado_states[i]`` by calling :meth:`extract`. + ``result.ado_states[i]`` by calling + :meth:`extract <.HierarchyADOsState.extract>`. * ``expect``: a list containing the values of each ``e_ops`` at time ``t``. @@ -559,7 +571,8 @@ class HEOMSolver(Solver): """ HEOM solver that supports multiple baths. - The baths must be all either bosonic or fermionic baths. + Each bath must be either bosonic or fermionic, but bosonic and fermionic + baths can be mixed. Parameters ---------- @@ -568,13 +581,17 @@ class HEOMSolver(Solver): QobjEvo. list of [:obj:`.Qobj`, :obj:`.Coefficient`] or callable that can be made into :obj:`.QobjEvo` are also accepted. - bath : Bath or list of Bath - A :obj:`Bath` containing the exponents of the expansion of the + bath : Bath specification or list of Bath specifications + A bath containing the exponents of the expansion of the bath correlation funcion and their associated coefficients and coupling operators, or a list of baths. - If multiple baths are given, they must all be either fermionic - or bosonic baths. + Each bath can be specified as *either* an object of type + :class:`.Bath`, :class:`.BosonicBath`, :class:`.FermionicBath`, or + their subtypes, *or* as a tuple ``(env, Q)``, where ``env`` is an + :class:`.ExponentialBosonicEnvironment` or an + :class:`.ExponentialFermionicEnvironment` and ``Q`` the system coupling + operator. odd_parity : Bool For fermionic baths only. Default is "False". "Parity" refers to the @@ -712,12 +729,15 @@ def _initialize_stats(self): def _combine_bath_exponents(self, bath): """ Combine the exponents for the specified baths. """ - if not isinstance(bath, (list, tuple)): - exponents = bath.exponents - else: - exponents = [] - for b in bath: - exponents.extend(b.exponents) + # Only one bath provided, not a list of baths + if (not isinstance(bath, (list, tuple)) + or self._is_environment_api(bath)): + bath = [bath] + bath = [self._to_bath(b) for b in bath] + exponents = [] + for b in bath: + exponents.extend(b.exponents) + if not all(exp.Q.dims == exponents[0].Q.dims for exp in exponents): raise ValueError( "All bath exponents must have system coupling operators" @@ -726,6 +746,40 @@ def _combine_bath_exponents(self, bath): ) return exponents + def _is_environment_api(self, bath_spec): + if not isinstance(bath_spec, (list, tuple)) or len(bath_spec) < 2: + return False + env, Q, *args = bath_spec + + if not isinstance(env, (BosonicEnvironment, FermionicEnvironment)): + return False + + if not isinstance(Q, (Qobj, QobjEvo)): + return False + + return True + + def _to_bath(self, bath_spec): + if isinstance(bath_spec, (Bath, BosonicBath, FermionicBath)): + return bath_spec + + if not self._is_environment_api(bath_spec): + raise ValueError( + "Environments must be passed as either Bath instances or" + " as a tuple or list corresponding to an environment and a" + " coupling operator, (env, Q)" + ) + env, Q, *args = bath_spec + + if isinstance(env, ExponentialBosonicEnvironment): + return BosonicBath.from_environment(env, Q, *args) + if isinstance(env, ExponentialFermionicEnvironment): + return FermionicBath.from_environment(env, Q, *args) + raise ValueError("The HEOM solver requires the environment to have" + " a multi-exponential correlation function. Use" + " one of the `approx_by_` functions to generate a" + " multi-exponential approximation.") + def _grad_n(self, he_n): """ Get the gradient for the hierarchy ADO at level n. """ vk = self.ados.vk @@ -1186,11 +1240,6 @@ class HSolverDL(HEOMSolver): See :class:`HEOMSolver` and :class:`DrudeLorentzBath` for more descriptions of the underlying solver and bath construction. - An exact copy of the QuTiP 4.6 HSolverDL is provided in - ``qutip.nonmarkov.dlheom_solver`` for cases where the functionality of - the older solver is required. The older solver will be completely - removed in QuTiP 5. - .. note:: Unlike the version of ``HSolverDL`` in QuTiP 4.6, this solver @@ -1222,10 +1271,12 @@ class HSolverDL(HEOMSolver): See parameter ``Q`` in :class:`BosonicBath` for a complete description. coup_strength : float - Coupling strength. Referred to as ``lam`` in :class:`DrudeLorentzBath`. + Coupling strength. Referred to as ``lam`` in + :class:`.DrudeLorentzEnvironment`. temperature : float - Bath temperature. Referred to as ``T`` in :class:`DrudeLorentzBath`. + Bath temperature. Referred to as ``T`` in + :class:`.DrudeLorentzEnvironment`. N_cut : int The maximum depth of the hierarchy. See ``max_depth`` in @@ -1238,7 +1289,7 @@ class HSolverDL(HEOMSolver): cut_freq : float Bath spectral density cutoff frequency. Referred to as ``gamma`` in - :class:`DrudeLorentzBath`. + :class:`.DrudeLorentzEnvironment`. bnd_cut_approx : bool Use boundary cut off approximation. If true, the Matsubara @@ -1253,9 +1304,11 @@ class HSolverDL(HEOMSolver): combine : bool, default: True Whether to combine exponents with the same frequency (and coupling - operator). See :meth:`BosonicBath.combine` for details. + operator). See :meth:`.ExponentialBosonicEnvironment.combine` for + details. Keyword only. Default: True. """ + def __init__( self, H_sys, coup_op, coup_strength, temperature, N_cut, N_exp, cut_freq, *, bnd_cut_approx=False, options=None, @@ -1304,6 +1357,7 @@ class _GatherHEOMRHS: nhe : int The number of ADOs in the hierarchy. """ + def __init__(self, f_idx, block, nhe): self._block_size = block self._n_blocks = nhe diff --git a/qutip/tests/core/test_environment.py b/qutip/tests/core/test_environment.py new file mode 100644 index 0000000000..61046c2246 --- /dev/null +++ b/qutip/tests/core/test_environment.py @@ -0,0 +1,1286 @@ +import pytest + +from numbers import Number + +import numpy as np +import mpmath as mp +from scipy.integrate import quad_vec +from qutip.utilities import n_thermal + +from qutip.core.environment import ( + BosonicEnvironment, + DrudeLorentzEnvironment, + UnderDampedEnvironment, + OhmicEnvironment, + CFExponent, + ExponentialBosonicEnvironment, + LorentzianEnvironment, + ExponentialFermionicEnvironment +) + + +def assert_guarantees(env, skip_sd=False, skip_cf=False, skip_ps=False): + """ + Checks the argument types accepted by the SD, CF and PS of the provided + environment, and that these functions satisfy certain symmetries + """ + functions = [] + if not skip_sd: + functions.append(env.spectral_density) + if not skip_cf: + functions.append(env.correlation_function) + if not skip_ps: + functions.append(env.power_spectrum) + + # SD, CF and PS can be called with number and array + # Result is required to have same type, ndim and len as argument + for fun in functions: + assert isinstance(fun(1), Number) + assert isinstance(fun(0.), Number) + assert isinstance(fun(-1.), Number) + + res = fun(np.array([])) + assert isinstance(res, np.ndarray) + assert res.ndim == 1 + assert len(res) == 0 + + res = fun(np.array([-1, 0, 1])) + assert isinstance(res, np.ndarray) + assert res.ndim == 1 + assert len(res) == 3 + + if not skip_sd: + # For negative frequencies, SD must return zero + res = env.spectral_density(np.linspace(-10, 10, 20)) + np.testing.assert_allclose(res[:10], np.zeros_like(res[:10])) + # SD must be real + np.testing.assert_allclose(np.imag(res), np.zeros_like(res)) + + if not skip_cf: + # CF symmetry + # Intentionally not testing at t=0 because of D-L spectral density + res = env.correlation_function(np.linspace(-10, 10, 20)) + res_reversed = res[::-1] + np.testing.assert_allclose(res, np.conjugate(res_reversed)) + + if not skip_ps: + # PS must be real + res = env.power_spectrum(np.linspace(-10, 10, 20)) + np.testing.assert_allclose(np.imag(res), np.zeros_like(res)) + +def assert_equivalent(env1, env2, *, tol, + skip_sd=False, skip_cf=False, skip_ps=False, + tMax=25, wMax=10): + """ + Checks that two environments have the same SD, CF and PS + (up to given tolerance) + """ + tlist = np.linspace(0, tMax, 100) + wlist = np.linspace(0, wMax, 100) + wlist2 = np.linspace(-wMax, wMax, 100) + + if not skip_sd: + assert_allclose(env1.spectral_density(wlist), + env2.spectral_density(wlist), tol) + if not skip_cf: + assert_allclose(env1.correlation_function(tlist), + env2.correlation_function(tlist), tol) + if not skip_ps: + assert_allclose(env1.power_spectrum(wlist2), + env2.power_spectrum(wlist2), tol) + +def assert_allclose(actual, desired, tol): + # We want to compare to arrays and provide both an abs and a rel tolerance + # We use the same parameter for both for simplicity + # However we reduce the abs tolerance if the numerical values are all small + atol = tol * min(np.max(np.abs(desired)), 1) + # If desired diverges somewhere, we allow actual to be anything + desired_finite = np.isfinite(desired) + np.testing.assert_allclose(actual[desired_finite], desired[desired_finite], + rtol=tol, atol=atol) + + +class UDReference: + def __init__(self, T, lam, gamma, w0): + self.T = T + self.lam = lam + self.gamma = gamma + self.w0 = w0 + + def spectral_density(self, w): + # only valid for w >= 0 + return self.lam**2 * self.gamma * w / ( + (w**2 - self.w0**2)**2 + (self.gamma * w)**2 + ) + + def power_spectrum(self, w, eps=1e-16): + # at zero frequency, take limit w->0 + w = np.array(w) + w[w == 0] = eps + return 2 * np.sign(w) * self.spectral_density(np.abs(w)) * ( + n_thermal(w, self.T) + 1 + ) + + def correlation_function(self, t, Nk=1000): + # only valid for t >= 0 + if self.T == 0: + return self._cf_zeroT(t) + + def coth(x): + return 1 / np.tanh(x) + + Om = np.sqrt(self.w0**2 - (self.gamma / 2)**2) + Gamma = self.gamma / 2 + beta = 1 / self.T + + ckAR = [ + (self.lam**2 / (4*Om)) * coth(beta * (Om + 1.0j * Gamma) / 2), + (self.lam**2 / (4*Om)) * coth(beta * (Om - 1.0j * Gamma) / 2), + ] + ckAR.extend( + (-2 * self.lam**2 * self.gamma / beta) * (2 * np.pi * k / beta) / + (((Om + 1.0j * Gamma)**2 + (2 * np.pi * k / beta)**2) * + ((Om - 1.0j * Gamma)**2 + (2 * np.pi * k / beta)**2)) + 0.j + for k in range(1, Nk + 1) + ) + vkAR = [ + -1.0j * Om + Gamma, + 1.0j * Om + Gamma, + ] + vkAR.extend( + 2 * np.pi * k * self.T + 0.j + for k in range(1, Nk + 1) + ) + + factor = 1j / 4 + ckAI = [-factor * self.lam**2 / Om, factor * self.lam**2 / Om] + vkAI = [-(-1j * Om - Gamma), -(1j * Om - Gamma)] + + result = np.zeros_like(t, dtype=complex) + for ck, vk in zip(ckAR, vkAR): + result += ck * np.exp(-vk * t) + for ck, vk in zip(ckAI, vkAI): + result += 1j * ck * np.exp(-vk * t) + return result + + def _cf_zeroT(self, t, eps=1e-3): + Om = np.sqrt(self.w0**2 - (self.gamma / 2)**2) + Gamma = self.gamma / 2 + + term1 = self.lam**2 / (2 * Om) * np.exp(-(1j * Om + self.gamma / 2) * t) + + def integrand(x, t): + return x * np.exp(-x * t) / ( + ((Om + 1j * Gamma)**2 + x**2) * ((Om - 1j * Gamma)**2 + x**2) + ) + integral = quad_vec( + lambda w: integrand(w, t), 0, np.inf, epsabs=eps, epsrel=eps + ) + result = term1 - self.gamma * self.lam**2 / np.pi * integral[0] + return result + + def _cf_by_integral(self, t, eps=1e-3): + # alternative implementation of correlation function + # currently unused because it makes tests slow + # but kept here in case we might need it in the future + def integrand(w, t): + return self.spectral_density(w) / np.pi * ( + (2 * n_thermal(w, self.T) + 1) * np.cos(w * t) + - 1j * np.sin(w * t) + ) + + result = quad_vec(lambda w: integrand(w, t), 0, np.inf, + epsabs=eps, epsrel=eps) + return result[0] + + +class DLReference: + def __init__(self, T, lam, gamma): + self.T = T + self.lam = lam + self.gamma = gamma + + def spectral_density(self, w): + # only valid for w >= 0 + return 2 * self.lam * self.gamma * w / (w**2 + self.gamma**2) + + def power_spectrum(self, w, eps=1e-16): + # at zero frequency, take limit w->0 + w = np.array(w) + w[w == 0] = eps + return 2 * np.sign(w) * self.spectral_density(np.abs(w)) * ( + n_thermal(w, self.T) + 1 + ) + + def correlation_function(self, t, Nk=1000): + def cot(x): + return 1 / np.tan(x) + + # C_real expansion terms: + ck_real = [self.lam * self.gamma / np.tan(self.gamma / (2 * self.T))] + ck_real.extend([ + (8 * self.lam * self.gamma * self.T * np.pi * k * self.T / + ((2 * np.pi * k * self.T)**2 - self.gamma**2)) + for k in range(1, Nk + 1) + ]) + vk_real = [self.gamma] + vk_real.extend([2 * np.pi * k * self.T for k in range(1, Nk + 1)]) + + # C_imag expansion terms (this is the full expansion): + ck_imag = [self.lam * self.gamma * (-1.0)] + vk_imag = [self.gamma] + + result = np.zeros_like(t, dtype=complex) + for ck, vk in zip(ck_real, vk_real): + result += ck * np.exp(-vk * np.abs(t)) + for ck, vk in zip(ck_imag, vk_imag): + result += self._sign(t) * 1j * ck * np.exp(-vk * np.abs(t)) + result[t == 0] += np.inf # real part of CF diverges at t=0 + return result + + def _sign(self, x): + return np.sign(x) + (x == 0) + + +class OhmicReference: + def __init__(self, T, alpha, wc, s): + self.T = T + self.alpha = alpha + self.wc = wc + self.s = s + + def spectral_density(self, w): + # only valid for w >= 0 + return (self.alpha * self.wc**(1 - self.s) * + w**self.s * np.exp(-w / self.wc)) + + def power_spectrum(self, w, eps=1e-16): + # at zero frequency, take limit w->0 + w = np.array(w) + w[w == 0] = eps + return 2 * np.sign(w) * self.spectral_density(np.abs(w)) * ( + n_thermal(w, self.T) + 1 + ) + + def correlation_function(self, t): + # only valid for t >= 0 + if self.T == 0: + return self._cf_zeroT(t) + + beta = 1 / self.T + corr = ( + (self.alpha / np.pi) * self.wc**(1 - self.s) * + beta**(-(self.s + 1)) * mp.gamma(self.s + 1) + ) + z1_u = (1 + beta * self.wc - 1.0j * self.wc * t) / (beta * self.wc) + z2_u = (1 + 1.0j * self.wc * t) / (beta * self.wc) + + return np.array([ + complex(corr * (mp.zeta(self.s + 1, u1) + mp.zeta(self.s + 1, u2))) + for u1, u2 in zip(z1_u, z2_u) + ]) + + def _cf_zeroT(self, t): + return ( + self.alpha / np.pi * self.wc**(self.s + 1) * + complex(mp.gamma(self.s + 1)) * + (1 + 1j * self.wc * t)**(-(self.s + 1)) + ) + + +class SingleExponentReference: + def __init__(self, coefficient, exponent, T): + self.coefficient = coefficient + self.exponent = exponent + self.T = T + + def spectral_density(self, w): + # only valid for w >= 0 + w = np.asarray(w) + result = self.power_spectrum(w) / 2 / (n_thermal(w, self.T) + 1) + result[w == 0] = 0 + return result + + def power_spectrum(self, w): + return 2 * np.real( + self.coefficient / (self.exponent - 1j * w) + ) + + def correlation_function(self, t): + # only valid for t >= 0 + return self.coefficient * np.exp(-self.exponent * t) + + +class TestBosonicEnvironment: + @pytest.mark.parametrize(["ref", "info"], [ + pytest.param(UDReference(gamma=.1, lam=.5, w0=1, T=.5), + {'tMax': 125, 'npoints': 300, 'tol': 5e-3}, + id="UD finite T"), + pytest.param(DLReference(gamma=.5, lam=.1, T=.5), + {'tMax': 15, 'npoints': 300, 'tol': 1e-2}, + id="DL finite T"), + pytest.param(OhmicReference(T=0, alpha=1, wc=10, s=1), + {'tMax': 10, 'npoints': 500, 'tol': 1e-2}, + id="Ohmic zero T"), + ]) + @pytest.mark.parametrize(["interpolate", "provide_tmax"], [ + [True, False], [False, False], [False, True] + ]) + @pytest.mark.parametrize("provide_temp", [True, False]) + def test_from_cf( + self, ref, info, interpolate, provide_tmax, provide_temp + ): + tMax = info['tMax'] + npoints = info['npoints'] + tol = info['tol'] + + # Collect arguments + if interpolate: + tlist = np.linspace(0, tMax, npoints) + clist = ref.correlation_function(tlist) + if np.isfinite(clist[0]): + args1 = {'tlist': tlist, 'C': clist} + else: + # CF may diverge at t=0, e.g. Drude-Lorentz + args1 = {'tlist': tlist[1:], 'C': clist[1:]} + else: + args1 = {'C': ref.correlation_function} + args2 = {'tMax': tMax} if provide_tmax else {} + args3 = {'T': ref.T} if provide_temp else {} + + env = BosonicEnvironment.from_correlation_function( + **args1, **args2, **args3, + ) + + # Determine which characteristic functions should be accessible + skip_ps = False + skip_sd = False + if not interpolate and not provide_tmax: + skip_ps = True + skip_sd = True + with pytest.raises(ValueError) as err: + env.power_spectrum(0) + assert str(err.value) == ( + 'The support of the correlation function (tMax) must be ' + 'specified for this operation.') + with pytest.raises(ValueError) as err: + env.spectral_density(0) + assert str(err.value) == ( + 'The support of the correlation function (tMax) must be ' + 'specified for this operation.') + elif not provide_temp: + skip_sd = True + with pytest.raises(ValueError) as err: + env.spectral_density(0) + assert str(err.value) == ( + 'The temperature must be specified for this operation.') + + assert_guarantees(env, skip_sd=skip_sd, skip_ps=skip_ps) + assert_equivalent(env, ref, skip_sd=skip_sd, skip_ps=skip_ps, + tol=tol, tMax=tMax) + + @pytest.mark.parametrize(["ref", "info"], [ + pytest.param(UDReference(gamma=.1, lam=.5, w0=1, T=.5), + {'wMax': 5, 'npoints': 200, 'tol': 5e-3}, + id="UD finite T"), + pytest.param(DLReference(gamma=.5, lam=.1, T=.5), + {'wMax': 200, 'npoints': 1500, 'tol': 5e-3}, + id="DL finite T"), + pytest.param(OhmicReference(T=0, alpha=1, wc=10, s=1), + {'wMax': 150, 'npoints': 150, 'tol': 5e-3}, + id="Ohmic zero T"), + ]) + @pytest.mark.parametrize(["interpolate", "provide_wmax"], [ + [True, False], [False, False], [False, True] + ]) + @pytest.mark.parametrize("provide_temp", [True, False]) + def test_from_sd( + self, ref, info, interpolate, provide_wmax, provide_temp + ): + wMax = info['wMax'] + npoints = info['npoints'] + tol = info['tol'] + + # Collect arguments + if interpolate: + wlist = np.linspace(0, wMax, npoints) + jlist = ref.spectral_density(wlist) + args1 = {'wlist': wlist, 'J': jlist} + else: + args1 = {'J': ref.spectral_density} + args2 = {'wMax': wMax} if provide_wmax else {} + args3 = {'T': ref.T} if provide_temp else {} + + env = BosonicEnvironment.from_spectral_density( + **args1, **args2, **args3, + ) + + # Determine which characteristic functions should be accessible + skip_ps = False + skip_cf = False + if not provide_temp: + skip_ps = True + skip_cf = True + with pytest.raises(ValueError) as err: + env.power_spectrum(0) + assert str(err.value) == ( + 'The temperature must be specified for this operation.') + with pytest.raises(ValueError) as err: + env.correlation_function(0) + assert str(err.value) == ( + 'The temperature must be specified for this operation.') + elif not interpolate and not provide_wmax: + skip_cf = True + with pytest.raises(ValueError) as err: + env.correlation_function(0) + assert str(err.value) == ( + 'The support of the spectral density (wMax) must be ' + 'specified for this operation.') + + assert_guarantees(env, skip_cf=skip_cf, skip_ps=skip_ps) + assert_equivalent(env, ref, skip_cf=skip_cf, skip_ps=skip_ps, + tol=tol, wMax=wMax) + + @pytest.mark.parametrize(["ref", "info"], [ + pytest.param(UDReference(gamma=.1, lam=.5, w0=1, T=.5), + {'wMax': 5, 'npoints': 200, 'tol': 5e-3}, + id="UD finite T"), + pytest.param(DLReference(gamma=.5, lam=.1, T=.5), + {'wMax': 200, 'npoints': 1500, 'tol': 5e-3}, + id="DL finite T"), + pytest.param(OhmicReference(T=0, alpha=1, wc=10, s=1), + {'wMax': 150, 'npoints': 2000, 'tol': 5e-3}, + id="Ohmic zero T"), + ]) + @pytest.mark.parametrize(["interpolate", "provide_wmax"], [ + [True, False], [False, False], [False, True] + ]) + @pytest.mark.parametrize("provide_temp", [True, False]) + def test_from_ps( + self, ref, info, interpolate, provide_wmax, provide_temp + ): + wMax = info['wMax'] + npoints = info['npoints'] + tol = info['tol'] + + # Collect arguments + if interpolate: + wlist = np.linspace(-wMax, wMax, 2 * npoints + 1) + slist = ref.power_spectrum(wlist) + args1 = {'wlist': wlist, 'S': slist} + else: + args1 = {'S': ref.power_spectrum} + args2 = {'wMax': wMax} if provide_wmax else {} + args3 = {'T': ref.T} if provide_temp else {} + + env = BosonicEnvironment.from_power_spectrum( + **args1, **args2, **args3, + ) + + # Determine which characteristic functions should be accessible + skip_sd = False + skip_cf = False + if not provide_temp: + skip_sd = True + with pytest.raises(ValueError) as err: + env.spectral_density(0) + assert str(err.value) == ( + 'The temperature must be specified for this operation.') + if not interpolate and not provide_wmax: + skip_cf = True + with pytest.raises(ValueError) as err: + env.correlation_function(0) + assert str(err.value) == ( + 'The support of the power spectrum (wMax) must be ' + 'specified for this operation.') + + assert_guarantees(env, skip_cf=skip_cf, skip_sd=skip_sd) + assert_equivalent(env, ref, skip_cf=skip_cf, skip_sd=skip_sd, + tol=tol, wMax=wMax) + + @pytest.mark.parametrize(["reference", "tMax"], [ + pytest.param(DLReference(.5, .1, .5), 15, id="DL Example"), + ]) + @pytest.mark.parametrize(["full_ansatz", "tol"], [ + pytest.param(True, 1e-3, id="Full ansatz"), + pytest.param(False, 5e-3, id="Not-full ansatz"), + ]) + def test_fixed_cf_fit(self, reference, tMax, full_ansatz, tol): + # fixed number of exponents + env = BosonicEnvironment.from_correlation_function( + reference.correlation_function, T=reference.T + ) + tlist = np.linspace(0, tMax, 100)[1:] # exclude t=0 + fit, info = env.approx_by_cf_fit( + tlist, target_rsme=None, Nr_max=2, Ni_max=2, + full_ansatz=full_ansatz, combine=False + ) + + assert isinstance(fit, ExponentialBosonicEnvironment) + assert len(fit.exponents) == 8 + assert fit.T == env.T + assert fit.tag is None + assert_equivalent( + fit, env, tol=tol, skip_sd=True, skip_ps=True, tMax=tMax + ) + + assert info["Nr"] == 2 + assert info["Ni"] == 2 + assert info["rmse_real"] < 5e-3 + assert info["rmse_imag"] < 5e-3 + for key in ["fit_time_real", "fit_time_imag", + "params_real", "params_imag", "summary"]: + assert key in info + + @pytest.mark.parametrize(["reference", "tMax", "tol"], [ + pytest.param(DLReference(.5, .1, .5), 15, 5e-2, id="DL Example"), + pytest.param(SingleExponentReference(1, 2 + .5j, None), 10, 1e-5, + id="Exponential function") + ]) + @pytest.mark.parametrize("full_ansatz", [True, False]) + def test_dynamic_cf_fit(self, reference, tMax, tol, full_ansatz): + # dynamic number of exponents + env = BosonicEnvironment.from_correlation_function( + reference.correlation_function, tag="test" + ) + tlist = np.linspace(0, tMax, 100)[1:] # exclude t=0 + fit, info = env.approx_by_cf_fit( + tlist, target_rsme=0.01, Nr_max=3, Ni_max=3, + full_ansatz=full_ansatz + ) + + assert isinstance(fit, ExponentialBosonicEnvironment) + assert fit.T == env.T + assert fit.tag == ("test", "CF Fit") + assert_equivalent( + fit, env, tol=tol, skip_sd=True, skip_ps=True, tMax=tMax + ) + + assert info["Nr"] == 1 + assert info["Ni"] == 1 + assert info["rmse_real"] <= 0.01 + assert info["rmse_imag"] <= 0.01 + for key in ["fit_time_real", "fit_time_imag", + "params_real", "params_imag", "summary"]: + assert key in info + + @pytest.mark.parametrize(["reference", "wMax", "tol"], [ + pytest.param(OhmicReference(3, .75, 10, 1), 15, 5e-2, id="DL Example"), + ]) + def test_fixed_sd_fit(self, reference, wMax, tol): + # fixed number of lorentzians + env = BosonicEnvironment.from_spectral_density( + reference.spectral_density, T=reference.T + ) + wlist = np.linspace(0, wMax, 100) + fit, info = env.approx_by_sd_fit( + wlist, Nk=1, target_rmse=None, Nmax=4, combine=False + ) + + assert isinstance(fit, ExponentialBosonicEnvironment) + assert len(fit.exponents) == 4 * (1 + 4) + assert fit.T == env.T + assert fit.tag is None + assert_equivalent( + fit, env, tol=tol, skip_cf=True, skip_ps=True, wMax=wMax + ) + + assert info["N"] == 4 + assert info["Nk"] == 1 + assert info["rmse"] < 5e-3 + for key in ["fit_time", "params", "summary"]: + assert key in info + + @pytest.mark.parametrize(["reference", "wMax", "tol", "params"], [ + pytest.param(OhmicReference(3, .75, 10, 1), 15, .2, {}, + id="DL Example"), + pytest.param(UDReference(1, .5, .1, 1), 2, 1e-4, + {'guess': [.2, .5, 1]}, id='UD Example'), + ]) + def test_dynamic_sd_fit(self, reference, wMax, tol, params): + # dynamic number of exponents + env = BosonicEnvironment.from_spectral_density( + reference.spectral_density, T=reference.T, tag="test" + ) + wlist = np.linspace(0, wMax, 100) + fit, info = env.approx_by_sd_fit( + wlist, Nk=1, target_rmse=0.01, Nmax=5, **params + ) + + assert isinstance(fit, ExponentialBosonicEnvironment) + assert fit.T == env.T + assert fit.tag == ("test", "SD Fit") + assert_equivalent( + fit, env, tol=tol, skip_cf=True, skip_ps=True, wMax=wMax + ) + + assert info["N"] < 5 + assert info["Nk"] == 1 + assert info["rmse"] < 0.01 + for key in ["fit_time", "params", "summary"]: + assert key in info + + +@pytest.mark.parametrize("params", [ + {'gamma': 2.5, 'lam': .75, 'T': 1.5} +]) +class TestDLEnvironment: + def test_matches_reference(self, params): + ref = DLReference(**params) + env = DrudeLorentzEnvironment(**params) + + assert_guarantees(env) + assert_equivalent(env, ref, tol=1e-8) + + def test_zero_temperature(self, params): + params_T0 = {**params, 'T': 0} + ref = DLReference(**params_T0) + env = DrudeLorentzEnvironment(**params_T0) + + with pytest.raises(ValueError) as err: + env.correlation_function(0) + assert str(err.value) == ( + 'The Drude-Lorentz correlation function ' + 'diverges at zero temperature.') + with pytest.raises(ValueError) as err: + env.approx_by_matsubara(10) + assert str(err.value) == ( + 'The Drude-Lorentz correlation function ' + 'diverges at zero temperature.') + with pytest.raises(ValueError) as err: + env.approx_by_pade(10) + assert str(err.value) == ( + 'The Drude-Lorentz correlation function ' + 'diverges at zero temperature.') + + assert_guarantees(env, skip_cf=True) + assert_equivalent(env, ref, tol=1e-6, skip_cf=True) + + @pytest.mark.parametrize("tag", [None, "test"]) + def test_matsubara_approx(self, params, tag): + Nk = 25 + original_tag = object() + env = DrudeLorentzEnvironment(**params, tag=original_tag) + + approx = env.approx_by_matsubara(Nk, combine=False, tag=tag) + assert isinstance(approx, ExponentialBosonicEnvironment) + assert len(approx.exponents) == Nk + 2 # (Nk+1) real + 1 imag + if tag is None: + assert approx.tag == (original_tag, "Matsubara Truncation") + else: + assert approx.tag == tag + assert approx.T == env.T + + # CF at t=0 might not match, which is okay + assert_equivalent(approx, env, tol=1e-2, skip_cf=True) + + approx_combine, delta = env.approx_by_matsubara( + Nk, compute_delta=True, combine=True + ) + assert isinstance(approx_combine, ExponentialBosonicEnvironment) + assert len(approx_combine.exponents) < Nk + 2 + assert approx_combine.T == env.T + + assert_equivalent(approx_combine, approx, tol=1e-8) + + # Check terminator amplitude + delta_ref = 2 * env.lam * env.T / env.gamma - 1j * env.lam + for exp in approx_combine.exponents: + delta_ref -= exp.coefficient / exp.exponent + assert_allclose(delta, delta_ref, tol=1e-8) + + @pytest.mark.parametrize("tag", [None, "test"]) + def test_pade_approx(self, params, tag): + Nk = 4 + original_tag = object() + env = DrudeLorentzEnvironment(**params, tag=original_tag) + + approx = env.approx_by_pade(Nk, combine=False, tag=tag) + assert isinstance(approx, ExponentialBosonicEnvironment) + assert len(approx.exponents) == Nk + 2 # (Nk+1) real + 1 imag + if tag is None: + assert approx.tag == (original_tag, "Pade Truncation") + else: + assert approx.tag == tag + assert approx.T == env.T + + # Wow, Pade is so much better + assert_equivalent(approx, env, tol=1e-8, skip_cf=True) + + approx_combine, delta = env.approx_by_pade( + Nk, combine=True, compute_delta=True + ) + assert isinstance(approx_combine, ExponentialBosonicEnvironment) + assert len(approx_combine.exponents) < Nk + 2 + assert approx_combine.T == env.T + + assert_equivalent(approx_combine, approx, tol=1e-8) + + # Check terminator amplitude + delta_ref = 2 * env.lam * env.T / env.gamma - 1j * env.lam + for exp in approx_combine.exponents: + delta_ref -= exp.coefficient / exp.exponent + assert_allclose(delta, delta_ref, tol=1e-8) + + +class TestUDEnvironment: + @pytest.mark.parametrize("params", [ + pytest.param({'gamma': 2.5, 'lam': .75, 'w0': 5, 'T': 1.5}, + id='finite T'), + pytest.param({'gamma': 2.5, 'lam': .75, 'w0': 5, 'T': 0}, + id='zero T'), + ]) + def test_matches_reference(self, params): + ref = UDReference(**params) + env = UnderDampedEnvironment(**params) + + assert_guarantees(env) + # a bit higher tolerance since we currently compute CF via FFT + assert_equivalent(env, ref, tol=1e-3) + + @pytest.mark.parametrize("params", [ + {'gamma': 2.5, 'lam': .75, 'w0': 5, 'T': 0}, + ]) + def test_zero_temperature(self, params): + # Attempting a Matsubara expansion at T=0 should raise a warning + # and return only the resonant (non-Matsubara) part of the expansion + env = UnderDampedEnvironment(**params) + + tlist = np.linspace(0, 5, 100) + resonant_approx = env.approx_by_matsubara(Nk=0) + + with pytest.warns(UserWarning) as record: + test_approx = env.approx_by_matsubara(Nk=3) + assert str(record[0].message) == ( + 'The Matsubara expansion cannot be performed at zero temperature. ' + 'Use other approaches such as fitting the correlation function.') + + assert_allclose(resonant_approx.correlation_function(tlist), + test_approx.correlation_function(tlist), + tol=1e-8) + + @pytest.mark.parametrize("params", [ + {'gamma': 2.5, 'lam': .75, 'w0': 5, 'T': 1.5}, + ]) + @pytest.mark.parametrize("tag", [None, "test"]) + def test_matsubara_approx(self, params, tag): + Nk = 25 + original_tag = object() + env = UnderDampedEnvironment(**params, tag=original_tag) + + approx = env.approx_by_matsubara(Nk, combine=False, tag=tag) + assert isinstance(approx, ExponentialBosonicEnvironment) + assert len(approx.exponents) == Nk + 4 # (Nk+2) real + 2 imag + if tag is None: + assert approx.tag == (original_tag, "Matsubara Truncation") + else: + assert approx.tag == tag + assert approx.T == env.T + + assert_equivalent(approx, env, tol=1e-3) + + approx_combine = env.approx_by_matsubara(Nk, combine=True) + assert isinstance(approx_combine, ExponentialBosonicEnvironment) + assert len(approx_combine.exponents) < Nk + 4 + assert approx_combine.T == env.T + + assert_equivalent(approx_combine, approx, tol=1e-8) + + +@pytest.mark.parametrize("params", [ + pytest.param({'alpha': .75, 'wc': 10, 's': 1, 'T': 3}, + id='finite T ohmic'), + pytest.param({'alpha': .75, 'wc': .5, 's': 1, 'T': 0}, + id='zero T ohmic'), + pytest.param({'alpha': .75, 'wc': 10, 's': .5, 'T': 3}, + id='finite T subohmic'), + pytest.param({'alpha': .75, 'wc': .5, 's': .5, 'T': 0}, + id='zero T subohmic'), + pytest.param({'alpha': .75, 'wc': 10, 's': 5, 'T': 3}, + id='finite T superohmic'), + pytest.param({'alpha': .75, 'wc': .5, 's': 5, 'T': 0}, + id='zero T superohmic'), +]) +class TestOhmicEnvironment: + def test_matches_reference(self, params): + ref = OhmicReference(**params) + env = OhmicEnvironment(**params) + + assert_guarantees(env) + assert_equivalent(env, ref, tol=1e-8) + + +class TestCFExponent: + def test_create(self): + exp_r = CFExponent("R", ck=1.0, vk=2.0) + check_exponent(exp_r, "R", 1.0, 2.0) + + exp_i = CFExponent("I", ck=1.0, vk=2.0) + check_exponent(exp_i, "I", 1.0j, 2.0) + + exp_ri = CFExponent("RI", ck=1.0, vk=2.0, ck2=3.0) + check_exponent(exp_ri, "RI", 1.0 + 3.0j, 2.0) + + exp_p = CFExponent("+", ck=1.0, vk=2.0) + check_exponent(exp_p, "+", 1.0, 2.0) + + exp_m = CFExponent("-", ck=1.0, vk=2.0) + check_exponent(exp_m, "-", 1.0, 2.0) + + exp_tag = CFExponent("R", ck=1.0, vk=2.0, tag="tag1") + check_exponent(exp_tag, "R", 1.0, 2.0, tag="tag1") + + for exp_type in ["R", "I", "+", "-"]: + with pytest.raises(ValueError) as err: + CFExponent(exp_type, ck=1.0, vk=2.0, ck2=3.0) + assert str(err.value) == ( + "Second co-efficient (ck2) should only be specified for RI" + " exponents" + ) + + with pytest.raises(ValueError) as err: + CFExponent('RI', ck=1.0, vk=2.0) + assert str(err.value) == "RI exponents require ck2" + + def test_repr(self): + exp1 = CFExponent("R", ck=1.0, vk=2.0) + assert repr(exp1) == ( + "" + ) + exp2 = CFExponent("+", ck=1.0, vk=2.0, tag="bath1") + assert repr(exp2) == ( + "" + ) + + +def check_exponent(exp, type, coefficient, exponent, tag=None): + assert exp.type is CFExponent.types[type] + assert exp.fermionic == (type in ["+", "-"]) + assert exp.coefficient == pytest.approx(coefficient) + assert exp.exponent == pytest.approx(exponent) + assert exp.tag == tag + + +class TestExpBosonicEnv: + def test_create(self): + env = ExponentialBosonicEnvironment([1.], [0.5], [2.], [0.6]) + exp_r, exp_i = env.exponents + check_exponent(exp_r, "R", 1.0, 0.5) + check_exponent(exp_i, "I", 2.0j, 0.6) + + env = ExponentialBosonicEnvironment([1.], [0.5], [2.], [0.5]) + [exp_ri] = env.exponents + check_exponent(exp_ri, "RI", 1.0 + 2.0j, 0.5) + + env = ExponentialBosonicEnvironment( + [1.], [0.5], [2.], [0.6], tag="bath1" + ) + exp_r, exp_i = env.exponents + check_exponent(exp_r, "R", 1.0, 0.5, tag="bath1") + check_exponent(exp_i, "I", 2.0j, 0.6, tag="bath1") + + exp1 = CFExponent("RI", ck=1.0, vk=2.0, ck2=3.0) + exp2 = CFExponent("I", ck=1.0, vk=2.0) + env = ExponentialBosonicEnvironment( + exponents=[exp1, exp2], combine=False) + assert env.exponents == [exp1, exp2] + + env = ExponentialBosonicEnvironment( + [], [], [1.0], [2.0], exponents=[exp1], combine=True) + [exp_ri] = env.exponents + check_exponent(exp_ri, "RI", 1.0 + 4.0j, 2) + + with pytest.raises(ValueError) as err: + ExponentialBosonicEnvironment([1.], [], [2.], [0.6]) + assert str(err.value) == ( + "The exponent lists ck_real and vk_real, and ck_imag and" + " vk_imag must be the same length." + ) + + with pytest.raises(ValueError) as err: + ExponentialBosonicEnvironment([1.], [0.5], [2.], []) + assert str(err.value) == ( + "The exponent lists ck_real and vk_real, and ck_imag and" + " vk_imag must be the same length." + ) + + with pytest.raises(ValueError) as err: + ExponentialBosonicEnvironment() + assert str(err.value) == ( + "Either the parameter `exponents` or the parameters " + "`ck_real`, `vk_real`, `ck_imag`, `vk_imag` must be provided." + ) + + with pytest.raises(ValueError) as err: + ExponentialBosonicEnvironment( + exponents=[exp1, CFExponent('+', 1.0, 2.0)] + ) + assert str(err.value) == ( + "Fermionic exponent passed to exponential bosonic environment." + ) + + @pytest.mark.parametrize("provide_temp", [True, False]) + def test_matches_reference(self, provide_temp): + if provide_temp: + args = {'T': .5} + skip_sd = False + else: + args = {} + skip_sd = True + + env = ExponentialBosonicEnvironment([1.], [0.5], [2.], [0.5], **args) + assert_guarantees(env, skip_sd=skip_sd) + ref = SingleExponentReference(coefficient=(1 + 2j), exponent=.5, T=.5) + assert_equivalent(env, ref, tol=1e-8, skip_sd=skip_sd) + + if not provide_temp: + with pytest.raises(ValueError) as err: + env.spectral_density(0) + assert str(err.value) == ( + 'The temperature must be specified for this operation.' + ) + + +# ----- Fermionic Environments ----- + + +def assert_guarantees_f(env, check_db=False, beta=None, mu=None): + """ + Checks the argument types accepted by the SD, CFs and PSs of the provided + fermionic environment, that the CFs satisfy their symmetries and that the + PSs are real. If `check_db` is True, `beta` and `mu` must be provided and + we also check the detailed balance condition S^- = e^{β(ω-μ)} S^+. + """ + + # SD, CF and PS can be called with number and array + # Result is required to have same type, ndim and len as argument + for fun in [env.spectral_density, + env.correlation_function_plus, + env.correlation_function_minus, + env.power_spectrum_plus, + env.power_spectrum_minus]: + assert isinstance(fun(1), Number) + assert isinstance(fun(0.), Number) + assert isinstance(fun(-1.), Number) + + res = fun(np.array([])) + assert isinstance(res, np.ndarray) + assert res.ndim == 1 + assert len(res) == 0 + + res = fun(np.array([-1, 0, 1])) + assert isinstance(res, np.ndarray) + assert res.ndim == 1 + assert len(res) == 3 + + # SD must be real + sd = env.spectral_density(np.linspace(-10, 10, 20)) + np.testing.assert_allclose(np.imag(sd), np.zeros_like(sd)) + + # CF symmetry + for fun in [env.correlation_function_plus, env.correlation_function_minus]: + cf = fun(np.linspace(-10, 10, 20)) + cf_reversed = cf[::-1] + np.testing.assert_allclose(cf, np.conjugate(cf_reversed)) + + # PS must be real and their sum must be SD + wlist = np.linspace(-10, 10, 20) + psp = env.power_spectrum_plus(wlist) + psm = env.power_spectrum_minus(wlist) + np.testing.assert_allclose(np.imag(psp), np.zeros_like(psp)) + np.testing.assert_allclose(np.imag(psm), np.zeros_like(psm)) + np.testing.assert_allclose(psp + psm, sd) + + # Detailed balance + if not check_db: + return + + if beta == np.inf: + foo = psp[wlist > mu] + np.testing.assert_allclose(foo, np.zeros_like(foo)) + bar = psp[wlist < mu] + np.testing.assert_allclose(bar, np.zeros_like(bar)) + return + + factor = np.exp(beta * (wlist - mu)) + np.testing.assert_allclose(psm, factor * psp) + +def assert_equivalent_f(env1, env2, *, tol, + tMax=25, wMax=10): + """ + Checks that two fermionic environments have the same SD, CFs and PSs + (up to given tolerance) + """ + tlist = np.linspace(0, tMax, 100) + wlist = np.linspace(-wMax, wMax, 100) + + assert_allclose(env1.spectral_density(wlist), + env2.spectral_density(wlist), tol) + assert_allclose(env1.correlation_function_plus(tlist), + env2.correlation_function_plus(tlist), tol) + assert_allclose(env1.correlation_function_minus(tlist), + env2.correlation_function_minus(tlist), tol) + assert_allclose(env1.power_spectrum_plus(wlist), + env2.power_spectrum_plus(wlist), tol) + assert_allclose(env1.power_spectrum_minus(wlist), + env2.power_spectrum_minus(wlist), tol) + + +class LorentzianReference: + def __init__(self, T, mu, gamma, W, omega0): + self.T = T + self.mu = mu + self.gamma = gamma + self.W = W + self.omega0 = omega0 + + def _f(self, x): + if self.T == 0: + return np.heaviside(-x, 0.5) + else: + return 1 / (np.exp(x / self.T) + 1) + + def spectral_density(self, w): + return self.gamma * self.W**2 / ( + (w - self.omega0)**2 + self.W**2 + ) + + def power_spectrum_plus(self, w): + return self.spectral_density(w) * self._f(w - self.mu) + + def power_spectrum_minus(self, w): + return self.spectral_density(w) * self._f(-w + self.mu) + + def correlation_function_plus(self, t, Nk=5000): + # only valid for t >= 0 + T = self.T + mu = self.mu + gamma = self.gamma + W = self.W + omega0 = self.omega0 + + # zero temperature CFs not implemented yet + if T == 0: + assert False + + ck = [(gamma * W / 2 * self._f(omega0 - mu + 1j * W))] + ck.extend( + 1j * gamma * W**2 * T / ( + (1j * (omega0 - mu) + (2 * k - 1) * np.pi * T)**2 - W**2 + ) for k in range(1, Nk + 1) + ) + vk = [W - 1j * omega0] + vk.extend( + (2 * k - 1) * np.pi * T - 1j * mu + for k in range(1, Nk + 1) + ) + + result = np.zeros_like(t, dtype=complex) + for c, v in zip(ck, vk): + result += c * np.exp(-v * t) + return result + + def correlation_function_minus(self, t, Nk=5000): + # only valid for t >= 0 + T = self.T + mu = self.mu + gamma = self.gamma + W = self.W + omega0 = self.omega0 + + # zero temperature CFs not implemented yet + if T == 0: + assert False + + ck = [(gamma * W / 2 * self._f(mu - omega0 + 1j * W))] + ck.extend( + 1j * gamma * W**2 * T / ( + (1j * (mu - omega0) + (2 * k - 1) * np.pi * T)**2 - W**2 + ) for k in range(1, Nk + 1) + ) + vk = [W + 1j * omega0] + vk.extend( + (2 * k - 1) * np.pi * T + 1j * mu + for k in range(1, Nk + 1) + ) + + result = np.zeros_like(t, dtype=complex) + for c, v in zip(ck, vk): + result += c * np.exp(-v * t) + return result + + +class SinglePlusExpReference: + def __init__(self, coefficient, exponent): + self.coefficient = coefficient + self.exponent = exponent + + def spectral_density(self, w): + return self.power_spectrum_plus(w) + + def power_spectrum_plus(self, w): + return 2 * np.real( + self.coefficient / (self.exponent + 1j * w) + ) + + def power_spectrum_minus(self, w): + return np.zeros_like(w) + + def correlation_function_plus(self, t): + # only valid for t >= 0 + return self.coefficient * np.exp(-self.exponent * t) + + def correlation_function_minus(self, t): + return np.zeros_like(t) + + +class SingleMinusExpReference: + def __init__(self, coefficient, exponent): + self.coefficient = coefficient + self.exponent = exponent + + def spectral_density(self, w): + return self.power_spectrum_minus(w) + + def power_spectrum_plus(self, w): + return np.zeros_like(w) + + def power_spectrum_minus(self, w): + return 2 * np.real( + self.coefficient / (self.exponent - 1j * w) + ) + + def correlation_function_plus(self, t): + return np.zeros_like(t) + + def correlation_function_minus(self, t): + # only valid for t >= 0 + return self.coefficient * np.exp(-self.exponent * t) + + +@pytest.mark.parametrize("params", [ + pytest.param({'T': 1.5, 'mu': 1, 'gamma': .75, 'W': .5, 'omega0': 1}, + id='finite T - mu=omega0'), + pytest.param({'T': 1.5, 'mu': 1, 'gamma': .75, 'W': .5, 'omega0': 3}, + id='finite T - mu can't be compared using '=='. + # Check that all properties of the CFExponents are retained. + for exp1, exp2 in zip(list1, list2): + if (exp1.type != exp2.type or + exp1.ck != exp2.ck or + exp1.vk != exp2.vk or + exp1.ck2 != exp2.ck2 or + exp1.tag != exp2.tag): + return False + return True + + def test_create_bosonic(self): + Q = sigmaz() + H = sigmax() + exponents = [ + CFExponent("R", ck=1.1, vk=2.1), + CFExponent("I", ck=1.2, vk=2.2), + CFExponent("RI", ck=1.3, vk=2.3, ck2=3.3), + ] + env = ExponentialBosonicEnvironment(exponents=exponents) + + hsolver = HEOMSolver(H, (env, Q), 2) + assert self._exponents_equal(hsolver.ados.exponents, exponents) + assert hsolver.ados.max_depth == 2 + + hsolver = HEOMSolver(H, [(env, Q)], 2) + assert self._exponents_equal(hsolver.ados.exponents, exponents) + assert hsolver.ados.max_depth == 2 + + hsolver = HEOMSolver(H, [(env, Q)] * 3, 2) + assert self._exponents_equal(hsolver.ados.exponents, exponents * 3) + assert hsolver.ados.max_depth == 2 + + def test_create_fermionic(self): + Q = sigmaz() + H = sigmax() + exponents = [ + CFExponent("+", ck=1.1, vk=2.1), + CFExponent("-", ck=1.2, vk=2.2), + ] + env = ExponentialFermionicEnvironment(exponents=exponents) + + hsolver = HEOMSolver(H, (env, Q), 2) + assert self._exponents_equal(hsolver.ados.exponents, exponents) + assert hsolver.ados.max_depth == 2 + + hsolver = HEOMSolver(H, [(env, Q)], 2) + assert self._exponents_equal(hsolver.ados.exponents, exponents) + assert hsolver.ados.max_depth == 2 + + hsolver = HEOMSolver(H, [(env, Q)] * 3, 2) + assert self._exponents_equal(hsolver.ados.exponents, exponents * 3) + assert hsolver.ados.max_depth == 2 + + def test_create_mixed_api(self): + Q = sigmaz() + H = sigmax() + exponent = CFExponent("R", ck=1.2, vk=2.2) + bath_exponent = BathExponent("R", 2, Q=Q, ck=1.2, vk=2.2) + env = ExponentialBosonicEnvironment(exponents=[exponent]) + bath = Bath([bath_exponent]) + + hsolver = HEOMSolver(H, [(env, Q), bath, (env, Q)], 2) + assert self._exponents_equal(hsolver.ados.exponents, [exponent] * 3) + assert hsolver.ados.max_depth == 2 + + hsolver = HEOMSolver(H, [bath, (env, Q), bath] * 3, 2) + assert self._exponents_equal(hsolver.ados.exponents, [exponent] * 3) + assert hsolver.ados.max_depth == 2 + + def test_create_bath_errors(self): + Q = sigmaz() + H = sigmax() + exponents = [ + CFExponent("I", ck=1.2, vk=2.2), + CFExponent("R", ck=1.2, vk=2.2), + ] + env = ExponentialBosonicEnvironment(exponents=exponents) + + with pytest.raises(ValueError) as err: + HEOMSolver(H, [(env, Q), (env, Q & Q)], 2) + assert str(err.value) == ( + "All bath exponents must have system coupling operators with the" + " same dimensions but a mixture of dimensions was given." + ) + + def test_create_h_sys_errors(self): + H = object() + Q = sigmax() + empty_env = ExponentialBosonicEnvironment(exponents=[]) + + with pytest.raises(TypeError) as err: + HEOMSolver(H, (empty_env, Q), 2) + assert str(err.value) == ( + "The Hamiltonian (H) must be a Qobj or QobjEvo" + ) + + H = [sigmax()] + with pytest.raises(TypeError) as err: + HEOMSolver([H], (empty_env, Q), 2) + assert str(err.value) == ( + "The Hamiltonian (H) must be a Qobj or QobjEvo" + ) + + @pytest.mark.parametrize(['method'], [ + pytest.param("run", id="run"), + pytest.param("start", id="start"), + ]) + def test_invalid_rho0_errors(self, method): + Q = sigmaz() + H = sigmax() + exponents = [ + CFExponent("R", ck=1.1, vk=2.1), + CFExponent("I", ck=1.2, vk=2.2), + CFExponent("RI", ck=1.3, vk=2.3, ck2=3.3), + ] + env = ExponentialBosonicEnvironment(exponents=exponents) + hsolver = HEOMSolver(H, (env, Q), 2) + + if method == "run": + def solve_method(rho0): + return hsolver.run(rho0, [0, 1]) + elif method == "start": + def solve_method(rho0): + return hsolver.start(rho0, 0) + else: + assert False, f"method {method} not supported by test" + + with pytest.raises(ValueError) as err: + solve_method(basis(3, 0)) + assert str(err.value) == ( + "Initial state rho has dims [[3], [1]]" + " but the system dims are [[2], [2]]" + ) + + with pytest.raises(TypeError) as err: + solve_method("batman") + assert str(err.value) == ( + "Initial ADOs passed have type but a " + "HierarchyADOsState or a numpy array-like instance was expected" + ) + + with pytest.raises(ValueError) as err: + solve_method(np.array([1, 2, 3])) + assert str(err.value) == ( + "Initial ADOs passed have shape (3,) but the solver hierarchy" + " shape is (10, 2, 2)" + ) + + @pytest.mark.parametrize(['evo'], [ + pytest.param("qobj", id="qobj"), + pytest.param("qobjevo_const", id="qobjevo_const"), + pytest.param("qobjevo_timedep", id="qobjevo_timedep"), + ]) + @pytest.mark.parametrize(['liouvillianize'], [ + pytest.param(False, id="hamiltonian"), + pytest.param(True, id="liouvillian"), + ]) + def test_pure_dephasing_model_bosonic_env( + self, evo, liouvillianize, atol=1e-3 + ): + dlm = DrudeLorentzPureDephasingModel( + lam=0.025, gamma=0.05, T=1/0.95, Nk=2, + ) + ck_real, vk_real, ck_imag, vk_imag = dlm.bath_coefficients() + H_sys = hamiltonian_to_sys(dlm.H, evo, liouvillianize) + + env = ExponentialBosonicEnvironment(ck_real, vk_real, ck_imag, vk_imag) + options = {"nsteps": 15000, "store_states": True} + hsolver = HEOMSolver(H_sys, (env, dlm.Q), 14, options=options) + + tlist = np.linspace(0, 10, 21) + result = hsolver.run(dlm.rho(), tlist) + + test = dlm.state_results(result.states) + expected = dlm.analytic_results(tlist) + np.testing.assert_allclose(test, expected, atol=atol) + + if evo != "qobjevo_timedep": + rho_final, ado_state = hsolver.steady_state() + test = dlm.state_results([rho_final]) + expected = dlm.analytic_results([100]) + np.testing.assert_allclose(test, expected, atol=atol) + assert rho_final == ado_state.extract(0) + else: + assert_raises_steady_state_time_dependent(hsolver) + + @pytest.mark.parametrize(['terminator'], [ + pytest.param(True, id="terminator"), + pytest.param(False, id="noterminator"), + ]) + @pytest.mark.parametrize('approx', ["matsubara", "pade"]) + def test_pure_dephasing_model_drude_lorentz_baths( + self, terminator, approx, atol=1e-3 + ): + dlm = DrudeLorentzPureDephasingModel( + lam=0.025, gamma=0.05, T=1/0.95, Nk=2, + ) + env = DrudeLorentzEnvironment(lam=dlm.lam, gamma=dlm.gamma, T=dlm.T) + if approx == 'matsubara': + approx_env, delta = env.approx_by_matsubara( + Nk=dlm.Nk, compute_delta=True) + else: + approx_env, delta = env.approx_by_pade( + Nk=dlm.Nk, compute_delta=True) + if terminator: + terminator_op = system_terminator(dlm.Q, delta) + H_sys = liouvillian(dlm.H) + terminator_op + else: + H_sys = dlm.H + + options = {"nsteps": 15000, "store_states": True} + hsolver = HEOMSolver(H_sys, (approx_env, dlm.Q), 14, options=options) + + tlist = np.linspace(0, 10, 21) + result = hsolver.run(dlm.rho(), tlist) + + test = dlm.state_results(result.states) + expected = dlm.analytic_results(tlist) + np.testing.assert_allclose(test, expected, atol=atol) + + rho_final, ado_state = hsolver.steady_state() + test = dlm.state_results([rho_final]) + expected = dlm.analytic_results([100]) + np.testing.assert_allclose(test, expected, atol=atol) + assert rho_final == ado_state.extract(0) + + def test_underdamped_pure_dephasing_model_underdamped_bath( + self, atol=1e-3 + ): + udm = UnderdampedPureDephasingModel( + lam=0.1, gamma=0.05, w0=1, T=1/0.95, Nk=2, + ) + env = UnderDampedEnvironment( + lam=udm.lam, T=udm.T, gamma=udm.gamma, w0=udm.w0 + ).approx_by_matsubara(Nk=udm.Nk) + + options = {"nsteps": 15000, "store_states": True} + hsolver = HEOMSolver(udm.H, (env, udm.Q), 14, options=options) + + tlist = np.linspace(0, 10, 21) + result = hsolver.run(udm.rho(), tlist) + + test = udm.state_results(result.states) + expected = udm.analytic_results(tlist) + np.testing.assert_allclose(test, expected, atol=atol) + + rho_final, ado_state = hsolver.steady_state() + test = udm.state_results([rho_final]) + expected = udm.analytic_results([5000]) + np.testing.assert_allclose(test, expected, atol=atol) + assert rho_final == ado_state.extract(0) + + @pytest.mark.parametrize(['evo'], [ + pytest.param("qobj"), + pytest.param("qobjevo_const"), + pytest.param("qobjevo_timedep"), + ]) + @pytest.mark.parametrize(['liouvillianize'], [ + pytest.param(False, id="hamiltonian"), + pytest.param(True, id="liouvillian"), + ]) + def test_discrete_level_model_fermionic_bath(self, evo, liouvillianize): + dlm = DiscreteLevelCurrentModel( + gamma=0.01, W=1, T=0.025851991, lmax=10, + ) + H_sys = hamiltonian_to_sys(dlm.H, evo, liouvillianize) + ck_plus, vk_plus, ck_minus, vk_minus = dlm.bath_coefficients() + + options = { + "store_states": True, + "store_ados": True, + "nsteps": 15_000, + "rtol": 1e-7, + "atol": 1e-7, + } + env = ExponentialFermionicEnvironment( + ck_plus, vk_plus, ck_minus, vk_minus + ) + # for a single impurity we converge with max_depth = 2 + hsolver = HEOMSolver(H_sys, (env, dlm.Q), 2, options=options) + + tlist = [0, 600] + result = hsolver.run(dlm.rho(), tlist) + current = dlm.state_current(result.ado_states[-1]) + analytic_current = dlm.analytic_current() + np.testing.assert_allclose(analytic_current, current, rtol=1e-3) + + if evo != "qobjevo_timedep": + rho_final, ado_state = hsolver.steady_state() + current = dlm.state_current(ado_state) + analytic_current = dlm.analytic_current() + np.testing.assert_allclose(analytic_current, current, rtol=1e-3) + else: + assert_raises_steady_state_time_dependent(hsolver) + + @pytest.mark.parametrize(['approx', 'analytic_current'], [ + pytest.param("matsubara", 0.001101, id="matsubara"), + pytest.param("pade", 0.000813, id="pade"), + ]) + def test_discrete_level_model_lorentzian_baths( + self, approx, analytic_current, atol=1e-3 + ): + dlm = DiscreteLevelCurrentModel( + gamma=0.01, W=1, T=0.025851991, lmax=10, + ) + + options = { + "nsteps": 15_000, "store_states": True, "store_ados": True, + "rtol": 1e-7, "atol": 1e-7, + } + env_l = LorentzianEnvironment( + gamma=dlm.gamma, W=dlm.W, T=dlm.T, mu=dlm.theta / 2 + ) + env_r = LorentzianEnvironment( + gamma=dlm.gamma, W=dlm.W, T=dlm.T, mu=-dlm.theta / 2 + ) + if approx == 'matsubara': + env_l = env_l.approx_by_matsubara(Nk=dlm.lmax) + env_r = env_r.approx_by_matsubara(Nk=dlm.lmax) + else: + env_l = env_l.approx_by_pade(Nk=dlm.lmax) + env_r = env_r.approx_by_pade(Nk=dlm.lmax) + # for a single impurity we converge with max_depth = 2 + hsolver = HEOMSolver( + dlm.H, [(env_r, dlm.Q), (env_l, dlm.Q)], 2, options=options + ) + + tlist = [0, 600] + result = hsolver.run(dlm.rho(), tlist) + current = dlm.state_current(result.ado_states[-1]) + # analytic_current = dlm.analytic_current() + np.testing.assert_allclose(analytic_current, current, rtol=1e-3) + + rho_final, ado_state = hsolver.steady_state() + current = dlm.state_current(ado_state) + # analytic_current = dlm.analytic_current() + np.testing.assert_allclose(analytic_current, current, rtol=1e-3) + + def test_parity(self): + depth = 2 + Nk = 2 + # system: two fermions + N = 2 + d_1 = fdestroy(N, 0) + d_2 = fdestroy(N, 1) + # bath params: + mu = 0. # chemical potential + Gamma = 1 # coupling strenght + W = 2.5 # bath width + # system params: + # coulomb repulsion + U = 3 * np.pi * Gamma + # impurity energy + w0 = - U / 2. + beta = 1 / (0.2 * Gamma) # Inverse Temperature + H = w0 * (d_1.dag() * d_1 + d_2.dag() + * d_2) + U * d_1.dag() * d_1 * d_2.dag() * d_2 + + L = liouvillian(H) + env1 = LorentzianEnvironment( + gamma=2 * Gamma, W=W, mu=mu, T=1 / beta, tag="Lead 1" + ).approx_by_pade(Nk=Nk) + env2 = LorentzianEnvironment( + gamma=2 * Gamma, W=W, mu=mu, T=1 / beta, tag="Lead 2" + ).approx_by_pade(Nk=Nk) + solver = HEOMSolver( + L, [(env1, d_1), (env2, d_2)], depth, odd_parity=True + ) + rhoss, _ = solver.steady_state(use_mkl=False) + rhoss = rhoss.full() + expected_odd = np.diag([-0.18472, 0.68472, 0.68472, -0.18472]) + np.testing.assert_allclose(rhoss, expected_odd, atol=1e-5) + + solver = HEOMSolver( + L, [(env1, d_1), (env2, d_2)], depth, odd_parity=False + ) + rhoss, _ = solver.steady_state(use_mkl=False) + rhoss = rhoss.full() + expected = np.diag([0.10623, 0.39376, 0.39376, 0.10623]) + np.testing.assert_allclose(rhoss, expected, atol=1e-5) + + +class TestHeomsolveFunctionWithEnv: + # Copied from TestHeomsolveFunction but uses "environment" API instead of + # "bath"s + @pytest.mark.parametrize(['evo'], [ + pytest.param("qobj", id="qobj"), + pytest.param("listevo_const", id="listevo_const"), + pytest.param("qobjevo_const", id="qobjevo_const"), + pytest.param("qobjevo_timedep", id="qobjevo_timedep"), + ]) + @pytest.mark.parametrize(['liouvillianize'], [ + pytest.param(False, id="hamiltonian"), + pytest.param(True, id="liouvillian"), + ]) + def test_heomsolve_with_pure_dephasing_model( + self, evo, liouvillianize, atol=1e-3 + ): + dlm = DrudeLorentzPureDephasingModel( + lam=0.025, gamma=0.05, T=1/0.95, Nk=2, + ) + ck_real, vk_real, ck_imag, vk_imag = dlm.bath_coefficients() + H_sys = hamiltonian_to_sys(dlm.H, evo, liouvillianize) + + env = ExponentialBosonicEnvironment(ck_real, vk_real, ck_imag, vk_imag) + options = {"nsteps": 15000, "store_states": True} + + e_ops = { + "11": basis(2, 0) * basis(2, 0).dag(), + "22": basis(2, 1) * basis(2, 1).dag(), + } + + tlist = np.linspace(0, 10, 21) + result = heomsolve( + H_sys, (env, dlm.Q), 14, dlm.rho(), tlist, + e_ops=e_ops, args={"foo": 1}, options=options) + + test = dlm.state_results(result.states) + expected = dlm.analytic_results(tlist) + np.testing.assert_allclose(test, expected, atol=atol) + + for label in ["11", "22"]: + np.testing.assert_allclose( + result.e_data[label], + [expect(rho, e_ops[label]) for rho in result.states], + atol=atol, + ) + + class TestHSolverDL: @pytest.mark.parametrize(['bnd_cut_approx', 'atol'], [ pytest.param(True, 1e-4, id="bnd_cut_approx"), diff --git a/qutip/tests/solver/heom/test_heom.py b/qutip/tests/solver/heom/test_heom.py index c8988a168a..afdb1289da 100644 --- a/qutip/tests/solver/heom/test_heom.py +++ b/qutip/tests/solver/heom/test_heom.py @@ -1,8 +1,7 @@ """ -Tests for qutip.nonmarkov.heom. +Tests for qutip.solver.heom. """ -import pytest from qutip.solver.heom import ( BathExponent, Bath, diff --git a/qutip/tests/solver/test_brmesolve.py b/qutip/tests/solver/test_brmesolve.py index 6e0aeb929d..a0d5677c07 100644 --- a/qutip/tests/solver/test_brmesolve.py +++ b/qutip/tests/solver/test_brmesolve.py @@ -2,6 +2,7 @@ import numpy as np import qutip from qutip.solver.brmesolve import brmesolve +from qutip.core.environment import DrudeLorentzEnvironment def pauli_spin_operators(): @@ -10,7 +11,7 @@ def pauli_spin_operators(): _simple_qubit_gamma = 0.25 coeff = qutip.coefficient(lambda t, w: _simple_qubit_gamma * (w >= 0), - args={'w':0}) + args={'w': 0}) _m_c_op = np.sqrt(_simple_qubit_gamma) * qutip.sigmam() _z_c_op = np.sqrt(_simple_qubit_gamma) * qutip.sigmaz() _x_a_op = [qutip.sigmax(), coeff] @@ -132,8 +133,8 @@ def test_tensor_system(): + w3/2. * qutip.tensor(id2, id2, qutip.sigmaz())) # White noise - S2 = qutip.coefficient(lambda t, w: gamma2, args={'w':0}) - S3 = qutip.coefficient(lambda t, w: gamma3, args={'w':0}) + S2 = qutip.coefficient(lambda t, w: gamma2, args={'w': 0}) + S3 = qutip.coefficient(lambda t, w: gamma3, args={'w': 0}) qubit_2_x = qutip.tensor(id2, qutip.sigmax(), id2) qubit_3_x = qutip.tensor(id2, id2, qutip.sigmax()) @@ -209,10 +210,10 @@ def _string_w_interpolating_t(kappa, times): @pytest.mark.slow @pytest.mark.parametrize("time_dependence_tuple", [ - _mixed_string, - _separate_strings, - _string_w_interpolating_t, - ]) + _mixed_string, + _separate_strings, + _string_w_interpolating_t, +]) def test_time_dependence_tuples(time_dependence_tuple): N = 10 a = qutip.destroy(N) @@ -312,3 +313,27 @@ def test_feedback(): args={"A": qutip.BRSolver.ExpectFeedback(qutip.num(N))} ) assert np.all(result.expect[0] > 4. - tol) + + +@pytest.mark.parametrize("lam,gamma,beta", [(0.05, 1, 1), (0.1, 5, 2)]) +def test_accept_environment(lam, gamma, beta): + DL = ( + "2 * pi * 2.0 * {lam} / (pi * {gamma} * {beta}) if (w==0) " + "else 2 * pi * (2.0 * {lam} * {gamma} * w / (pi * (w**2 + {gamma}**2))) " + "* ((1 / (exp(w * {beta}) - 1)) + 1)" + ).format(gamma=gamma, beta=beta, lam=lam) + H = 0.5 * qutip.sigmax() + psi0 = (2 * qutip.basis(2, 0) + qutip.basis(2, 1)).unit() + times = np.linspace(0, 10, 100) + resultBR_str = brmesolve( + H, psi0, times, + a_ops=[[qutip.sigmaz(), DL]], + e_ops=[qutip.sigmaz()] + ) + env = DrudeLorentzEnvironment(T=1/beta, lam=lam, gamma=gamma) + resultBR_env = brmesolve( + H, psi0, times, + a_ops=[[qutip.sigmaz(), env]], + e_ops=[qutip.sigmaz()] + ) + assert np.allclose(resultBR_env.expect[0], resultBR_str.expect[0]) diff --git a/qutip/tests/test_utilities.py b/qutip/tests/test_utilities.py index f29d8224c6..aaab14dfb2 100644 --- a/qutip/tests/test_utilities.py +++ b/qutip/tests/test_utilities.py @@ -10,8 +10,10 @@ pytest.param(np.log(2)*5, 5, 1, id='5*log(2)'), pytest.param(0, 1, 0, id="0_energy"), pytest.param(1, -1, 0, id="neg_temp"), - pytest.param(np.array([np.log(2), np.log(3), np.log(4)]), 1, - np.array([1, 1/2, 1/3]), id="array"), + pytest.param(np.array([0, np.log(2), np.log(3), np.log(4)]), 1, + np.array([0, 1, 1/2, 1/3]), id="array"), + pytest.param(np.array([0, np.log(2), np.log(3), np.log(4)]), 0, + np.array([0, 0, 0, 0]), id="array_0temp") ]) def test_n_thermal(w, w_th, expected): np.testing.assert_allclose(n_thermal(w, w_th), expected) @@ -119,3 +121,25 @@ def test_cpu_count(monkeypatch): monkeypatch.setenv("QUTIP_NUM_PROCESSES", str(0)) new_ncpus = available_cpu_count() assert new_ncpus >= 1 + + +class TestFitting: + def model(self, x, a, b, c): + return np.real(a * np.exp(-(b + 1j * c) * x)) + + def test_same_model(self): + x = np.linspace(0, 10, 100) + + # Generate data to fit + fparams1 = [1, .5, 0] + fparams2 = [3, 2, .5] + y = self.model(x, *fparams1) + self.model(x, *fparams2) + + rmse, params = utils.iterated_fit( + self.model, num_params=3, xdata=x, ydata=y, + lower=[-np.inf, -np.inf, 0], target_rmse=1e-8, Nmax=2) + + assert rmse < 1e-8 + print(params) + assert (np.all(np.isclose(params, [fparams1, fparams2], atol=1e-3)) or + np.all(np.isclose(params, [fparams2, fparams1], atol=1e-3))) \ No newline at end of file diff --git a/qutip/utilities.py b/qutip/utilities.py index 46a028a5da..db7d9ec4a6 100644 --- a/qutip/utilities.py +++ b/qutip/utilities.py @@ -3,9 +3,17 @@ qutip modules. """ -__all__ = ['n_thermal', 'clebsch', 'convert_unit'] +# Required for Sphinx to follow autodoc_type_aliases +from __future__ import annotations + +__all__ = ['n_thermal', 'clebsch', 'convert_unit', 'iterated_fit'] + +from typing import Callable import numpy as np +from numpy.typing import ArrayLike +from scipy.optimize import curve_fit + def n_thermal(w, w_th): """ @@ -34,14 +42,17 @@ def n_thermal(w, w_th): """ - if isinstance(w, np.ndarray): - return 1.0 / (np.exp(w / w_th) - 1.0) + w = np.array(w, dtype=float) + result = np.zeros_like(w) - else: - if (w_th > 0) and np.exp(w / w_th) != 1.0: - return 1.0 / (np.exp(w / w_th) - 1.0) - else: - return 0.0 + if w_th <= 0: + result[w < 0] = -1 + return result.item() if w.ndim == 0 else result + + non_zero = w != 0 + result[non_zero] = 1 / (np.exp(w[non_zero] / w_th) - 1) + + return result.item() if w.ndim == 0 else result def _factorial_prod(N, arr): @@ -110,8 +121,8 @@ def clebsch(j1, j2, j3, m1, m2, m3): s_factors = np.zeros(((vmax + 1 - vmin), (int(j1 + j2 + j3))), np.int32) # `S` and `C` are large integer,s if `sign` is a np.int32 it could oveflow sign = int((-1) ** (vmin + j2 + m2)) - for i,v in enumerate(range(vmin, vmax + 1)): - factor = s_factors[i,:] + for i, v in enumerate(range(vmin, vmax + 1)): + factor = s_factors[i, :] _factorial_prod(j2 + j3 + m1 - v, factor) _factorial_prod(j1 - m1 + v, factor) _factorial_div(j3 - j1 + j2 - v, factor) @@ -120,7 +131,7 @@ def clebsch(j1, j2, j3, m1, m2, m3): _factorial_div(v, factor) common_denominator = -np.min(s_factors, axis=0) numerators = s_factors + common_denominator - S = sum([(-1)**i * _to_long(vec) for i,vec in enumerate(numerators)]) * \ + S = sum([(-1)**i * _to_long(vec) for i, vec in enumerate(numerators)]) * \ sign / _to_long(common_denominator) return C * S @@ -331,3 +342,184 @@ def _version2int(version_string): "post")[0].split('.') return sum([int(d if len(d) > 0 else 0) * (100 ** (3 - n)) for n, d in enumerate(str_list[:3])]) + + +# ----------------------------------------------------------------------------- +# Fitting utilities +# + +def iterated_fit( + fun: Callable[..., complex], + num_params: int, + xdata: ArrayLike, + ydata: ArrayLike, + target_rmse: float = 1e-5, + Nmin: int = 1, + Nmax: int = 10, + guess: ArrayLike | Callable[[int], ArrayLike] = None, + lower: ArrayLike = None, + upper: ArrayLike = None, + sigma: float | ArrayLike = None, + maxfev: int = None +) -> tuple[float, ArrayLike]: + r""" + Iteratively tries to fit the given data with a model of the form + + .. math:: + y = \sum_{k=1}^N f(x; p_{k,1}, \dots, p_{k,n}) + + where `f` is a model function depending on `n` parameters, and the number + `N` of terms is increased until the normalized rmse (root mean square + error) falls below the target value. + + Parameters + ---------- + fun : callable + The model function. Its first argument is the array `xdata`, its other + arguments are the fitting parameters. + num_params : int + The number of fitting parameters per term (`n` in the equation above). + The function `fun` must take `num_params+1` arguments. + xdata : array_like + The independent variable. + ydata : array_like + The dependent data. + target_rmse : optional, float + Desired normalized root mean squared error (default `1e-5`). + Nmin : optional, int + The minimum number of terms to be used for the fit (default 1). + Nmax : optional, int + The maximum number of terms to be used for the fit (default 10). + If the number `Nmax` of terms is reached, the function returns even if + the target rmse has not been reached yet. + guess : optional, array_like or callable + This can be either a list of length `n`, with the i-th entry being the + guess for the parameter :math:`p_{k,i}` (for all terms :math:`k`), or a + function that provides different initial guesses for each term. + Specifically, given a number `N` of terms, the function returns an + array `[[p11, ..., p1n], [p21, ..., p2n], ..., [pN1, ..., pNn]]` of + initial guesses. + lower : optional, list of length `num_params` + Lower bounds on the parameters for the fit. + upper : optional, list of length `num_params` + Upper bounds on the parameters for the fit. + sigma : optional, float or array_like + The uncertainty in the dependent data, see the documentation of + ``scipy.optimize.curve_fit``. + maxfev : optional, int + The maximum number of function evaluations (per value of ``N``). + + Returns + ------- + rmse : float + The normalized mean squared error of the fit + params : array_like + The model parameters in the form + `[[p11, ..., p1n], [p21, ..., p2n], ..., [pN1, ..., pNn]]`. + """ + + if len(xdata) != len(ydata): + raise ValueError( + "The shape of the provided fit data is not consistent") + + if lower is None: + lower = np.full(num_params, -np.inf) + if upper is None: + upper = np.full(num_params, np.inf) + if not (len(lower) == num_params and len(upper) == num_params): + raise ValueError( + "The shape of the provided fit bounds is not consistent") + + N = Nmin + rmse1 = np.inf + + while rmse1 > target_rmse and N <= Nmax: + if guess is None: + guesses = np.ones((N, num_params), dtype=float) + elif callable(guess): + guesses = np.array(guess(N)) + if guesses.shape != (N, num_params): + raise ValueError( + "The shape of the provided fit guesses is not consistent") + else: + guesses = np.tile(guess, (N, 1)) + + lower_repeat = np.tile(lower, N) + upper_repeat = np.tile(upper, N) + rmse1, params = _fit(fun, num_params, xdata, ydata, + guesses, lower_repeat, + upper_repeat, sigma, maxfev) + N += 1 + + return rmse1, params + + +def _pack(params): + # Pack parameter lists for fitting. + # Input: array of parameters like `[[p11, ..., p1n], ..., [pN1, ..., pNn]]` + # Output: packed parameters like `[p11, ..., p1n, p21, ..., p2n, ...]` + return params.ravel() # like flatten, but doesn't copy data + + +def _unpack(params, num_params): + # Inverse of _pack, `num_params` is "n" + N = len(params) // num_params + return np.reshape(params, (N, num_params)) + + +def _evaluate(fun, x, params): + result = 0 + for term_params in params: + result += fun(x, *term_params) + return result + + +def _rmse(fun, xdata, ydata, params): + """ + The normalized root mean squared error for the fit with the given + parameters. (The closer to zero = the better the fit.) + """ + yhat = _evaluate(fun, xdata, params) + if (yhat == ydata).all(): + return 0 + return ( + np.sqrt(np.mean((yhat - ydata) ** 2) / len(ydata)) + / (np.max(ydata) - np.min(ydata)) + ) + + +def _fit(fun, num_params, xdata, ydata, guesses, lower, upper, sigma, + maxfev, method='trf'): + # fun: model function + # num_params: number of parameters in fun + # xdata, ydata: data to be fit + # N: number of terms + # guesses: initial guesses [[p11, ..., p1n],..., [pN1, ..., pNn]] + # lower, upper: parameter bounds + # sigma: data uncertainty (useful to control when values are small) + # maxfev: how many times the parameters can be altered, lower is faster but + # less accurate + if (upper <= lower).all(): + return _rmse(fun, xdata, ydata, guesses), guesses + + # Depending on the method, scipy uses leastsq or least_squares, and the + # `maxfev` parameter has different names in the two functions + if method == 'lm': + maxfev_arg = {'maxfev': maxfev} + else: + maxfev_arg = {'max_nfev': maxfev} + + # Scipy only supports scalar sigma since 1.12 + if sigma is not None and not hasattr(sigma, "__len__"): + sigma = [sigma] * len(xdata) + + packed_params, _ = curve_fit( + lambda x, *packed_params: _evaluate( + fun, x, _unpack(packed_params, num_params) + ), + xdata, ydata, p0=_pack(guesses), bounds=(lower, upper), + method=method, sigma=sigma, **maxfev_arg + ) + params = _unpack(packed_params, num_params) + rmse = _rmse(fun, xdata, ydata, params) + return rmse, params diff --git a/setup.cfg b/setup.cfg index dfd09c1964..9fdc2261a8 100644 --- a/setup.cfg +++ b/setup.cfg @@ -61,6 +61,7 @@ ipython = extras = loky tqdm + mpmath mpi = mpi4py ; This uses ConfigParser's string interpolation to include all the above