Hi,
I am studying the Bose-Hubbard model using Boson site:
sites = Boson(L, {"MaxOcc=", n, "ConserveQNs", true, "ConserveNb", false});
The reason why I set "ConserveNb" as false is that there is a term in Hamiltonian which breaks the particle number conservation: $$H_e=\sum_i{\left( \hat{a}_{j}^{\dagger 2}+\hat{a}_{j}^{2} \right)}$$. At this case, the U(1) symmetry of the system is broken and so it is very slow to get the ground state MPS using DMRG. Without symmetry the DMRG cannot spped up and calculate slowly. In my personal computer I can simply calculate @@L=12, n=3@@, which is a really small system when I want to study the correlation function. I am wondering if there is any efficient way to deal with this case while U(1) symmetry is broken.
Looking forward to the help. Thank you.
My Hamiltonian is:
Hamiltonian:@@\widehat{H}=\sum_{j=1}^L{\left[ -\mu \widehat{n}_j+\frac{U}{2}\widehat{n}_j\left( \widehat{n}_j-1 \right) -J\left( \widehat{a}_j\widehat{a}_{j+1}^{\dagger}+\widehat{a}_{j}^{\dagger}\widehat{a}_j \right) +V\widehat{n}_j\widehat{n}_{j+1} \right]}+em\sum_i{\left( \hat{a}_{j}^{\dagger 2}+\hat{a}_{j}^{2} \right)}@@
And my parameters are as follows:
auto sweeps = Sweeps(8);
sweeps.maxdim() = 40, 80, 400, 400, 800, 800, 1000, 1000;
sweeps.cutoff() = 1E-16;
Real mu = 1.0;
Real U = 1.0;
Real J = 1.0;
Real V = 1.0;
Real em = 1.0;
It takes me around 3 hours to calculate the ground state as L = 12, N = 3.
When I set the "ConserveNb" as true, it just take around 1 hour to calculate as L = 24, N = 3.
The total code is as follow:
auto sweeps = Sweeps(8);
sweeps.maxdim() = 40, 80, 400, 400, 800, 800, 1000, 1000;
sweeps.cutoff() = 1E-16;
Real mu = 1.0;
Real U = 1.0;
Real J = 1.0;
Real V = 1.0;
Real em = 1.0;
int L = 12;
int n = 3;
auto sites = Boson(L, {"MaxOcc=", n, "ConserveQNs", true, "ConserveNb", false});
auto ampo = AutoMPO(sites);
for (int i = 1; i < L; i += 1)
{
ampo += -miu - U / 2, "N", i;
ampo += U / 2, "N", i, "N", i;
ampo += -J, "A", i, "Adag", i + 1;
ampo += -J, "Adag", i, "A", i + 1;
ampo += V, "N", i, "N", i + 1;
ampo += -em / 2, "Adag", i, "Adag", i;
ampo += -em / 2, "A", i, "A", i;
}
ampo += -miu - U / 2, "N", L;
ampo += U / 2, "N", L, "N", L;
ampo += -em / 2, "Adag", L, "Adag", L;
ampo += -em / 2, "A", L, "A", L;
ampo += -J, "A", 1, "Adag", L;
ampo += -J, "Adag", 1, "A", L;
ampo += V, "N", 1, "N", L;
auto H = toMPO(ampo);
auto state = InitState(sites);
for (int i : range1(L))
{
if (i % 2 == 1)
state.set(i, "1");
else
state.set(i, "1");
}
auto psi = randomMPS(state);
auto [energy, psi0] = dmrg(H, psi, sweeps, "Quiet");