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.
Python: 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
