Ground state wavefunction

[Faridun]

Dear all,

I try to obtain ground state wavefunction using dmrg:

energy, psi = dmrg(H, psi0, sweeps)
println(psi)

Is there any way to obtain array like wavefunction for further calculations (for example to calculate fractal dimensions)?
Also, I tried psi.data This gives NDTensors.BlockSparse and I do not know how to manage with this data.

Thank you in advance.

[kmp5]

The returned value psi is the DMRG optimized wavefunction in an MPS form. If you would like to turn the wavefunction into a single order-N tensor (where N is the number of sites) you have to contract over all of the bond dimensions. This can be done like this

ψ = contract(psi)

however it will be very large (if site size = 2, 30 sites will create a tensor of dimension 2^30 which is 8Gb of storage in Float64)
Note also you are finding that psi.data is a BlockSparse because you have the variable conserve_qns=true in the siteinds function.

[mtfishman]

@Faridun

Also, I tried psi.data This gives NDTensors.BlockSparse and I do not know how to manage with this data.

I’m confused why you say psi.data returns NDTensors.BlockSparse, it should be Vector{ITensor} (a vector of the tensors of the MPS).

[Faridun]

Thank you so much for your reply!
Yes, a data is very large. Is there any way to obtain the value psi, when I have the variable conserve_qns=true in the siteinds function?

[kmp5]

@mtfishman is right I actually didn’t catch that psi.data does return a Vector{ITensors}

julia> typeof(psi.data)
Vector{ITensor} (alias for Array{ITensor, 1})

potentially you are looking at the output of psi.data and seeing something like this

ITensor ord=2
Dim 1: (dim=3|id=146|"S=1,Site,n=5") <Out>
 1: QN("Sz",2) => 1
 2: QN("Sz",0) => 1
 3: QN("Sz",-2) => 1
Dim 2: (dim=3|id=479|"Link,l=4") <In>
 1: QN("Sz",-2) => 1
 2: QN("Sz",0) => 1
 3: QN("Sz",2) => 1
NDTensors.BlockSparse{Float64, Vector{Float64}, 2}
 3×3
Block(3, 1)
 [3:3, 1:1]
 1.0

Block(2, 2)
 [2:2, 2:2]
 1.0

Block(1, 3)
 [1:1, 3:3]
 1.0]

which is an ITensor with a BlokcSparse storage structure.
Would you please show more lines of code so we can know what is going on.
Also if psi is in fact a Vector{ITensor} you can create the full non-mps wavefunction with the code
contract(psi)

[Faridun]

Sure!
Here is full code

using ITensors
N = 12
Npart = 6
t1 = 1.0
U = 0.5
h = 1.0
sites = siteinds("Electron", N; )

ampo = AutoMPO()
for f = 1:N-1
    ampo += -t1, "Cdagup", f, "Cup", f+1
    ampo += -t1, "Cdagup", f+1, "Cup", f
    ampo += -t1, "Cdagdn", f, "Cdn", f+1
    ampo += -t1, "Cdagdn", f+1, "Cdn", f
end

for i = 1:N
    ampo += U, "Nupdn", i
end
    # Magnetic field (Zeeman) term
for i = 1:N
    hz = 0.1 * (rand() )  # [-1, 1]
        #println(hz)
    ampo += hz, "Sz", i
end

    # Potential field term
for i = 1:N
    hE = rand() - 0.5  # [-0.5, 0.5]
        #println(hE)
    ampo += hE, "Ntot", i
end


H = MPO(ampo, sites)
state = [isodd(n) ? "Up" : "Dn" for n = 1:N]
psi0 = randomMPS(sites, state)
@show flux(psi0)

sweeps = Sweeps(6)
maxdim!(sweeps, 50, 100, 200, 400, 800, 800)
cutoff!(sweeps, 1E-12)

@show sweeps
energy, psi = dmrg(H, psi0, sweeps)

[Faridun]

Here you are

[kmp5]

I don’t see anything related to BlockSparse in this code. And adding something like
psi_dense = contract(psi)
does successful construct the order-12 representation of psi I just tested.

[mtfishman]

@Faridun

Is there any way to obtain array like wavefunction for further calculations (for example to calculate fractal dimensions)?

Could you explain what might be involved in computing fractal dimensions? It’s possible that you could do it with an MPS, I think you just need to give us some more details about what kind of calculation you are trying to perform and what operations you might need.

[Faridun]

Thank you for your reply. I mean something like:

S_q = \ln(\sum_s^{N_H}|\psi(s)|^2)

where N_H is dimensional Hilbert space. It means the logarithm of the square of the sum of all elements of wavefunction.

[miles]

If I understand your formula correctly, the value inside the log is the squared norm of the wavefunction, which as you know of course is usually taken to equal 1.0 by convention (and certainly our DMRG code and other codes will generally output a state normalized to 1). So then the above quantity would be ln(1) = 0 unless I am missing something such as about what the sum runs over.

What I would suggest for you is to test your calculation idea out on a small system just by finding the ground state (if that’s what you’re after) in a direct (‘exact diagonalization’ or ‘full state’) manner, to test your observable and whether it gives the result you want. Then if it does you can next scale up by using MPS and DMRG.

[Faridun]

That’s right. What if I want to calculate:

Sq=\ln(\sum_s^{N_H}|\psi(s)|^4)

In this case we need, what we should to do?

[Faridun]

In case of degree equal 4, is any way to use MPS to calculate?

[miles]

I believe this is possible, but it is not straightforward to do. The way I could think to do it would be to use a kind of machine learning or tomography method to “learn” the state with amplitudes that are the square of the original ones. Then by squaring that new state again one can get the sums of the fourth powers. One such tomography method is called “tensor cross interpolation”.

Another option is to square the MPS using the network for squaring functions described in this paper:
https://journals.aps.org/prx/abstract/10.1103/PhysRevX.13.021015

[Faridun]

Thank you so much for you reply. It was really useful. I found another property which may be calculated. It can be written as follows:
D_1 = \sum_{i}^N |c_i|^2 \ln |c_i|^2

where c_i is an expansion coefficient of a given state \ket{\psi} in some (finite) orthonormal basis {\ket{\psi_i}}, i.e. \ket{\psi} = \sum_{i}^N c_i {\ket{i}}.

Any ideas how to perform this for the system with number of sites L>100 ?

[miles]

I have an answer to this you might find interesting, as I think it highlights the basis dependence of the D_1 quantity.

Any MPS, when you bring it into “orthogonal form” and then if you SVD the ‘center’ tensor to expose the singular values matrix on a bond, is equivalent to a Schmidt decomposition:

\ket{\psi} = \sum_i s_i \ket{L_i} \ket{R_i}

You can think of defining basis states \ket{i} = \ket{L_i} \ket{R_i} and then the expansion is of the form you are looking for.

Then the D_1 quantity you defined, when computed in this basis, would be exactly the von Neumann entanglement entropy, perhaps up to a minus sign. And then while this quantity is equal for a large system on any bond deep in the “bulk”, near the edges it depends strongly on which bond you compute it across.

If you are interested in computing D_1 in this basis, since it is equivalent to entanglement entropy, please see our Code Example on how to compute entanglement:
https://itensor.github.io/ITensors.jl/dev/examples/MPSandMPO.html#Computing-the-Entanglement-Entropy-of-an-MPS