auto sites = Hubbard(Nspins); // Conserving all QNs

//The code below generates the Hubbard Hamiltonian for different U's
void generate_H()
{

  for (int ii=0; ii<=15; ii++)
  {
    double U = UU + ii*0.5;
    auto ampo = AutoMPO(sites);

        for(int b = 1; b < Nspins; ++b)
            {
                ampo += -1.0,"Cdagup",b,"Cup",b+1; //Hopping for spin up
                ampo += -1.0,"Cdagup",b+1,"Cup",b; // <-
                ampo += -1.0,"Cdagdn",b,"Cdn",b+1; //Hopping for spin down
                ampo += -1.0,"Cdagdn",b+1,"Cdn",b;// <-
                ampo += U, "Nup", b, "Ndn", b;
            }
                ampo += U, "Nup", Nspins, "Ndn", Nspins;
                ampo += -1.0,"Cdagup",1,"Cup",Nspins; //Hopping for spin up
                ampo += -1.0,"Cdagup",Nspins,"Cup",1; // <-
                ampo += -1.0,"Cdagdn",1,"Cdn",Nspins; //Hopping for spin down
                ampo += -1.0,"Cdagdn",Nspins,"Cdn",1;// <-
    MPO H_ = toMPO(ampo);
    H.at(ii) = H_;
  }
}

//Here the Phi_0 states on which I want to act upon with c_p!c_q. It is a product of dimers.
double create_state(const gsl_vector *cicoeff, void *params)
{
    int N = ((params_t *) params)->N;
    int index = ((params_t *) params)->index;

    auto psi = MPS(sites);
    double normalization;

    for(int n = 1; n <= N; n += 2)
        {

        auto s1 = sites(n);
        auto s2 = sites(n+1);

        normalization = 0;
        int pos = n/2*ms;
        for (int i=pos; i<=pos+ms-1; i++)
          normalization += pow(gsl_vector_get(cicoeff,i),2);

        auto wf = ITensor(s1,s2);
        wf.set(s1(2),s2(3), gsl_vector_get(cicoeff, pos)/sqrt(normalization));
        wf.set(s1(3),s2(2), gsl_vector_get(cicoeff, pos+1)/sqrt(normalization));
        wf.set(s1(1),s2(4), gsl_vector_get(cicoeff, pos+2)/sqrt(normalization));
        wf.set(s1(4),s2(1), gsl_vector_get(cicoeff, pos+3)/sqrt(normalization));

        psi.ref(n) = ITensor(s1);
        psi.ref(n+1) = ITensor(s2);
        auto [U,S,V] = svd(wf,{inds(psi.ref(n))},{"Cutoff=",1E-8});
        psi.ref(n) = U;
        psi.ref(n+1) = S*V;

      }

//The code below is for swapping the indices of the dimers in case they are not nearest-neighbors.
      int swaps = 1;
      vector<pair <int,int>> swaps_;
      vector <int> dimer;
      for (int i=0; i<N; i++)
        dimer.push_back(dimers[index][i]);

      while (swaps > 0)
      {
        swaps = 0;
        pair <int,int> pair1;
        for (int i=1; i<N; i++)
        {
          int pos1, pos2;
          pos1 = find_index(dimer, i);
          pos2 = find_index(dimer, i+1);
          if (pos1 > pos2)
          {
            pair1.first = dimer[pos1];
            pair1.second = dimer[pos2];
            swaps_.push_back(pair1);
            swaps += 1;
            int temp = dimer[pos1];
            dimer[pos1] = dimer[pos2];
            dimer[pos2] = temp;
            continue;
          }
        }
      }
      auto s1 = vector<Index>(1);
      auto s2 = vector<Index>(1);
      int ar1, ar2;
      ITensor wf, S,V,U;
      for (int i=swaps_.size()-1; i>=0; i--){
        ar1 = swaps_[i].first;
        ar2 = swaps_[i].second;
        s1.clear(); s1.push_back(siteIndex(psi,ar1));
        s2.clear(); s2.push_back(siteIndex(psi,ar2));
        auto is1 = IndexSet(s1);
        auto is2 = IndexSet(s2);

        psi.position(ar1);
        wf = psi.A(ar1)*psi.A(ar2);
        wf.swapInds(is1,is2);
        wf.noPrime();
        U = psi.A(ar1);
        svd(wf,U,S,V,{"Cutoff=",1E-8});
        psi.setA(ar1,U);
        psi.setA(ar2,S*V);
      }
      states.at(index) = psi;
      int point = ((params_t *) params)->point;
      return inner(states.at(index),H[point],states.at(index));
}

//Here is the code where all the problems arise.
void computeG(gsl_matrix *G, int point, int ind)
{
//The first few lines are for generating C!_pC_q with up's and down's.
        for (int i=1; i<=Nspins; i++)
        {
                for (int j=1; j<=Nspins; j++)
                {
                        for (int sgi=0; sgi<=1; sgi+=1) {
                                for (int sgj=0; sgj<=1; sgj+=1) {
                                         char Cd[7] = "Cdag"; char C[4] = "C";
                                        char Sgi[3], Sgj[3]; char up[3]="up"; char dn[3]="dn";
        if (sgi == 0) {
            strcpy(Sgi,"dn");
        }
        if (sgj == 0) {
            strcpy(Sgj,"dn");
        }
        if (sgi == 1) {
            strcpy(Sgi,"up");
        }
        if (sgj == 1) {
            strcpy(Sgj, "up");
        }
                                        MPS Phi_0 = removeQNs(states.at(ind)); //This is one of the state that is a product
                                                                                of dimers.
                                        MPO h = H[point]; //This is the Hubbard Hamiltonian for a specific U.
                                        auto ampo = AutoMPO(sites);
                                        if (sgi == sgj) {
                                        ampo += convertToString(Cd)+convertToString(Sgi), i,                ConvertToString(C)+convertToString(Sgj), j; //Generating C!_pC_q
                                        MPO cc = toMPO(ampo);
                                        MPS h_mps = applyMPO(h,Phi_0); MPS cc_mps = applyMPO(cc,Phi_0); //Acting with H and C!C
                                     cout << " " << inner(Phi_0,cc,h_mps) << " " << inner(Phi_0, h, cc_mps) << " " << i << j << endl;
                                        double elem = -(inner(Phi_0, h, cc_mps) - inner(Phi_0,cc,h_mps));
                                        gsl_matrix_set(G, 2*(i-1)+sgi,2*(j-1)+sgj,elem);}
                                }
                        }
                }
        }
}