Groundstate and GS energy from DMRG

[skothegh]

Dear Sir or Madam,

Currently I try to understand your library as I need it for work and I have some issues understanding some parts of it. I was advised to post my questions on your GitHub page and the "issues" tab is the closest thing to a forum I found.

My current problem is to understand how your DMRG class works exactly. I could not find a function that allows me to simply extract the groundstate enegery or the groundstate. The closest thing I could find is "run_two_site()", which should return the energy after the optimization. It returns a different energy after each run, though. So I assume it is not returning the groundstate energy.

Furthermore, if I measure the local spin states of the chain it is as random after the dmrg run as before. I would assume that the state would be either in the ferromagnetic or anti-ferromagnetic state.

Now to the questions: Which function returns the groundstate (energy)? Or did I completely misunderstand how the dmrg class works? Do you plan on writing a documentation and maybe a tutorial on MPOs and DMRG? And last but not least, is this the correct place to post such questions? Given the lack of a documentation on these things, I would think that I will come up with many more questions in the future.

All the best sk

My code:

mps = tn.FiniteMPS(initial_state)

mpo = tn.matrixproductstates.mpo.FiniteMPO(op_list)

test_dmrg = tn.matrixproductstates.dmrg.FiniteDMRG(mps,mpo)

print(np.real(mps.measure_local_operator(loc_op,site_list)))

test_dmrg.run_two_site(max_bond_dim = 50, num_sweeps=10,verbose=0)

print(np.real(mps.measure_local_operator(loc_op,site_list)))

[mganahl]

Hi @skothegh, thanks for the issue. Yes, this is indeed the right place to ask such questions. The energy can be obtained using the FiniteDMRG.compute_energy() method.

Unfortunately, the code above is incomplete, so I can't test it. Please post a complete, self-contained working example, and I will have a look.

[skothegh]

Sorry mganahl,

I thought that was sufficient. This here should do the trick. I did use compute_energy() but it spit out the same value before and after running run_two_site(), so I thought it was either something different or I made another mistake.

Best, sk

import tensornetwork as tn
import numpy as np

#### PARAMETERS ####
J = -12.0
h = 0.0

time = 10.0
tsteps = 1000
dtime = time/tsteps

Nsites = 50

#### OPERATORS ####
sz = np.array([[1,0],[0,-1]])

sx = np.array([[0,1],[1,0]])

sp = np.array([[0,1],[0,0]])

sm = np.array([[0,0],[1,0]])

Id = np.array([[1,0],[0,1]])

spsp = np.zeros([2,2,2,2])
smsm = np.zeros([2,2,2,2])
smsm[1][1][0][0] = 1.0
spsp[0][0][1][1] = 1.0

szsz = np.zeros([2,2,2,2])
szsz[0][0][0][0] = 1.0
szsz[0][1][0][1] = -1.0
szsz[1][0][1][0] = -1.0
szsz[1][1][1][1] = 1.0

sxI = np.zeros([2,2,2,2])
sxI[0][0][1][0] = 1.0
sxI[0][1][1][1] = -1.0
sxI[1][0][0][0] = -1.0
sxI[1][1][0][1] = 1.0

spsm = np.zeros([2,2,2,2])
spsm[0][1][1][0] = 1.0
spsm[1][0][0][1] = 1.0

Isx = np.zeros([2,2,2,2])
Isx[0][0][0][1] = 1.0
Isx[0][1][0][0] = -1.0
Isx[1][0][1][1] = -1.0
Isx[1][1][1][0] = 1.0

I = np.zeros([2,2,2,2])
I[0][0][0][0] = 1.0
I[0][1][0][1] = 1.0
I[1][0][1][0] = 1.0
I[1][1][1][1] = 1.0

#### INITIAL STATE ####
up = np.zeros([1,2,1])
up[0,0,0] = 1.0

down = np.zeros([1,2,1])
down[0,1,0] = 1.0

initial_state = []
for site in range(Nsites):
    if(np.random.rand() < 0.5):
        initial_state.append(up)
    else:
        initial_state.append(down)

op_list = []

op1 = np.zeros([1,3,2,2])
op1[0][1][:][:] = -h*sx
op1[0][1][:][:] = -J*sz
op1[0][2][:][:] = Id
op_list.append(op1)

op2 = np.zeros([3,3,2,2])
op2[0][0][:][:] = Id
op2[1][0][:][:] = sz
op2[1][1][:][:] = -h*sx
op2[2][1][:][:] = -J*sz
op2[2][2][:][:] = Id

for site in range(1,Nsites-1,1):
    op_list.append(op2)

op3 = np.zeros([3,1,2,2])
op3[0][0][:][:] = Id
op3[1][0][:][:] = sz
op3[2][0][:][:] = -h*sx
op_list.append(op3)

loc_op = []
for site in range(Nsites):
    loc_op.append(sz)

site_list = []
for site in range(Nsites):
    site_list.append(site) 

mps = tn.FiniteMPS(initial_state)

mpo = tn.matrixproductstates.mpo.FiniteMPO(op_list)

test_dmrg = tn.matrixproductstates.dmrg.FiniteDMRG(mps,mpo)

print(np.real(mps.measure_local_operator(loc_op,site_list)))

test_dmrg.run_two_site(max_bond_dim = 50, num_sweeps=10,verbose=0)

print(np.real(mps.measure_local_operator(loc_op,site_list)))

[mganahl]

Thanks! It seems this is an overcomplete example, but fair enough :). Looking at your MPO, what model would you actually like to implement?

[mganahl]

Hi @skothegh, I had a look at the code! First, there is a bug in the definition of op1:

op1[0][1][:][:] = -h*sx
op1[0][1][:][:] = -J*sz

You are overwriting the first assignment in the second line there. Also, your definition of op2 is strange. Are you simulating the Ising model? If so, the MPO should look different. That said, the above doesn't affect the h=0 case you are considering here.

The reason the simulation is going nowhere is the initial state. For h=0, your initial state is already an eigenstate of the Hamiltonian. In more technical terms, the sparse solver used to minimize the local optimization problem within DMRG immediately returns because the local initial state is already an eigenstate of the local Hamiltonian at each site. To remedy this problem, you should initialize your simulation with an entangled state. Here is minimal working example:

import tensornetwork as tn
import numpy as np
import matplotlib.pyplot as plt

#### PARAMETERS ####
J = -12.0
Nsites = 50
backend='numpy'

#### OPERATORS ####
sz = np.array([[1,0],[0,-1]])
sx = np.array([[0,1],[1,0]])
Id = np.array([[1,0],[0,1]])

#### INITIAL STATE ####
up = np.zeros([1,2,1])
up[0,0,0] = 1.0

down = np.zeros([1,2,1])
down[0,1,0] = 1.0

initial_state = []
for site in range(Nsites):
    if(np.random.rand() < 0.5):
        initial_state.append(up)
    else:
        initial_state.append(down)

op_list = []

op1 = np.zeros([1,3,2,2])
op1[0][1][:][:] = -J*sz
op1[0][2][:][:] = Id
op_list.append(op1)

op2 = np.zeros([3,3,2,2])
op2[0][0][:][:] = Id
op2[1][0][:][:] = sz
op2[2][1][:][:] = -J*sz
op2[2][2][:][:] = Id

for site in range(1,Nsites-1,1):
    op_list.append(op2)

op3 = np.zeros([3,1,2,2])
op3[0][0][:][:] = Id
op3[1][0][:][:] = sz
op_list.append(op3)

loc_op = []
for site in range(Nsites):
    loc_op.append(sz)

site_list = []
for site in range(Nsites):
    site_list.append(site) 
D=2
mps = tn.FiniteMPS.random(d = [2] * Nsites, D = [D] * (Nsites-1),
                          dtype=np.float64,backend=backend)
plt.plot(np.real(mps.measure_local_operator(loc_op,site_list)))
plt.xlabel('site')
plt.ylabel('<Sz>')
plt.draw()
plt.show()
mpo = tn.matrixproductstates.mpo.FiniteMPO(op_list, backend=backend)
dmrg = tn.matrixproductstates.dmrg.FiniteDMRG(mps,mpo)
dmrg.run_two_site(max_bond_dim=10,num_sweeps=3,verbose=0,num_krylov_vecs=10)

plt.plot(np.real(mps.measure_local_operator(loc_op,site_list)))
plt.xlabel('site')
plt.ylabel('<Sz>')
plt.draw()
plt.show()

hope this helps!

[skothegh]

Hi mganahl,

you are right in all points. God knows where my head was yesterday and I am sorry to have bothered you with such a bad code and stupid questions. The first bug completely slipped past me and of course its an eigenstate. I even thought about that just to have it slip out of my mind a second later. And yes, it is supposed to be the MPO of the TFI model. And yes, its wrong. I faultily adapted the simple spin chain. If I am not mistaken, it does actually reproduce the TFI in the case of only two sites. It becomes wrong once you add more, though.

Thanks for your kind help and your patience in the face of such sloppy work. sk