About bondgate and beyond nearest neigbor

Hi,
I’m trying to do time evolution for a 1D chain with next nearest neighbor interactions. The following code worked for me, but I want to make sure that this is indeed doing what I am intending. This term is supposed to be N_j N_{j+2} interaction term (boson site).

auto s1 = BondGate(sites,j+1,j+2); // swap j+1 and j+2
gates.push_back(s1);
hterm = op(sites,"N",j)*op(sites,"N",j+1);
auto g3 = BondGate(sites,j,j+1,BondGate::tReal,tstep/2.,hterm);
gates.push_back(g3);	
gates.push_back(s1);// swap back.

Since I swapped j+1 and j+2 in the beginning, I thought I should write hterm as it is written above, between j and j+1 not j and j+2. After pushing g3, I swap back j+1 and j+2.
Is this correct?

Also, is this the recommended way of coding TEBD for chains beyond nearest neighbor interactions? Is there any other method to achieve this?

Thanks for the comments!
Ceren.

EDIT: Based on a check with ED, this code seems to be giving the right physics indeed.