ABACUS/src/TESTS/LiebLin_Rho_vs_Psidag_Psi.cc

/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  ln_Density_ME.cc

Purpose:  Computes the density operator \rho(x = 0) matrix element

***********************************************************/

#include "ABACUS.h"

using namespace std;
using namespace ABACUS;


/**
   The matrix elements of \f$\hat{\rho}(x=0)\f$ and \f$\Psi^\dagger (x=0) \Psi(x=0)\f$
   are compared between two random bra and ket states. The density matrix element is
   directly computed from the determinant representation. The field operator product is
   computed by introducing a resolution of the identity. The function outputs the degree
   of accuracy achieved.
 */
int main()
{

  DP c_int = 4.344;
  DP L = 8.0;
  int N = 8;
  DP kBT = 8.0;

  int DiKmax = 10*N;

  // Use a thermal state to get some nontrivial distribution
  //LiebLin_Bethe_State spstate = Canonical_Saddle_Point_State (c_int, L, N, kBT);
  LiebLin_Bethe_State spstate (c_int, L, N);
  // spstate.Ix2[0] -= 8;
  // spstate.Ix2[1] -= 4;
  // spstate.Ix2[2] -= 2;
  // spstate.Ix2[N-3] += 2;
  // spstate.Ix2[N-2] += 4;
  // spstate.Ix2[N-1] += 6; // keep asymmetry
  spstate.Compute_All(true);
  //cout << setprecision(16) << spstate << endl;


  // Now define the anchor for scanning
  LiebLin_Bethe_State SeedScanState = Remove_Particle_at_Center (spstate);
  SeedScanState.Compute_All(true);
  //cout << setprecision(16) << SeedScanState << endl;


  // Map out the available unoccupied quantum nrs
  // We scan from Ix2min - DiKmax to Ix2max + DiKmax, and there are N-1 occupied, so there are
  int nr_par_positions = DiKmax + (SeedScanState.Ix2[N-2] - SeedScanState.Ix2[0])/2 + 1 - N + 1;
  Vect<int> Ix2_available (nr_par_positions);
  int nfound = 0;
  for (int i = SeedScanState.Ix2[0] - DiKmax; i <= SeedScanState.Ix2[N-2] + DiKmax; i += 2) {
    bool available = true;

    for (int j = 0; j < N-1; ++j)
      if (i == SeedScanState.Ix2[j]) {
	available = false;
	break;
      }
    if (available)
      Ix2_available[nfound++] = i;
  }
  if (nfound != nr_par_positions) {
    cout << "nround = " << nfound << "\tnr_par_positions = " << nr_par_positions << endl;
    ABACUSerror("Wrong number of particle positions found.");
  }


  // Define bra and ket states
  LiebLin_Bethe_State lstate = spstate;
  lstate.Ix2[0] -= 4;
  lstate.Ix2[1] -= 2;
  lstate.Compute_All(false);
  cout << setprecision(16) << lstate << endl;


  LiebLin_Bethe_State rstate = spstate;
  rstate.Ix2[N-1] += 8;
  rstate.Ix2[N-2] += 2;
  rstate.Compute_All(false);
  cout << setprecision(16) << rstate << endl;



  //lstate = rstate;

  // Matrix element to reproduce
  complex<DP> ln_rho_ME = ln_Density_ME(lstate, rstate);
  cout << "ln_Density_ME = " << ln_rho_ME << endl;

  // We now do hard-coded scanning of up to fixed nr of particle-hole pairs
  complex<DP> Psidag_Psi = 0.0;

  LiebLin_Bethe_State scanstate = SeedScanState;


  // Zero particle-hole:
  scanstate = SeedScanState;
  Psidag_Psi += exp(conj(ln_Psi_ME(scanstate, lstate)) + ln_Psi_ME(scanstate, rstate));

  cout << "Ratio Rho/PsidagPsi After scanning 0 p-h: "
       << setprecision(16) << real(Psidag_Psi/exp(ln_rho_ME))
       << "\t" << imag(Psidag_Psi/exp(ln_rho_ME)) << endl;


  // One particle-hole:
  char a;
  for (int ih1 = 0; ih1 < N-1; ++ih1)
    for (int ipar1 = 0; ipar1 < nr_par_positions - 1; ++ipar1) {

      scanstate = SeedScanState;
      scanstate.Ix2[ih1] = Ix2_available[ipar1];
      scanstate.Ix2.QuickSort();
      scanstate.Compute_All(true);
      Psidag_Psi += exp(conj(ln_Psi_ME(scanstate, lstate)) + ln_Psi_ME(scanstate, rstate));

      // cout << "ih1 = " << ih1 << "\tipar1 = " << ipar1
      // 	   << "\tPsidag_Psi = " << setprecision(16) << Psidag_Psi << endl;
      //cout << scanstate << endl;
      //cin >> a;
    }

  cout << "Ratio Rho/PsidagPsi After scanning 1 p-h: "
       << setprecision(16) << real(Psidag_Psi/exp(ln_rho_ME))
       << "\t" << imag(Psidag_Psi/exp(ln_rho_ME)) << endl;

  // Two particle-hole:
  for (int ih1 = 0; ih1 < N-2; ++ih1)
    for (int ih2 = ih1+1; ih2 < N-1; ++ih2)
      for (int ipar1 = 0; ipar1 < nr_par_positions - 2; ++ipar1) {
  	for (int ipar2 = ipar1+1; ipar2 < nr_par_positions - 1; ++ipar2) {

  	  scanstate = SeedScanState;
  	  scanstate.Ix2[ih1] = Ix2_available[ipar1];
  	  scanstate.Ix2[ih2] = Ix2_available[ipar2];
  	  scanstate.Ix2.QuickSort();
  	  scanstate.Compute_All(true);
  	  Psidag_Psi += exp(conj(ln_Psi_ME(scanstate, lstate)) + ln_Psi_ME(scanstate, rstate));

  	}
      }

  cout << "Ratio Rho/PsidagPsi After scanning 2 p-h: "
       << setprecision(16) << real(Psidag_Psi/exp(ln_rho_ME))
       << "\t" << imag(Psidag_Psi/exp(ln_rho_ME)) << endl;


  // Three particle-hole:
  for (int ih1 = 0; ih1 < N-3; ++ih1)
    for (int ih2 = ih1+1; ih2 < N-2; ++ih2)
      for (int ih3 = ih2+1; ih3 < N-1; ++ih3)
  	for (int ipar1 = 0; ipar1 < nr_par_positions - 3; ++ipar1) {
  	  for (int ipar2 = ipar1+1; ipar2 < nr_par_positions - 2; ++ipar2) {
  	    for (int ipar3 = ipar2+1; ipar3 < nr_par_positions - 1; ++ipar3) {

  	      scanstate = SeedScanState;
  	      scanstate.Ix2[ih1] = Ix2_available[ipar1];
  	      scanstate.Ix2[ih2] = Ix2_available[ipar2];
  	      scanstate.Ix2[ih3] = Ix2_available[ipar3];
  	      scanstate.Ix2.QuickSort();
  	      scanstate.Compute_All(true);
  	      Psidag_Psi += exp(conj(ln_Psi_ME(scanstate, lstate)) + ln_Psi_ME(scanstate, rstate));
  	    }
  	  }
  	  // cout << "Ratio Rho/PsidagPsi while scanning 3 p-h: "
  	  //      << setprecision(16) << real(Psidag_Psi/exp(ln_rho_ME))
  	  //      << "\t" << imag(Psidag_Psi/exp(ln_rho_ME)) << endl;
  	}

  cout << "Ratio Rho/PsidagPsi After scanning 3 p-h: "
       << setprecision(16) << real(Psidag_Psi/exp(ln_rho_ME))
       << "\t" << imag(Psidag_Psi/exp(ln_rho_ME)) << endl;


  return(0);
}