187 lines
5.8 KiB
C++
187 lines
5.8 KiB
C++
/**********************************************************
|
|
|
|
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);
|
|
}
|