Reproducing results of the purification code snippet for Heisenberg chain

[spiralsnip]

Hi,

I have a very low level question.

Having in mind to use the purification method to compute the energies of the Gibbs states of Fermi-Hubbard chains at different temperatures, I thought I’d start by understanding the purification code snippet I could find in ITensorMPS on the example of the Heisenberg chain (ITensors.jl/src/lib/ITensorMPS/examples/finite_temperature/purification.jl at main · ITensor/ITensors.jl · GitHub).

The problem is I cannot obtain results which match my understanding of what the code is doing (namely by comparison with brute force diagonalization, implemented with Python, see below).

Let us take the example of 3 sites up to beta=1.

Using the snippet (replacing N=10 with N=3), I get the following energies on the one hand:
β = 0.00 energy = 0.00000000
β = 0.10 energy = -0.03841640
β = 0.20 energy = -0.07854052
β = 0.30 energy = -0.12016190
β = 0.40 energy = -0.16302627
β = 0.50 energy = -0.20684070
β = 0.60 energy = -0.25128101
β = 0.70 energy = -0.29600131
β = 0.80 energy = -0.34064501
β = 0.90 energy = -0.38485656
β = 1.00 energy = -0.42829320

whereas with the Python code, I obtain:
beta=0.00 energy=0.00000000
beta=0.10 energy=-0.65201966
beta=0.20 energy=-1.36238215
beta=0.30 energy=-2.04609685
beta=0.40 energy=-2.62925102
beta=0.50 energy=-3.07783848
beta=0.60 energy=-3.39691169
beta=0.70 energy=-3.61204006
beta=0.80 energy=-3.75236435
beta=0.90 energy=-3.84224174
beta=1.00 energy=-3.89933520

Python code:

import numpy as np
import functools as ft
import scipy.sparse as sp

sigma_x = np.zeros((2,2), dtype='complex128')
sigma_x[0, 1] = 1
sigma_x[1, 0] = 1

sigma_y = np.zeros((2,2), dtype='complex128')
sigma_y[0, 1] = -1j
sigma_y[1, 0] = 1j

sigma_z = np.zeros((2,2), dtype='complex128')
sigma_z[0, 0] = 1
sigma_z[1, 1] = -1

id2 = np.eye(2)

sigma_plus = sigma_x + 1j*sigma_y
sigma_minus = sigma_x - 1j*sigma_y

def construct_heisenberg_model_mat(nsites):
    
    heisenberg_mat = np.zeros((2**nsites, 2**nsites), dtype='complex128')
    
    for i in range(nsites - 1):
        ops_list = [id2]*i + [sigma_plus] + [id2]*(nsites - i - 1)
        mat1 = ft.reduce(np.kron, ops_list)  
        ops_list = [id2]*(i+1) + [sigma_minus] + [id2]*(nsites - i - 2)
        mat2 = ft.reduce(np.kron, ops_list)
        
        sp_sm_term = 0.5*mat1@mat2                      
        heisenberg_mat += sp_sm_term
        
        ops_list = [id2]*i + [sigma_minus] + [id2]*(nsites - i - 1)
        mat1 = ft.reduce(np.kron, ops_list)  
        ops_list = [id2]*(i+1) + [sigma_plus] + [id2]*(nsites - i - 2)
        mat2 = ft.reduce(np.kron, ops_list)
        
        sm_sp_term = 0.5*mat1@mat2                      
        heisenberg_mat += sm_sp_term
        
        ops_list = [id2]*i + [sigma_z] + [id2]*(nsites - i - 1)
        mat1 = ft.reduce(np.kron, ops_list)  
        ops_list = [id2]*(i+1) + [sigma_z] + [id2]*(nsites - i - 2)
        mat2 = ft.reduce(np.kron, ops_list)
        
        sz_sz_term = mat1@mat2                      
        heisenberg_mat += sz_sz_term

    return heisenberg_mat

H_mat = construct_heisenberg_model_mat(3)

betas = np.arange(0, 1.1, 0.1)

for beta in betas:
    sig_beta = sp.linalg.expm(-beta*H_mat)
    sig_beta /= np.trace(sig_beta)
    energy = np.trace(sig_beta@H_mat)
    print('beta=%.2f energy=%.8f '%(beta, energy.real))

I also have a version of the code with imposes periodic boundary condition but still get a mismatch.

Could you please explain to me why the Python code does not yield similar results? Am I missing something?

Best regards,

SpSn

[sandipanmanna]

Hi! Thanks for posting. Just a suggestion. Sp,Sn are defined with a factor of (1/2) for the specific definition of your Pauli sigma matrices. Specifically, in the following segment.

sigma_plus = sigma_x + 1j*sigma_y
sigma_minus = sigma_x - 1j*sigma_y

Can you try defining it as follows.

sigma_plus =0.5*( sigma_x + 1j*sigma_y)
sigma_minus =0.5*( sigma_x - 1j*sigma_y)

You have to take care of Sz operator too (with a 0.5 factor). In ITensor, Sx,Sy,Sz include 0.5. Your definition of Pauli matrices correspond to using X,Y,Z operators respectively in ITensor ( the later is more common in QI texts).

[spiralsnip]

Hi,

Thank you very much! Indeed I hadn’t paid attention to the fact those were spin operators instead of the usual Pauli operators and their related ladder operators.

Adding the 0.5 coefficient and not imposing periodic boundary conditions enabled me to get the same results as ITensor.

Best regards,

SpSn