## Asimmetry in trivial 2D expansion using TDVP

How do I use this algorithm? What does that parameter do?
tanahy
Posts: 4
Joined: 17 Feb 2021, 15:24

### Asimmetry in trivial 2D expansion using TDVP

Hello,

I am finding some issues in the time evolution calculation I am doing with TDVP, surely this is a matter of me using the library the wrong way, but I can't catch where. In principle, I have a Bose Hubbard model with hard-core bosons in a square lattice where we also include dipolar interactions up to next nearest neighbours. So the Hamiltonian looks like

$$\hat{H} = -t \sum_j (\hat{b}_j^\dagger \hat{b}_{j+1} + \mathrm{h.c.}) + V \sum_{<<i,j>>} \frac{1-3\cos\theta_{ij}}{r_{i,j}^3} \hat{n}_i \hat{n}_j$$.

Then I define my initial state as a cluster of particles of size n x n in the center of the lattice. For example, a 3 x 3 cluster in the centre of a 9 x 9 lattice. I then do time evolution using TDVP and check the expansion in different directions in the lattice of this cloud. Turns out that with my current parameters, I see a preference for the system to expand faster in Y than in X, no matter what I choose for the dipole direction or the magnitude of V.

So I simplified the model and set V = 0 so that only the tunnelling term remains. I also set my initial state as a single particle in a 5x5 lattice. Turns out that if I set my trunc_cut to 10^-6, I see a discrepancy in the expansion radius and expansion velocity of the cloud. This can be perceived in the following animation in the first time steps. If I set trunc_cut to 10^-12, the discrepancy is much smaller.

The issue reappears for a larger system though. With the same parameters but a 3x3 cluster in an 11x11 lattice, I get a discrepancy in the expansion on X and Y. What am I doing wrong?

I attach some figures corresponding to the latest example and the script I'm using to calculate the results.

Code: Select all

import sys
import pickle

import numpy as np

import tenpy

from tenpy.models import lattice
from tenpy.models.model import CouplingMPOModel
from tenpy.models.hubbard import BoseHubbardChain, BoseHubbardModel
import tenpy.linalg.np_conserved as npc
from tenpy.networks.mps import MPS
from tenpy.tools.params import asConfig

from tenpy.algorithms import tebd, tdvp

from scipy.linalg import eig, eigh
from scipy.constants import h

J = 27
V = 0

t_f = 10/J

N_steps = 50

L_x = 11#7
L_y = 11#7
l = 3

M = L_x * L_y

direction = [2, 1, 0]
cutoff = 2

if len(sys.argv) > 1:
print(sys.argv)
V = int(sys.argv[1])
direction = [float(value) for value in sys.argv[2:5]]
if len(sys.argv) > 5:
cutoff = int(sys.argv[5])
print('cutoff:', cutoff)

dipole = np.array(direction)
dipole = dipole/np.linalg.norm(dipole)

def cluster(L_x, L_y, shape):
x, y = shape
loc = []
for j in range(y):
for i in range(x):
pos = i + L_x//2 - x//2
pos += (j + L_y//2 - y//2) * L_y
loc.append(pos)

return loc

locations = cluster(L_x, L_y, (l,l))

N = len(locations)

def dipolar_int_coeff(V, s_1, s_2, r):
#a = 266 * 1e-9 #[m]
if len(r) == 2:
r = np.array((r[0], r[1], 0))
#r = [a * item for item in r]
assert len(s_1) == len(s_2) == len(r)
return V/np.linalg.norm(r)**3 * (3 * s_1.dot(r) * s_2.dot(r)/np.linalg.norm(r)**2 - s_1.dot(s_2))

class DipolarModel2D(BoseHubbardModel):
def init_terms(self, model_params):
V = model_params.get('V_dip', 0)
dip = model_params.get('dir_dip', np.zeros(2))
cutoff = model_params.get('cutoff_dip', 1)
for u1, u2, dx in self.lat.pairs['nearest_neighbors']:
# TO-DO: define dipolar interaction with function based on vector direction
coeff = dipolar_int_coeff(V, dip, dip, dx)
print('nn', dx, coeff)
self.add_coupling(coeff, u1, 'N', u2, 'N', dx)
if cutoff > 1:
for u1, u2, dx in self.lat.pairs['next_nearest_neighbors']:
coeff = dipolar_int_coeff(V, dip, dip, dx)
print('nnn', dx, coeff)
self.add_coupling(coeff, u1, 'N', u2, 'N', dx)
super().init_terms(model_params)

model_params = {'lattice': lattice.Square,
'bc_x': 'open',
'bc_y': 'open',
'Lx': L_x,
'Ly': L_y,
'order': 'snake',
't': J,
'U': 0,
'V': 0,
'mu': 0,
'n_max': 1,
'conserve': 'None', #'parity', # 'N', #'None',
'V_dip': V,
'dir_dip': dipole,
'cutoff_dip': cutoff,
'bc_MPS': 'finite'}

model = DipolarModel2D(model_params)

p_state = [0 for _ in range(M)]
for loc in locations:
p_state[loc] = 1

psi = MPS.from_product_state(model.lat.mps_sites(), p_state, bc=model.lat.bc_MPS)

tdvp_params = {
'N_steps': 1,
'dt': t_f/N_steps,
'active_sites': 2,
#'trunc_params': {'chi_max': 250, 'svd_min': 1.e-12},
#'trunc_params': {'chi_max': 200, 'svd_min': 1.e-12, 'trunc_cut': 1e-7},
'trunc_params': {'chi_max': 500, 'svd_min': 1.e-12, 'trunc_cut': 1e-12},
}

eng = tdvp.TwoSiteTDVPEngine(psi, model, tdvp_params)

corr = eng.psi.correlation_function('N', 'N')

def calc_locs(lat, flag = 'nearest_neighbors', double_count = False):
pos = [tuple(item[:2]) for item in lat.order]
res = []
for i in range(len(pos)):
res_i = []
for nn in lat.pairs[flag]:
_, _, vec = nn
for c in range(double_count+1):
new_pos = (pos[i][0] + (1-2*c) * vec[0],
pos[i][1] + (1-2*c) * vec[1])
if new_pos in pos:
res_i.append(pos.index(new_pos))
res.append(res_i)

return res

def measurement(eng, data):
keys = ['t', 'N', 'NN', 'NNN', 'corr']
if data is None:
data = dict([(k, []) for k in keys])
data['t'].append(eng.evolved_time)
corr = eng.psi.correlation_function('N', 'N')
n = eng.psi.expectation_value('N')

#occ = [np.sqrt(corr[i,i]) for i in range(corr.shape[0])]
data['N'].append(n)

NN_indexes = calc_locs(model.lat)
NN = [corr[i, j].sum() for i, j in enumerate(NN_indexes)]
data['NN'].append(NN)

NNN_indexes = calc_locs(model.lat, 'next_nearest_neighbors')
NNN = [corr[i, j].sum() for i, j in enumerate(NNN_indexes)]
data['NNN'].append(NNN)

data['corr'].append(corr)

return data

data = measurement(eng, None)

for _ in range(N_steps):
eng.run()
measurement(eng, data)

with open(f'output/occ_corr_J{J}_V{V}_Lx{L_x}_Ly{L_y}_dipole{direction}_cutoff{cutoff}.pickle', 'wb') as handle:
pickle.dump(data, handle, protocol=pickle.HIGHEST_PROTOCOL)

Tana
Attachments
Screenshot from 2024-01-30 10-34-53.png (77.54 KiB) Viewed 42014 times
lattice_2D_V0_Lx11_Ly11.png (60.28 KiB) Viewed 42014 times
tanahy
Posts: 4
Joined: 17 Feb 2021, 15:24

### Re: Asimmetry in trivial 2D expansion using TDVP

Ok, so I got confused and thought two-side TDVP was as resilient to long time-steps as one-site TDVP, which is not the case according to the documentation. Thus, I heavily reduced the time step and the expansion seems isotropic now!
The thing is, there is a peak (or peaks) in the total particle number at the first time steps still, which I guess is correlated with this issue of the system somehow preferring the bonds along the y-axis.
Is this an expected result in any way?

Best,
Tana
Attachments
N_total_L10_V0_dir[0.0, 0.0, 2.0].png (18.7 KiB) Viewed 41983 times
Johannes
Posts: 442
Joined: 21 Jul 2018, 12:52
Location: TU Munich

### Re: Asimmetry in trivial 2D expansion using TDVP

This reminds me of one of my first work on tensor networks, arXiv:1509.00696 - using back then new W_II method that we now have in the ExpMPOEvolution in TeNPy

While the Hamiltonian is symmetric under exchange of x/y, the MPS snaking through the square lattice explicitly breaks it - hopping in y-direction is nearest-neighbor in the MPS, while x-direction hopping in the square lattice is long-range in the MPS (Ly-th nearest neighbor in the center).
Both 1- and 2-site TDVP approximate the exact time evolution by sequential left-right sweeps in the MPS, introducing an O(dt^2) error (similar to the trotter error in TEBD). The hopping in x-direction needs to happen by successively transporting a particle over the Ly-distance, which can be problematic if you truncate along the way.

Back then, we found that we can only reach a few hopping times before the truncation error became too large anyways. Thus, a smaller lattice was sufficient to capture the expansion; the particles couldn't reach the boundaries anyways. I would expect that a smaller Ly probably translates into a larger dt to still get accurate results, maybe you want to try that.

The particle number should be conserved by the evolution (since H commutes with sum_i n_i) - I see that you already tried to put conserve='N' in your model - I'd recommend that as well. if you calculate sum(psi.expectation_value('N')) you might still see it differ from the expected total N, as in you plot in the second post - that should come from the state not being in canonical form perfectly, but it's basically unaffected by the truncation error (at least if you conserve N). In other words, if you only check N(t) - N(t=0), you might vastly underestimate the error along the evolution...