API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Sep 2nd, 2020
70
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 22.94 KB | None | 0 0
raw download clone embed print report
  1. #include "itensor/all.h"
  2. #include <iostream>
  3. #include <chrono>
  4. #include <cstdlib>
  5. #include <stdio.h>
  6. #include <ctime>
  7. #include <gsl/gsl_blas.h>
  8. #include <gsl/gsl_linalg.h>
  9. #include <gsl/gsl_eigen.h>
  10. #include <gsl/gsl_multimin.h>
  11. #include <gsl/gsl_complex_math.h>
  12. #include <string>
  13. #include <set>
  14. #include <vector>
  15. #include <functional>
  16. #include <algorithm>
  17.  
  18. using namespace itensor;
  19. using namespace std;
  20. using namespace std::chrono;
  21.  
  22. #define Ndimers 1
  23. #define Nspins 4
  24. #define UU 1.0
  25. #define optimization true
  26. #define ms 4
  27. #define groups 2
  28.  
  29. typedef std::chrono::high_resolution_clock Clock;
  30.  
  31. typedef struct params_t {
  32.   int N;
  33.   double U;
  34.   double h;
  35.   int index;
  36.   int point;
  37. } params_t;
  38.  
  39. vector<vector<int>> dimers;
  40. auto states = vector<MPS>(Ndimers);
  41. auto sites_dmrg = Hubbard(Nspins);
  42. auto sites = Hubbard(Nspins);//, {"ConserveSz", false});
  43. auto H = vector<MPO>(16);
  44.  
  45. int find_index(vector <int> numbers, int num)
  46. {
  47.   vector<int>:: iterator it;
  48.   int pos=0;
  49.   for (it=numbers.begin(); it!=numbers.end(); it++) {
  50.       if (*it == num)
  51.         return pos;
  52.       ++pos;
  53.     }
  54.     return -1;
  55. }
  56.  
  57. double create_state(const gsl_vector *cicoeff, void *params)
  58. {
  59.     int N = ((params_t *) params)->N;
  60.     int index = ((params_t *) params)->index;
  61.  
  62.     auto psi = MPS(sites);
  63.     double normalization;
  64.  
  65.     for(int n = 1; n <= N; n += 2)
  66.         {
  67.  
  68.         auto s1 = sites(n);
  69.         auto s2 = sites(n+1);
  70.  
  71.         normalization = 0;
  72.         int pos = n/2*ms;
  73.         for (int i=pos; i<=pos+ms-1; i++)
  74.           normalization += pow(gsl_vector_get(cicoeff,i),2);
  75.  
  76.         auto wf = ITensor(s1,s2);
  77.         wf.set(s1(2),s2(3), gsl_vector_get(cicoeff, pos)/sqrt(normalization));
  78.         wf.set(s1(3),s2(2), gsl_vector_get(cicoeff, pos+1)/sqrt(normalization));
  79.         wf.set(s1(1),s2(4), gsl_vector_get(cicoeff, pos+2)/sqrt(normalization));
  80.         wf.set(s1(4),s2(1), gsl_vector_get(cicoeff, pos+3)/sqrt(normalization));
  81.         //wf.set(s1(2),s2(2), gsl_vector_get(cicoeff, pos+4)/sqrt(normalization));
  82.         //wf.set(s1(3),s2(3), gsl_vector_get(cicoeff, pos+5)/sqrt(normalization));
  83.  
  84.  
  85.        /* wf.set(s1(4),s2(1), gsl_vector_get(cicoeff, pos)/sqrt(normalization));
  86.         wf.set(s1(2),s2(3), gsl_vector_get(cicoeff, pos+1)/sqrt(normalization));
  87.         wf.set(s1(3),s2(2), gsl_vector_get(cicoeff, pos+2)/sqrt(normalization));
  88.         wf.set(s1(1),s2(4), gsl_vector_get(cicoeff, pos+3)/sqrt(normalization));
  89.         wf.set(s1(2),s2(2), gsl_vector_get(cicoeff, pos+4)/sqrt(normalization));
  90.         wf.set(s1(3),s2(3), gsl_vector_get(cicoeff, pos+5)/sqrt(normalization));
  91.  
  92.         wf.set(s1(2),s2(1), gsl_vector_get(cicoeff, pos+6)/sqrt(normalization));
  93.         wf.set(s1(3),s2(1), gsl_vector_get(cicoeff, pos+7)/sqrt(normalization));
  94.         wf.set(s1(1),s2(2), gsl_vector_get(cicoeff, pos+8)/sqrt(normalization));
  95.         wf.set(s1(1),s2(3), gsl_vector_get(cicoeff, pos+9)/sqrt(normalization));
  96.  
  97.     wf.set(s1(2),s2(4), gsl_vector_get(cicoeff, pos+10)/sqrt(normalization));
  98.         wf.set(s1(3),s2(4), gsl_vector_get(cicoeff, pos+11)/sqrt(normalization));
  99.         wf.set(s1(4),s2(2), gsl_vector_get(cicoeff, pos+12)/sqrt(normalization));
  100.         wf.set(s1(4),s2(3), gsl_vector_get(cicoeff, pos+13)/sqrt(normalization));
  101.        
  102.     wf.set(s1(4),s2(4), gsl_vector_get(cicoeff, pos+14)/sqrt(normalization));
  103.         wf.set(s1(1),s2(1), gsl_vector_get(cicoeff, pos+15)/sqrt(normalization));
  104. */
  105.         psi.ref(n) = ITensor(s1);
  106.         psi.ref(n+1) = ITensor(s2);
  107.         auto [U,S,V] = svd(wf,{inds(psi.ref(n))},{"Cutoff=",1E-8});
  108.         psi.ref(n) = U;
  109.         psi.ref(n+1) = S*V;
  110.  
  111.       }
  112.       int swaps = 1;
  113.       vector<pair <int,int>> swaps_;
  114.       vector <int> dimer;
  115.       for (int i=0; i<N; i++)
  116.         dimer.push_back(dimers[index][i]);
  117.  
  118.       while (swaps > 0)
  119.       {
  120.         swaps = 0;
  121.         pair <int,int> pair1;
  122.         for (int i=1; i<N; i++)
  123.         {
  124.           int pos1, pos2;
  125.           pos1 = find_index(dimer, i);
  126.           pos2 = find_index(dimer, i+1);
  127.           if (pos1 > pos2)
  128.           {
  129.             pair1.first = dimer[pos1];
  130.             pair1.second = dimer[pos2];
  131.             swaps_.push_back(pair1);
  132.             swaps += 1;
  133.             int temp = dimer[pos1];
  134.             dimer[pos1] = dimer[pos2];
  135.             dimer[pos2] = temp;
  136.             continue;
  137.           }
  138.         }
  139.       }
  140.  
  141.       auto s1 = vector<Index>(1);
  142.       auto s2 = vector<Index>(1);
  143.       int ar1, ar2;
  144.       ITensor wf, S,V,U;
  145.       for (int i=swaps_.size()-1; i>=0; i--){
  146.         ar1 = swaps_[i].first;
  147.         ar2 = swaps_[i].second;
  148.         s1.clear(); s1.push_back(siteIndex(psi,ar1));
  149.         s2.clear(); s2.push_back(siteIndex(psi,ar2));
  150.         auto is1 = IndexSet(s1);
  151.         auto is2 = IndexSet(s2);
  152.  
  153.         psi.position(ar1);
  154.         wf = psi.A(ar1)*psi.A(ar2);
  155.         wf.swapInds(is1,is2);
  156.         wf.noPrime();
  157.         U = psi.A(ar1);
  158.         svd(wf,U,S,V,{"Cutoff=",1E-8});
  159.         psi.setA(ar1,U);
  160.         psi.setA(ar2,S*V);
  161.       }
  162.       states.at(index) = psi;
  163.       int point = ((params_t *) params)->point;
  164.       return inner(states.at(index),H[point],states.at(index));
  165.  
  166. }
  167.  
  168. void find_grad(const gsl_vector *cicoeff, void *params, gsl_vector *grad)
  169. {
  170.  
  171.   /*  for (int i=0; i<N; i++) {
  172.     etacoeff[i] += 2*h;
  173.     double energy = 0;
  174.         energy -= eval_energy(N,M,etacoeff,delta);
  175.     etacoeff[i] -= h;
  176.     energy += 8*eval_energy(N,M,etacoeff,delta);
  177.     etacoeff[i] -= 2*h;
  178.     energy -= 8*eval_energy(N,M,etacoeff,delta);
  179.     etacoeff[i] -= h;
  180.     energy += eval_energy(N,M,etacoeff,delta);
  181.     grad[i] = energy/(12*h);
  182.     etacoeff[i] += 2*h;
  183. }*/
  184.    double h = ((params_t *) params)->h;
  185.    int N = ((params_t *) params)->N;
  186.  
  187.    gsl_vector *test_cicoeff;
  188.    test_cicoeff = gsl_vector_alloc(ms*groups);
  189.    for (int i=0; i<ms*groups; i++)
  190.       gsl_vector_set(test_cicoeff, i, gsl_vector_get(cicoeff, i));
  191.  
  192.    for (int i=0; i<ms*groups; i++) {
  193.         gsl_vector_set(test_cicoeff,i, gsl_vector_get(test_cicoeff,i)+ h);
  194.         double energy = 0;
  195.       energy += create_state(test_cicoeff, params);
  196.       gsl_vector_set(test_cicoeff,i, gsl_vector_get(test_cicoeff,i)- h);
  197.       energy -= create_state(test_cicoeff, params);
  198.       gsl_vector_set(grad,i, energy/h);
  199.     }
  200. }
  201.  
  202. void fdf(const gsl_vector *cicoeff, void *params, double *energy, gsl_vector *grad)
  203. {
  204.     *energy = create_state(cicoeff, params);
  205.     find_grad(cicoeff, params, grad);
  206. }
  207.  
  208. void debug_print(const gsl_multimin_fdfminimizer *s, const int k, int N)
  209. {
  210.     const gsl_vector   *x          = NULL;
  211.     const gsl_vector   *g          = NULL;
  212.  
  213.     /**
  214.      * Display current iteration.
  215.      */
  216.     //cout << "\nAt iteration k = " << k << endl;
  217.  
  218.    
  219.     /*g = gsl_multimin_fdfminimizer_gradient(s);
  220.     cout << "g = [ ";
  221.     for (int i=0; i<N; i++)
  222.       cout << gsl_vector_get(g, i) << " ";
  223.     cout << endl;
  224.  
  225.     x = gsl_multimin_fdfminimizer_x(s);
  226.     cout << "x = [ ";
  227.     for (int i=0; i<N; i++)
  228.       cout << gsl_vector_get(x, i) << " ";
  229.     cout << endl;*/
  230.  
  231.     cout << "f(x) = " << gsl_multimin_fdfminimizer_minimum(s) << endl;
  232. } /** end function print_result() */
  233.  
  234. string convertToString(char* a)
  235. {
  236.     int i, size=strlen(a);
  237.     string s = "";
  238.     for (i = 0; i < size; i++) {
  239.         s = s + a[i];
  240.     }
  241.     return s;
  242. }
  243.  
  244. void generate_H()
  245. {
  246.  
  247.   for (int ii=0; ii<=15; ii++)
  248.   {
  249.     double U = UU + ii*0.5;
  250.     auto ampo = AutoMPO(sites);
  251.  
  252.         for(int b = 1; b < Nspins; ++b)
  253.             {
  254.                 ampo += -1.0,"Cdagup",b,"Cup",b+1; //Hopping for spin up
  255.                 ampo += -1.0,"Cdagup",b+1,"Cup",b; // <-
  256.                 ampo += -1.0,"Cdagdn",b,"Cdn",b+1; //Hopping for spin down
  257.                 ampo += -1.0,"Cdagdn",b+1,"Cdn",b;// <-
  258.                 ampo += U, "Nup", b, "Ndn", b;
  259.             }
  260.                 ampo += U, "Nup", Nspins, "Ndn", Nspins;
  261.                 ampo += -1.0,"Cdagup",1,"Cup",Nspins; //Hopping for spin up
  262.                 ampo += -1.0,"Cdagup",Nspins,"Cup",1; // <-
  263.                 ampo += -1.0,"Cdagdn",1,"Cdn",Nspins; //Hopping for spin down
  264.                 ampo += -1.0,"Cdagdn",Nspins,"Cdn",1;// <-
  265.     MPO H_ = toMPO(ampo);
  266.     H.at(ii) = H_;
  267.   }
  268. }
  269.  
  270. void orth_states()
  271. {
  272.     auto states_temp = vector<MPS>(Ndimers);
  273.     states_temp.at(0) = states[0];
  274.     for (int i=1; i<Ndimers; i++)
  275.     {
  276.       states_temp.at(i) = states.at(i);
  277.       states_temp.at(i).plusEq(-inner(states.at(0),states.at(i))*states.at(0));
  278.     }
  279.  
  280.     for (int i=1; i<Ndimers; i++)
  281.     {
  282.       states[i] = states_temp[i];
  283.       states[i] /= norm(states_temp[i]); //Check if necessary
  284.     }
  285. }
  286.  
  287. double fRand(double fMin, double fMax)
  288. {
  289.     double f = (double)rand() / RAND_MAX;
  290.     return fMin + f * (fMax - fMin);
  291. }
  292.  
  293. void computeG(gsl_matrix *G, int point, int ind)
  294. {
  295.     for (int i=1; i<=Nspins; i++)
  296.     {
  297.         for (int j=1; j<=Nspins; j++)
  298.         {
  299.             for (int sgi=0; sgi<=1; sgi+=1) {
  300.                 for (int sgj=0; sgj<=1; sgj+=1) {
  301.                      char Cd[7] = "Cdag"; char C[4] = "C";
  302.                         char Sgi[3], Sgj[3]; char up[3]="up"; char dn[3]="dn";
  303.         if (sgi == 0) {
  304.             strcpy(Sgi,"dn");
  305.         }
  306.         if (sgj == 0) {
  307.             strcpy(Sgj,"dn");
  308.         }
  309.         if (sgi == 1) {
  310.             strcpy(Sgi,"up");
  311.         }
  312.         if (sgj == 1) {
  313.             strcpy(Sgj, "up");
  314.         }
  315.                     MPS Phi_0 = removeQNs(states.at(ind));
  316.                     MPO h = H[point];
  317.                             auto ampo = AutoMPO(sites);
  318.                     if (sgi == sgj) {
  319.                     ampo += convertToString(Cd)+convertToString(Sgi), i, convertToString(C)+convertToString(Sgj), j;
  320.                     MPO cc = toMPO(ampo);
  321.                                         MPS h_mps = applyMPO(h,Phi_0); MPS cc_mps = applyMPO(cc,Phi_0);
  322.                      cout << " " << inner(Phi_0,cc,h_mps) << " " << inner(Phi_0, h, cc_mps) << " " << i << j << endl;
  323.                     double elem = -(inner(Phi_0, h, cc_mps) - inner(Phi_0,cc,h_mps));
  324.                             gsl_matrix_set(G, 2*(i-1)+sgi,2*(j-1)+sgj,elem);}
  325.                                 }
  326.                         }
  327.         }
  328.     }
  329. }
  330.  
  331. int main()
  332. {
  333.     srand (time(NULL));
  334.     //srand(2);
  335.     int N = Nspins;
  336.     params_t   *params = NULL;
  337.     params = (params_t *) malloc(sizeof(params_t));
  338.     generate_H();
  339.     params->N = N;
  340.     params->U = UU;
  341.     params->h = 1e-4;
  342.  
  343.     gsl_vector *cicoeff;
  344.  
  345.     cicoeff = gsl_vector_alloc(ms*groups);
  346.  
  347.     FILE *out = fopen("results.txt", "w");
  348.  
  349.     vector<int> dimer;
  350.  
  351.     FILE *fp = fopen("input.dat", "r");
  352.     for (int i=0; i<Ndimers; i++)
  353.     {
  354.       int temp;
  355.       for (int j=0; j<N; j++)
  356.       {
  357.         fscanf(fp, "%d ", &temp);
  358.         dimer.push_back(temp);
  359.       }
  360.       dimers.push_back(dimer);
  361.       dimer.clear();
  362.     }
  363.     fclose(fp);
  364.  
  365.     params->point = 0;
  366.     while (params->point <= 15)
  367.     {
  368.  
  369.       for (int i=0; i<Ndimers; i++)
  370.       {
  371.  
  372.  
  373.            for (int i=0; i<ms*groups; i++)
  374.              gsl_vector_set(cicoeff,i,fRand(0,1));
  375.  
  376.  
  377.  
  378. /*    for (int i=ms; i<ms*groups; i+=ms)
  379.     {
  380.       gsl_vector_set(cicoeff, i, gsl_vector_get(cicoeff, 0));
  381.       gsl_vector_set(cicoeff, i+1, gsl_vector_get(cicoeff, 1));
  382.       gsl_vector_set(cicoeff, i+2, gsl_vector_get(cicoeff, 2));
  383.       gsl_vector_set(cicoeff, i+3, gsl_vector_get(cicoeff, 3));
  384.       gsl_vector_set(cicoeff, i+4, gsl_vector_get(cicoeff, 4));
  385.       gsl_vector_set(cicoeff, i+5, gsl_vector_get(cicoeff, 5));
  386.       gsl_vector_set(cicoeff, i+6, gsl_vector_get(cicoeff, 6));
  387.       gsl_vector_set(cicoeff, i+7, gsl_vector_get(cicoeff, 7));
  388.       gsl_vector_set(cicoeff, i+8, gsl_vector_get(cicoeff, 8));
  389.       gsl_vector_set(cicoeff, i+9, gsl_vector_get(cicoeff, 9));
  390.       gsl_vector_set(cicoeff, i+10, gsl_vector_get(cicoeff, 10));
  391.       gsl_vector_set(cicoeff, i+11, gsl_vector_get(cicoeff, 11));
  392.       gsl_vector_set(cicoeff, i+12, gsl_vector_get(cicoeff, 12));
  393.       gsl_vector_set(cicoeff, i+13, gsl_vector_get(cicoeff, 13));
  394.       gsl_vector_set(cicoeff, i+14, gsl_vector_get(cicoeff, 14));
  395.       gsl_vector_set(cicoeff, i+15, gsl_vector_get(cicoeff, 15));
  396. }*/
  397.  
  398.             params->index = i;
  399.  
  400.             if (optimization)
  401.             {
  402.  
  403.               int status;
  404.  
  405.               const gsl_multimin_fdfminimizer_type *T;
  406.               gsl_multimin_fdfminimizer *s;
  407.  
  408.               gsl_multimin_function_fdf my_func;
  409.  
  410.               my_func.n = ms*groups;
  411.               my_func.f = &create_state;
  412.               my_func.df = &find_grad;
  413.               my_func.fdf = &fdf;
  414.               my_func.params = params;
  415.  
  416.             //  T = gsl_multimin_fdfminimizer_vector_bfgs;
  417.             // T = gsl_multimin_fdfminimizer_conjugate_pr;
  418.                 T = gsl_multimin_fdfminimizer_conjugate_fr;
  419.            //   T = gsl_multimin_fdfminimizer_steepest_descent;
  420.  
  421.               s = gsl_multimin_fdfminimizer_alloc(T, ms*groups);
  422.  
  423.               gsl_multimin_fdfminimizer_set(s,
  424.                                             &my_func,
  425.                                             cicoeff,
  426.                                             0.1,      /** step size */
  427.                                             1e-4); /** tol */
  428.  
  429.               int  k = 0;
  430.               double old_val = 0;
  431.               double new_val;
  432.               do {
  433.  
  434.                   status = gsl_multimin_fdfminimizer_iterate(s);
  435.                   /*new_val = gsl_multimin_fdfminimizer_minimum(s);
  436.                   if (abs(old_val-new_val) > 1e-6)
  437.                     status = -2;
  438.                   old_val = new_val;*/
  439.                   status = gsl_multimin_test_gradient(s->gradient, 1e-4);
  440.                   //debug_print(s, k, ms*groups);
  441.                   ++k;
  442.               } while(status == GSL_CONTINUE && k <= 200);
  443.               if (k >= 200)
  444.                 cout << "Didn't converge properly.\n";
  445.               cout << "New dimer: " << " " << gsl_multimin_fdfminimizer_minimum(s) << "\n";
  446.               const gsl_vector   *x          = NULL;
  447.               x = gsl_multimin_fdfminimizer_x(s);
  448.               for (int i=0; i<ms*groups; i++) {
  449.                 gsl_vector_set(cicoeff,i,gsl_vector_get(x,i));
  450.           }
  451.    
  452.           }
  453.           else
  454.           {
  455.             create_state(cicoeff, params);
  456.           }
  457.  
  458.     for (int gg=0; gg<groups; gg++) {
  459.     double normalization=0;
  460.         for (int i=0; i<ms; i++) {
  461.             normalization += pow(gsl_vector_get(cicoeff,gg*ms+i),2);
  462.         }
  463.  
  464.         for (int i=0; i<ms; i++) {
  465.                 //cout << "cicoeff[" << gg*ms+i+1 << "] = " << gsl_vector_get(cicoeff, gg*ms+i)/sqrt(normalization) << "\t" << endl;
  466.         }
  467.            
  468.     /*SpinHalf sites4 = SpinHalf(4, {"ConserveQNs=",false});
  469.         auto psi4 = vector<MPS>(ms);
  470.  
  471.         auto init4 = InitState(sites4);
  472.  
  473.     init4.set(1, "Up"); init4.set(2, "Up"); init4.set(3, "Dn"); init4.set(4, "Dn");
  474.         psi4.at(0) = MPS(init4); init4 = InitState(sites4);
  475.         init4.set(1, "Up"); init4.set(2, "Dn"); init4.set(3, "Up"); init4.set(4, "Dn");
  476.         psi4.at(1) = MPS(init4); init4 = InitState(sites4);
  477.         init4.set(1, "Up"); init4.set(2, "Dn"); init4.set(3, "Dn"); init4.set(4, "Up");
  478.         psi4.at(2) = MPS(init4); init4 = InitState(sites4);
  479.         init4.set(1, "Dn"); init4.set(2, "Up"); init4.set(3, "Up"); init4.set(4, "Dn");
  480.         psi4.at(3) = MPS(init4); init4 = InitState(sites4);
  481.         init4.set(1, "Dn"); init4.set(2, "Up"); init4.set(3, "Dn"); init4.set(4, "Up");
  482.         psi4.at(4) = MPS(init4); init4 = InitState(sites4);
  483.         init4.set(1, "Dn"); init4.set(2, "Dn"); init4.set(3, "Up"); init4.set(4, "Up");
  484.         psi4.at(5) = MPS(init4); init4 = InitState(sites4);
  485.         init4.set(1, "Dn"); init4.set(2, "Up"); init4.set(3, "Up"); init4.set(4, "Up");
  486.         psi4.at(6) = MPS(init4);init4 = InitState(sites4);
  487.         init4.set(1, "Up"); init4.set(2, "Dn"); init4.set(3, "Up"); init4.set(4, "Up");
  488.         psi4.at(7) = MPS(init4);init4 = InitState(sites4);
  489.         init4.set(1, "Up"); init4.set(2, "Up"); init4.set(3, "Dn"); init4.set(4, "Up");
  490.         psi4.at(8) = MPS(init4); init4 = InitState(sites4);
  491.         init4.set(1, "Up"); init4.set(2, "Up"); init4.set(3, "Up"); init4.set(4, "Dn");
  492.         psi4.at(9) = MPS(init4); init4 = InitState(sites4);
  493.         init4.set(1, "Dn"); init4.set(2, "Dn"); init4.set(3, "Dn"); init4.set(4, "Up");
  494.         psi4.at(10) = MPS(init4);init4 = InitState(sites4);
  495.         init4.set(1, "Dn"); init4.set(2, "Dn"); init4.set(3, "Up"); init4.set(4, "Dn");
  496.         psi4.at(11) = MPS(init4);init4 = InitState(sites4);
  497.         init4.set(1, "Dn"); init4.set(2, "Up"); init4.set(3, "Dn"); init4.set(4, "Dn");
  498.         psi4.at(12) = MPS(init4);init4 = InitState(sites4);
  499.         init4.set(1, "Up"); init4.set(2, "Dn"); init4.set(3, "Dn"); init4.set(4, "Dn");
  500.         psi4.at(13) = MPS(init4); init4 = InitState(sites4);
  501.         init4.set(1, "Dn"); init4.set(2, "Dn"); init4.set(3, "Dn"); init4.set(4, "Dn");
  502.         psi4.at(14) = MPS(init4); init4 = InitState(sites4);
  503.         init4.set(1, "Up"); init4.set(2, "Up"); init4.set(3, "Up"); init4.set(4, "Up");
  504.         psi4.at(15) = MPS(init4); init4 = InitState(sites4);
  505.  
  506.  
  507.         auto sz4 = AutoMPO(sites4);
  508.         for(int i=1; i<=4; i++)
  509.         {
  510.                 sz4 += "Sz", i;
  511.         }
  512.          MPO Sz4 = toMPO(sz4);
  513.  
  514.         for (int i=0; i<ms; i++)
  515.       psi4.at(i) *= gsl_vector_get(cicoeff,gg*ms+i)/sqrt(normalization);
  516.    
  517.         for (int i=1; i<ms; i++)
  518.         {
  519.              auto temp = psi4.at(i);
  520.              psi4.at(0).plusEq(temp);
  521.         }
  522.         cout << "Cluster " << gg+1 << ", Sz = " << inner(psi4.at(0),Sz4,psi4.at(0))/inner(psi4.at(0),psi4.at(0)) << endl;
  523.  
  524.       */}
  525.     }  
  526.        //orth_states();
  527.         gsl_matrix *metric;
  528.         gsl_matrix *Hmat;
  529.         Hmat = gsl_matrix_alloc(Ndimers,Ndimers);
  530.         metric = gsl_matrix_alloc(Ndimers,Ndimers);
  531.  
  532.         for (int i=0; i<Ndimers; i++)
  533.           gsl_matrix_set(metric, i, i, 1.0);
  534.  
  535.         for (int i=0; i<Ndimers; i++){
  536.           for (int j=i+1; j<Ndimers; j++){
  537.             double val = inner(states.at(i), states.at(j));
  538.             gsl_matrix_set(metric, i, j, val);
  539.             gsl_matrix_set(metric, j, i, val);}}
  540.  
  541.         gsl_vector *evals;
  542.         gsl_matrix *evec;
  543.         evals = gsl_vector_alloc(Ndimers);
  544.         evec = gsl_matrix_alloc(Ndimers,Ndimers);
  545.         gsl_eigen_gensymmv_workspace *w;
  546.         w = gsl_eigen_gensymmv_alloc (Ndimers);
  547.  
  548.         for (int i=0; i<Ndimers; i++){
  549.       for (int j=i; j<Ndimers; j++){
  550.         double val =  inner(states.at(i), H[params->point], states.at(j));
  551.         gsl_matrix_set(Hmat, j, i, val);
  552.         gsl_matrix_set(Hmat, i, j, val);}}
  553.  
  554.       gsl_eigen_gensymmv (Hmat, metric, evals, evec, w);
  555.  
  556.       cout << params->U << " ";
  557.       fprintf(out, "%lf ", params->U);
  558.  
  559.       double lowest=1e6;
  560.       for (int i=0; i<Ndimers; i++)
  561.       {
  562.         if (gsl_vector_get(evals, i) < lowest)
  563.         {
  564.           lowest = gsl_vector_get(evals, i);
  565.         }
  566.       }
  567.       cout << lowest << " ";
  568.       fprintf(out, "%.12lf\n", lowest);
  569.  
  570.  
  571.  auto dmrg_state = InitState(sites_dmrg);
  572.      for(auto i : range1(N))
  573.             {
  574.             if(i%2 == 1) dmrg_state.set(i,"Up");
  575.             else         dmrg_state.set(i,"Dn");
  576.             }
  577.  
  578.         auto psi0 = randomMPS(dmrg_state);
  579.  
  580.         auto ampo_dmrg = AutoMPO(sites_dmrg);
  581.         for(int b = 1; b < N; ++b)
  582.             {
  583.         ampo_dmrg += -1.0,"Cdagup",b,"Cup",b+1; //Hopping for spin up
  584.         ampo_dmrg += -1.0,"Cdagup",b+1,"Cup",b; // <-
  585.         ampo_dmrg += -1.0,"Cdagdn",b,"Cdn",b+1; //Hopping for spin down
  586.         ampo_dmrg += -1.0,"Cdagdn",b+1,"Cdn",b;// <-
  587.             ampo_dmrg += params->U, "Nup", b, "Ndn", b;
  588.             }
  589.                 ampo_dmrg += params->U, "Nup", N, "Ndn", N;
  590.             ampo_dmrg += -1.0,"Cdagup",1,"Cup",N; //Hopping for spin up
  591.                 ampo_dmrg += -1.0,"Cdagup",N,"Cup",1; // <-
  592.                 ampo_dmrg += -1.0,"Cdagdn",1,"Cdn",N; //Hopping for spin down
  593.                 ampo_dmrg += -1.0,"Cdagdn",N,"Cdn",1;// <-
  594.     MPO H_dmrg = toMPO(ampo_dmrg);
  595.         auto sweeps = Sweeps(10);
  596.         sweeps.maxdim() = 10,40,100,200,200;
  597.         sweeps.cutoff() = 1E-8;
  598.  
  599.         auto [energy_dmrg, psi] = dmrg(H_dmrg,psi0,sweeps,{"Silent",true});
  600.  
  601.  
  602. auto ampo = AutoMPO(sites);
  603.         for(int i = 1; i <= N; ++i)
  604.             {
  605.             ampo += "Ntot", i;
  606.             }
  607.         auto Ntot = toMPO(ampo);
  608.         cout << energy_dmrg << " " << endl; //<< inner(states.at(0), Ntot, states.at(0)) << endl;
  609.  
  610.  
  611. /*
  612.     auto s2 = AutoMPO(sites);
  613.     auto lattice = squareNextNeighbor(Nx,Ny,{"YPeriodic=", yperiodic});
  614.     for(int i=1; i<=Nx*Ny; i++)
  615.         {
  616.            
  617.                 //s2 += 0.5,"S+",i,"S-",i;
  618.                 //s2 += 0.5,"S-",i,"S+",i;
  619.                 //s2 +=     "Sz",i,"Sz",i;
  620.                 s2 += "Sz", i;
  621.         }
  622.     MPO S2 = toMPO(s2);
  623.     cout << "Sz = " << inner(states.at(0),S2,states.at(0)) << endl;    
  624.  
  625.     auto sz = AutoMPO(sites);
  626.     for(int i=1; i<=Nx*Ny; i++)
  627.         {
  628.  
  629.                 sz += 0.5,"S+",i,"S-",i;
  630.                 sz += 0.5,"S-",i,"S+",i;
  631.                 sz +=     "Sz",i,"Sz",i;
  632.         }
  633.         MPO SZ = toMPO(sz);
  634.         cout << "S^2 = " << inner(states.at(0),SZ, states.at(0)) << endl;
  635.   */    
  636. gsl_matrix *G; G = gsl_matrix_alloc(Nspins*2,Nspins*2);
  637. for (int p=0; p<Nspins; p++){
  638.         for (int q=0; q<Nspins; q++)
  639.         {for (int pp=0; pp<=1; pp++) { for (int qq=0;qq<=1;qq++) {
  640. gsl_matrix_set(G,2*p+pp,2*q+qq,0);}}}}
  641.  
  642. computeG(G, params->point,0);
  643. for (int p=0; p<Nspins; p++){
  644.         for (int q=0; q<Nspins; q++)
  645.         {for (int pp=0; pp<=1; pp++)  {
  646.                 double elem =gsl_matrix_get(G,2*p+pp,2*q+pp);
  647.                 cout << elem << " " << p << " " << q << " " << pp << endl;
  648.         }}}
  649.  
  650.  
  651.         params->U += 0.5;
  652.         params->point += 1;
  653.     }
  654.  
  655.    
  656. /*
  657. //Beginning DMRG
  658.      auto dmrg_state = InitState(sites);
  659.      for(auto i : range1(N))
  660.             {
  661.             if(i%2 == 1) dmrg_state.set(i,"Up");
  662.             else         dmrg_state.set(i,"Dn");
  663.             }
  664.  
  665.         auto psi0 = randomMPS(dmrg_state);
  666.  
  667.         auto ampo_dmrg = AutoMPO(sites);
  668.         for(int i = 1; i < N; ++i)
  669.             {
  670.             ampo_dmrg += -1, "Cdagup", i, "Cup", i+1;
  671.             ampo_dmrg += -1, "Cdagup", i+1, "Cup", i;
  672.             ampo_dmrg += -1, "Cdagdn", i, "Cdn", i+1;
  673.             ampo_dmrg += -1, "Cdagdn", i+1, "Cdn", i;
  674.             ampo_dmrg += params->U, "Nup", i, "Ndn", i;
  675.             }
  676.         ampo_dmrg += params->U, "Nup", N, "Ndn", N;
  677.             ampo_dmrg += -1, "Cdagup", N, "Cup", 1;
  678.             ampo_dmrg += -1, "Cdagup", 1, "Cup", N;
  679.             ampo_dmrg += -1, "Cdagdn", N, "Cdn", 1;
  680.             ampo_dmrg += -1, "Cdagdn", 1, "Cdn", N;
  681.         auto H_dmrg = toMPO(ampo_dmrg);
  682.  
  683.         auto sweeps = Sweeps(10);
  684.         sweeps.maxdim() = 10,40,100,200,200;
  685.         sweeps.cutoff() = 1E-8;
  686.  
  687.         auto [energy_dmrg, psi] = dmrg(H_dmrg,psi0,sweeps,{"Silent",true});
  688.  
  689.         cout << energy_dmrg << endl;*/
  690.  
  691.  
  692.  
  693.     /*double pt2 = eval_pt2(params);
  694.     cout << pt2+eval_element(params,0,0)  << endl;
  695.     fprintf(out, "%.12lf\n", pt2+eval_element(params,0,0));*/
  696.  
  697.  
  698.     fclose(out);
  699.   return 0;
  700.  
  701. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!