From bd1ff09da38574e20e741d31d8e5c0d39464a2db Mon Sep 17 00:00:00 2001 From: Simon Cross Date: Thu, 2 Dec 2021 20:17:27 +0100 Subject: [PATCH] Replace terminator=True/False argument with a terminator method on the baths that support it. --- qutip/nonmarkov/bofin_baths.py | 160 ++++++++++++++++-------------- qutip/nonmarkov/bofin_solvers.py | 4 +- qutip/tests/test_bofin_baths.py | 41 +++++--- qutip/tests/test_bofin_solvers.py | 4 +- 4 files changed, 117 insertions(+), 92 deletions(-) diff --git a/qutip/nonmarkov/bofin_baths.py b/qutip/nonmarkov/bofin_baths.py index 559ae27547..faded10c63 100644 --- a/qutip/nonmarkov/bofin_baths.py +++ b/qutip/nonmarkov/bofin_baths.py @@ -335,14 +335,6 @@ class DrudeLorentzBath(BosonicBath): Number of exponential terms used to approximate the bath correlation functions. - terminator : bool, default False - Whether to calculate the Matsubara terminator for the bath. If - true, the calculated terminator and the calculated approximation - discrepancy are provided via the ``.terminator`` and ``.delta`` - attributes, respectively. Otherwise, ``.terminator`` and ``.delta`` are - both set to ``None``. The terminator is a Liouvillian term, so it's - dimensions are those of a superoperator of ``Q``. - combine : bool, default True Whether to combine exponents with the same frequency (and coupling operator). See :meth:`BosonicBath.combine` for details. @@ -351,23 +343,9 @@ class DrudeLorentzBath(BosonicBath): 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. - - Attributes - ---------- - terminator : Qobj - The Matsubara 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 : float - The approximation discrepancy. That is, the difference between the true - correlation function of the Drude-Lorentz bath and the sum of the - ``Nk`` exponential terms is approximately ``2 * delta * dirac(t)``, - where ``dirac(t)`` denotes the Dirac delta function. """ def __init__( - self, Q, lam, gamma, T, Nk, terminator=False, combine=True, tag=None, + self, Q, lam, gamma, T, Nk, combine=True, tag=None, ): ck_real, vk_real, ck_imag, vk_imag = self._matsubara_params( lam=lam, @@ -380,32 +358,33 @@ def __init__( Q, ck_real, vk_real, ck_imag, vk_imag, combine=combine, tag=tag, ) - if terminator: - self.delta, self.terminator = DrudeLorentzBath._terminator( - Q=Q, lam=lam, gamma=gamma, T=T, exponents=self.exponents, - ) - else: - self.delta = None - self.terminator = None + self._dl_terminator = _DrudeLorentzTerminator( + Q=Q, lam=lam, gamma=gamma, T=T, + ) - @staticmethod - def _terminator(Q, lam, gamma, T, exponents): - """ Calculate the hierarchy terminator term for the Liouvillian. """ - beta = 1 / T + def terminator(self): + """ + Return the Matsubara terminator for the bath and the calculated + approximation discrepancy. - op = -2*spre(Q)*spost(Q.dag()) + spre(Q.dag()*Q) + spost(Q.dag()*Q) - delta = 2 * lam / (beta * gamma) - 1j * lam + Returns + ------- + delta: float - 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 + The approximation discrepancy. That is, the difference between the + true correlation function of the Drude-Lorentz bath and the sum of + the ``Nk`` exponential terms is approximately ``2 * delta * + dirac(t)``, where ``dirac(t)`` denotes the Dirac delta function. - L_bnd = -delta * op - return delta, L_bnd + terminator : Qobj + + The Matsubara 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. """ @@ -462,14 +441,6 @@ class DrudeLorentzPadeBath(BosonicBath): Number of Padé exponentials terms used to approximate the bath correlation functions. - terminator : bool, default False - Whether to calculate the Padé terminator for the bath. If true, the - calculated terminator and the calculated approximation discrepancy are - provided via the ``.terminator`` and ``.delta`` attributes, - respectively. Otherwise, ``.terminator`` and ``.delta`` are both set to - ``None``. The terminator is a Liouvillian term, so it's dimensions are - those of a superoperator of ``Q``. - combine : bool, default True Whether to combine exponents with the same frequency (and coupling operator). See :meth:`BosonicBath.combine` for details. @@ -478,23 +449,9 @@ class DrudeLorentzPadeBath(BosonicBath): 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. - - Attributes - ---------- - terminator : Qobj - 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 : float - The approximation discrepancy. That is, the difference between the true - correlation function of the Drude bath and the correlation function - represented by this object is approximately ``2 delta * dirac(t)``, - where ``dirac(t)`` denotes the Dirac delta function. """ def __init__( - self, Q, lam, gamma, T, Nk, terminator=False, combine=True, tag=None + self, Q, lam, gamma, T, Nk, combine=True, tag=None ): eta_p, gamma_p = self._corr(lam=lam, gamma=gamma, T=T, Nk=Nk) @@ -509,13 +466,33 @@ def __init__( Q, ck_real, vk_real, ck_imag, vk_imag, combine=combine, tag=tag, ) - if terminator: - self.delta, self.terminator = DrudeLorentzBath._terminator( - Q=Q, lam=lam, gamma=gamma, T=T, exponents=self.exponents, - ) - else: - self.delta = None - self.terminator = None + self._dl_terminator = _DrudeLorentzTerminator( + Q=Q, lam=lam, gamma=gamma, T=T, + ) + + def terminator(self): + """ + Return the Padé terminator for the bath and the calculated + approximation discrepancy. + + Returns + ------- + delta: float + + The approximation discrepancy. That is, the difference between the + true correlation function of the Drude-Lorentz bath and the sum of + the ``Nk`` exponential terms is approximately ``2 * delta * + dirac(t)``, where ``dirac(t)`` denotes the Dirac delta function. + + terminator : Qobj + + 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 _corr(self, lam, gamma, T, Nk): beta = 1. / T @@ -581,6 +558,39 @@ def _calc_chi(self, Nk): 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 + + class UnderDampedBath(BosonicBath): """ A helper class for constructing an under-damped bosonic bath from the diff --git a/qutip/nonmarkov/bofin_solvers.py b/qutip/nonmarkov/bofin_solvers.py index d7d8907467..877e99ccae 100644 --- a/qutip/nonmarkov/bofin_solvers.py +++ b/qutip/nonmarkov/bofin_solvers.py @@ -1037,7 +1037,6 @@ def __init__( Nk=N_exp - 1, T=temperature, combine=combine, - terminator=bnd_cut_approx, ) if bnd_cut_approx: @@ -1049,7 +1048,8 @@ def __init__( is_hamiltonian = H0.type == "oper" if is_hamiltonian: H_sys = liouvillian(H_sys) - H_sys = H_sys + bath.terminator + _, terminator = bath.terminator() + H_sys = H_sys + terminator super().__init__( H_sys, bath=bath, max_depth=N_cut, options=options, diff --git a/qutip/tests/test_bofin_baths.py b/qutip/tests/test_bofin_baths.py index 7dbc861b05..32fe5589ef 100644 --- a/qutip/tests/test_bofin_baths.py +++ b/qutip/tests/test_bofin_baths.py @@ -202,8 +202,6 @@ def test_create(self): ck=ck_real[1], vk=vk_real[1], tag="bath1", ) - assert bath.delta is None - assert bath.terminator is None bath = DrudeLorentzBath( Q=Q, lam=0.025, T=1 / 0.95, Nk=1, gamma=0.05, combine=False, @@ -213,20 +211,21 @@ def test_create(self): check_exponent(exp2, "R", dim=None, Q=Q, ck=ck_real[1], vk=vk_real[1]) check_exponent(exp3, "I", dim=None, Q=Q, ck=ck_imag[0], vk=vk_imag[0]) - bath = DrudeLorentzBath( - Q=Q, lam=0.025, T=1 / 0.95, Nk=1, gamma=0.05, terminator=True, - ) - bath2 = DrudeLorentzBath( - Q=Q, lam=0.025, T=1 / 0.95, Nk=1, gamma=0.05, terminator=True, - combine=False, - ) + @pytest.mark.parametrize(['combine'], [ + pytest.param(True, id="combine"), + pytest.param(False, id="no-combine"), + ]) + def test_terminator(self, combine): + Q = sigmaz() op = -2*spre(Q)*spost(Q.dag()) + spre(Q.dag()*Q) + spost(Q.dag()*Q) - assert np.abs(bath.delta - (0.00031039 / 4.0)) < 1e-8 - assert np.abs(bath2.delta - (0.00031039 / 4.0)) < 1e-8 + bath = DrudeLorentzBath( + Q=Q, lam=0.025, T=1 / 0.95, Nk=1, gamma=0.05, combine=combine, + ) + delta, terminator = bath.terminator() - assert isequal(bath.terminator, - (0.00031039 / 4.0) * op, tol=1e-8) - assert isequal(bath2.terminator, - (0.00031039 / 4.0) * op, tol=1e-8) + assert np.abs(delta - (0.00031039 / 4.0)) < 1e-8 + assert isequal(terminator, - (0.00031039 / 4.0) * op, tol=1e-8) class TestDrudeLorentzPadeBath: @@ -260,6 +259,22 @@ def test_create(self): check_exponent(exp2, "R", dim=None, Q=Q, ck=ck_real[1], vk=vk_real[1]) check_exponent(exp3, "I", dim=None, Q=Q, ck=ck_imag[0], vk=vk_imag[0]) + @pytest.mark.parametrize(['combine'], [ + pytest.param(True, id="combine"), + pytest.param(False, id="no-combine"), + ]) + def test_terminator(self, combine): + Q = sigmaz() + op = -2*spre(Q)*spost(Q.dag()) + spre(Q.dag()*Q) + spost(Q.dag()*Q) + + bath = DrudeLorentzPadeBath( + Q=Q, lam=0.025, T=1 / 0.95, Nk=1, gamma=0.05, combine=combine, + ) + delta, terminator = bath.terminator() + + assert np.abs(delta - (0.0 / 4.0)) < 1e-8 + assert isequal(terminator, - (0.0 / 4.0) * op, tol=1e-8) + class TestUnderDampedBath: def test_create(self): diff --git a/qutip/tests/test_bofin_solvers.py b/qutip/tests/test_bofin_solvers.py index 9252899271..647b3b5d9c 100644 --- a/qutip/tests/test_bofin_solvers.py +++ b/qutip/tests/test_bofin_solvers.py @@ -622,10 +622,10 @@ def test_pure_dephasing_model_drude_lorentz_baths( ) bath = bath_cls( Q=dlm.Q, lam=dlm.lam, gamma=dlm.gamma, T=dlm.T, Nk=dlm.Nk, - terminator=terminator, ) if terminator: - H_sys = liouvillian(dlm.H) + bath.terminator + _, terminator_op = bath.terminator() + H_sys = liouvillian(dlm.H) + terminator_op else: H_sys = dlm.H