Menu

Heisenberg chain with symmetries

Technical
2018-03-27
2018-05-30
  • Rafael Barfknecht

    Hello everyone,

    I am currently trying to calculate ground state properties in the Heisenberg chain. I am starting with a small chain (L=3,4) and calculating the expected value of Sz (see code below). I have run into two problems:

    1. Upon running the code twice, even as I clear the folder /TMP and /OUTPUTS, I get different results. They seem to be internally consistent (the ratio between the projections at different sites seem correct) but the numbers differ a lot.
    2. I then try to add symmetries (I want to fix the magnetization by choosing a total value for Sz - ). For a case of L=3, I fix Sz=1/2, which in my understanding would mean two up spins and one down spin. Then I get the following message:

    Exception: The shifted and scaled quantum number 1.5 computed from (0--1.5)/1.0
    is not an integer and so does not represent a valid quantum number!
    If you are certain this is an error file a bug report!

    Even when I try to add symmetries in the L=2 case, in which I would fix Sz as an integer (Sz=0 or Sz=1), I get

    STOP Infinite system algorithm cannot be used for systems with less than 4 sites!
    Traceback (most recent call last):
    File "Heisenberg.py", line 75, in <module>
    main(PostProcess=Post)
    File "Heisenberg.py", line 53, in main
    Outputs = mps.ReadStaticObservables(parameters)
    File "/usr/local/lib/python2.7/dist-packages/MPSPyLib/Obs.py", line 1204, in ReadStaticObservables
    suppress_errors)
    File "/usr/local/lib/python2.7/dist-packages/MPSPyLib/Obs.py", line 1287, in openfile_getunitobservable
    Obsfilebase + '_' + str(whichconv) + '.dat')
    IOError: Result file was missing: OUTPUTS/Heisenberg_L_2ObsOut0_1.dat

    I don't see where I am stating that I want to simulate an infinite system.

    Any help?

    Cheers,

    Rafael Barfknecht

    ~~~
    import numpy as np
    import sys

    def main(PostProcess=False):

    L=3
    
    # Build operators
    Operators = mps.BuildSpinOperators(0.5)
    
    Operators['sx']=0.5*(Operators['splus']+Operators['sminus'])
    Operators['sy']=(-0.5j)*(Operators['splus']-Operators['sminus'])
    
    # Define Hamiltonian MPO
    H = mps.MPO(Operators)
    H.AddMPOTerm('bond', ['sx', 'sx'], hparam='J', weight=1.0)
    H.AddMPOTerm('bond', ['sy', 'sy'], hparam='J', weight=1.0)
    H.AddMPOTerm('bond', ['sz', 'sz'], hparam='J', weight=1.0)
    
    # Observables
    myObservables = mps.Observables(Operators)
    myObservables.AddObservable('site', 'sz', 'zproj')
    
    myConv = mps.MPSConvParam()
    
    parameters = [{ 
        # Directories
        'simtype':'Finite',
        'jobID':'Heisenberg',
        'uniqueID':'L'+str(L),
        'WriteDirectory':'TMP/', 
        'OutputDirectory':'OUTPUTS/', 
        # System size and Hamiltonian parameters
        'L':L,
        'J':1,
        'verbose':2,
        'logfile':True,
        # Specification of symmetries and good quantum numbers
        #'Abeliangenerators':['sz'],
        #'Abelianquantumnumbers':[1/2],
        'MPSObservables':myObservables,
        'MPSConvergenceParameters':myConv
    }]
    
    # Write Fortran-readable main files
    MainFiles = mps.WriteFiles(parameters, Operators, H, PostProcess=PostProcess)
    
    #if(not PostProcess):
    # Run the simulations
    mps.runMPS(MainFiles, RunDir='')
    
    Outputs = mps.ReadStaticObservables(parameters)
    
    energ = Outputs[0]['energy']
    zproj = Outputs[0]['zproj'] 
    convg = Outputs[0]['converged']
    
    #state = Outputs[0]['totalspin']
    
    print('total energy',energ)
    print('z projection',zproj)
    print('convergence',convg)
    
    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')

    # Run main function
    main(PostProcess=Post)**
    

    ~~~*

     
  • Daniel Jaschke

    Daniel Jaschke - 2018-03-28

    Hello Rafael,

    I'll start with the IMPS error and then move to the symmetry:

    1) Yes indeed, you are simulating a finite-size system. But the IMPS algorithm is used even for finite-size system to grow the system to the target size and then start the variational algorithm for fixed L. This initial guess from IMPS is better than just filling tensors with random numbers. The minimal size of four sites is difficult to overcome and was usually not a problem for the many-body applications. Either start with 4 sites, or extract the ground state with the built-in ED methods.

    2) Symmetry: If you are using python2, 1/2 = 0 and zero is indeed not a valid quantum number for three sites. Try [0.5] written out or [1. / 2]. From your output I assume that is the problem.

    3) I get different results when rerunning, but they are all numerical zero. Do you have different values above 1e-14?

    You can also represent the y-y-interaction with real operators. That allows MPS to use real tensors throughout the ground state calculations at the cost of increasing the bond dimension of the MPO by one. Might give a speed-up for larger systems, but I never compared the run times.

    Kind regards,

    Daniel

    Real y-y-interaction (bond term takes care of hermitian conjugate, see also weight=0.5)

    # Build operators
    Operators = mps.BuildSpinOperators(0.5)
    
    Operators['sx'] = 0.5 * (Operators['splus']+Operators['sminus'])
    Operators['sy1'] = 0.5 * np.array([[0, -1], [1, 0.]])
    Operators['sy2'] = 0.5 * np.array([[0, 1], [-1, 0.]])
    
    # Define Hamiltonian MPO
    H = mps.MPO(Operators)
    H.AddMPOTerm('bond', ['sx', 'sx'], hparam='J', weight=1.0)
    H.AddMPOTerm('bond', ['sy1', 'sy2'], hparam='J', weight=0.5)
    H.AddMPOTerm('bond', ['sz', 'sz'], hparam='J', weight=1.0)
    
     
  • Rafael Barfknecht

    Hi Daniel,

    Thank you for the reply.

    1) OK, my idea was to do some cases with just a few spins, but I can check with larger numbers as well.

    2) and 3) I do get different numerical values which are different from zero, if the total number of sites is odd (if L is even, then it's all numerically zero). For instance, upon running the code twice for L=5 (no symmetries enforced), I get

    ('total energy', -1.92788625331799)
    ('z projection', array([-0.13529511, 0.07786379, -0.14955828, 0.07786379, -0.13529511]))
    ('convergence', 'T')

    and then

    ('total energy', -1.927886253318)
    ('z projection', array([-0.25570358, 0.14716017, -0.28266054, 0.14716017, -0.25570358]))
    ('convergence', 'T')

    The results are internally consistent in that they seem to represent an antiferromagnetic ground state, but I cannot figure out why the values of 'sz' are different.

    If I add symmetries, by including

    'Abelian_generators':['sz'],
    'Abelian_quantum_numbers':[0.5],

    I get convergence problems:

    ('total energy', -11.3522592142027)
    ('z projection', array([ 0.28724895, -0.02540869, 0.07286293, 0.27230301, -0.29197794]))
    ('convergence', 'F')

    In these attempts I am using the real y-y operator as you described above.

    Best regards,

    Rafael

     
  • Daniel Jaschke

    Daniel Jaschke - 2018-04-09

    Hello Rafael,

    whenever you have (nearly) degenerate ground states, the ground state found by MPS can be a superposition of these (nearly) degererate ground states. This leads to different values for the observables, each consistend within one simulations.

    Now, the symmetries can prevent this effect, if the degenerate states are in two different symmetry sectors. I checked on the non-converging simulation. The convergence flag is only a very quick check, and the variance tells you a more detailed story. In your case, the simulation was not close to being converged and I will look into the problem.

    Best regards,

    Daniel

     
  • Daniel Jaschke

    Daniel Jaschke - 2018-05-30

    Hello Rafael,

    Sorry, it took me some time to get back to you. Symmetries can get quite technical which was necessary to check in all directions here. I was trying to find the reason why the current code fails, but it is difficult.

    Instead, I have another improved version of the Hamiltonian, i.e., less terms. The key is to express both sx_i sx_j + sy_i sy_j as spin raising and lowing operators and j = i + 1. This step leads to a bond term with splus_i + sminus_j + h.c.; thus, the bond dimension of our MPO is reduced by one in comparison to the previous approach and still keeping it real-valued. My test case has no degeneracy and the magnetization for the ground state stays constant over runs.

    Let me know if you have any problems with this approach.

    Best regards,

    Daniel

    P.S. Some more details for developers just in case we have to get back to this or a similar case: By adding sx_i sx_j and sy_i sy_j individually, there are terms with +/- splus_i splus_j and +/- sminus_i sminus_j. The +/- case conserves the symmetry together, but not on its own. So OSMPS has to rely on the fact that they cancel out during the calculation, while creating intermediate states which violate the symmetry. This intermediate violation of the symmetry is most likely causing the problems.

     

Log in to post a comment.

Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.