2nd neighbor simple update results

[ogauthe]

I try to reproduce my results on the J1-J2 model at finite temperature using PEPSKit. I use the 3 site simple update code from #171. An example of script can be found at https://github.com/ogauthe/BenchmarkPEPS/blob/main/run/run_j1j2.jl.

The simple update Hamiltonian is defined as

    terms = [CartesianIndex.(t) => h_pa for t in nearest_neighbours(lattice)]  # 1st order SU, 1st neighbor
    append!(
            terms,
            [CartesianIndex.(t) => J2 * h_pa for t in next_nearest_neighbours(lattice)],
        )
    su_hamilt = LocalOperator(physical_spaces, terms...)

The rest of the code is nearly the same for the first neighbor.

For J2=0.0, the results look correct and match those obtained with both frostspin and run_heisenberg.jl. However, for J2 > 0, the results do not match frostspin. I am pretty confident in frostspin results as they match high temperature series expansion. It is unclear to me where the difference comes from. It is pretty easy to get a factor 2 wrong, so I computed several data points. There is no match outside of J2=0:

D=7, chi=49, beta=1.0

J2 / energy	PEPSKit 1sntei	frostspsin 1stnei	PEPSKit 2nd	frostspin 2nd
J2 = 0.0	-0.1927	-0.1928	0.07598	0.07599
J2 = 0.15	-0.0941	-0.1742	0.00836	0.04712
J2 = 0.30	-0.0877	-0.1561	-0.00695	0.01630
J2 = 0.60	-0.075	-0.1196	-0.0386	-0.04736

D=7, chi=49, beta=2.0

J2 / energy	PEPSKit 1sntei	frostspsin 1stnei	PEPSKit 2nd	frostspin 2nd
J2 = 0.0	-0.29060	-0.29060	0.1576	0.1576
J2 = 0.15	-0.175	-0.267	0.04728	0.125
J2 = 0.30	-0.156	-0.2371	0.0158	0.0769
J2 = 0.60	-0.117	-0.166	-0.0527	-0.03637

There are some similarities when matching PEPSKit at beta=2 with frostspin at beta=1, so I thought about a double counting problem, but then considering J2/2 or J2*2 should fix it. Also the case J2=0 matches precisely for both values of beta.

Comparing different values for the bond energies in the unit cell, the results are much more asymmetric (breaking rotation/translation) with PEPSKit. This may be explained by 1st order Trotter vs 2nd order in frostspin, but maybe there is something more.
@Yue-Zhengyuan do you have any clue?

[Yue-Zhengyuan]

Thank you for reporting this! I'm new to finite temperature calculations, and need some time to figure it out. Would you like to also share the frostspin benchmark code, and the reference that performs the high-T expansion as well?

Can you clarify if 1stnei and 2nd in the energy table mean the average value of $\langle S_i \cdot S_j \rangle$ on each nearest neighbor and next-nearest neighbor bond?

[ogauthe]

Would you like to also share the frostspin benchmark code, and the reference that performs the high-T expansion as well?

The HTSE data come from https://doi.org/10.1103/PhysRevB.67.014416
frostspin run script is available at https://framagit.org/ogauthe/run_thermal_TN. You first need to clone frostspin and add it in your PYTHONPATH (I did not packaged it). You can then use call run_J2.py with the input. To be honest I do not really know to what extend it is usable by someone else.

Can you clarify if 1stnei and 2nd in the energy table mean the average value of ⟨ S i ⋅ S j ⟩ on each nearest neighbor and next-nearest neighbor bond?

Yes, I measure the average value of <S.S> over respectively first neighbor and second neighbor bonds (one should average on all bonds of the unit cell, but I was too lazy for that)

Here is the per site energy for J2=0.0 and J2=0.50. This really looks like some double counting for J2. However setting J2=0.25 or J2=1.00 still does not match HTSE data.

logs_benchmark_J1-J2_J2_0.00_SU2_D7.log
logs_benchmark_J1-J2_J2_0.50_SU2_D7.log
logs_J2_0.00.log
logs_J2_0.50.log

[ogauthe]

I am not sure where this is coming from. I could imagine there is a factor 2 hidden in the conventions, but then I do not understand why the case J2=0.0 works. Does SimpleUpdate detect that the 2nd neighbor term is actually zero and uses a different code path in such a case?

For the second neighbor interaction, it is common to implement each update in 2 steps, with 2 different intermediate sites to reduce lattice symmetry breaking. Did you implement that?

Slightly unrelated, I am also interested in benchmarking the full update. I did not implement it and I am curious about its accuracy. Am I correct that _, a2, b2 = PEPSKit.fixgauge_benv(Z,a,b) gives me the 2 new effective tensors that I can reunify with X and Y from PEPSKit._qr_bond or do I need to run it several time to converge a2 and b2?

[Yue-Zhengyuan]

Thank you for the additional information! I'll try to reproduce your findings first and then see what was wrong.

Slightly unrelated, I am also interested in benchmarking the full update. I did not implement it and I am curious about its accuracy.

I'll soon add the complete full update code (hopefully in this week; I'll notify you when I open the PR). It's actually the "fast" full update (PRB 92, 035142 (2015)); "fast" in the sense that the CTMRG environment is updated together with the iPEPS (instead of first updating the iPEPS only, and then recalculate all CTMRG tensors), so we don't have to fully re-converge the environment after every Trotter step (every 10 to 20 steps is enough). There will be a few limitations, though:

• I didn't implement FU for bipartite iPEPS, which needs some special care that (currently) does not fit well with the general algorithm for arbitrary $N_r \times N_c$ unit cell.
• It can only handle Hamiltonians with only nearest neighbor interactions. I don't think there's an easy generalization for the "fast" part when there's NNN interactions.

[Yue-Zhengyuan]

Does SimpleUpdate detect that the 2nd neighbor term is actually zero and uses a different code path in such a case?

Yes; for an NN Hamiltonian it switches to the usual 2-site SU. See the following part of code. I have a test on Hubbard model at half-filling to check that the 2-site and 3-site versions give (physically) the same results. But you can set force_3site = true to always use the 3-site version.

PEPSKit.jl/src/algorithms/time_evolution/simpleupdate.jl

Lines 176 to 193 in 9aca427

	function simpleupdate(
	peps::InfiniteWeightPEPS,
	ham::LocalOperator,
	alg::SimpleUpdate;
	bipartite::Bool=false,
	force_3site::Bool=false,
	check_interval::Int=500,
	)
	# determine if Hamiltonian contains nearest neighbor terms only
	nnonly = is_nearest_neighbour(ham)
	use_3site = force_3site || !nnonly
	@assert !(bipartite && use_3site) "3-site simple update is incompatible with bipartite lattice."
	if use_3site
	return _simpleupdate3site(peps, ham, alg; check_interval)
	else
	return _simpleupdate2site(peps, ham, alg; bipartite, check_interval)
	end
	end

For the second neighbor interaction, it is common to implement each update in 2 steps, with 2 different intermediate sites to reduce lattice symmetry breaking. Did you implement that?

I didn't do it this way. In my approach, each update is applied to an ⅃-shaped 3-site cluster:

PEPSKit.jl/src/algorithms/time_evolution/simpleupdate3site.jl

Lines 346 to 352 in 9aca427

	```
	r-1 M3
	|
	↓
	r M1 -←- M2
	c c+1
	```

I didn't simply use SVD to decompose gate * cluster into 3-sites, but first decompose the gate into a 3-site MPO,

PEPSKit.jl/src/algorithms/time_evolution/evoltools.jl

Lines 164 to 172 in 9aca427

	function gate_to_mpo3(
	gate::AbstractTensorMap{T,S,3,3}, trunc=truncbelow(MPSKit.Defaults.tol)
	) where {T<:Number,S<:ElementarySpace}
	Os = MPSKit.decompose_localmpo(MPSKit.add_util_leg(gate), trunc)
	g1 = removeunit(Os[1], 1)
	g2 = Os[2]
	g3 = removeunit(Os[3], 4)
	return [g1, g2, g3]
	end

regard the 3-site cluster as an MPS, and then used an approximate way to compress the MPO-MPS product.

PEPSKit.jl/src/algorithms/time_evolution/simpleupdate3site.jl

Lines 324 to 333 in 9aca427

	function apply_gatempo!(
	Ms::Vector{T1}, gs::Vector{T2}; trunc::TensorKit.TruncationScheme
	) where {T1<:PEPSTensor,T2<:AbstractTensorMap}
	@assert length(Ms) == length(gs)
	revs = [isdual(space(M, 1)) for M in Ms[2:end]]
	@assert !all(revs)
	_apply_gatempo!(Ms, gs)
	wts, ϵs, = _cluster_truncate!(Ms, trunc, revs)
	return wts, ϵs
	end

Maybe we can make the decomposition of the gate into MPO (gate_to_mpo3) more symmetric.

I tried to reduce symmetry breaking by rotating the iPEPS: I update the L cluster using the same code for ⅃, but with the iPEPS rotated by 90 degrees.

PEPSKit.jl/src/algorithms/time_evolution/simpleupdate3site.jl

Lines 417 to 424 in 9aca427

	for i in 1:2
	for site in CartesianIndices(peps2.vertices)
	r, c = site[1], site[2]
	gs = gatempos[i][r, c]
	_su3site_se!(r, c, gs, peps2, alg)
	end
	peps2 = (i == 1) ? rotl90(peps2) : rotr90(peps2)
	end

Oh I may now understand the mismatch with frostspin: looking at this part of code, you can figure out each bond in the unit cell is actually updated twice, hence the match between 3-site $\beta / 2$ with 2-site $\beta$. It doesn't matter for $T = 0$ ground state search, but for finite-temperature calculation this will be an issue.

[ogauthe]

Ok so I understand what is happening. LocalOperator is too clever, it automatically removed the empty terms for J2=0.0 and is_nearest_neighbour is true. As I did not set force_3site, I fell back to the 2 site, bipartite case. So the case J2=0.0 was actually a first neighbor algorithm.

I then checked that setting J2=0.0001 and there is indeed a factor 2 in the energy scale.

So I now think with this factor 2 included, the code gives correct values. I think this factor 2 should be fixed such that very small J2 match J2=0.0. I also think the interface should avoid ideally avoid such confusion, I know that is a point you already discussed with @lkdvos . Maybe the best solution is to wait for #207 to define a tiling interface and use it here too.

[ogauthe]

Thank you for the code explanations.

I tried to reduce symmetry breaking by rotating the iPEPS: I update the L cluster using the same code for ⅃, but with the iPEPS rotated by 90 degrees.

I still observe the algorithm to be quite asymmetric, which I can observe from the CTMRG corners having very different spaces and dimensions. It is not necessarily a problem, but it will bias the optimized wavefunction towards a specific kind of C4v breaking. For finite temperature applications, I try to avoid C4v breaking so I apply twice sqrt(gate) = exp(-tau/2*h) over both ┘ and ┌ (I use directly the gate, no MPO) and I get much more symmetric CTMRG corners, but also bond energies (I also use 2nd order Trotter, it helps)

Oh I may now understand the mismatch with frostspin: looking at this part of code, you can figure out each bond in the unit cell is actually updated twice, hence the match between 3-site β / 2 with 2-site β . It doesn't matter for T = 0 ground state search, but for finite-temperature calculation this will be an issue.

Precisely, when applying twice the operator it makes sense to use exp(-tau/2*h) (to be fair I had the same issue in frostspin).

[ogauthe]

There is still something wrong: applying the update twice in 2 different configurations should double beta, while here we observe the second neighbor code to end up applying half beta.

I think I found out the issue. In the MPO decomposition inside _get_gatempo_se at PEPSKit.jl/src/algorithms/time_evolution/evoltools.jl line 196, the 3 site operator is written as

op = (1 / 2) * (nb1x ⊗ unit + unit ⊗ nb1y) + permute(nb2 ⊗ unit, ((1, 3, 2), (4, 6, 5)))

This is correct for operators, there will be a sum of 2 half to recover the full operator. But for gates, this is not correct:
$$\frac{1}{2}\exp(-\tau h_1) + \frac{1}{2}\exp(-\tau h_2) = 1 -\frac{1}{2}\tau h_1 -\frac{1}{2}\tau h_2 + O(\tau^2) = \exp(-\frac{\tau}{2}H) + O(\tau^2) $$
When considering more terms on the plaquette, including the identity for a next nearest neighbor in the limit $J_2\rightarrow0$, one still gets a factor 1/2 in tau inside the exponential (instead of a scaling factor).

I think this construction is not correct. Instead of computing the MPO of the gates (=exponentiated Hamiltonian), the code should compute the MPO of the Hamiltonian, then exponentiate this MPO (there are of course other solutions)

[Yue-Zhengyuan]

Thanks for pointing this error in Trotter decomposition! I have opened #219 which addresses both the dt and the $C_{4v}$ breaking issues. Could you please help test it?

Now that dt is actually equal to the time evolution step of each SU iteration, I think in the following line of your code:
https://github.com/ogauthe/BenchmarkPEPS/blob/86c0d57e9bafc58fcae2a30ed69d259ec36c68d0/run/run_j1j2.jl#L152
you can change su_niter back to Int(round(Δβ / dt; digits=10)).

[Yue-Zhengyuan]

One thing to keep in mind is that the improved C4v symmetry might also be due to the reduced dt for each gate MPO. You may also need to change dt in your frostspin benchmark code.

[ogauthe]

Thanks for pointing this error in Trotter decomposition! I have opened #219 which addresses both the dt and the C 4 v breaking issues. Could you please help test it?

Indeed I observe the new 2nd neighbor to be much more symmetric now. I would be curious to know if this also improves ground state energies.

Now that dt is actually equal to the time evolution step of each SU iteration, I think in the following line of your code: https://github.com/ogauthe/BenchmarkPEPS/blob/86c0d57e9bafc58fcae2a30ed69d259ec36c68d0/run/run_j1j2.jl#L152 you can change su_niter back to Int(round(Δβ / dt; digits=10)).

This factor 2 is coming from the definition of the thermal density matrix. It is defined as
$$\rho(\beta) = \text{Tr}_{\text{ancilla}} \ket{\Psi(\beta/2)}\bra{\Psi(\beta/2)}$$
basically you contract $\exp(-\beta H/2)$ over the ket, then contracting the ket with the associated bra gives $\beta$. It is also present in the first neighbor case.

One thing to keep in mind is that the improved C4v symmetry might also be due to the reduced dt for each gate MPO. You may also need to change dt in your frostspin benchmark code.

Smaller dt help of course, but the asymmetry was there even for pretty small ones. It is now 2 orders of magnitude smaller, even when doubling dt between the 2 versions.

[Yue-Zhengyuan]

This factor 2 is coming from the definition of the thermal density matrix.

Sorry I was ignorant of this detail 😅 Also I'm glad to know that the C4v symmetry of the result is improved.

I would be curious to know if this also improves ground state energies.

I run the SU part of our J1-J2 example (source here, with spin U(1) symmetry). There J1 = 1, and we gradually increase J2 from 0 to 0.5. The current master branch gives

Dict{String, Any} with 4 entries:
  "corH" => [-0.338079 -0.338085; -0.338079 -0.338085]
  "mag" => [[0.112882 -0.112497; -0.112881 0.112496]]
  "corV" => [-0.291535 -0.291582; -0.291541 -0.291589]
  "e_site" => -0.491239
(E - E_ref) / abs(E_ref) = 0.006091732078066108

while the updated 3-site SU gives,

Dict{String, Any} with 4 entries:
  "corH" => [-0.318609 -0.318623; -0.318609 -0.318623]
  "mag" => [[0.103574 -0.103575; -0.103574 0.103575]]
  "corV" => [-0.318623 -0.318623; -0.318609 -0.318609]
  "e_site" => -0.490848
(E - E_ref) / abs(E_ref) = 0.00688255949657097

The energy is slightly higher but still acceptable in my opinion... corH, corV are $\langle S_i \cdot S_j\rangle$ on the horizontal, vertical NN bonds in the 2 x 2 unit cell. For the updated SU the C4v breaking is much less significant, as seen from the spin correlations. At the level of SU, it seems that some C4v breaking can slightly lower the energy.

Actually I encountered this issue of C4v breaking before, but obviously I didn't fully solve it...

[Yue-Zhengyuan]

BTW, for nearest neighbor Hamiltonians, I also have the "Neighborhood Tensor Update" almost ready (PhysRevB.104.094411; #144). In its simplest version, the nearest neighbor tensors around a bond (instead of SU weights or the full CTMRG environment) are used as the environment when updating and truncating a bond. If you think it's promising, I can create an NTU branch for you to test its performance.