core dumped issue on DMRG calculation with Itenso (c++ version3)

Dear ITensor

I am working on a frustrated antifermagnetic XX model on the 2d honeycomb lattice with nearest- and next-nearest neighbour interaction. Using ITensor, I was able to study lattice up to N=60 sites, but going to a bigger lattice size, N=90, I got a “core dumped” message on my pc with 96GB ram. I asked a friend to run the code with 500GB accessible memory, but again we received a “core dumped” error. There are 370 bound links( total of nn + nnn). Does the program need such a big ram, or I have done something wrong. The attached code is working correctly for lattice size N=60, and the changes for bigger system size related to the variables “bnd1” and “bnd1”. Any help would be appreciated.

Regards
Javad

P.S: the code attached below:

using namespace std;
using namespace itensor;

int
main()
    {

    auto N = 90;

    //
    // Initialize the site degrees of freedom.
    //
    auto sites = SpinHalf(N,{"ConserveQNs=",true});

    // Set the initial wavefunction matrix product state
    // to be a Neel state.
    //
    // This choice implicitly sets the global Sz quantum number
    // of the wavefunction to zero. Since it is an MPS
    // it will remain in this quantum number sector.
    //
    auto state = InitState(sites);
    for(int i = 1; i <= N; ++i)
    {
   		if(i%2 == 1)
    		     state.set(i,"Up");
	     else
     		     state.set(i,"Dn");
	}
    auto psi0 = MPS(state);




    auto  bnd1 = {1, 2, 3, 3, 4, 5, 5, 6, 7, 7, 8, 9, 9, 11,
      12, 13, 13, 14, 15, 15, 16, 17, 17, 18, 19, 19, 21, 22,
      23, 23, 24, 25, 25, 26, 27, 27, 28, 29, 29, 31, 32, 33,
      33, 34, 35, 35, 36, 37, 37, 38, 39, 39, 41, 42, 43, 43,
      44, 45, 45, 46, 47, 47, 48, 49, 49, 51, 52, 53, 53, 54,
      55, 55, 56, 57, 57, 58, 59, 59, 61, 62, 63, 63, 64, 65,
      65, 66, 67, 67, 68, 69, 69, 71, 72, 73, 73, 74, 75, 75,
      76, 77, 77, 78, 79, 79, 81, 82, 83, 84, 85, 86, 87, 88,
      89, 10, 20, 20, 30, 30, 40, 40, 50, 50, 60, 60, 70, 70,
      80, 80, 90, 90,
      1, 1, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7,
      8, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 13, 14, 14,
      14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19,
      20, 20, 21, 21, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25,
      26, 26, 26, 27, 27, 27, 28, 28, 28, 29, 29, 30, 30, 31, 31,
      32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 37,
      37, 37, 38, 38, 38, 39, 39, 40, 40, 41, 41, 42, 42, 43, 43,
      43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48,
      48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 53, 54, 54, 54,
      55, 55, 55, 56, 56, 56, 57, 57, 57, 58, 58, 58, 59, 59, 60,
      60, 61, 61, 62, 62, 63, 63, 63, 64, 64, 64, 65, 65, 65, 66,
      66, 66, 67, 67, 67, 68, 68, 68, 69, 69, 70, 70, 71, 71, 72,
      72, 73, 73, 73, 74, 74, 74, 75, 75, 75, 76, 76, 76, 77, 77,
      77, 78, 78, 78, 79, 79, 80, 80, 81, 82, 83, 84, 85, 86, 87,
      88, 10, 20, 20, 30, 30, 40, 40, 50, 50, 60, 60, 70, 70, 80,
      80, 90, 90, 1, 1, 11, 11, 21, 21, 31, 31, 41, 41, 51, 51, 61,
      61, 71, 71, 81};

    auto  bnd2 = {2, 3, 4, 12, 5, 6, 14, 7, 8, 16, 9, 10, 18, 12, 13,
      14, 22, 15, 16, 24, 17, 18, 26, 19, 20, 28, 22, 23, 24, 32, 25,
      26, 34, 27, 28, 36, 29, 30, 38, 32, 33, 34, 42, 35, 36, 44, 37,
      38, 46, 39, 40, 48, 42, 43, 44, 52, 45, 46, 54, 47, 48, 56, 49,
      50, 58, 52, 53, 54, 62, 55, 56, 64, 57, 58, 66, 59, 60, 68, 62,
      63, 64, 72, 65, 66, 74, 67, 68, 76, 69, 70, 78, 72, 73, 74, 82,
      75, 76, 84, 77, 78, 86, 79, 80, 88, 82, 83, 84, 85, 86, 87, 88,
      89, 90, 1, 1, 11, 11, 21, 21, 31, 31, 41, 41, 51, 51, 61, 61, 71,
      71, 81, 3, 11, 4, 12, 5, 11, 13, 6, 12, 14, 7, 13, 15, 8, 14, 16,
      9, 15, 17, 10, 16, 18, 17, 19, 18, 20, 13, 21, 14, 22, 15, 21, 23,
      16, 22, 24, 17, 23, 25, 18, 24, 26, 19, 25, 27, 20, 26, 28, 27, 29,
      28, 30, 23, 31, 24, 32, 25, 31, 33, 26, 32, 34, 27, 33, 35, 28, 34,
      36, 29, 35, 37, 30, 36, 38, 37, 39, 38, 40, 33, 41, 34, 42, 35, 41,
      43, 36, 42, 44, 37, 43, 45, 38, 44, 46, 39, 45, 47, 40, 46, 48, 47,
      49, 48, 50, 43, 51, 44, 52, 45, 51, 53, 46, 52, 54, 47, 53, 55, 48,
      54, 56, 49, 55, 57, 50, 56, 58, 57, 59, 58, 60, 53, 61, 54, 62, 55,
      61, 63, 56, 62, 64, 57, 63, 65, 58, 64, 66, 59, 65, 67, 60, 66, 68,
      67, 69, 68, 70, 63, 71, 64, 72, 65, 71, 73, 66, 72, 74, 67, 73, 75,
      68, 74, 76, 69, 75, 77, 70, 76, 78, 77, 79, 78, 80, 73, 81, 74, 82,
      75, 81, 83, 76, 82, 84, 77, 83, 85, 78, 84, 86, 79, 85, 87, 80, 86,
      88, 87, 89, 88, 90, 83, 84, 85, 86, 87, 88, 89, 90, 2, 2, 12, 12, 22,
      22, 32, 32, 42, 42, 52, 52, 62, 62, 72, 72, 82, 9, 19, 19, 29, 29, 39,
      49, 49, 59, 59, 69, 69, 79, 79, 89, 89};

    auto Jb=0.0;
    int  n1=130;
    std::ofstream enf("dmrg_45_46_47_svn_new.txt");
    enf  << "Jb,     E0,       Sz45,        Sz46,        Sz47,         Sz45Sz46,        Sz45Sz47,     Sp45Sn46,      Sp45Sn47,    Svn,    SvN_norm\n"  << std::flush;


    for (int jj=0; jj<=20; ++jj)

    {
            auto ampo = AutoMPO(sites);
            int cte=0;
            for(auto bnd : zip(bnd1,bnd2)){
                auto [s1,s2] = bnd;
                //printfln(" s1 : %i , s2 : %i ",bnd.s1,bnd.s2);
                if (cte < n1)
                {
                  //printfln(" cte : %i,  s1 : %i , s2 : %i ",cte, s1,s2);
                  ampo += 0.5,"S+",s1,"S-",s2;
                  ampo += 0.5,"S-",s1,"S+",s2;
                }
                else
                {
                  Jb=jj*0.02;
                  ampo += 0.5*Jb,"S+",s1,"S-",s2;
                  ampo += 0.5*Jb,"S-",s1,"S+",s2;
                }
                cte +=1;
            }

            auto H = toMPO(ampo);

            //
            // overlap calculates matrix elements of MPO's with respect to MPS's
            // inner(psi0,H,psi0) = <psi0|H|psi0>
            //
            printfln("Initial energy = %.5f", inner(psi0,H,psi0) );


            //
            // Set the parameters controlling the accuracy of the DMRG
            // calculation for each DMRG sweep.
            // Here less than 5 cutoff values are provided, for example,
            // so all remaining sweeps will use the last one given (= 1E-10).
            //
            auto sweeps = Sweeps(7);
            sweeps.maxdim() = 10,20,200,400,600,800,1000;
            sweeps.cutoff() = 1E-12;
            sweeps.niter() = 2;
            sweeps.noise() = 1E-7,1E-8,0.0;
            //println(sweeps);

            //
            // Begin the DMRG calculation
            // for the ground state
            //
            auto [en,psi] = dmrg(H,psi0,sweeps,{"Quiet=",true});
            auto psi0 = psi;

            //
            // Measuring Sz
            //
            //println("==============================");
            //println("======= Measuring Sz ==========");

            //auto Sz     = Vector(3);
            float Sz  [3]  = {};
            int cte0 = 0;
            for( auto m : {45,46,47} )
                {
                //re-gauge psi to get ready to measure at position j
                psi.position(m);
                auto ket = psi(m);
                auto bra = dag(prime(ket,"Site"));
                auto szop = op(sites,"Sz",m);
                //take an inner product
                Sz[cte0] = elt(bra*szop*ket);
                //printfln("%d %.12f %.12f",j,szj,sxj);
                cte0+=1;
                }



            println("==============================");
            println("=======    SziSzj   ==========");

            int i    = 45;
            int cte1 = 0;
            //auto SziSzj  = Vector(2);
            float SziSzj [2]  = {};
            for(auto j : {46,47})
               {
                auto Szi = op(sites,"Sz",i);
                auto Szj = op(sites,"Sz",j);
                psi.position(i);
                auto C = psi(i);
                C *=Szi;
                auto ir = commonIndex(psi(i),psi(i+1),"Link");
                C *= dag(prime(prime(psi(i),"Site"),ir));
                for (int k=i+1; k<j; ++k)
                {
                    C *= psi(k);
                    C *= dag(prime(psi(k),"Link"));
                }
                C *= psi(j);
                C *= Szj;
                auto il = commonIndex(psi(j),psi(j-1),"Link");
                C *=  dag(prime(prime(psi(j),"Site"),il));
                SziSzj[cte1]=elt(C);
                //printfln("\n Jb : %.4f,  SzSz  : %.10f",Jb, elt(C));
                cte1+=1;
            }


            println("==============================");
            println("=======    SpiSnj   ==========");
            //auto SpiSnj=Vector(2);
            float SpiSnj [2]  = {};
            int n    = 45;
            int cte2 = 0;
            for(auto j : {46,47})
               {
                auto Spi = op(sites,"S+",n);
                auto Snj = op(sites,"S-",j);
                psi.position(n);
                auto C = psi(n);
                C *=Spi;
                auto ir = commonIndex(psi(n),psi(n+1),"Link");
                C *= dag(prime(prime(psi(n),"Site"),ir));
                for (int k=n+1; k<j; ++k)
                {
                    C *= psi(k);
                    C *= dag(prime(psi(k),"Link"));
                }
                C *= psi(j);
                C *=  Snj;
                auto il = commonIndex(psi(j),psi(j-1),"Link");
                C *=  dag(prime(prime(psi(j),"Site"),il));
                SpiSnj[cte2]=elt(C);
                //printfln("\n Jb : %.4f,  S+S-  : %.10f",Jb, elt(C));
                cte2+=1;
            }

            //compute EE
	    auto b = 45;
            psi.position(b);
            ITensor wf = psi.A(b)*psi.A(b+1);
            auto U = psi.A(b);
            ITensor S,V;
            auto spectrum = svd(wf,U,S,V);
            Real SvN = 0.;
            for(auto p : spectrum.eigs())
	        {
		    if(p > 1E-12) SvN += -p*log(p);
		 }
            auto normSvN = SvN /log(2);

	    printfln("\n Jb : %.4f,  E0  : %.10f,  Sz1  : %.10f,  Sz2  : %.10f ,  Sz3  : %.10f,  Sz1Sz2  : %.10f,  Sz1Sz3  : %.10f, Sp1Sn2  : %.10f,  Sp1Sn3  : %.10f, SvN  : %.10f,  normSvN  : %.10f",
			                             Jb,  en,   Sz[0],   Sz[1],   Sz[2],  SziSzj[0],   SziSzj[1],   SpiSnj[0],   SpiSnj[1], SvN, normSvN);

	    enf  << format("%.4f   %.10f   %.10f    %.10f   %.10f   %.10f   %.10f   %.10f   %.10f   %.10f   %.10f\n",
			                                 Jb,  en,    Sz[0],  Sz[1],   Sz[2],   SziSzj[0],   SziSzj[1],   SpiSnj[0],   SpiSnj[1], SvN, normSvN)<< std::flush;

          }

        enf.close();

    return 0;
    }

and her is the lattice figure

The lattice figure it helpful to answer your question. When you say that you go to N=90 sites, are you growing the system in both the x and y directions? Or just making it longer in the x direction say?

The reason I ask this is that let’s say there are N=N_x N_y sites on a quasi-2D lattice like yours. Then DMRG scales only linearly in N_x but exponentially in N_y. So the relevant thing to understand for how costly in terms of time and memory a calculation will be is not the total size N but the y-size N_y.

Usually the way people have to work for 2D DMRG is to fix N_y to a small size like 4,6,8 then study various sizes in N_x. Basically it means studying laders or cylinders and not square shaped patches of 2D systems.

A good reference about best practices for 2D DMRG is this one:
Studying Two Dimensional Systems With the Density Matrix Renormalization Group

1 Like

Dear Miles

Many thanks for the hint and also bringing my attention to your useful paper.

Javad

1 Like