diff --git a/src/SCAN/General_Scan.cc b/src/SCAN/General_Scan.cc index 88d878a..150c05b 100644 --- a/src/SCAN/General_Scan.cc +++ b/src/SCAN/General_Scan.cc @@ -34,11 +34,11 @@ namespace ABACUS { // Types of descendents: // 14 == iK stepped up, leading exc one step further (can lead to ph recombination) - // 13 == iK stepped up, next exc (nr of ph preserved, not taking possible ph recombination into account) - // 12 == iK stepped up, next ext (nr ph increased, not taking possible ph recombination into account) - // 11 == iK stepped down, leading exc one step further (can lead to ph recombination) - // 10 == iK stepped down, next exc (nr of ph preserved, not taking possible ph recombination into account) - // 9 == iK stepped down, next exc (nr ph increased, not taking possible ph recombination into account) + // 13 == iK up, next exc (nr of ph preserved, not taking possible ph recombination into account) + // 12 == iK up, next ext (nr ph increased, not taking possible ph recombination into account) + // 11 == iK down, leading exc one step further (can lead to ph recombination) + // 10 == iK down, next exc (nr of ph preserved, not taking possible ph recombination into account) + // 9 == iK down, next exc (nr ph increased, not taking possible ph recombination into account) // 8 == iK preserved, 14 and 11 (up to 2 ph recombinations) // 7 == iK preserved, 14 and 10 // 6 == iK preserved, 14 and 9 @@ -50,8 +50,10 @@ namespace ABACUS { // 0 == iK preserved, 12 and 9 // For scanning over symmetric states, the interpretation is slightly different. - // Types 14, 13 and 12 step iK up using the Ix2 on the right only, and mirrors the change on the left Ix2. - // Types 11, 10 and 9 step iK down using the Ix2 on the right only, and mirrors the change on the left Ix2. + // Types 14, 13 and 12 step iK up using the Ix2 on the right only, + // and mirrors the change on the left Ix2. + // Types 11, 10 and 9 step iK down using the Ix2 on the right only, + // and mirrors the change on the left Ix2. // There is then no need of scanning over types 0 - 8. // By convention, types 9, 10 and 11 can call types 9 - 14; types 12-14 can only call types 12-14. @@ -61,7 +63,8 @@ namespace ABACUS { // This function returns true if descending further can lead to a particle-hole recombination. // The criteria which are used are: // - the active excitation has moved at least one step (so it has already created its p-h pair) - // - there exists an OriginIx2 between the active Ix2 and the next Ix2 (to right or left depending on type of descendent) + // - there exists an OriginIx2 between the active Ix2 and the next Ix2 + // (to right or left depending on type of descendent) Vect > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2); @@ -71,7 +74,7 @@ namespace ABACUS { bool excfound = false; do { exclevel++; - if (exclevel == ScanIx2.size()) { // there isn't a single right-moving quantum number in ScanIx2 + if (exclevel == ScanIx2.size()) { // there is no right-moving quantum number in ScanIx2 break; } for (int alpha = 0; alpha < ScanIx2[exclevel].size(); ++alpha) @@ -125,7 +128,8 @@ namespace ABACUS { // This function returns true if descending further can lead to a particle-hole recombination. // The criteria which are used are: // - the active excitation has moved at least one step (so it has already created its p-h pair) - // - there exists an OriginIx2 between the active Ix2 and the next Ix2 (to right or left depending on type of descendent) + // - there exists an OriginIx2 between the active Ix2 and the next Ix2 + // (to right or left depending on type of descendent) Vect > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2); @@ -153,7 +157,8 @@ namespace ABACUS { if (excfound && !BaseScanIx2[exclevel].includes(ScanIx2[exclevel][excindex])) { // there exists an already dispersing excitation which isn't in Origin // Is there a possible recombination? - if (excindex > 0) { // a particle to the left of excitation has already moved left, so there is a hole + if (excindex > 0) { + // a particle to the left of excitation has already moved left, so there is a hole // check that there exists an occupied Ix2 in Origin sitting between the excitation // and the next Ix2 to its left in ScanIx2 for (int alpha = 0; alpha < BaseScanIx2[exclevel].size(); ++alpha) @@ -227,13 +232,18 @@ namespace ABACUS { int Max_Secs_alert = int(0.95 * Max_Secs); // we break any ongoing ithread loop beyond this point stringstream filenameprefix; - Data_File_Name (filenameprefix, whichDSF, iKmin, iKmax, kBT, AveragingState, SeedScanState, defaultScanStatename); + Data_File_Name (filenameprefix, whichDSF, iKmin, iKmax, kBT, + AveragingState, SeedScanState, defaultScanStatename); - if (in_parallel) for (int r = 0; r < paralevel; ++r) filenameprefix << "_" << rank[r] << "_" << nr_processors[r]; + if (in_parallel) + for (int r = 0; r < paralevel; ++r) + filenameprefix << "_" << rank[r] << "_" << nr_processors[r]; string prefix = filenameprefix.str(); - stringstream filenameprefix_prevparalevel; // without the rank and nr_processors of the highest paralevel + stringstream filenameprefix_prevparalevel; + // without the rank and nr_processors of the highest paralevel + Data_File_Name (filenameprefix_prevparalevel, whichDSF, iKmin, iKmax, kBT, AveragingState, SeedScanState, defaultScanStatename); if (in_parallel) for (int r = 0; r < paralevel - 1; ++r) @@ -330,7 +340,7 @@ namespace ABACUS { ScanStateList.Populate_List(whichDSF, SeedScanState); if (refine && !in_parallel) ScanStateList.Load_Info (SUM_Cstr); - else if (in_parallel && rank.sum() == 0) {}; // do nothing, keep the info in the higher .sum file! + else if (in_parallel && rank.sum() == 0) {}; // do nothing, keep info in the higher .sum file! DP Chem_Pot = Chemical_Potential (AveragingState); DP sumrule_factor = Sumrule_Factor (whichDSF, AveragingState, Chem_Pot, iKmin, iKmax); @@ -342,8 +352,8 @@ namespace ABACUS { int Ndata_previous_cycle = 0; - int ninadm = 0; // number of inadmissible states for which we save some info in .inadm file. Save first 1000. - int nconv0 = 0; // number of unconverged states for which we save some info in .conv0 file. Save first 1000. + int ninadm = 0; // nr of inadm states for which we save info in .inadm file. Save first 1000. + int nconv0 = 0; // nr of unconv states for which we save info in .conv0 file. Save first 1000. double start_time_omp = omp_get_wtime(); double current_time_omp = omp_get_wtime(); @@ -373,12 +383,14 @@ namespace ABACUS { double start_time_flags = omp_get_wtime(); // First flag the new base/type 's that we need to include: - ScanStateList.Raise_Scanning_Flags (exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0)); + ScanStateList.Raise_Scanning_Flags + (exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0)); // Get these base/type started: for (int i = 0; i < ScanStateList.ndef; ++i) { - if (ScanStateList.flag_for_scan[i] && ScanStateList.info[i].Ndata == 0 && !ScanStateList.scan_attempted[i]) { + if (ScanStateList.flag_for_scan[i] && ScanStateList.info[i].Ndata == 0 + && !ScanStateList.scan_attempted[i]) { Scan_Info scan_info_flags; diff --git a/src/TESTS/LiebLin_Rho_vs_Psidag_Psi.cc b/src/TESTS/LiebLin_Rho_vs_Psidag_Psi.cc new file mode 100644 index 0000000..be97765 --- /dev/null +++ b/src/TESTS/LiebLin_Rho_vs_Psidag_Psi.cc @@ -0,0 +1,186 @@ +/********************************************************** + +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 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 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 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); +}