+1 vote
asked by (190 points)

I'm calculating time evolution in the Bose Hubbard model (using my own version of the hubbard site set). I would like to obtain the full state of my system in the Fock (number) basis at a given point to work further with it using different (Python) code. Is there an easy way to do this?

commented by (52.6k points)
Hi Andreas,
So just to understand: you aren't looking to export the MPS form of your wavefunction to your python code, correct? But instead the full wavefunction in an exponentially large many-body basis (presumably for a pretty small basis size)?
commented by (190 points)
edited by
Yes, exactly. For example, for a system of 6 lattice sites and 8 particles.

More precisely, I  simulate the time evolution of  a larger system via TEBD (e.g. 16 sites and 8 particles) and would like to export the reduced density matrix of a small subsystem (e.g. 6 sites and up to 8 particles, dimension 3003x3003) in the Fock basis.

1 Answer

0 votes
answered by (52.6k points)

Assuming that you are using IQTensors (IQMPS) for your calculations, I would recommend the following:

(1) convert whichever tensors you want to export to just ITensors. So for an IQMPS convert it to an MPS. Or for an IQTensor (like a reduced density matrix) convert it to an ITensor.

(2) then call the "ordered" function on this ITensor to get a simple, plain tensor with a well-defined memory access pattern and contiguous storage in memory. For more information see the documentation of the ordered function on this page: http://itensor.org/docs.cgi?page=classes/itensor

auto t = ordered(T,i1,i2,i3); // <- here T is an ITensor
t(1,2,3); // <- how to access an element of the plain tensor t

(3) either use the .data() and .size() methods of the tensor object returned by the ordered function to write the tensor data to a file (binary stream) or explicitly loop over the indices of the returned tensor object and write them one at a time to a file in the way you wish (ascii or binary)

commented by (680 points)
Hi Miles,

I tried the method and it works perfectly well for IQTensor. However, for QDiagReal I got the error message
    doTask not defined for task ToITensor and storage type QDiagReal
What is the correct way of doing this? Thanks!
commented by (52.6k points)
Thanks for pointing this out; it was a missing feature. But it was easy to add so I just added it. Please pull the latest version of ITensor, compile, and let me know if it still doesn't work for you.
commented by (680 points)
edited by
Thanks very much for the quick reply. Using the new version, the toITensor works well, but the "orderedC" function gives an error message:
    doTask not defined for task [unknown] and storage type DiagReal
The function "ordered" gives another error message:
     applyFunc: function object has no operator() method for storage type
commented by (52.6k points)
I see. So do you really need to read a diagonal tensor using the ordered method? For example could you first multiply the diagonal tensor into a dense one and call ordered on that? Otherwise you will need to use the .real or .cplx methods to get the elements of the diagonal tensor. The reason ordered can't be defined for diagonal tensors is that we don't have sparse storage types for the basic tensors or we'd have to convert the diagonal tensors to dense to call ordered on them which could lead to accidental inefficient behavior.
commented by (680 points)
Thanks for the information. I am currently using a for loop to get all the elements and it works well. I'm just wondering whether a consistent method can be used so as to make the code more clear.
commented by (52.6k points)
Hm it is worth thinking about. One thing that could be done is to make a method specifically for getting the data of a diagonal tensor. There is actually a nice system for making on-the-fly functions that can easily be applied to "unwrap" and then operate on an ITensor's storage object.

Here's a quick example for DiagReal storage (the storage type of an diagonal-sparse ITensor):

//func is a "lambda" function
auto func = [](DiagReal const& d)
for(auto& el : d.store) print(" ",el);

applyFunc(func,T.store()); //here T is an ITensor with DiagReal storage
//If T has a different kind of storage that func doesn't know how to act on,
//then an exception will be thrown by the above call
commented by (680 points)
Thanks for the nice function.
commented by (52.6k points)
edited by
Glad it is helpful. Eventually I want to have more documentation on the storage system (it's a really interesting and powerful system!). But for now if you want to understand a bit about each storage type you can look at their type definitions in the itdata folder (like itdata/diag.h for the DiagReal and DiagCplx storage types).

Also instead of "func" taking a single type of argument, it can actually be an object that overloads operator() to take multiple different storage types. So this can be used to have different behaviors for various storage types or to distinguish one storage type from another.

The most powerful part of all this is that the storage type is determined dynamically i.e. at run time. So it's similar to Julia's "multiple dispatch" feature but in C++. And it lets the code compile even if the function applied only accepts a subset of the storage types; for the missing ones code is injected at compile time to create run-time exceptions for those cases.
commented by (190 points)
Hi Miles,
I am using your function to obtain the elements of a QDiagReal storage

auto func = [](QDiagReal const& d)
        std::vector<double> spectrum;
        for(auto& el : d.store) {
        return spectrum;


auto spectrum = applyFunc(func,D.store());

I get the elements of the QDiagReal D, but information about the QN is lost.  Is there a simply way to keep also track of the quantum numbers associated with the eigenvalues? For exampl, to get for each quantum number a corresponding list with its values.

For example, D has a structure like this

1<In>QN({"Sz",2}), 1<Out>QN({"Sz",2})
(1,1) 0.1
2<In>QN({"Sz",0}), 2<Out>QN({"Sz",0})
(2,2) 0.2
(3,3) 0.3
1<In>QN({"Sz",-2}), 1<Out>QN({"Sz",-2})
(4,4) 0.4

and I would like to obtain the following list (vectors):

spectrum= [ 0.1, 0.2, 0.3, 0.4]
SZ_QN  = [ -2 , 0, 0 , 2 ]

Thank you!
commented by (52.6k points)
Hi Andreas,
For different questions like this which aren't quite the same as the original one, or are spaced out very far in time, please post a new question. It helps us to find which questions have not yet been answered.

To briefly answer your question, I think there are different ways of getting the information you want but the simplest would be to obtain one of the indices of the ITensor which has this QDiagReal storage inside it, then loop over its subspaces (you can use i.nblock() to find out how many and i.qn(n) to get the QN of the nth one) and get the QN of each, then print these out. Hope that answers your question!

Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.