I am currentlly using openMPS to try and obtain results for mixtures of atoms in a harmonic trap. I noticed that, while I can get the correct results for the spatial distributions (e.g the number of particles per site), the energies are always coming out wrong (or so I believe).
For instance, I am trying to reproduce the results for a single particle in a harmonic trap, with m=1 and \omega=1. So essentially I have one particle (a boson) distributed over a large number of sites. My goal is to reproduce the gaussian ground state expected for the harmonic oscillator.
I fix the lattice potential strength so that the energy in each site corresponds to the energy in the continuum. But somehow the energies come out negative, instead of E=(1/2)\hbar\omega. Can someone help (see code below)?
import MPSPyLib as mps
import numpy as np
import sys
def main(PostProcess=False):
L = 101
t = 0.5
N = 1
k = 12.5/(L-(L+1)/2.)**2
grid=np.linspace(1,L,L)
trappot=k*(grid-(L+1)/2.)**2
Operators = mps.BuildBoseOperators(N)
H = mps.MPO(Operators)
H.AddMPOTerm('bond', ['bdagger', 'b'], hparam='t', weight=-1.0)
H.AddMPOTerm('bond', ['b', 'bdagger'], hparam='t', weight=-1.0)
H.AddMPOTerm('site', 'nbtotal', hparam='trap', weight=1.0)
myObservables = mps.Observables(Operators)
myObservables.AddObservable('site', 'nbtotal', 'number')
myConv = mps.MPSConvParam(max_bond_dimension=30, max_num_sweeps=10)
parameters = [{
# Directories
'simtype':'Finite',
'job_ID':'BoseHubbard_',
'unique_ID':'L_'+str(L)+'N'+str(N),
'Write_Directory':'TMP/',
'Output_Directory':'OUTPUTS/',
# System size and Hamiltonian parameters
'L':L,
't':t,
'trap':trappot,
'verbose':2,
'logfile':True,
# Specification of symmetries and good quantum numbers
'Abelian_generators':['nbtotal'],
'Abelian_quantum_numbers':[N],
'MPSObservables':myObservables,
'MPSConvergenceParameters':myConv
}]
MainFiles = mps.WriteFiles(parameters, Operators, H, PostProcess=PostProcess)
if(not PostProcess):
mps.runMPS(MainFiles, RunDir='')
Outputs = mps.ReadStaticObservables(parameters)
bosenum = Outputs[0]['number']
en = Outputs[0]['energy']
conv = Outputs[0]['converged']
print('pot',trappot)
print('frequency',k)
print('energy',en)
print('convergence',conv)
file=open('data.txt','w')
for j in bosenum:
print>>file, j
file.close()
return
if(name == 'main'):
# Check for command line arguments
Post = False
for arg in sys.argv[1:]:
key, val = arg.split('=')
if(key == '--PostProcess'): Post = (val == 'T') or (val == 'True')
main(PostProcess=Post)
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
before I get to the problem, I realized one thing in your script. You add the tunneling terms [b, bdagger] and [bdagger, b]. openMPS adds the hermitian conjugate of such bond terms automatically. Thus, your tunneling strength is (most likely) twice as big than you expect.
This might contribute actually to your problem. The energy you get depends on the tunneling and the potential term. The tunneling had a negative contribution and outweights the positive contribution from the potential in my opinion. If you want to check what the tunneling contributes, you can measure the correlation [bdagger, b] and extract the first off-diagonal or work with two-site density matrices. You already have access to the energy from the potential via the nbtotal measurment. This way you can track your contributions.
I assume you will continue to get the excited states for this type of system. Once you have those energies, I would look at the actual spectrum and if they follow the spacing you would expect from the harmonic oscillator. Then one can start to shift (add another site term with the identity) and scale (multiply t and trappot with scalar) the energy from OSMPS. I would have to think longer about the problem why the energies do not match from the start, but the spatial basis and discretization are two points coming to my mind right away. And how do you choose the ratio between tunneling strength and potential?
I hope that helps you to some extend and let me know if this is not enough to resolve the issue.
Kind regards,
Daniel
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
You are correct in that an extra term has to be included in the Hamiltonian to reproduce the energies of the continuum. I find that this term has to be proportional to the lattice spacing, that is, the total length of the system divided by the number of sites.
However, how do I obtain the eigenstates of the system? I am able to calculate the static observables, such as total number of particles per site, and also the energy, but I'd like to look at the actual vectors for the ground and excited states.
Best regards,
Rafael
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
the only case where you could get the complete state vector from the MPS is for small systems below 20 sites, if you stretch it, you can go up to around 30. For this reason, such conversion is not implemented at the moment. You have 101 sites and a local dimension of 2 (0 or 1 boson per site) in your example above and that results in a state vector of 2.5e+30 entries. If each entry is a 8 byte entry (64bit), you need 2.02e+12 Terrabyte of memory.
I think the closest you get is the spatial distribution, i.e., measurement of nbtotal, where you said before that it represents the spatial representation of some continuum wave function psi(x). Most likely, up to some normalization since nbtotal is already the probability. Now, you may argue that in psi(x) you can have complex entries and some phase relation between two points in space once you do time evolution, which is not captured by nbtotal. (For the ground/excited state, a completely real Hamiltonian yields a ground state with real entries only.) I don't know the exact answer how to capture that from MPS. If I had to come up with an algorithm right now, I could point out one idea, but involves iterations over an enormous number of configurations; fortunately, it would be bounded by the bond dimension. I can think about it and see if it is feasible. I don't think it will be different from nbtotal for ground/excited states of real Hamiltonians.
Let me know if that information helps you in any way.
Kind regards,
Daniel
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
I see. But you say that even if I chose a small system, with just a few sites (or a system spanning a smaller Hilbert space, such as a short spin chain) I cannot look at the eigenvectors?
Rafael
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
actually, the Hilbert space in your example would not be the size I said; I did not consider the symmetry which reduces it, of course, to 101 basis states. And you could solve the same system with exact diagonalization. Hopefully, I will find a script for the ED setup by Monday, most of the syntax is the same as for MPS.
Due to the limitation without symmetry, the subroutines to get the state vector from an MPS are not implemented. I will put it on the development to-do list, you convinced me that it has some applications beyond testing and debugging.
Kind regards,
Daniel
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The answer is now going from MPS methods to ED (exact diagonalization) methods, but if it helps you to solve your problem, that's fine to me.
I attached below the script to run with ED, which is slower, but the current size can run on my laptop. I highlighted the lines you have to change/add with ###. I had to fix one issue and updated the library via svn.
With ED, you have a much better hold of the state vectors if you need them for the single particle problem on a discretized space. Saving out the states is no default method. You can go to the function def statics_simulation in MPSPyLib/EDLib/exactdiag.py. You could write out elem[0] in the for loop with np.save().
Maybe the optimal solution as long as you do not want to time evolution: you can use the EDLib in openMPS to construct the Hamiltonian matrix directly from your MPO, e.g., Hmat = H.build_hamiltonian(parameters[0], SymmSec, False, maxdim) without calling EDLib or MPS directly. (False for no sparse matrix, maxdim is the maximal dimension of your Hilbert space.) Then, you have all the eigenvectors after an eigenvalue decomposition with numpy/scipy.
Let me know if there are any follow-up questions.
Kind regards,
Daniel
P.S. If you are not familiar with svn to get the newest revision, I can dump the fixed exactdiag.py here.
Thank you for the code. That will be really helpful. I assume I can use it for other hamiltonians as well (I am also working with the Heisenberg chain).
However, I am also interested in the time evolution. Obtaining the single particle states is just a step for me to study quench dynamics in a system with a few interacting particles (then I'd be using the Bose Hubbard model to simulate the continuum). I have recently started looking at the time evolution of the single particle after a quench in the trapping potential. I have run into some issues there, so I'll probably open another topic regarding that.
Thanks again,
Rafael
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Hello there,
I am currentlly using openMPS to try and obtain results for mixtures of atoms in a harmonic trap. I noticed that, while I can get the correct results for the spatial distributions (e.g the number of particles per site), the energies are always coming out wrong (or so I believe).
For instance, I am trying to reproduce the results for a single particle in a harmonic trap, with m=1 and \omega=1. So essentially I have one particle (a boson) distributed over a large number of sites. My goal is to reproduce the gaussian ground state expected for the harmonic oscillator.
I fix the lattice potential strength so that the energy in each site corresponds to the energy in the continuum. But somehow the energies come out negative, instead of E=(1/2)\hbar\omega. Can someone help (see code below)?
import MPSPyLib as mps
import numpy as np
import sys
def main(PostProcess=False):
L = 101 t = 0.5 N = 1 k = 12.5/(L-(L+1)/2.)**2 grid=np.linspace(1,L,L) trappot=k*(grid-(L+1)/2.)**2 Operators = mps.BuildBoseOperators(N) H = mps.MPO(Operators) H.AddMPOTerm('bond', ['bdagger', 'b'], hparam='t', weight=-1.0) H.AddMPOTerm('bond', ['b', 'bdagger'], hparam='t', weight=-1.0) H.AddMPOTerm('site', 'nbtotal', hparam='trap', weight=1.0) myObservables = mps.Observables(Operators) myObservables.AddObservable('site', 'nbtotal', 'number') myConv = mps.MPSConvParam(max_bond_dimension=30, max_num_sweeps=10) parameters = [{ # Directories 'simtype':'Finite', 'job_ID':'BoseHubbard_', 'unique_ID':'L_'+str(L)+'N'+str(N), 'Write_Directory':'TMP/', 'Output_Directory':'OUTPUTS/', # System size and Hamiltonian parameters 'L':L, 't':t, 'trap':trappot, 'verbose':2, 'logfile':True, # Specification of symmetries and good quantum numbers 'Abelian_generators':['nbtotal'], 'Abelian_quantum_numbers':[N], 'MPSObservables':myObservables, 'MPSConvergenceParameters':myConv }] MainFiles = mps.WriteFiles(parameters, Operators, H, PostProcess=PostProcess) if(not PostProcess): mps.runMPS(MainFiles, RunDir='') Outputs = mps.ReadStaticObservables(parameters) bosenum = Outputs[0]['number'] en = Outputs[0]['energy'] conv = Outputs[0]['converged'] print('pot',trappot) print('frequency',k) print('energy',en) print('convergence',conv) file=open('data.txt','w') for j in bosenum: print>>file, j file.close() returnif(name == 'main'):
# Check for command line arguments
Post = False
for arg in sys.argv[1:]:
key, val = arg.split('=')
if(key == '--PostProcess'): Post = (val == 'T') or (val == 'True')
Hello Rafael,
before I get to the problem, I realized one thing in your script. You add the tunneling terms [b, bdagger] and [bdagger, b]. openMPS adds the hermitian conjugate of such bond terms automatically. Thus, your tunneling strength is (most likely) twice as big than you expect.
This might contribute actually to your problem. The energy you get depends on the tunneling and the potential term. The tunneling had a negative contribution and outweights the positive contribution from the potential in my opinion. If you want to check what the tunneling contributes, you can measure the correlation [bdagger, b] and extract the first off-diagonal or work with two-site density matrices. You already have access to the energy from the potential via the nbtotal measurment. This way you can track your contributions.
I assume you will continue to get the excited states for this type of system. Once you have those energies, I would look at the actual spectrum and if they follow the spacing you would expect from the harmonic oscillator. Then one can start to shift (add another site term with the identity) and scale (multiply t and trappot with scalar) the energy from OSMPS. I would have to think longer about the problem why the energies do not match from the start, but the spatial basis and discretization are two points coming to my mind right away. And how do you choose the ratio between tunneling strength and potential?
I hope that helps you to some extend and let me know if this is not enough to resolve the issue.
Kind regards,
Daniel
Hello Daniel,
Thank you again for the reply.
You are correct in that an extra term has to be included in the Hamiltonian to reproduce the energies of the continuum. I find that this term has to be proportional to the lattice spacing, that is, the total length of the system divided by the number of sites.
However, how do I obtain the eigenstates of the system? I am able to calculate the static observables, such as total number of particles per site, and also the energy, but I'd like to look at the actual vectors for the ground and excited states.
Best regards,
Rafael
Hello Rafael,
the only case where you could get the complete state vector from the MPS is for small systems below 20 sites, if you stretch it, you can go up to around 30. For this reason, such conversion is not implemented at the moment. You have 101 sites and a local dimension of 2 (0 or 1 boson per site) in your example above and that results in a state vector of 2.5e+30 entries. If each entry is a 8 byte entry (64bit), you need 2.02e+12 Terrabyte of memory.
I think the closest you get is the spatial distribution, i.e., measurement of nbtotal, where you said before that it represents the spatial representation of some continuum wave function psi(x). Most likely, up to some normalization since nbtotal is already the probability. Now, you may argue that in psi(x) you can have complex entries and some phase relation between two points in space once you do time evolution, which is not captured by nbtotal. (For the ground/excited state, a completely real Hamiltonian yields a ground state with real entries only.) I don't know the exact answer how to capture that from MPS. If I had to come up with an algorithm right now, I could point out one idea, but involves iterations over an enormous number of configurations; fortunately, it would be bounded by the bond dimension. I can think about it and see if it is feasible. I don't think it will be different from nbtotal for ground/excited states of real Hamiltonians.
Let me know if that information helps you in any way.
Kind regards,
Daniel
Hi Daniel,
I see. But you say that even if I chose a small system, with just a few sites (or a system spanning a smaller Hilbert space, such as a short spin chain) I cannot look at the eigenvectors?
Rafael
Hello Rafael,
actually, the Hilbert space in your example would not be the size I said; I did not consider the symmetry which reduces it, of course, to 101 basis states. And you could solve the same system with exact diagonalization. Hopefully, I will find a script for the ED setup by Monday, most of the syntax is the same as for MPS.
Due to the limitation without symmetry, the subroutines to get the state vector from an MPS are not implemented. I will put it on the development to-do list, you convinced me that it has some applications beyond testing and debugging.
Kind regards,
Daniel
Hello Rafael,
The answer is now going from MPS methods to ED (exact diagonalization) methods, but if it helps you to solve your problem, that's fine to me.
I attached below the script to run with ED, which is slower, but the current size can run on my laptop. I highlighted the lines you have to change/add with ###. I had to fix one issue and updated the library via svn.
With ED, you have a much better hold of the state vectors if you need them for the single particle problem on a discretized space. Saving out the states is no default method. You can go to the function def statics_simulation in MPSPyLib/EDLib/exactdiag.py. You could write out elem[0] in the for loop with np.save().
Maybe the optimal solution as long as you do not want to time evolution: you can use the EDLib in openMPS to construct the Hamiltonian matrix directly from your MPO, e.g., Hmat = H.build_hamiltonian(parameters[0], SymmSec, False, maxdim) without calling EDLib or MPS directly. (False for no sparse matrix, maxdim is the maximal dimension of your Hilbert space.) Then, you have all the eigenvectors after an eigenvalue decomposition with numpy/scipy.
Let me know if there are any follow-up questions.
Kind regards,
Daniel
P.S. If you are not familiar with svn to get the newest revision, I can dump the fixed exactdiag.py here.
Hi Daniel,
Thank you for the code. That will be really helpful. I assume I can use it for other hamiltonians as well (I am also working with the Heisenberg chain).
However, I am also interested in the time evolution. Obtaining the single particle states is just a step for me to study quench dynamics in a system with a few interacting particles (then I'd be using the Bose Hubbard model to simulate the continuum). I have recently started looking at the time evolution of the single particle after a quench in the trapping potential. I have run into some issues there, so I'll probably open another topic regarding that.
Thanks again,
Rafael