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