6. Time Evolution¶
Time evolutions in quimb
are handled by the class Evolution
, which is initialized with a starting state and hamiltonian.
6.1. Basic Usage¶
Set up the Evolution
object with a initial state and hamiltonian.
from quimb import *
p0 = rand_ket(2**10)
h = ham_heis(10, sparse=True)
evo = Evolution(p0, h)
Update it in a single shot to a new time and get the state,
evo.update_to(1)
evo.pt
[[0.006074+0.011257j]
[0.028495+0.016557j]
[0.01617 -0.015705j]
...
[0.018271+0.026071j]
[0.018582-0.016482j]
[0.027511+0.048377j]]
Lazily generate the state at multiple times:
for pt in evo.at_times([2, 3, 4]):
print(expec(pt, p0))
0.003244693043705576
0.008360765055266119
0.009413166876677937
6.2. Methods of Updating¶
There are three methods of updating the state:
Evolution(..., method='integrate')
: use definite integration. Get system at each time step, only need action of Hamiltonian on state. Generally efficient. For pure and mixed states. The additional optionint_small_step={False, True}
determines whether a low or high order adaptive stepping scheme is used, giving naturally smaller or larger times steps. Seescipy.integrate.ode
for details,False
corresponds to"dop853"
,True
to"dopri5"
.
Evolution(..., method='solve')
. Diagonalize the hamiltonian, which once done, allows quickly updating to arbitrary times. Supports pure and mixed states, recomended for small systems.
Evolution(..., method='expm')
: compute the evolved state using the action of the matrix exponential in a ‘single shot’ style. Only needs action of Hamiltonian, for very large systems can use distributed MPI. Only for pure states.
6.3. Computing on the fly¶
Sometimes, if integrating, it is best to just query the state at time-steps chosen dynamically by the adaptive scheme. This is achieved using the compute
keyword supplied to Evolution
. It can also just be a convenient way to set up calculations as well:
p0 = rand_product_state(10)
h = ham_heis(10, sparse=True)
dims = [2] * 10
sysa, sysb = (0, 1), (2, 3)
def calc_t_and_logneg(t, pt):
ln = logneg_subsys(pt, dims, sysa, sysb)
return t, ln
evo = Evolution(p0, h, compute=calc_t_and_logneg, progbar=True)
evo.update_to(1)
ts, lns = zip(*evo.results)
100%|##########| 100/100 [00:00<00:00, 20144.58%/s]
ts
(0.0,
0.2494153162899183,
0.4631400312076022,
0.6768647461252861,
0.9475557600771725,
1.0)
lns
(0.0,
0.2589846231770502,
0.4423181326621942,
0.592697104829058,
0.7379099150663264,
0.7601643146969251)
If a dict of callables is supplied to compute
, (each should take two arguments, the time, and the state, as above), Evolution.results
will itself be a dictionary containing the results of each function at each time step, under the respective key. This can be more convenient:
def calc_t(t, _):
return t
def calc_logneg(_, pt):
return logneg_subsys(pt, [2]*10, 0, 1)
evo = Evolution(p0, h, compute={'t': calc_t, 'ln': calc_logneg}, progbar=True)
evo.update_to(1)
100%|##########| 100/100 [00:00<00:00, 18631.41%/s]
evo.results
{'t': [0.0,
0.2494153162899183,
0.4631400312076022,
0.6768647461252861,
0.9475557600771725,
1.0],
'ln': [0.0,
0.150158101589196,
0.23732495484813082,
0.2922696608578911,
0.3242059171688551,
0.3263498611123527]}
6.4. Time-Dependent Evolutions¶
If you are using method='integrate'
you can supply a callable to ham
to evolve the state with a time dependent Hamiltonian. It should take a single argument t
and return the Hamiltonian at the time. It probably makes sense to use a custom class here to avoid reconstructing as much of the Hamiltonian as possible at each step.
Here we’ll evolve the Neel state:
with the Hamiltonian:
class MyTimeDepIsingHam:
def __init__(self, L):
self.h_interaction = ham_ising(L, sparse=True, jz=1.0, bx=0.0, cyclic=False)
self.h_field = ham_ising(L, sparse=True, jz=0.0, bx=1.0, cyclic=False)
def __call__(self, t):
return self.h_interaction + cos(t) * self.h_field
L = 16
# our initial state
psi0 = neel_state(L)
# instantiate the ham object, it's __call__ method will be used by Evolution
fn_ham_t = MyTimeDepIsingHam(L)
We still want to compute some properties during the evolution:
compute = {
'time': lambda t, p: t,
'entropy': lambda t, p: entropy_subsys(p, dims=[2] * L, sysa=range(L // 2))
}
Now we set up the evolution object again:
evo = Evolution(psi0, fn_ham_t, progbar=True, compute=compute)
evo.update_to(10)
100%|##########| 100/100 [00:07<00:00, 12.53%/s]
We can plot the half chain entropy that we computed on the fly:
%matplotlib inline
import matplotlib.pyplot as plt
with plt.style.context(NEUTRAL_STYLE):
plt.plot(evo.results['time'], evo.results['entropy'])
Or we can use the final state:
fidelity(psi0, evo.pt)
0.003302180752068547