sc.exponential() raises error and other questions

[sajjan02purdue]

Hi Jose, This is an amazing package and I have made it a point to let all my collaborators know about it. Kudos for making such nicely constructed !! I have a couple of questions which are appended below. I must say that i do not have much experience in DMRG but I am just beginning to find my way through.

1) I am interested in computing String correlation function for S=1 Haldane chains. I am using the exponential attribute for an sc object as given in example (https://github.com/joselado/dmrgpy/blob/master/examples/exponential_EV/main.py). My code is as follows

#String correlator
#******************
def String_correlation_fn(self):
         String_correlator=[]
         item='z'
         Op_S={'x':sc.Sx, 'y':sc.Sy, 'z':sc.Sz}
         ground_st = (sc.get_gs(mode=self.mode)).normalize()   

         def String_sum(index, item):
               op_string = Op_S[item][1]
               for k in np.arange(2, index):
                   op_string += Op_S[item][k]
               return op_string 

          if item in Op_S.keys():

                 for i in np.arange(self.no_spins):
                      if i in [0,1]:
                          String_correlator.append(ground_st.dot(Op_S[item][0]*Op_S[item][i]*ground_st))
                      else:    
                          String_correlator.append(ground_st.dot(Op_S[item][0]*sc.exponential(1j*np.pi*String_sum(i, item), Op_S[item][i]*ground_st, mode=self.mode))) 

I can produce the results for N=4, 6 spins but the code returns an error saying " No active exception to re-raise" for N>6 onwards. It has to do with maxsize in algebra.py. Do you have any suggestions?

2) Is it possible to get excited state wavefunctions (say the first or second excited state) like we have access to ground state through the attribute sc.get_gs()? I dont yet see any examples doing that hence the question.

3) Is it possible to get the overlap of the wavefunction (both ground and excited state of choice) over the constituent basis states. Say for example for N=4 , S=1 case I would want the contribution of the basis (0,+1, -1, 0)^T in the ground singlet state etc from DMRG and ED? What is the way to do that. One option would be to use the sc.overlap () attribute but the basis vector of choice ([0, +1, -1, 0)^T here) would then have to be defined as an sc object externally by the user . What is protocol to do it ?

4) For smaller N (<10) it is possible to get the Hamiltonian matrix ? Right now sc.set_hamiltonian() produces an object with an address pointer . Can it be converted to a dense matrix object ? For example the way one does with scipy etc ? like scipy. sparse.csr_matrix.todense() for a sparse csr_matrix

5) I havent yet tried time evolution but I am going to soon with superposition of DMRG eigenstates . Do the imaginary_evolution() attribute as given in example https://github.com/joselado/dmrgpy/blob/master/examples/evolve_wavefunction/main.py has any restriction on size of N (spins) as the sc.exponential() attribute in point 1). I am talking about S=1 case for close to N=10-25 spins

It would be extremely beneficial if you have any suggestions about the aforementioned queries ? Thanks a lot in advance and also for the excellently written package.

[joselado]

Dear Manas Sajjan,

thank you for the kind feedback! It is wonderful to hear that the software is useful!

I address your questions below

1) I think that the issue is that you are running the calculations with exact diagonalization (mode="ED"), which becomes too computationally expensive after N>6 (the code raises error at that point). If you use mode="DMRG" the calculation should run without issues for any N

2) Yes, you can use the function sc.get_excited_states(), which returns both the energies and wavefunctions. You can see an example in https://github.com/joselado/dmrgpy/tree/master/examples/excited_states

3) Yes, to compute the overlap you would need to build the wavefunction onto which you want to project. A simple way of doing it would be to define a Hamiltonian whose ground state is exactly that one, and obtain it with sc.get_gs() (for example, by defining a dummy Hamiltonian that only has magnetic fields)

4) Yes, it would be possible, although currently it is not implemented in the main object. If this is something that would be important for you let me know and I can include it.

5) If you do time evolution with the DMRG mode in principle there is no limitation on the number of sites (apart from the intrinsic accuracy of DMRG, of course). It is certainly possible to reach N=25 sites with DMRG, but of course the calculations are more expensive that just computing the ground state.

Thanks again for the nice feedback!

Best regards, Jose

[sajjan02purdue]

Hi Jose, Thanks a lot for your quick response. With your answers I have done the following:

1) Thanks for pointing it out. Yes I do see, if the mode is 'DMRG', there are no errors . However just to clarify this issue a bit , for mode='ED' one can obtain good results (energies of ground and excited states) upto N=10 for S=1 case. For the computation of the String correlation (code given above) I see the error is raised with mode='ED' even for N=8 even though the energies and other things are obtained fine. Is this normal ?

2)Thanks for this. I tried this and It works perfectly fine. This question is resolved

3) Thanks for this. I need to try this . I will do so and may be let you know if any major hiccups arises.

4) Thanks for this. It would be great if this is available as there is some other analysis on the hamiltonian we can do. Let me know whenever the feature is available in future.

5) Thanks. I am trying this now. Will keep you updated.

Also there is a 6th point which I forgot to mention. Can you point me to one or two papers from which the basic structure of the algorithm as used in the package is adapted.? The reason I ask is that there is likely to be many variants of DMRG now (I am no expert but that is usually the case for most techniques) so it would help to know what is the specific algorithm that is being used at the backend for mode='DMRG' in this package ? Also what is the format you would prefer for citing the package in academic work ?

Thanks once again for the responses and all the help.

[joselado]

Dear Manas Sajjan,

thank you for your response! I address some of the points below (I will also address the others later)

4 I have now added this feature, you can see an example in https://github.com/joselado/dmrgpy/tree/master/examples/full_hamiltonian

6. The details of the algorithm are very nicely explained in https://arxiv.org/abs/2007.14822

for citing the package, you could refer both to the paper above (which is the library in which dmrgpy is built on), and if you would like to also cite dmrgpy you can just refer to the github webpage (https://github.com/joselado/dmrgpy)

Best regards, Jose

[sajjan02purdue]

Dear Jose, Thanks once again. The responses are very helpful especially point 4) above.

With Regards, Manas

[sajjan02purdue]

Hi Jose, I have noticed the following while working with your package lately :

1) The functions site_entropy() and pair_entropy() inside entropy.py dont quite work for S=1 spinchains ? Besides even for S=1/2 for wf mode='ED', the function call fails . It only works for mode='DMRG'. I am trying to access the 1 -RDM (one-particle density matrix) and 2-RDM (two-particle density matrix) for S=1 spinchains . So I was wondering is there a way to bypass this and/or what needs to be done to retrieve the reduced density matrices? I can compute any body or two properties after that using those like (but not limited to) site_entropy, pair_entropy. Thus accessing 1-RDM, 2-RDM for both S=1 and S=1/2 chains specifically would be immensely beneficial to me

2) Also I was trying to compute the dynamical properties of the system following a Hamiltonian ramp of the form H(s) = (1-s)H_1 + s H_2 where s= f(t) for S=1 spinchains in DMRG and ED mode. However the evolve_wf method could not handle it ? It works fine for time independent hamiltonians, however in my case H(s) is explicitly time dependant . I see there is timedependant.py module which computes overlap with initial state and dynamical correlators , however i cant see any function yielding the full wf object directly to access all properties at any time . What needs to be done in this case to access the wavefunction (wf) in time so as to compute various properties thereafter ? Is there any other way of doing it that i am missing ?

It will be really helpful if you can shed some light on the above two questions .

Many thanks in advance

With Regards, Manas

[joselado]

Hi Manas,

thanks a lot for the two comments, I address them below.

1- You are totally right, that function still does not work in all the cases. So far it only works for S=1/2 with DMRG. I should have put a warning in the function that is still not fully functional, and that it only works in specific cases. 2- For this case I would say that the way to proceed is that you recompute the ground state for each "s", and compute your dynamical correlator with that computed ground state for each H(s). Time evolution currently only works for time-independent Hamiltonians. I would think that, even it worked for time dependent, depending on task you have in mind it may be more convenient to recompute the ground state for each s, as otherwise if you perform time evolution you will likely push you state away from the ground state manifold (especially if your time evolution is faster that the lowest excited state).

Thanks again for the comments! Once I finish with the full implementation of your point 1) I will let you know!

Best regards, Jose

[sajjan02purdue]

Hi Jose, Thanks a lot for the answers. They are very helpful. I will look forward to hearing from you when point 1) (accessing 1-RDM, 2-RDM for S=1) is implemented.

With Regards, Manas

[sajjan02purdue]

Hi Jose, I wanted to touch base on point 1) (1-RDM and 2-RDM for S=1 in DMRG, ED) above on Feb for which you mentioned that you would get back to me once its fully implemented. I was wondering if the current version of the package support the functionality ? I will wait in anticipation to hear from your end. It would really benefit my research immensely

With Regards, Manas

[joselado]

Hi Manas,

thanks for the follow-up. That functionality is not implemented in code at this stage, yet of course it is the todo list. I will let you know once it is implemented.

Best regards, Jose

[sajjan02purdue]

Hi Jose, I wanted to follow-up on your last message. Even if all the entropic calculations are not fully implemented, just getting access to the 1-RDM and 2-RDM (1st and 2nd order reduced density matrices) for S=1 systems in DMRG/ED would be good enough. Will you please let me know if there is a way to retrieve that from the codebase by some means. It would really benefit me if this functionality is present. I will wait to hear from your end. Thanks a lot once again

With Regards, Manas

[joselado]

Hi Manas,

thank you for the follow up message. Formally, it would be possible to extract such reduced density matrices by explicitly building MPO projectors to each local sector and applying then to the MPS of interest. This can be directly done at the Python level with the MPO-MPS algebra, so it could be perhaps the easiest way to get that functionality with the current version of the code.

Best regards, Jose