Creating random state with sign definite quantum numbers

Hi everyone,

I am working on symmetry broken phi^4 scalar field theory. I expect the field expectation to pick up a (near) constant nonzero value after the DMRG calculation because the mass squared term in the Hamiltonian is <0 while the coupling is >0. However, when working in ITensor, I can’t seem to get this configuration, I always get numerous kinks in an otherwise almost flat field configuration, where, between the kinks, the field expectation is of the same magnitude but opposite sign. I do not know what is going wrong exactly, my matrix elements of field and momenta operators should be correct, as well as my written Hamiltonian. I suspect the issue lies in the initialized random state corresponding to a highly fluctuating field expectation, that gets stuck with kinks when the DMRG tries to minimize the energy. Is there any functional control over the random MPS function where one can get better control over the initial state such that the initialized state would correspond to a positive or negative definite field expectation. Running the code below should give you similar results to what I’m getting. Any input would be greatly appreciated. Thank you.


using ITensors
using Plots

function kd(x,y)
    if x == y
        return 1
    else
        return 0
    end
end


let

    d= 5
    N= 1000
    mu2= -5
    ld = -3*mu2
    tau = 10
    sites = siteinds("Boson", N, dim=d , conserve_qns= false   )
    #return sites
   


    function ITensors.op(::OpName"phi",::SiteType"Boson", d::Int )
        mat = zeros(d,d)
        for i in 1:d
            for j in 1:d
                mat[i ,j] = (1/ sqrt(2))*( kd(i-1, j)*sqrt(j+1)+ kd(i, j-1)*sqrt(i+1) )
            end
        end
        return mat
    end
    function ITensors.op(::OpName"cpi",::SiteType"Boson", d::Int )
        mat = zeros(d,d)
        for i in 1:d
            for j in 1:d
                mat[i ,j] = (im/sqrt(2))* ( kd( i, j+1)* sqrt(j+1)- kd(i+1,j)*sqrt(i+1))
            end
        end
        return mat
    end
    function ITensors.op(::OpName"phi2",::SiteType"Boson", d::Int )
        mat = zeros(d,d)
        for i in 1:d
            for j in 1:d
                mat[i ,j] =(1/2)*( kd(i-1, j+1)*sqrt((j+1)*(j+2))+ kd(i,j)*(2*i +1)+ kd( i+1,j-1)*sqrt((i+1)*(i+2)))
            end
        end
        return mat
    end
    function ITensors.op(::OpName"cpi2",::SiteType"Boson", d::Int )
        mat = zeros(d,d)
        for i in 1:d
            for j in 1:d
                mat[i ,j] = (1/2)*( -kd(i-1, j+1)*sqrt((j+1)*(j+2))+ kd(i,j)*(2*i+1)- kd(i+1,j-1)*sqrt((i+1)*(i+2)))
            end
        end
        return mat
    end


    op = OpSum()
    for l in 1: N-1 #bulk terms 
        op += (-1),"phi",l,"phi",l+1
        op += (1/2)*(1+mu2), "phi",l,"phi",l
        op += (1/2), "phi",l+1,"phi",l+1
        op += (1/2), "cpi2",l
        op +=  (ld/24), "phi2",l,"phi2",l
    end 
# boundary terms for PBC interaction
        op += (1/2)*(1+ mu2), "phi",N,"phi",N
        op += (1/2), "phi",1,"phi",1
        op += (1/2), "cpi2",N
        op +=  (ld/24), "phi2",N,"phi2",N
        op += (-1),"phi",N,"phi",1

    H = MPO( op, sites)

    inis= randomMPS(sites)

    sweeps = Sweeps(100)
    setmaxdim!( sweeps, 100 )
    setcutoff!(sweeps,1E-9)
   #@show sweeps

    energy, state= dmrg(H, inis, sweeps)
    #println(" Final energy = $energy")


    x= collect(1:N)
    y = expect(state, "phi", sites= x)
    plot(x,y, label=["exp phi"])
end

If you’re running into issues with DMRG convergence, there’s a good chance the advice in this FAQ listing could help you:
https://itensor.github.io/ITensors.jl/dev/faq/DMRG.html#Preventing-DMRG-from-getting-stuck-in-a-local-minimum
Specifically, I noticed that you aren’t using a non-zero “noise” value in your code.

Also please see the listing just above it, about ensuring a DMRG calculation is converged.

I hope one of those tips help!

Hi Jake, It will be very tough to converge this model from a random state. As you have noticed, getting to the ground state requires all domain wall pairs to meet and effectively annihilate each other. As N grows and mu goes more negative the number of sweeps required probably grows nonlinearly. Without using the noise feature you are effectively quenching a T=inf state to T<T_c which tends to freeze in many domain wall boundaries. I found for N=100 I could easily get to a GS using the noise feature (which emulatues simulated annealing), but for N=1000 I usually got two domain walls separated by 100s of lattice sites. With enough sweeps they will eventually meet and cancel.
You can get a GS very quickly by avoiding a random state, using

state= MPS(sites,1) #<phi>=0.0, state[n]=[1.0  0.0  0.0  0.0  0.0]
energy, state= dmrg(H, state;nsweeps=10, maxdim=1, cutoff=1e-6, noise=0) #iterate to <phi>~=1.15

If you increase maxdim all it does is fine tune the state for some tiny edge effects.
I hope this helps.
JR

Thank you both for your feedback. I got convergence by simply not using a random MPS but just one of my construction. Thanks again.