Is it possible to export a state (density matrix) in matrix form in a given basis?

Question

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?

Comments

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)?

---

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.

Answer 1

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)

Comments

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!

---

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.

---

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

---

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.

---

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.

---

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);
println();
};

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

---

Thanks for the nice function.

---

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.

---

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) {
            spectrum.push_back(el);
        }
        return spectrum;
};

with

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

QDiagReal
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!

---

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!

Miles