Different coupling terms on different sites

[on0001st]

Hi all,

I am trying to implement a DMRG calculation using the spin_nnn model with the coupling term being different on different sites. To my understanding, we can do this for magnetic fields by defining a list with the different magnetic fields for different sites and passing the list as the model parameters. However, when I do the same for the coupling, I get a vector length issue. It would seem this is because the array length for the onsite terms are N and N-1 for the coupling terms. Hence, I would like to ask what is the way to go about implementing non-homogeneous couplings.

Sorry if this is the wrong place to post this. Please move it to the correct section if it is wrong.

Edit: I would also like to ask what is the reason that there is onsite terms in spin_nnn.py but not in spin.py for the couplings.

[Johannes]

To answer first in general (That's the perfect place to ask this ):

If you give the strength parameter to add_coupling, it is tiled to the "correct" size of the couplings.
You are right, for onsite terms in add_onsite, this is simply the Ls of the lattice.

Now what is the "correct" size in add_coupling?
That depends on the dx and the bc_coupling.
For open boundary conditions, it's (Ls[0] - |dx[0]|, Ls[1] - |dx[1]|, ...), for periodic boundary conditions it's simply the Ls again,
for a mix of open and periodic it's the straightforward generalization.
(For the MultiCouplingModel, it's a bit more complicated, but again the answer to how much you can shift the "box" of couplings around without violating the boundary conditions.)

Now in your case of the SpinChainNNN, it actually doesn't work for open/finite boundaries
You have a simple Chain with `L` DoubleSite, so `2*L` usual sites in total (!).
If you look into the model code, you see that it uses the same parameters for couplings with dx=0 and dx=1. For that reason, your array shape is always wrong - except if you just give a single number...
The DoubleSite is anyways just a kind of hack to allow a time evolution of a next-nearest neighbor model with TEBD, as long as we don't have time evolution from an MPO. This is only a temporary solution and also makes code run slowlier.
If you're interested in DMRG and not in time evolution, I strongly advise you to take the code of the usual SpinChain instead and extend it by a next-nearest neighbor coupling instead, without making it a NearestNeighborModel (which a next-nearest neighbor model usualy isn't...).
I've just added the according code in the spin_nnn.py as the SpinChainNNN2 in the newest git commit.
With that model, it should be clear which array sizes you can give to the model parameters, right?
For open boundary conditions, hz gets `L`, Jz gets `L-1` and Jzp gets `L-2` entries

[on0001st]

Thank you so much for the swift reply!

I have pulled the latest git commit and am able to run the code using spinChainNNN2. However, as not it is not a nearest neighbour model, I am not able to to use the method of

Python: Select all

E = np.sum(psi.expectation_value(M.H_bond[1:]))

to calculate the energies. I have tried fiddling around the code and one way I can come up with is to directly read the energies from the dmrg run by

Python: Select all

E, shelve, bond, sweep = dmrg.run(psi, M, dmrg_params)

but I think that it is quite "brute force" like and not really the way it should be done. I was thinking of directly calculating the expectation value by using the MPS psi and the MPO of the Hamiltonian but I cannot really figure it out. Sorry if this is a really trivial question but how should I calculate energies now?

[Johannes]

Taking the energy as returned by the DMRG code is a perfectly reasonable thing And it's what I would recommend you for now.

In general, you could still try to get the energy from local expectation values and/or correlation functions, but you would need to carefully figure all the terms out and sum them up.
The best way would be to have a method MPO.expectation_value(psi).
Unfortunately, such a method is not yet implemented, since calculating the expectation values for a general MPO with respect to a general MPS is very non-trivial for infinite boundary conditions. But using the MPOEnvironment, it should not be too hard to implement it for a finite-range Hamiltonian/MPO and/or finite systems, which would do the job for most cases...
In fact, for finite systems, you can just use MPOEnvironment(psi, H, psi).full_contraction().

Do you consider finite or infinite systems?

[on0001st]

In fact, for finite systems, you can just use MPOEnvironment(psi, H, psi).full_contraction().

Yes, this was what I was wondering how it can be done. I see in the method of full contraction that there is a need to specify a site index. As we are contracting the entire MPS to calculate for energy, is my understanding correct in saying that it does not matter which site we choose in this case?

Do you consider finite or infinite systems?

For now, we are considering a finite system. Based on what I understand, iDMRG can still be calculated for such a system but the only way to obtain the energy is directly from the iDMRG output?

[Johannes]

As we are contracting the entire MPS to calculate for energy, is my understanding correct in saying that it does not matter which site we choose in this case?

Yes, that's correct!

Based on what I understand, iDMRG can still be calculated for such a system but the only way to obtain the energy is directly from the iDMRG output?

For the moment, yes, until someone implements a function for it.

[orkiss]

Hi!

How can we include site-dependent coefficients for 2d lattices, using a square lattice? The size of coupling seems to be (L-1), but i would like to have 2*(L-1)**2 different coefficients (so all different).

Thank you!!

[Johannes]

As documented, the strength argumente of add_coupling can be a numpy array - for a square lattice with Lx by Ly sites, for nearest-neighbor hopping/interaction in y -direction, it could be of shape (Lx, Ly-1), or for the one in x-direction (Lx-1, Ly), respectively.