MPO for an exponential operator

I am trying to define the MPO for the operator exp(i theta n_j) for the Anyon Hubbard model(mapped to an interacting boson model). I have tried to define an operator using the op function for the matrix of exp(i theta n_j) and then use it in the operator sum for the Hamiltonian but am recieving a method error in trying to do so.

image

(ImageSource:DOI: 10.1038/ncomms1353)

function Ham(θ,N,d)
    s = siteind("Qudit",dim=d)
    function expnmat(θ, d)
        return [(i == j) ?  exp(1im*θ*(i-1)) : 0.0 for i in 1:d, j in 1:d]
    end
    
    statphase = op(expnmat(θ,d), s)
    
    
    sites = siteinds("Qudit",N,conserve_qns=true,dim=d)
    os = OpSum()
    for j=1:N-1
    os += ("adag",j,"a",j+1,"statphase",j)
    os += ("statphase",j,"adag",j+1,"a",j)
    os += 1/2,"n",j, "n",j
    os -= 1/2,"n",j    
    end
    H = MPO(os,sites)
    return H 
end

To make a new operator definition that OpSum will recognize, you will need to make a new overload of the ITensors.op method. Please see the code examples and walkthrough in the documentation here:
https://itensor.github.io/ITensors.jl/dev/examples/Physics.html#custom_op

And let us know if you have any questions about that process.