HELP!!! Ising 1D model not working

Hello everybody, I’m a newbie to the ITensor library and in general to coding in Julia and I’ve been tasked for an exam to develop a code able to analyze basic properties (ground state energy, entanglement entropy etc.) of quantum many-body systems using ITensor so here I am desperately asking for help.

Looking through the documentation given to us by the devs of this wonderful library I found some code that was supposed to return ground state energy of a given model which i replaced with the infamous hamiltonian of 1d ising’s model.
I looked online and found a bunch of examples and I’ve been finding contrasting outcomes since.

In a 1d ising model considering h as our external magnetic field (transverse, so on the x axis) we’re supposed to have a critical point in h = 1 and in general with our field “shut off” as a function of the spin coupling J (the first term in our hamiltonian in front of the sum) the energy of our gs should look like a straight line with negative angular coefficient (feel free to underline any mistakes in my understanding of the model cause i’m still trying to figure it out theoretically).

What I’m getting from 3 different codes is a totally different behaviour.
What happens once i plot all the values returned by the code is that i get with my field shut off an upside down “module function” (like a triangle with no base) and once i turn my field on it turns out to be a parabola pointing downwards, i tried changing signs in my ampos but nothing changes, what am i doing wrong??

Here’s the code for the kind-hearted person that wants to dig deeper into this “Ising mistery”:

	# In this example we show how to pass a DMRGObserver to
	# the dmrg function which allows tracking energy convergence and
	# convergence of local operators.
	using ITensors
	  Get MPO of transverse field Ising model Hamiltonian with field strength h
	function tfimMPO(sites, h::Float64)
	  # Input operator terms which define a Hamiltonian
	  N = length(sites)
	  os = OpSum()
	  for j in 1:(N - 1)
	    os += 1, "Z", j, "Z", j + 1
	  for j in 1:N
	    os += h, "X", j
	  # Convert these terms to an MPO tensor network
	  return MPO(os, sites)
	  N = 100
	  sites = siteinds("S=1/2", N)
	  psi0 = randomMPS(sites; linkdims=10)
	  # define parameters for DMRG sweeps
	  nsweeps = 15
	  maxdim = [10, 20, 100, 100, 200]
	  cutoff = [1E-10]
	  create observer which will measure Sᶻ at each
	  site during the dmrg sweeps and track energies after each sweep.
	  in addition it will stop the computation if energy converges within
	  1E-7 tolerance
	    Sz_observer = DMRGObserver(["Sz"], sites; energy_tol=1E-7)
	    # we will now run DMRG calculation for different values
	    # of the transverse field and check how local observables
	    # converge to their ground state values
	    println("Running DMRG for TFIM with h=0.0")
	    H = tfimMPO(sites, 0.1)
	    energy, psi = dmrg(H, psi0; nsweeps, maxdim, cutoff, observer=Sz_observer)

and here’s my attempt of writing a code:

using ITensors

function main()
    # Definire i parametri del modello di Ising
    L = 10  # Numero di siti nella catena
    J = 1.0  # Accoppiamento di spin
    h = 0.5  # Campo magnetico esterno

    # Creare gli indici di sito per una catena di spin-1/2
    sites = siteinds("S=1/2", L)

    # Costruire l'Hamiltoniano utilizzando AutoMPO
    ampo = AutoMPO()
    for j=1:L-1
        ampo += J, "Sz", j, "Sz", j+1
        ampo += h, "Sx", j
    ampo += h, "Sx", L
    H = MPO(ampo, sites)

    # Inizializzare lo stato MPS casuale
    psi0 = randomMPS(sites)

    # Impostare i parametri per l'algoritmo DMRG
    sweeps = Sweeps(5)  # Numero di iterazioni DMRG
    maxdim!(sweeps, 100, 200, 400, 800, 1000)
    cutoff!(sweeps, 1E-10)

    # Eseguire l'algoritmo DMRG per trovare l'energia dello stato fondamentale
    energy, psi = dmrg(H, psi0, sweeps)

    println("Energia dello stato fondamentale: ", energy)


PS I’m sorry for my English, it’s been a lot since I last had to use it properly and it’s a bit “rusty”, add to it the fact that specifical math/physics terms are different in my native language and here you have this “melting pot” of weird terms tryin to make sense out of what I was tryin to say.
Thanks for the attention, hopefully someone wiser and more experienced than me can help me uncover the mistery! <3 <3 <3

While patiently refreshing the page every femtosecond I actually thought it could be useful to plot myself the expected theoretical behaviour using the equations is was given, so for completeness’ sake I used the following hamiltonian:

H = - J * \sum_{<i,j>} Sz_j * Sz_{j+1} - h * \sum_{k} Sx_j

and I calculated the ground state energy with the following equation:

E_k = J * \sqrt{1 + \frac{h^2}{J^2} + 2\frac{h}{J}*cos({\frac{2pik}{n}})}
with n, number of sites, k: [-n/2+1,n/2], and:
E_{ground state} = - \sum_{k} E_k

from which I obtained the following plot (i actually put h=0.05 but it’s nearly the same as the case with h=0):
ising 1d plot
where you can see the behaviour of the energy as a function of the spin coupling, which on the other hand looks like this when I plot my numerical results from the script:
ising J variabile

I even took a step further and plotted the difference between the ground state and first excited but i think it’s already a step further from what I’m tryin to solve, so let me know if they might be helpful!

Once again, thanks for the attention, hopefully this helps! fingers crossed!

I hope I can help, but I’m not really sure what is the main question you are asking? I think the question is that you are getting different energies than you expected when comparing to a theoretical prediction, is that correct?

Some points that may be helpful:

  • the "Sz" and "Sx" operators in ITensor are defined as have a relative factor of 1/2 compared to the operators "Z" and "X", which are the Pauli z and x matrices. So please do make sure you are using these consistently and comparing Hamiltonians consistently between theory and code. I believe most transverse-field Ising literature uses the Pauli matrices where the critical point is located at a field value of h=1.
  • the theoretical expression you posted may only be perturbative thus only valid in a certain range of h values. I’m not sure if this is true or not, but thought it could be good for you to double check.
  • it could be important to carefully take into account both finite-size and especially open-boundary effects when computing the energy. I would expect most theoretical expressions for the TFIM are either in the infinite-size limit or for systems with periodic boundary conditions. Whereas ITensor DMRG works on finite-size, open boundary systems. There is generically a large contribution to the energy from the boundaries. This contribution can be subtracted off but requires a few steps to do that correctly, which I would be happy to discuss with you.

Really appreciate the reply :pray: :pray:

Since I posted the question I had the opportunity to chat with my professors regarding the problems I’ve been facing for the last week and it came up that the code was actually working fine, it was the theoretical model that had a small, yet crucial problem: it was missing a module, which changed everything.

Sorry for the confusion, I was overwhelmed by the theory and coding and totally took for granted that the model was right.

Anyways, since you have been so kind to help me out, if you have any advices on how to deal with these kind of matters or how to actually learn how to use properly this awesome library I’m all ears!

Thank you once again for taking the time to reply, having an active community always ready to help is something that for a newbie like me is truly valued and apprecciated.