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

+1 vote

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 May 17 by (16,560 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 May 18 by (160 points)
edited May 18 by Andreas
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.

answered May 18 by (16,560 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 May 20 by (650 points)
Hi Miles,

I tried the method and it works perfectly well for IQTensor. However, for QDiagReal I got the error message
What is the correct way of doing this? Thanks!
commented May 20 by (16,560 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 May 20 by (650 points)
edited May 20 by chengshu
Thanks very much for the quick reply. Using the new version, the toITensor works well, but the "orderedC" function gives an error message:
The function "ordered" gives another error message:
applyFunc: function object has no operator() method for storage type
commented May 21 by (16,560 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 May 21 by (650 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 May 21 by (16,560 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);
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
commented May 22 by (650 points)
Thanks for the nice function.
commented May 22 by (16,560 points)
edited May 22 by miles
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.