Here is an example as simple as I could make it:
Everything is just (normalized) identities. The Hamiltonian is too (which is stupid for an actual simulation), but it's completely irrelevant since I'm only looking at dt=0. I get the same results below for any H.
Here is my code implementing single steps of TDVP on this state with dt=0:
Code: Select all
import tenpy
import numpy as np
p=2 # phys dim
e=2 # external dim (replace with 1 and set d=2 and this works)
d=4 # bond dim
delta_t = 0.0
N_steps = 1
sites = [tenpy.networks.site.Site(tenpy.linalg.charges.LegCharge.from_trivial(p)),
tenpy.networks.site.Site(tenpy.linalg.charges.LegCharge.from_trivial(p))]
H_MPO = tenpy.networks.mpo.MPO(sites,
[tenpy.linalg.np_conserved.Array.from_ndarray_trivial(np.expand_dims(np.expand_dims(np.eye(p),0),0), labels = ['wL','wR','p','p*']),
tenpy.linalg.np_conserved.Array.from_ndarray_trivial(np.expand_dims(np.expand_dims(np.eye(p),0),0), labels = ['wR','wL','p','p*']),
],
IdL = 0,
IdR = 2, bc = 'segment')
MPOModel = tenpy.models.model.MPOModel(tenpy.models.lattice.TrivialLattice(sites),H_MPO)
psi = tenpy.networks.mps.MPS([tenpy.networks.site.Site(tenpy.linalg.charges.LegCharge.from_trivial(p))]*2,
[tenpy.linalg.np_conserved.Array.from_ndarray_trivial(np.eye(d).reshape((e,p,d)), labels = ['vL','p','vR']),
tenpy.linalg.np_conserved.Array.from_ndarray_trivial(np.eye(d).reshape((d,p,e)), labels = ['vL','p','vR'])],
[np.ones(e), np.ones(d)/np.sqrt(d), np.ones(e)],
bc = 'segment',
form = 'G')
print(f'psi_0 = \n{psi.get_theta(0,2)}')
#Initialize TDVP
tdvp_params = {
'start_time': 0,
'dt': delta_t,
}
LP = tenpy.linalg.np_conserved.Array.from_ndarray_trivial(np.expand_dims(np.eye(e),1), labels = ['vR','wR','vR*'])
RP = tenpy.linalg.np_conserved.Array.from_ndarray_trivial(np.expand_dims(np.eye(e),1), labels = ['vL','wL','vL*'])
env = tenpy.networks.mpo.MPOEnvironment(psi, H_MPO, psi, init_LP=LP, init_RP=RP)
tdvp_engine = tenpy.algorithms.tdvp.TDVPEngine(psi, MPOModel, tdvp_params, env )
print('\nRunning TDVP:')
for i in range(4):
tdvp_engine.run_one_site(N_steps=1)
print(f'psi_{i+1} = \n{psi.get_theta(0,2)}')
And here are the outputs:
Code: Select all
psi_0 =
<npc.Array shape=(2, 2, 2, 2) labels=['vL', 'p0', 'p1', 'vR']
charge=ChargeInfo([], [])
+1  +1  +1  +1
0 []0 []0 []0 []
2 2 2 2
[[[[0.5 0. ]
[0. 0. ]]
[[0. 0.5]
[0. 0. ]]]
[[[0. 0. ]
[0.5 0. ]]
[[0. 0. ]
[0. 0.5]]]]
>
Running TDVP:
psi_1 =
<npc.Array shape=(2, 2, 2, 2) labels=['vL', 'p0', 'p1', 'vR']
charge=ChargeInfo([], [])
+1  +1  +1  1
0 []0 []0 []0 []
2 2 2 2
[[[[ 0. +0.j 0. +0.j]
[0.5+0.j 0. +0.j]]
[[ 0. +0.j 0. +0.j]
[ 0. +0.j 0.5+0.j]]]
[[[0.5+0.j 0. +0.j]
[ 0. +0.j 0. +0.j]]
[[ 0. +0.j 0.5+0.j]
[ 0. +0.j 0. +0.j]]]]
>
psi_2 =
<npc.Array shape=(2, 2, 2, 2) labels=['vL', 'p0', 'p1', 'vR']
charge=ChargeInfo([], [])
+1  +1  +1  1
0 []0 []0 []0 []
2 2 2 2
[[[[0. +0.j 0. +0.j]
[0. +0.j 0.5+0.j]]
[[0. +0.j 0. +0.j]
[0.5+0.j 0. +0.j]]]
[[[0. +0.j 0.5+0.j]
[0. +0.j 0. +0.j]]
[[0.5+0.j 0. +0.j]
[0. +0.j 0. +0.j]]]]
>
psi_3 =
<npc.Array shape=(2, 2, 2, 2) labels=['vL', 'p0', 'p1', 'vR']
charge=ChargeInfo([], [])
+1  +1  +1  1
0 []0 []0 []0 []
2 2 2 2
[[[[ 0. +0.j 0. +0.j]
[0.5+0.j 0. +0.j]]
[[ 0. +0.j 0. +0.j]
[ 0. +0.j 0.5+0.j]]]
[[[0.5+0.j 0. +0.j]
[ 0. +0.j 0. +0.j]]
[[ 0. +0.j 0.5+0.j]
[ 0. +0.j 0. +0.j]]]]
>
psi_4 =
<npc.Array shape=(2, 2, 2, 2) labels=['vL', 'p0', 'p1', 'vR']
charge=ChargeInfo([], [])
+1  +1  +1  1
0 []0 []0 []0 []
2 2 2 2
[[[[0. +0.j 0. +0.j]
[0.5+0.j 0. +0.j]]
[[0. +0.j 0. +0.j]
[0. +0.j 0.5+0.j]]]
[[[0.5+0.j 0. +0.j]
[0. +0.j 0. +0.j]]
[[0. +0.j 0.5+0.j]
[0. +0.j 0. +0.j]]]]
>
So after one dt=0 timestep, the state has become (I think)

which is obviously not the original state.
I don't think this behaviour is just because of the noninjectivity of my example, or because of the stupidity of an identity Hamiltonian: An even worse scrambling happens when my initial state is a random normalized MPS, and my initial Hamiltonian is a random Hermitian MPO, so long as the exterior bond dimension is greater than one.
What am I missing?