tensor_builder does not support 3-body+ terms; is there any workaround?

[lolazon22]

Hello all,
I am trying to find the ground state of the generalized cluster hamiltonian using DMRG. This hamiltonian contains a 3-body interaction term. This is the code I'm using to construct it and find the ground state:

from quimb import *
from quimb.tensor import *
builder = SpinHam1D(S=1 / 2)
builder -= 1.0, 'X', 'X'
builder -= 1.0, 'X', 'Z', 'X'
builder += 1.0, 'Z'
mpo_ham = builder.build_mpo(30)
dmrg = DMRG2(mpo_ham)
dmrg.solve(max_sweeps=6, verbosity=1, cutoffs=1e-6)
ground_state = dmrg.state

The following error appears: NotImplementedError: 3-body+ terms are not supported yet.

Is there any possible workaround? I appreciate any help. Thank you in advance.

[jcmgray]

Hi @lolazon22. There are few options for you, in order of implementation effort:

1. Build each product as an MPO, then simply add and compress them.
2. Work out how to manually construct the W-matrices in order to build the MPO directly - the discussion here might also be of use
3. Wait a little bit for me to add the arbitrary hamiltonian MPO builder to quimb.experimental.operatorbuilder.py. I currently have a basic method that should be much better than adding and compressing, but is not optimal in all cases. It would take a little bit of time to polish into shape even for the 'experimental' module, and I am a bit short of time, so wouldn't want to give any guarantees here...

[jcmgray]

I've pushed an initial general MPO builder here. See if it works for you, though it currently comes with no guarantees at all - you should test what it produces thoroughly!

[lolazon22]

Thank you @jcmgray, I really appreciate your answer! Given my current expertise (almost none) I decided to follow your discussion on how to manually construct the W-matrices. It was the most intuitive option for me. I already managed to construct the MPO and everything is running perfectly. I'll put the code below just in case you observe any possible error in it:

import quimb as qu
import quimb.tensor as qtn

L = 32
j = -2.0
I = qu.eye(2)
X = qu.pauli('X')
Z = qu.pauli('Z')

# we define the MPO tensor
W = np.zeros([4, 4, 2, 2], dtype=complex)
W[0, 0, :, :] = I
W[0, 1, :, :] = X
W[1, 2, :, :] = Z
W[2, 3, :, :] = j*X
W[3, 3, :, :] = I

# define the end pieces
Wl = W[0, :, :, :]
Wr = W[:, 3, :, :]

# build the MPO
H = qtn.MatrixProductOperator([Wl] + [W] * (L - 2) + [Wr])

# add both MPOs
builder = qtn.SpinHam1D(S=1 / 2)
builder -= 1.0, 'X', 'X'
builder += 1.0, 'Z'
mpo_ham = builder.build_mpo(L)

ham = H + mpo_ham
ham.compress()

# get ground state
dmrg = qtn.DMRG2(ham)
dmrg.solve(max_sweeps=6, verbosity=1, cutoffs=1e-6)
ground_state = dmrg.state

[jcmgray]

Yes all looks good to me, the only thing I'd note is you might want to explicitly construct the MPO in real form. I = qu.eye(2, dtype='float64'); W = zeros(..., dtype='float'64') etc or just cast ham.astype_('float64') (assuming the hamiltonian is going to be real of course).