correlation_2op_1t for harmonic oscillator ground state?
272 views
Skip to first unread message
gste...@gmail.com
Oct 7, 2019, 7:05:27 PM10/7/19
to QuTiP: Quantum Toolbox in Python
Hi all,
Maybe a bit of a noob question
It seems that when I use correlation_2op_1t() to calculate the two time correlation function for the ground state of the harmonic oscillator, I get a real valued oscillation signal (with no imaginary part)
Conventional knowledge that this should be a negative frequency, complex-valued exponential (representing the ability of the ground state to only absorb, not emit, photons...)
Am I missing something?
Thanks,
Gary
Keith Fratus
Nov 6, 2019, 11:30:12 AM11/6/19
to QuTiP: Quantum Toolbox in Python
Hi Gary,
I've recently been having the same problem, and so I came here to see if anyone else was seeing this effect. I think I might know what the problem is. Long story short, from my very quick look through the source code, I think it's possible that correlation_2op_1t() is calling mesolve() with a non-Hermitian initial density matrix, and mesolve() is then (wrongly) making it Hermitian somehow, thus losing the imaginary part of the correlation function.
For the model systems I was studying, I wrote my own simple solver in QuTIP that computes the two-time correlation functions "by hand," using a very simple time evolution of the density matrix, along with the Quantum Regression Theorem. Via the QRT, to calculate <A(t+T)B(t)>, one first time evolves the density matrix to get p(t), then does a second time evolution for time T with R = Bp(t) as the new initial density matrix (which may not be Hermitian), and then computes the final expectation value of A. When I do this, I successfully recover the imaginary part of the correlation function that I expect to be there.
However, correlation_2op_1t() always returns zero for the imaginary part. To try to reproduce this effect, I took my own QRT solver, but before time-evolving R = Bp(t) in the second step, I modified R -> 0.5*(R + R+), to get a Hermitian density matrix. This then recovered the exact problem I was having with correlation_2op_1t().
As a further experiment, I tried calling mesolve directly with an explicitly non-Hermitian initial density matrix, and observed that after the time evolution, it returned a Hermitian matrix to me, even though action under the Liouville operator I was using should have preserved the non-Hermitian nature of the density matrix.
So, I think correlation_2op_1t() is calling mesolve and passing a non-Hermitian density matrix to it, which is accidentally getting converted to a Hermitian matrix (perhaps only the upper half triangle of the matrix is being used, or something like that? It would of course "normally" make sense to assume the initial density matrix is Hermitian when doing time evolution).
Perhaps someone more knowledgeable than me can confirm this? I'm also quite new to QuTIP and the study of open quantum systems in general, so it's possible I've made a mistake somewhere. I also haven't looked through the source code too deeply.
- Keith
I've recently been having the same problem, and so I came here to see if anyone else was seeing this effect. I think I might know what the problem is. Long story short, from my very quick look through the source code, I think it's possible that correlation_2op_1t() is calling mesolve() with a non-Hermitian initial density matrix, and mesolve() is then (wrongly) making it Hermitian somehow, thus losing the imaginary part of the correlation function.
For the model systems I was studying, I wrote my own simple solver in QuTIP that computes the two-time correlation functions "by hand," using a very simple time evolution of the density matrix, along with the Quantum Regression Theorem. Via the QRT, to calculate <A(t+T)B(t)>, one first time evolves the density matrix to get p(t), then does a second time evolution for time T with R = Bp(t) as the new initial density matrix (which may not be Hermitian), and then computes the final expectation value of A. When I do this, I successfully recover the imaginary part of the correlation function that I expect to be there.
However, correlation_2op_1t() always returns zero for the imaginary part. To try to reproduce this effect, I took my own QRT solver, but before time-evolving R = Bp(t) in the second step, I modified R -> 0.5*(R + R+), to get a Hermitian density matrix. This then recovered the exact problem I was having with correlation_2op_1t().
As a further experiment, I tried calling mesolve directly with an explicitly non-Hermitian initial density matrix, and observed that after the time evolution, it returned a Hermitian matrix to me, even though action under the Liouville operator I was using should have preserved the non-Hermitian nature of the density matrix.
So, I think correlation_2op_1t() is calling mesolve and passing a non-Hermitian density matrix to it, which is accidentally getting converted to a Hermitian matrix (perhaps only the upper half triangle of the matrix is being used, or something like that? It would of course "normally" make sense to assume the initial density matrix is Hermitian when doing time evolution).
Perhaps someone more knowledgeable than me can confirm this? I'm also quite new to QuTIP and the study of open quantum systems in general, so it's possible I've made a mistake somewhere. I also haven't looked through the source code too deeply.
- Keith
Andrew M. C. Dawes
Nov 6, 2019, 1:56:42 PM11/6/19
to qu...@googlegroups.com
You can specify the solver that is used. Gary, changing your correlation line to:
corr = correlation_2op_1t(H, psi0, tlist, c_ops, x, x, solver="es")
which uses the Exponential Series solver and gives the results you expect (complex-valued oscillation and negative freq).
There are many options for the solvers, this question likely indicates a need to document some of assumptions that are being made within the solvers themselves.
Best,
Andy
.jpeg?part=0.2&view=1)
--
You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/766351fd-1f79-4a71-801b-3eaf0062e316%40googlegroups.com.
Gary Steele
Nov 6, 2019, 3:25:45 PM11/6/19
to qu...@googlegroups.com
Thanks Andrew!
Curiously enough, when I make the change to solver='es', I get consistently a dead kernel....on the terminal from which I start the notebook server, I get the following:
python(44025,0x106f7d5c0) malloc: Incorrect checksum for freed object 0x7f94bad7db90: probably modified after being freed.
Corrupt value: 0x3ff0000000000000
python(44025,0x106f7d5c0) malloc: *** set a breakpoint in malloc_error_break to debug
Sounds like a underlying library having some memory allocation issues...or perhaps some qutip-generated cython? (I'm not so familiar with Cython, I don't know if you need to do manual memory allocation...?)
Curious: I wonder why you also don't see the same? It could depend on what platform you are running on: I'm running MacOS 10.14.6. Maybe something about how strict the OS is on memory allocation...? In any case, here is my qutip version info:
| Software | Version |
|---|---|
| QuTiP | 4.4.1 |
| Numpy | 1.16.4 |
| SciPy | 1.3.1 |
| matplotlib | 3.1.1 |
| Cython | 0.29.13 |
| Number of CPUs | 4 |
| BLAS Info | INTEL MKL |
| IPython | 7.8.0 |
| Python | 3.7.4 (default, Aug 13 2019, 15:17:50) [Clang 4.0.1 (tags/RELEASE_401/final)] |
| OS | posix [darwin] |
| Wed Nov 06 21:11:12 2019 CET | |
Does anybody else run into the malloc() fatal error when using the 'es' solver?
As it happens, I have also access to a Jupyterhub server that is running linux. On this server the code runs fine with no errors...so clearly a platform-dependent problem:
| Software | Version |
|---|---|
| QuTiP | 4.4.1 |
| Numpy | 1.17.2 |
| SciPy | 1.3.1 |
| matplotlib | 3.1.1 |
| Cython | 0.29.13 |
| Number of CPUs | 15 |
| BLAS Info | Generic |
| IPython | 7.8.0 |
| Python | 3.7.3 | packaged by conda-forge | (default, Jul 1 2019, 21:52:21) [GCC 7.3.0] |
| OS | posix [linux] |
| Wed Nov 06 20:22:23 2019 UTC | |
Interestingly, the same version of Cython...
Does anyone else have a similar problem with solver='es'? Does it work on darwin (macos) for anyone? Or anyone know where I can chase the problem down further?
Thanks!
Gary
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/CAJBFORDnA7LzP8TNkyd0O2bU3gdg6b8hFjTxREzRNhfB2wWXXA%40mail.gmail.com.
Prof. Gary Steele
Antoni van Leeuwenhoek Professor
Quantum Nanoscience Department, Room D103
Kavli Institute of Nanoscience
Delft University of Technology
Lorentzweg 1, 2628 CJ Delft
The Netherlands
Office: +31-15-278-3402
Mobile: +31-64-667-4938
Email: g.a.s...@tudelft.nl
Website: http://steelelab.tudelft.nl
Quantum Nanoscience Department, Room D103
Kavli Institute of Nanoscience
Delft University of Technology
Lorentzweg 1, 2628 CJ Delft
The Netherlands
Office: +31-15-278-3402
Mobile: +31-64-667-4938
Email: g.a.s...@tudelft.nl
Website: http://steelelab.tudelft.nl
Gary Steele
Nov 6, 2019, 3:34:15 PM11/6/19
to qu...@googlegroups.com
And also maybe a point:
Does anyone know why 'mesolve' gives the wrong answer?
Note that we are calculating the expectation value of a non-hermitian operator when you calculate the two-time correlation function. In particular, in the schroedinger picture:
<x(0)x(t)> = <psi| x e^{iHt} x e^{iHt} | psi>
The unless x happens to commute with H, then the operator x e^{iHt} x e^{iHt} is NOT HERMITIAN
(and, also, curiously, the two-time correlator is not a physical observable...)
So I wonder if this has something to do with the earlier slightly cryptic question about the expectation value of non-hermitian operators being real-valued when using mesolve?
(while one would expect, generally, that the expectation value of non-hermitian operators could and sometimes should have an imaginary component)
And also:
I don't think it makes sense for it even to be possible to use the 2-time correlation function with the me solver if it gives sometimes, and even for the most basic problem one could imagine, the physically incorrect answer...
ie. I see it more as a bug in the code than in the documentation of the solvers
But I am happy to try to help fix it! I should have a bit more free time in a few weeks.
Cheers,
Gary
Andrew M. C. Dawes
Nov 6, 2019, 4:11:42 PM11/6/19
to qu...@googlegroups.com
I'm running on a mac, os 10.13.6 and slightly older versions of everything:
QuTiP Version: 4.4.0
Numpy Version: 1.14.2
Scipy Version: 1.3.0
Cython Version: 0.28.5
Matplotlib Version: 3.0.3
Python Version: 3.6.7
Number of CPUs: 2
BLAS Info: INTEL MKL
OPENMP Installed: False
INTEL MKL Ext: True
Platform Info: Darwin (x86_64)
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/CANZnjsWugEY5mhpW7%3DYpzjhmeFkRZS644wXgHuk6Ba_kGqiTOA%40mail.gmail.com.
Keith Fratus
Nov 7, 2019, 3:47:14 AM11/7/19
to QuTiP: Quantum Toolbox in Python
Hi Gary and Andrew,
Thanks for your thoughts and input on this! I used correlation_2op_1t with solver="es", and it works fine for me, reproducing the correct physics. I'm on Linux Mint, with qutip version 4.4.1, numpy version 1.17.1, and Cython version 0.29.13.
Gary: I actually don't have any issues computing the expectation values of non-Hermitian operators. I just tested that I can compute the expectation values of things like Sx*Sy for a qubit in the spin-up state, and I get a value which is (correctly) imaginary. In fact, I even discovered that I can pass an explicitly non-Hermitian density matrix p to the expect() method, and it still correctly returns Tr[pA] for expect(A, p), for some other operator A.
I think the issue is that internally, correlation_2op_1t is constructing a density matrix which is non-Hermitian, and then it's passing that non-Hermitian density matrix as an initial state to a time-evolution solver. By default it seems to use mesolve, which is not equipped to handle such a case. This seems to be supported by the source code, and what the user manual itself states:
"We therefore first calculate ρ(t)=V(t,0){ρ(0)} using one of the QuTiP evolution solvers with ρ(0) as initial state, and then again use the same solver to calculate V(t+τ,t){Bρ(t)} using Bρ(t) as initial state."
Of course it makes sense to me that mesolve() might not be automatically equipped to handle this - to be fair, a non-Hermitian initial density matrix doesn't "usually" make sense. But it does seem to be what correlation_2op_1t is doing internally by default.
I think I agree with Gary that this seems like it should maybe be classified as a proper bug - if a user treats correlation_2op_1t as a black box, and passes it a very simple problem (two Hermitian operators under the unitary time evolution of a Hamiltonian), one shouldn't expect the method to return incorrect physics by default. Presumably the fix is straightforward though - just change the default internal solver to always be essolve (since it seems to me that mesolve() should quite generally fail at this problem). Or, perhaps it could be possible to modify mesolve() so that it can handle the time-evolution of non-Hermitian density matrices (perhaps there is an assumption in the solver that could be relaxed?)
In any event, thank you both for your help!
- Keith
Thanks for your thoughts and input on this! I used correlation_2op_1t with solver="es", and it works fine for me, reproducing the correct physics. I'm on Linux Mint, with qutip version 4.4.1, numpy version 1.17.1, and Cython version 0.29.13.
Gary: I actually don't have any issues computing the expectation values of non-Hermitian operators. I just tested that I can compute the expectation values of things like Sx*Sy for a qubit in the spin-up state, and I get a value which is (correctly) imaginary. In fact, I even discovered that I can pass an explicitly non-Hermitian density matrix p to the expect() method, and it still correctly returns Tr[pA] for expect(A, p), for some other operator A.
I think the issue is that internally, correlation_2op_1t is constructing a density matrix which is non-Hermitian, and then it's passing that non-Hermitian density matrix as an initial state to a time-evolution solver. By default it seems to use mesolve, which is not equipped to handle such a case. This seems to be supported by the source code, and what the user manual itself states:
"We therefore first calculate ρ(t)=V(t,0){ρ(0)} using one of the QuTiP evolution solvers with ρ(0) as initial state, and then again use the same solver to calculate V(t+τ,t){Bρ(t)} using Bρ(t) as initial state."
Of course it makes sense to me that mesolve() might not be automatically equipped to handle this - to be fair, a non-Hermitian initial density matrix doesn't "usually" make sense. But it does seem to be what correlation_2op_1t is doing internally by default.
I think I agree with Gary that this seems like it should maybe be classified as a proper bug - if a user treats correlation_2op_1t as a black box, and passes it a very simple problem (two Hermitian operators under the unitary time evolution of a Hamiltonian), one shouldn't expect the method to return incorrect physics by default. Presumably the fix is straightforward though - just change the default internal solver to always be essolve (since it seems to me that mesolve() should quite generally fail at this problem). Or, perhaps it could be possible to modify mesolve() so that it can handle the time-evolution of non-Hermitian density matrices (perhaps there is an assumption in the solver that could be relaxed?)
In any event, thank you both for your help!
- Keith
To unsubscribe from this group and stop receiving emails from it, send an email to qu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/766351fd-1f79-4a71-801b-3eaf0062e316%40googlegroups.com.
--
You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qu...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/CAJBFORDnA7LzP8TNkyd0O2bU3gdg6b8hFjTxREzRNhfB2wWXXA%40mail.gmail.com.
--Prof. Gary Steele
Antoni van Leeuwenhoek Professor
Quantum Nanoscience Department, Room D103
Kavli Institute of Nanoscience
Delft University of Technology
Lorentzweg 1, 2628 CJ Delft
The Netherlands
Office: +31-15-278-3402
Mobile: +31-64-667-4938
--Prof. Gary Steele
Antoni van Leeuwenhoek Professor
Quantum Nanoscience Department, Room D103
Kavli Institute of Nanoscience
Delft University of Technology
Lorentzweg 1, 2628 CJ Delft
The Netherlands
Office: +31-15-278-3402
Mobile: +31-64-667-4938
Andrew M. C. Dawes
Nov 7, 2019, 4:55:41 PM11/7/19
to qu...@googlegroups.com
To add some more information to this situation, I do now see the malloc error after upgrading to python 3.7:
python(46653,0x7fffad7f0380) malloc: *** error for object 0x7fdbdec51610: incorrect checksum for freed object - object was probably modified after being freed.
*** set a breakpoint in malloc_error_break to debug
The "es" solver works fine in py36 on my mac.
This, at the very least, seems like a candidate for a bug report. Anyone have insight on to where to poke around for this?
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/2be40193-2e5f-4089-b565-94c89358f3be%40googlegroups.com.
Nathan Shammah
Nov 7, 2019, 9:39:40 PM11/7/19
to qu...@googlegroups.com
Dear Andrew,
Thank you, opening an issue may have this tracked more easily.
Bests,
Nathan
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/CAJBFORBjoP3faNGq-qb2n0-kjNCU5X2vTU7HkwyMpnti0aqZfQ%40mail.gmail.com.
keith...@gmail.com
Nov 8, 2019, 10:17:13 AM11/8/19
to QuTiP: Quantum Toolbox in Python
Hi Gary,
I've recently been having the same problem, and so I came here to see if anyone else was seeing this effect. I think I might know what the problem is. Long story short, from my very quick look through the source code, I think it's possible that correlation_2op_1t() is calling mesolve() with a non-Hermitian initial density matrix, and mesolve() is then (wrongly) making it Hermitian somehow, thus losing the imaginary part of the correlation function.
For the model systems I was studying, I wrote my own simple solver in QuTIP that computes the two-time correlation functions "by hand," using a very simple time evolution of the density matrix, along with the Quantum Regression Theorem. Via the QRT, to calculate <A(t+T)B(t)>, one first time evolves the density matrix to get p(t), then does a second time evolution for time T with R = Bp(t) as the new initial density matrix (which may not be Hermitian), and then computes the final expectation value of A. When I do this, I successfully recover the imaginary part of the correlation function that I expect to be there.
However, correlation_2op_1t() always returns zero for the imaginary part. To try to reproduce this effect, I took my own QRT solver, but before time-evolving R = Bp(t) in the second step, I modified R -> 0.5*(R + R+), to get a Hermitian density matrix. This then recovered the exact problem I was having with correlation_2op_1t().
As a further experiment, I tried calling mesolve directly with an explicitly non-Hermitian initial density matrix, and observed that after the time evolution, it returned a Hermitian matrix to me, even though action under the Liouville operator I was using should have preserved the non-Hermitian nature of the density matrix.
So, I think correlation_2op_1t() is calling mesolve and passing a non-Hermitian density matrix to it, which is accidentally getting converted to a Hermitian matrix (perhaps only the upper half triangle of the matrix is being used, or something like that? It would of course "normally" make sense to assume the initial density matrix is Hermitian when doing time evolution).
Perhaps someone more knowledgeable than me can confirm this? I'm also quite new to QuTIP and the study of open quantum systems in general, so it's possible I've made a mistake somewhere. I also haven't looked through the source code too deeply.
- Keith
On Tuesday, October 8, 2019 at 1:05:27 AM UTC+2, gste...@gmail.com wrote:
I've recently been having the same problem, and so I came here to see if anyone else was seeing this effect. I think I might know what the problem is. Long story short, from my very quick look through the source code, I think it's possible that correlation_2op_1t() is calling mesolve() with a non-Hermitian initial density matrix, and mesolve() is then (wrongly) making it Hermitian somehow, thus losing the imaginary part of the correlation function.
For the model systems I was studying, I wrote my own simple solver in QuTIP that computes the two-time correlation functions "by hand," using a very simple time evolution of the density matrix, along with the Quantum Regression Theorem. Via the QRT, to calculate <A(t+T)B(t)>, one first time evolves the density matrix to get p(t), then does a second time evolution for time T with R = Bp(t) as the new initial density matrix (which may not be Hermitian), and then computes the final expectation value of A. When I do this, I successfully recover the imaginary part of the correlation function that I expect to be there.
However, correlation_2op_1t() always returns zero for the imaginary part. To try to reproduce this effect, I took my own QRT solver, but before time-evolving R = Bp(t) in the second step, I modified R -> 0.5*(R + R+), to get a Hermitian density matrix. This then recovered the exact problem I was having with correlation_2op_1t().
As a further experiment, I tried calling mesolve directly with an explicitly non-Hermitian initial density matrix, and observed that after the time evolution, it returned a Hermitian matrix to me, even though action under the Liouville operator I was using should have preserved the non-Hermitian nature of the density matrix.
So, I think correlation_2op_1t() is calling mesolve and passing a non-Hermitian density matrix to it, which is accidentally getting converted to a Hermitian matrix (perhaps only the upper half triangle of the matrix is being used, or something like that? It would of course "normally" make sense to assume the initial density matrix is Hermitian when doing time evolution).
Perhaps someone more knowledgeable than me can confirm this? I'm also quite new to QuTIP and the study of open quantum systems in general, so it's possible I've made a mistake somewhere. I also haven't looked through the source code too deeply.
- Keith
On Tuesday, October 8, 2019 at 1:05:27 AM UTC+2, gste...@gmail.com wrote:
Andrew M. C. Dawes
Nov 8, 2019, 12:47:47 PM11/8/19
to qu...@googlegroups.com
Nathan,
I have added some comments to issue #1120 as this appears to be related. I suspect a larger problem that is new with py3.7 as that seems to be the primary distinction for when I see the issue or not.
Best,
Andy
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/CAKcZDfQrx4MMpntkxaTMWsbXwXU-KwCh9vovP8YnhBGGzDHA%3Dg%40mail.gmail.com.
Reply all
Reply to author
Forward
0 new messages