HEOM Parity support #2261
HEOM Parity support #2261
Conversation
… created an utils class
|
@mcditoos Thanks more making the PR. Would you mind reverting all the unrelated changes? We don't mind small clean ups, but big clean ups should rather go into separate PRs and it should be clear why they are being made. Once you're ready for this to be reviewed, remove the draft status and ping me on the PR. |
There was a problem hiding this comment.
Thank you for removing all the unnecessary changes. The parity support changes are much clearer now.
I've request some small changes. In addition we will need a test for the new feature and ideally for the feature to be mentioned in the HEOM guide in the documentation (perhaps also with a small example).
Co-authored-by: Simon Cross <hodgestar+github@gmail.com>
Co-authored-by: Simon Cross <hodgestar+github@gmail.com>
qutip/solver/heom/bofin_solvers.py
Outdated
| sign1 = (-1) ** (n_excite + 1 - self.odd_parity) | ||
| sign2 = (-1)**(n_excite_before_m + self.odd_parity) |
There was a problem hiding this comment.
Could you put these lines back where they were so that there is a small block of code for the calculation of each sign instead of interleaving the two calculations?
| n_excite_before_m = sum(he_fermionic_n[:k]) | ||
| sign2 = (-1) ** (n_excite_before_m) | ||
| sign1 = (-1) ** (n_excite + 1 - self.odd_parity) |
qutip/solver/heom/bofin_solvers.py
Outdated
| @@ -1281,6 +1293,7 @@ class _GatherHEOMRHS: | |||
| nhe : int | |||
| The number of ADOs in the hierarchy. | |||
| """ | |||
|
|
|||
There was a problem hiding this comment.
Could you again undo these unrelated changes?
| @@ -54,8 +54,8 @@ def assert_raises_steady_state_time_dependent(hsolver): | |||
| with pytest.raises(ValueError) as err: | |||
| hsolver.steady_state() | |||
| assert str(err.value) == ( | |||
| "A steady state cannot be determined for a time-dependent" | |||
| " system" | |||
| "A steady state cannot be determined for a time-dependent" | |||
There was a problem hiding this comment.
This can stay -- it's a valid indentation fix.
| @@ -232,6 +232,7 @@ class DrudeLorentzPureDephasingModel: | |||
| """ Analytic Drude-Lorentz pure-dephasing model for testing the HEOM | |||
| solver. | |||
| """ | |||
|
|
|||
| @@ -9,7 +9,7 @@ | |||
|
|
|||
| from qutip import ( | |||
| basis, destroy, expect, liouvillian, qeye, sigmax, sigmaz, | |||
| tensor, Qobj, QobjEvo | |||
| tensor, Qobj, QobjEvo, fcreate, fdestroy | |||
There was a problem hiding this comment.
fcreate appears to be unused. Could you remove it?
There was a problem hiding this comment.
Sure, I shouldn't have included it
| resultHEOMPade = HEOMSolver(L, [bath1, bath2], depth, odd_parity=True) | ||
| rhoss, _ = resultHEOMPade.steady_state(use_mkl=False) | ||
| rhoss = rhoss.full() | ||
| expected_odd = np.diag([-0.18472, 0.68472, 0.68472, -0.18472]) | ||
| expected = np.diag([0.10623, 0.39376, 0.39376, 0.10623]) | ||
| assert np.isclose(rhoss, expected_odd,atol=1e-5).all() | ||
| resultHEOMPade = HEOMSolver(L, [bath1, bath2], depth, odd_parity=False) | ||
| rhoss, _ = resultHEOMPade.steady_state(use_mkl=False) | ||
| rhoss = rhoss.full() | ||
| assert np.isclose(rhoss, expected,atol=1e-5).all() |
There was a problem hiding this comment.
| resultHEOMPade = HEOMSolver(L, [bath1, bath2], depth, odd_parity=True) | |
| rhoss, _ = resultHEOMPade.steady_state(use_mkl=False) | |
| rhoss = rhoss.full() | |
| expected_odd = np.diag([-0.18472, 0.68472, 0.68472, -0.18472]) | |
| expected = np.diag([0.10623, 0.39376, 0.39376, 0.10623]) | |
| assert np.isclose(rhoss, expected_odd,atol=1e-5).all() | |
| resultHEOMPade = HEOMSolver(L, [bath1, bath2], depth, odd_parity=False) | |
| rhoss, _ = resultHEOMPade.steady_state(use_mkl=False) | |
| rhoss = rhoss.full() | |
| assert np.isclose(rhoss, expected,atol=1e-5).all() | |
| solver = HEOMSolver(L, [bath1, bath2], 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, [bath1, bath2], 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) |
|
|
||
|
|
||
|
|
There was a problem hiding this comment.
There appear to be too many blank lines here.
| Q=d_2, gamma=2 * Gamma, w=W, mu=mu, T=1 / beta, Nk=Nk, | ||
| tag="Lead 2") | ||
| resultHEOMPade = HEOMSolver(L, [bath1, bath2], depth, odd_parity=True) | ||
| rhoss, _ = resultHEOMPade.steady_state(use_mkl=False) |
There was a problem hiding this comment.
Could you explain why you needed to pass use_mkl=False here? I think we would expect it to automatically do something sensible for users in all cases and we don't appear to have had to use it in any of the other HEOM solver tests yet.
There was a problem hiding this comment.
I used used_mkl=False, because if I don't for the test fails for python 3.10 with no cython (only in that setup it passes in all others), the test fails with Exception: Zero pivot, numerical factorization or iterative refinement problem. I found this old thread #791 suggesting that sometimes mkl solvers give some trouble so I decided to set it to false
| This is normally prohibited in 'normal' electronic systems, where | ||
| charge is conserved, but is useful for the construction of | ||
| electron spectrum (under linear response assumptions). | ||
|
|
There was a problem hiding this comment.
Simon asked us to help distill and clarify this a little bit; how about:
" For fermionic baths only. Default is "False". "Parity" refers to the parity of the initial system state used with the HEOM. An example of an odd parity state is one made from applying an odd number of fermionic creation operators to a physical density operator. Physical systems have even parity, but allowing the generalization to odd-parity states allows one to calculate useful physical quantities like the system power spectrum or density of states. The form of the HEOM differs depending on the parity of the initial system state, so if this option is set to "True", a different RHS is constructed, which can then be used with a system state of odd parity. "
Description
This PR adds odd parity support to HEOM's fermionic baths, the user can choose whether to use odd or even parity according to his needs when construncting the right hand side of HEOM's equations
Related issues or PRs
Heom tutorial for this feature qutip/qutip-tutorials#74