12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229 |
- /**********************************************************
-
- This software is part of J.-S. Caux's ABACUS library.
-
- Copyright (c) J.-S. Caux.
-
- -----------------------------------------------------------
-
- File: src/SCAN/General_Scan.cc
-
- Purpose: universal implementation of state scanning:
- functions to descend down hierarchy of intermediate states.
-
- NOTE: since templated functions have to be in the same file,
- we put all scanning functions here. The externally-accessible
- functions are defined at the end of this file.
-
- ***********************************************************/
-
- #include <omp.h>
- #include "ABACUS.h"
-
-
- using namespace std;
- using namespace ABACUS;
-
-
- string LABEL_TO_CHECK = "bla";
- //string LABEL_TO_CHECK = "6_0_";
- //string LABEL_TO_CHECK = "6_2_22y32";
-
-
- namespace ABACUS {
-
- // Types of descendents:
- // 14 == iK stepped up, leading exc one step further (can lead to ph recombination)
- // 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
- // 5 == iK preserved, 13 and 11
- // 4 == iK preserved, 13 and 10
- // 3 == iK preserved, 13 and 9
- // 2 == iK preserved, 12 and 11
- // 1 == iK preserved, 12 and 10
- // 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.
- // 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.
-
- bool Expect_ph_Recombination_iK_Up (string ScanIx2_label, const Vect<Vect<int> >& OriginIx2,
- const Vect<Vect<int> >& BaseScanIx2)
- {
- // 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)
-
- Vect<Vect<int> > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2);
-
- // Determine the level and index of the bottom-most left-most right-moving quantum number sits:
- int exclevel = -1;
- int excindex = 0;
- bool excfound = false;
- do {
- exclevel++;
- if (exclevel == ScanIx2.size()) { // there is no right-moving quantum number in ScanIx2
- break;
- }
- for (int alpha = 0; alpha < ScanIx2[exclevel].size(); ++alpha)
- if (ScanIx2[exclevel][alpha] > BaseScanIx2[exclevel][alpha]) {
- excindex = alpha;
- excfound = true;
- break;
- }
- } while (!excfound);
- // If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0;
-
- 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 < ScanIx2[exclevel].size() - 1) {
- // a particle to the right of excitation has already move right, so there is a hole
- // check that there exists an occupied Ix2 in Origin sitting between the excitation
- // and the next Ix2 to its right in ScanIx2
- for (int alpha = BaseScanIx2[exclevel].size() - 1; alpha >= 0; --alpha)
- if (BaseScanIx2[exclevel][alpha] > ScanIx2[exclevel][excindex]
- && BaseScanIx2[exclevel][alpha] < ScanIx2[exclevel][excindex + 1]) {
- return(true);
- }
- }
- } // if (excfound)
-
- return(false);
- }
-
- // Specialization for Lieb-Liniger:
- bool Expect_ph_Recombination_iK_Up (string ScanIx2_label, const LiebLin_Bethe_State& OriginState)
- {
- Vect<Vect<int> > OriginIx2here(1);
- OriginIx2here[0] = OriginState.Ix2;
- Vect<Vect<int> > BaseScanIx2here(1);
- BaseScanIx2here[0] = OriginState.Ix2;
- return(Expect_ph_Recombination_iK_Up (ScanIx2_label, OriginIx2here, BaseScanIx2here));
- }
-
- // Specialization for Heis
- bool Expect_ph_Recombination_iK_Up (string ScanIx2_label, const Heis_Bethe_State& OriginState)
- {
- return(Expect_ph_Recombination_iK_Up (ScanIx2_label, OriginState.Ix2, OriginState.Ix2));
- }
-
-
-
- bool Expect_ph_Recombination_iK_Down (string ScanIx2_label, const Vect<Vect<int> >& OriginIx2,
- const Vect<Vect<int> >& BaseScanIx2)
- {
- // 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)
-
- Vect<Vect<int> > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2);
-
- // Determine the level and index of the bottom-most right-most left-moving quantum number sits:
- int exclevel = -1;
- int excindex = 0;
- bool excfound = false;
-
- do {
- exclevel++;
- if (exclevel == ScanIx2.size()) { // there isn't a single left-moving quantum number in ScanIx2
- break;
- }
- for (int alpha = ScanIx2[exclevel].size() - 1; alpha >= 0; --alpha) {
- if (ScanIx2[exclevel][alpha] < BaseScanIx2[exclevel][alpha]) {
- excindex = alpha;
- excfound = true;
- break;
- }
- }
- } while (!excfound);
- // If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0;
- if (!excfound) excindex = ScanIx2[exclevel].size() - 1;
-
- 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
- // 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)
- if (BaseScanIx2[exclevel][alpha] > ScanIx2[exclevel][excindex - 1]
- && BaseScanIx2[exclevel][alpha] < ScanIx2[exclevel][excindex]) {
- return(true);
- }
- }
- } // if (excfound)
-
- return(false);
- }
-
- // Specialization for Lieb-Liniger:
- bool Expect_ph_Recombination_iK_Down (string ScanIx2_label, const LiebLin_Bethe_State& OriginState)
- {
- Vect<Vect<int> > OriginIx2here(1);
- OriginIx2here[0] = OriginState.Ix2;
- Vect<Vect<int> > BaseScanIx2here(1);
- BaseScanIx2here[0] = OriginState.Ix2;
- return(Expect_ph_Recombination_iK_Down (ScanIx2_label, OriginIx2here, BaseScanIx2here));
- }
-
- // Specialization for Heis
- bool Expect_ph_Recombination_iK_Down (string ScanIx2_label, const Heis_Bethe_State& OriginState)
- {
- return(Expect_ph_Recombination_iK_Down (ScanIx2_label, OriginState.Ix2, OriginState.Ix2));
- }
-
-
-
-
- template<class Tstate>
- Scan_Info General_Scan (char whichDSF, int iKmin, int iKmax, int iKmod, DP kBT,
- Tstate& AveragingState, Tstate& SeedScanState, string defaultScanStatename,
- int Max_Secs, DP target_sumrule, bool refine,
- int paralevel, Vect<int> rank, Vect<int> nr_processors)
- {
- // Performs the scan over excited states, writing data to file.
-
- // AveragingState is the state on which the correlations are calculated.
- // SeedScanState is the originator of all scan states.
- // This distinction is kept to allow for quenches and finite temperatures.
-
- // This function is also called by the parallel implementation of ABACUS.
- // In this case, file names carry a rank and nr_processors suffix.
- // In fact, the parallelization can be done in incremental levels.
- // If paralevel == 0, the run is serial.
- // If paralevel == n, the run is parallelized in a tree with n levels of branching.
- // A paralevel == 1 branching's files have a suffix of the form "_3_8", meaning that this
- // is the rank 3 out of 8 processors.
- // A paralevel == 2 branching's files have a suffix of the form "_3_8_2_8", meaning that this
- // is the rank 2 out of 8 subscan of the _3_8 scan.
-
- bool in_parallel = (paralevel > 0);
- if (in_parallel && (rank.size() != paralevel || nr_processors.size() != paralevel)) {
- cout << "paralevel = " << paralevel << "\trank.size() = " << rank.size()
- << "\tnr_processors.size() = " << nr_processors.size() << endl;
- cout << "rank = " << rank << endl;
- cout << "nr_processors = " << nr_processors << endl;
- ABACUSerror("Inconsistent paralevel, rank or nr_processors in General_Scan.");
- }
-
- if (in_parallel && !refine) ABACUSerror("Must refine when using parallel ABACUS");
-
- // expected cost on data_value of adding a particle-hole excitation.
- DP ph_cost = Particle_Hole_Excitation_Cost (whichDSF, AveragingState);
-
- int Max_Secs_used = int(0.9 * Max_Secs); // we don't start any new ithread loop beyond this point
- 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);
-
- 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
-
- Data_File_Name (filenameprefix_prevparalevel, whichDSF, iKmin, iKmax, kBT,
- AveragingState, SeedScanState, defaultScanStatename);
- if (in_parallel) for (int r = 0; r < paralevel - 1; ++r)
- filenameprefix << "_" << rank[r] << "_" << nr_processors[r];
-
- string prefix_prevparalevel = filenameprefix_prevparalevel.str();
-
- stringstream RAW_stringstream; string RAW_string;
- stringstream INADM_stringstream; string INADM_string;
- stringstream CONV0_stringstream; string CONV0_string;
- stringstream STAT_stringstream; string STAT_string;
- stringstream LOG_stringstream; string LOG_string;
- stringstream THR_stringstream; string THR_string;
- stringstream THRDIR_stringstream; string THRDIR_string;
- stringstream SRC_stringstream; string SRC_string;
- stringstream SUM_stringstream; string SUM_string;
-
- RAW_stringstream << prefix << ".raw";
- INADM_stringstream << prefix << ".inadm";
- CONV0_stringstream << prefix << ".conv0";
- STAT_stringstream << prefix << ".stat";
- LOG_stringstream << prefix << ".log";
- THR_stringstream << prefix << ".thr";
- THRDIR_stringstream << prefix << "_thrdir";
- SRC_stringstream << prefix << ".src";
- SUM_stringstream << prefix << ".sum";
-
- RAW_string = RAW_stringstream.str(); const char* RAW_Cstr = RAW_string.c_str();
- INADM_string = INADM_stringstream.str(); const char* INADM_Cstr = INADM_string.c_str();
- CONV0_string = CONV0_stringstream.str(); const char* CONV0_Cstr = CONV0_string.c_str();
- STAT_string = STAT_stringstream.str(); const char* STAT_Cstr = STAT_string.c_str();
- LOG_string = LOG_stringstream.str(); const char* LOG_Cstr = LOG_string.c_str();
- THR_string = THR_stringstream.str(); const char* THR_Cstr = THR_string.c_str();
- SRC_string = SRC_stringstream.str(); const char* SRC_Cstr = SRC_string.c_str();
- SUM_string = SUM_stringstream.str(); const char* SUM_Cstr = SUM_string.c_str();
-
- THRDIR_string = THRDIR_stringstream.str();
-
- fstream RAW_outfile;
- if (!refine || in_parallel) RAW_outfile.open(RAW_Cstr, fstream::out | fstream::trunc);
- else RAW_outfile.open(RAW_Cstr, fstream::out | fstream::app);
- if (RAW_outfile.fail()) {
- cout << RAW_Cstr << endl;
- ABACUSerror("Could not open RAW_outfile... ");
- }
- RAW_outfile.precision(16);
-
- fstream INADM_outfile;
- if (!refine || in_parallel) INADM_outfile.open(INADM_Cstr, fstream::out | fstream::trunc);
- else INADM_outfile.open(INADM_Cstr, fstream::out | fstream::app);
- if (INADM_outfile.fail()) ABACUSerror("Could not open INADM_outfile... ");
- INADM_outfile.precision(16);
-
- fstream CONV0_outfile;
- if (!refine || in_parallel) CONV0_outfile.open(CONV0_Cstr, fstream::out | fstream::trunc);
- else CONV0_outfile.open(CONV0_Cstr, fstream::out | fstream::app);
- if (CONV0_outfile.fail()) ABACUSerror("Could not open CONV0_outfile... ");
- CONV0_outfile.precision(16);
-
- fstream STAT_outfile;
- if (!refine || in_parallel) STAT_outfile.open(STAT_Cstr, fstream::out | fstream::trunc);
- else STAT_outfile.open(STAT_Cstr, fstream::out | fstream::app);
- if (STAT_outfile.fail()) ABACUSerror("Could not open STAT_outfile... ");
- STAT_outfile.precision(8);
-
- ofstream LOG_outfile;
- if (!in_parallel) {
- if (!refine) LOG_outfile.open(LOG_Cstr, fstream::out | fstream::trunc);
- else LOG_outfile.open(LOG_Cstr, fstream::out | fstream::app);
- if (LOG_outfile.fail()) ABACUSerror("Could not open LOG_outfile... ");
- LOG_outfile.precision(16);
- }
- else { // in_parallel
- LOG_outfile.open(LOG_Cstr, fstream::out | fstream::trunc);
- if (LOG_outfile.fail()) ABACUSerror("Could not open LOG_outfile... ");
- LOG_outfile.precision(16);
- }
-
- Scan_Info scan_info;
-
- if (!refine) mkdir(THRDIR_string.c_str(), S_IRWXU | S_IRWXG | S_IRWXO);
- Scan_Thread_Data paused_thread_data (THRDIR_string, refine);
-
- if (refine) {
- paused_thread_data.Load();
- if (!in_parallel) scan_info.Load(SRC_Cstr);
- }
-
- Scan_Info scan_info_before = scan_info; // for LOG file
- Scan_Info scan_info_before_descent;
- Scan_Info scan_info_obtained_in_descent;
-
- Scan_State_List<Tstate> ScanStateList (whichDSF, SeedScanState);
- ScanStateList.Populate_List(whichDSF, SeedScanState);
-
- if (refine && !in_parallel) ScanStateList.Load_Info (SUM_Cstr);
- 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);
-
-
- // Now go for it !
-
- bool at_least_one_new_flag_raised = false;
-
- int Ndata_previous_cycle = 0;
-
- 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();
-
- #pragma omp parallel
- do {
-
- int omp_thread_nr = omp_get_thread_num();
-
- if ((paused_thread_data.lowest_il_with_nthreads_neq_0 == paused_thread_data.nlists - 1)
- && omp_thread_nr > 0) {
- double start_time_wait = omp_get_wtime();
- double stop_time_wait;
- do {
- for (int i = 0; i < 100000; ++i) { }
- stop_time_wait = omp_get_wtime();
- } while (stop_time_wait - start_time_wait < 5.0);
- }
-
- double start_time_cycle_omp = omp_get_wtime();
-
- at_least_one_new_flag_raised = false;
-
-
- #pragma omp master
- {
- 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));
-
-
- // 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]) {
-
- Scan_Info scan_info_flags;
-
- at_least_one_new_flag_raised = true;
-
- ScanStateList.scan_attempted[i] = true;
-
- Tstate ScanState;
-
- ScanState = ScanStateList.State[i];
-
- DP data_value = -1.0;
-
- bool admissible = ScanState.Check_Admissibility(whichDSF);
-
- if (admissible) {
-
- ScanState.Compute_All(true);
-
- if (ScanState.conv) {
-
- // Put momentum in fundamental window, if possible:
- int iKexc = ScanState.iK - AveragingState.iK;
- while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod;
- while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod;
-
- stringstream rawfile_entry;
- data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState, Chem_Pot, rawfile_entry);
- {
- #pragma omp critical
- RAW_outfile << rawfile_entry.str();
- }
-
- {
- #pragma omp critical
- if (iKexc >= iKmin && iKexc <= iKmax) {
- scan_info_flags.Ndata++;
- scan_info_flags.Ndata_conv++;
- scan_info_flags.sumrule_obtained += data_value*sumrule_factor;
- }
- }
-
- // If we force descent: modify data_value by hand so that descent is forced on next scanning pass
- for (int itype = 0; itype < 15; ++itype) {
- DP data_value_used = 0.1* exp(-paused_thread_data.logscale * ABACUS::min(0, paused_thread_data.lowest_il_with_nthreads_neq_0));
- if (Force_Descent(whichDSF, ScanState, AveragingState, itype, iKmod, Chem_Pot))
- data_value = data_value_used;
- }
-
- Vect<bool> allowed(false, 15);
- if (whichDSF == 'B') { // symmetric state scanning
- allowed[9] = true; allowed[10] = true; allowed[11] = true;
- allowed[12] = true; allowed[13] = true; allowed[14] = true;
- }
- else {
- allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
- allowed[1] = allowed[0]; allowed[2] = allowed[0];
- allowed[3] = allowed[0]; allowed[4] = allowed[0];
- allowed[5] = allowed[0]; allowed[6] = allowed[0];
- allowed[7] = allowed[0]; allowed[8] = allowed[0];
- allowed[9] = (iKexc > iKmin);
- allowed[10] = allowed[9]; allowed[11] = allowed[9];
- allowed[12] = (iKexc < iKmax);
- allowed[13] = allowed[12]; allowed[14] = allowed[12];
- }
- for (int type_required_here = 0; type_required_here < 15; ++type_required_here) {
- if (!allowed[type_required_here]) continue;
- // All cases here are such that the ScanState hasn't been descended yet,
- // so we simply use data_value as expected data value:
- {
- #pragma omp critical
- paused_thread_data.Include_Thread (abs(data_value), ScanState.label, type_required_here);
- }
- }
- } // if (ScanState.conv)
-
- else {
- if (nconv0++ < 1000)
- CONV0_outfile << setw(25) << ScanState.label << setw(25) << ScanState.diffsq
- << setw(5) << ScanState.Check_Rapidities()
- << setw(25) << ScanState.String_delta() << endl;
- scan_info_flags.Ndata++;
- scan_info_flags.Ndata_conv0++;
- }
- } // if admissible
-
- else { // if inadmissible, modify data_value by hand so that descent is forced on next scanning pass
-
- if (ninadm++ < 10000000) INADM_outfile << ScanState.label << endl;
- scan_info_flags.Ndata++;
- scan_info_flags.Ninadm++;
-
- // Put momentum in fundamental window, if possible:
- int iKexc = ScanState.iK - AveragingState.iK;
- while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod;
- while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod;
-
- DP data_value = 1.0e-32;
- for (int itype = 0; itype < 15; ++itype)
- if (Force_Descent(whichDSF, ScanState, AveragingState, itype, iKmod, Chem_Pot))
- data_value = 0.1* exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0);
-
- Vect<bool> allowed(false, 15);
- if (whichDSF == 'B') {
- // We scan over symmetric states. Only types 14 down to 9 are allowed.
- allowed[9] = true; allowed[10] = true; allowed[11] = true;
- allowed[12] = true; allowed[13] = true; allowed[14] = true;
- }
- else {
- allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
- allowed[1] = allowed[0]; allowed[2] = allowed[0];
- allowed[3] = allowed[0]; allowed[4] = allowed[0];
- allowed[5] = allowed[0]; allowed[6] = allowed[0];
- allowed[7] = allowed[0]; allowed[8] = allowed[0];
- allowed[9] = (iKexc > iKmin);
- allowed[10] = allowed[9]; allowed[11] = allowed[9];
- allowed[12] = (iKexc < iKmax);
- allowed[13] = allowed[12]; allowed[14] = allowed[12];
- }
- for (int type_required_here = 0; type_required_here < 15; ++type_required_here) {
- if (!allowed[type_required_here]) continue;
- {
- #pragma omp critical
- paused_thread_data.Include_Thread (abs(data_value), ScanState.label, type_required_here);
- }
- }
- } // inadmissible
-
- scan_info_flags.TT += omp_get_wtime() - start_time_flags;
-
- // Put this info into the appropriate ScanStateList.info
- {
- #pragma omp critical
- ScanStateList.Include_Info(scan_info_flags, ScanStateList.base_label[i]);
- scan_info += scan_info_flags;
- }
-
- } // if flag_for_scan
- } // for i
-
- } // #pragma omp master
-
-
- // Now we deal with the previously existing paused threads:
-
- Vect<Scan_Thread> threads_to_do;
- int il_to_do = paused_thread_data.lowest_il_with_nthreads_neq_0; // for resaving threads in case we're out of time
- DP current_threshold = exp(-paused_thread_data.logscale * il_to_do);
- {
- #pragma omp critical
- threads_to_do = paused_thread_data.Extract_Next_Scan_Threads();
- }
-
- int ithread;
-
- {
- for (ithread = 0; ithread < threads_to_do.size(); ++ithread) {
-
- Scan_Info scan_info_this_ithread;
- double start_time_this_ithread = omp_get_wtime();
-
- // If we don't have time anymore, resave the threads instead of computing them:
- if (start_time_this_ithread - start_time_omp > Max_Secs_alert) {
- for (int ith = ithread; ith < threads_to_do.size(); ++ith) {
- #pragma omp critical
- paused_thread_data.Include_Thread (il_to_do, threads_to_do[ith].label, threads_to_do[ith].type);
- }
- break; // jump out of ithread loop
- }
-
- Tstate ScanState;
- {
- #pragma omp critical
- ScanState = ScanStateList.Return_State(Extract_Base_Label(threads_to_do[ithread].label));
- }
- Tstate BaseScanState; BaseScanState = ScanState;
-
- ScanState.Set_to_Label(threads_to_do[ithread].label, BaseScanState.Ix2);
-
- Tstate ScanStateBeingDescended = ScanState;
-
- // STARTING Descend_and_Compute block:
- int type_required = threads_to_do[ithread].type;
-
- ScanState.Compute_Momentum();
- Vect<string> desc_label;
-
- bool disperse_only_current_exc_up = false;
- if (type_required == 14 || type_required == 8 || type_required == 7 || type_required == 6)
- disperse_only_current_exc_up = true;
- bool preserve_nexc_up = false;
- if (type_required == 13 || type_required == 5 || type_required == 4 || type_required == 3)
- preserve_nexc_up = true;
- bool disperse_only_current_exc_down = false;
- if (type_required == 11 || type_required == 8 || type_required == 5 || type_required == 2)
- disperse_only_current_exc_down = true;
- bool preserve_nexc_down = false;
- if (type_required == 10 || type_required == 7 || type_required == 4 || type_required == 1)
- preserve_nexc_down = true;
-
- if (whichDSF == 'B') { // symmetric state scanning
- if (type_required >= 9 && type_required <= 11)
- desc_label = Descendent_States_with_iK_Stepped_Down_rightIx2only
- (ScanState.label, BaseScanState, disperse_only_current_exc_down, preserve_nexc_down);
- else if (type_required >= 12 && type_required <= 14)
- desc_label = Descendent_States_with_iK_Stepped_Up_rightIx2only
- (ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up);
- }
- else {
- if (type_required >= 0 && type_required <= 8) {
- desc_label = Descendent_States_with_iK_Preserved
- (ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up,
- disperse_only_current_exc_down, preserve_nexc_down);
- }
- else if (type_required >= 9 && type_required <= 11)
- desc_label = Descendent_States_with_iK_Stepped_Down
- (ScanState.label, BaseScanState, disperse_only_current_exc_down, preserve_nexc_down);
- else if (type_required >= 12 && type_required <= 14)
- desc_label = Descendent_States_with_iK_Stepped_Up
- (ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up);
- }
-
- string label_here = ScanState.label;
-
- for (int idesc = 0; idesc < desc_label.size(); ++idesc) {
-
- ScanState.Set_to_Label (desc_label[idesc], BaseScanState.Ix2);
-
- bool admissible = ScanState.Check_Admissibility(whichDSF);
-
- DP data_value = 0.0;
-
- ScanState.conv = false;
- ScanState.Compute_Momentum(); // since momentum is used as forced descent criterion
-
-
- if (admissible) {
-
- ScanState.Compute_All (idesc == 0);
-
- if (ScanState.conv) {
-
- // Put momentum in fundamental window, if possible:
- int iKexc = ScanState.iK - AveragingState.iK;
- while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod;
- while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod;
-
- stringstream rawfile_entry;
- data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState,
- Chem_Pot, rawfile_entry);
-
- {
- #pragma omp critical
- RAW_outfile << rawfile_entry.str();
- if (iKexc >= iKmin && iKexc <= iKmax) {
- scan_info_this_ithread.Ndata++;
- scan_info_this_ithread.Ndata_conv++;
- scan_info_this_ithread.sumrule_obtained += data_value*sumrule_factor;
- }
- }
-
- // Uncomment line below if .stat file is desired:
- // STAT_outfile << setw(20) << label_here << "\t" << setw(5) << type_required
- // << "\t" << setw(16) << std::scientific
- // << exp(-paused_thread_data.logscale * il_to_do)
- // << "\t" << setw(20) << ScanState.label << "\t" << setw(16) << data_value
- // << "\t" << setw(16) << std::fixed << setprecision(8)
- // << data_value/exp(-paused_thread_data.logscale * il_to_do) << endl;
-
- // Uncomment below if alerts for unexpectedly high data_value (as compared to threshold) are desired
- // if (fabs(data_value) > 10.0* current_threshold) {
- // cout << "\nAlert: data_value > 10* threshold, " << data_value << "\t" << current_threshold << endl;
- // cout << " for state " << ScanState.label << " descendent of type " << type_required
- // << " of state " << ScanStateBeingDescended.label << endl;
- // cout << AveragingState.Ix2 << endl;
- // cout << ScanStateBeingDescended.Ix2 << endl;
- // cout << ScanState.Ix2 << endl;
- // }
-
- } // if (ScanState.conv)
- else {
- if (nconv0++ < 1000)
- CONV0_outfile << setw(25) << ScanState.label << setw(25)
- << ScanState.diffsq << setw(5) << ScanState.Check_Rapidities()
- << setw(25) << ScanState.String_delta() << endl;
- scan_info_this_ithread.Ndata++;
- scan_info_this_ithread.Ndata_conv0++;
- }
- } // if (admissible)
-
- else {
- if (ninadm++ < 1000000) INADM_outfile << ScanState.label << endl;
- scan_info_this_ithread.Ndata++;
- scan_info_this_ithread.Ninadm++;
- }
-
- Tstate state_to_descend; state_to_descend = ScanState; // for checking
-
- ScanState.Compute_Momentum();
- // Put momentum in fundamental window, if possible:
- int iKexc = ScanState.iK - AveragingState.iK;
- while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod;
- while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod;
-
- // Momentum-preserving are only descended to momentum-preserving.
- // Momentum-increasing are only descended to momentum-preserving and momentum-increasing.
- // Momentum-decreasing are only descended to momentum-preserving and momentum-decreasing.
- Vect<bool> allowed(false, 15);
- if (whichDSF == 'B') {
- // We scan over symmetric states. Only types 14 down to 9 are allowed.
- if (type_required >= 9 && type_required <= 11) { // iK stepped down on rightIx2; step further up or down
- allowed[9] = true; allowed[10] = true; allowed[11] = true;
- allowed[12] = true; allowed[13] = true; allowed[14] = true;
- }
- else if (type_required >= 12 && type_required <= 14) { // iK stepped up on rightIx2; only step further up
- allowed[12] = true; allowed[13] = true; allowed[14] = true;
- }
- }
- else {
- if (type_required >= 0 && type_required <= 8) { // momentum-preserving
- allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
- allowed[9] = false;
- allowed[12] = false;
- }
- if (type_required >= 9 && type_required <= 11) { // momentum-decreasing
- allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
- allowed[9] = (iKexc > iKmin);
- allowed[12] = false;
- }
- if (type_required >= 12 && type_required <= 14) { // momentum-increasing
- allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
- allowed[9] = false;
- allowed[12] = (iKexc < iKmax);
- }
- // The others are just copies of the ones above:
- allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0];
- allowed[10] = allowed[9]; allowed[11] = allowed[9];
- allowed[13] = allowed[12]; allowed[14] = allowed[12];
- }
-
-
- for (int type_required_here = 0; type_required_here < 15; ++type_required_here) {
-
- if (!allowed[type_required_here]) continue;
-
- // Reset ScanState to what it was, if change on first pass
- if (type_required_here > 0) ScanState = state_to_descend;
-
- // We determine if we carry on scanning based on the data_value obtained, or forcing conditions:
- DP expected_abs_data_value = abs(data_value);
-
- //++G_7 logic
- if ((type_required_here == 14 || type_required_here == 8
- || type_required_here == 7 || type_required_here == 6)
- && Expect_ph_Recombination_iK_Up (ScanState.label, BaseScanState))
- expected_abs_data_value /= ph_cost;
- if (type_required_here == 12 || type_required_here == 2
- || type_required_here == 1 || type_required_here == 0)
- expected_abs_data_value *= ph_cost;
- if ((type_required_here == 11 || type_required_here == 8
- || type_required_here == 5 || type_required_here == 2)
- && Expect_ph_Recombination_iK_Down (ScanState.label, BaseScanState))
- expected_abs_data_value /= ph_cost;
- if (type_required_here == 9 || type_required_here == 6
- || type_required_here == 3 || type_required_here == 0)
- expected_abs_data_value *= ph_cost;
-
- {
- #pragma omp critical
- paused_thread_data.Include_Thread (expected_abs_data_value, ScanState.label, type_required_here);
- }
- } // for type_required_here
-
- } // for idesc
-
- // FINISHED Descend_and_Compute block
-
-
- scan_info_this_ithread.TT += omp_get_wtime() - start_time_this_ithread;
-
- #pragma omp critical
- {
- scan_info += scan_info_this_ithread;
- ScanStateList.Include_Info(scan_info_this_ithread, Extract_Base_Label(threads_to_do[ithread].label));
- }
- } // for ithread
-
- } // omp parallel region
-
-
- #pragma omp master
- {
- if (!in_parallel)
- LOG_outfile << "Master cycling. Ndata_conv " << scan_info.Ndata_conv
- << ". Threshold " << paused_thread_data.lowest_il_with_nthreads_neq_0 << " "
- << setw(9) << setprecision(3)
- << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0)
- << ". " << setw(12) << scan_info.Ndata - Ndata_previous_cycle << " new data. Nr of threads: "
- << setw(14) << paused_thread_data.nthreads_total.sum()
- << ". Saturation: " << setprecision(12) << scan_info.sumrule_obtained << endl;
-
- Ndata_previous_cycle = scan_info.Ndata;
- }
-
-
- current_time_omp = omp_get_wtime();
-
- } while (current_time_omp - start_time_omp < Max_Secs_used
- && scan_info.sumrule_obtained < target_sumrule
- );
- // This closes the #pragram omp parallel block
-
- RAW_outfile.close();
- INADM_outfile.close();
- CONV0_outfile.close();
- STAT_outfile.close();
-
- scan_info.Save(SRC_Cstr);
-
- Scan_Info scan_info_refine = scan_info;
- scan_info_refine -= scan_info_before;
-
- if (!in_parallel) {
-
- if (scan_info.sumrule_obtained >= target_sumrule)
- LOG_outfile << endl << "Achieved sumrule saturation of " << scan_info.sumrule_obtained
- << "\t(target was " << target_sumrule << ")." << endl << endl;
-
- if (!refine) {
- LOG_outfile << "Main run info: " << scan_info << endl;
- LOG_outfile << "Latest threshold level " << paused_thread_data.lowest_il_with_nthreads_neq_0
- << " " << std::scientific << setprecision(3)
- << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl;
- }
- else if (refine) {
- LOG_outfile << "Refining info: " << scan_info_refine << endl;
- LOG_outfile << "Latest threshold level " << paused_thread_data.lowest_il_with_nthreads_neq_0
- << " " << std::scientific << setprecision(3)
- << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl;
- LOG_outfile << "Resulting info: " << scan_info << endl;
- }
- time_t current_time = time(nullptr);
- char timestr[100];
- strftime(timestr, sizeof(timestr), "%Y-%m-%d %H:%M:%S", gmtime(¤t_time));
- LOG_outfile << "Run completion timestamp: " << timestr << " UTC" << endl;
- LOG_outfile << "ABACUS version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl << endl;
- LOG_outfile.close();
- }
-
- else { // in_parallel
-
- LOG_outfile << "rank " << rank << " out of " << nr_processors << " processors: "
- << "run info: " << scan_info << endl << "Latest threshold = "
- << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl;
- }
-
- paused_thread_data.Save();
-
- ScanStateList.Order_in_SRC ();
-
- ScanStateList.Save_Info (SUM_Cstr);
-
-
- // Evaluate f-sumrule:
- if (!in_parallel ) if (whichDSF != 'q')
- Evaluate_F_Sumrule (prefix_prevparalevel, whichDSF, AveragingState, Chem_Pot, iKmin, iKmax);
-
- return(scan_info);
- }
-
-
-
- //******************************************************
-
- // Functions to initiate scans:
-
-
- // General version for equilibrium correlators at generic (possibly finite) temperature:
-
- Scan_Info Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
- int Max_Secs, DP target_sumrule, bool refine,
- int paralevel, Vect<int> rank, Vect<int> nr_processors)
- {
-
- // This function scans the Hilbert space of the LiebLin gas,
- // for the function identified by whichDSF.
-
- // whichDSF == 'Z': canonical partition function
- // whichDSF == 'd': density-density correlation function
- // whichDSF == 'g': Green's function < \Psi \Psi^{\dagger}>
- // whichDSF == 'o': one-body function < \Psi^{\dagger} \Psi >
-
- // Delta is the number of sites involved in the smoothing of the entropy
- //int Delta = int(sqrt(N))/2;//6;//N/20;
- //DP epsilon = log(L)/L; // using Gaussian for density in entropy.
- //DP epsilon = 1.0/L; // using Lorentzian for density in entropy.
-
- // Construct the finite-size saddle-point state:
- // if we refine, read the quantum numbers of the saddle point state (and seed sps) from the sps file:
-
- stringstream SPS_stringstream; string SPS_string;
- Data_File_Name (SPS_stringstream, whichDSF, c_int, L, N, iKmin, iKmax, kBT, 0.0, "");
- SPS_stringstream << ".sps";
- SPS_string = SPS_stringstream.str(); const char* SPS_Cstr = SPS_string.c_str();
-
- fstream spsfile;
- if (refine) spsfile.open(SPS_Cstr, fstream::in);
- else spsfile.open(SPS_Cstr, fstream::out | fstream::trunc);
- if (spsfile.fail()) {
- cout << SPS_Cstr << endl; ABACUSerror("Could not open spsfile.");
- }
-
- LiebLin_Bethe_State spstate;
-
- if (!refine) { // obtain the sps from discretized TBA
- //spstate = Canonical_Saddle_Point_State (c_int, L, N, whichDSF == 'Z' ? 0.0 : kBT);
- spstate = Canonical_Saddle_Point_State (c_int, L, N, kBT);
- }
- else { // read it from the sps file
- // Check that the sps has the right number of Ix2:
- int Nspsread;
- spsfile >> Nspsread;
- if (Nspsread != N) {
- cout << Nspsread << "\t" << N << endl;
- ABACUSerror("Wrong number of Ix2 in saddle-point state.");
- }
- spstate = LiebLin_Bethe_State (c_int, L, N);
- for (int i = 0; i < N; ++i) spsfile >> spstate.Ix2[i];
- }
- spstate.Compute_All(true);
-
- int Nscan = N;
- if (whichDSF == 'o') Nscan = N - 1;
- if (whichDSF == 'g') Nscan = N + 1;
-
- // Now construct or read off the seed scan state:
-
- // TO MODIFY: this is not a good idea, since this might construct a state with many p-h w/r to the AveragingState.
- LiebLin_Bethe_State SeedScanState;
-
- if (whichDSF != 'o' && whichDSF != 'g') SeedScanState = spstate;
-
- else if (whichDSF == 'o' || whichDSF == 'g') {
- if (!refine) {
- if (whichDSF == 'o') SeedScanState = Remove_Particle_at_Center (spstate);
- else SeedScanState = Add_Particle_at_Center (spstate);
- }
- else { // read it from the sps file
- // Check that the sps has the right number of Ix2:
- int Nsspsread;
- spsfile >> Nsspsread;
- if (Nsspsread != Nscan) {
- cout << Nsspsread << "\t" << Nscan << endl;
- ABACUSerror("Wrong number of Ix2 in scan saddle-point state.");
- }
- SeedScanState = LiebLin_Bethe_State (c_int, L, Nscan);
- for (int i = 0; i < Nscan; ++i) spsfile >> SeedScanState.Ix2[i];
- }
- } // if one-body or Green's function
-
- SeedScanState.Compute_All(true);
-
- LiebLin_Bethe_State ScanState = SeedScanState;
-
- DP delta = sqrt(DP(N)) * (spstate.lambdaoc[N-1] - spstate.lambdaoc[0])/N;
-
- if (!refine) { // we write data to the sps file
-
- spsfile << N << endl;
- spsfile << spstate.Ix2 << endl;
- spsfile << Nscan << endl;
- spsfile << SeedScanState.Ix2 << endl;
-
- spsfile << endl << spstate << endl << endl;
- for (int i = 1; i < spstate.N - 2; ++i)
- spsfile << 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1])
- << "\t" << 1.0/spstate.L * (0.25/(spstate.lambdaoc[i] - spstate.lambdaoc[i-1])
- + 0.5/(spstate.lambdaoc[i+1] - spstate.lambdaoc[i])
- + 0.25/(spstate.lambdaoc[i+2] - spstate.lambdaoc[i+1]))
- << "\t" << rho_of_lambdaoc_1 (spstate, 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1]), delta)
- << "\t" << rho_of_lambdaoc_2 (spstate, 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1]), delta)
- << endl;
- }
- spsfile.close();
-
- // Perform the scan:
- return General_Scan (whichDSF, iKmin, iKmax, 100000000, kBT, spstate, SeedScanState, "",
- Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
- }
-
- Scan_Info Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
- int Max_Secs, DP target_sumrule, bool refine)
- {
- int paralevel = 0;
- Vect<int> rank(0,1);
- Vect<int> nr_processors(0,1);
-
- return Scan_LiebLin (whichDSF, c_int, L, N, iKmin, iKmax, kBT, Max_Secs, target_sumrule,
- refine, paralevel, rank, nr_processors);
- }
-
-
- // Scanning on an excited state defined by a set of Ix2:
- Scan_Info Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename,
- int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine,
- int paralevel, Vect<int> rank, Vect<int> nr_processors)
- {
- // This function is as Scan_LiebLin for generic T defined above, except that the
- // averaging is now done on a state defined by AveragingStateIx2
-
- // PRECONDITIONS:
- // - the Ix2 of AveragingState are properly set.
-
- DP c_int = AveragingState.c_int;
- DP L = AveragingState.L;
- int N = AveragingState.N;
-
- // The label of the Averaging State is by definition the `empty' label
- AveragingState.Set_Label_from_Ix2 (AveragingState.Ix2);
- AveragingState.Compute_All(true);
-
- int Nscan = N;
- if (whichDSF == 'o') Nscan = N - 1;
- if (whichDSF == 'g') Nscan = N + 1;
-
- LiebLin_Bethe_State SeedScanState (c_int, L, Nscan);
- if (whichDSF == 'd' || whichDSF == 'B') SeedScanState.Ix2 = AveragingState.Ix2;
- // If 'o', remove midmost and shift quantum numbers by half-integer towards removed one:
- if (whichDSF == 'o') {
- for (int i = 0; i < N-1; ++i)
- SeedScanState.Ix2[i] = AveragingState.Ix2[i + (i >= N/2)] + 1 - 2*(i >= N/2);
- }
- // If 'g', add a quantum number in middle (explicitly: to right of index N/2)
- // and shift quantum numbers by half-integer away from added one:
- if (whichDSF == 'g') {
- SeedScanState.Ix2[N/2] = AveragingState.Ix2[N/2] - 1;
- for (int i = 0; i < N+1; ++i)
- SeedScanState.Ix2[i + (i >= N/2)] = AveragingState.Ix2[i] - 1 + 2*(i >= N/2);
- }
- SeedScanState.Compute_All(true);
-
- SeedScanState.Set_Label_from_Ix2 (SeedScanState.Ix2);
-
- DP kBT = 0.0;
-
- // Perform the scan:
- return General_Scan (whichDSF, iKmin, iKmax, 100000000, kBT,
- AveragingState, SeedScanState, defaultScanStatename,
- Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
- }
-
- // Simplified function call of the above:
- Scan_Info Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename,
- int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine)
- {
- int paralevel = 0;
- Vect<int> rank(0,1);
- Vect<int> nr_processors(0,1);
-
- return Scan_LiebLin (whichDSF, AveragingState, defaultScanStatename, iKmin, iKmax, Max_Secs,
- target_sumrule, refine, paralevel, rank, nr_processors);
- }
-
- // Scanning on a previously-defined AveragingState
- Scan_Info Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, string defaultScanStatename,
- int iKmin, int iKmax,
- int Max_Secs, DP target_sumrule, bool refine,
- int paralevel, Vect<int> rank, Vect<int> nr_processors)
- {
- // General state scanning for Heisenberg chains
-
- // PRECONDITIONS:
- // - the Ix2 of AveragingState are properly set.
-
- // Prepare the AveragingState:
- AveragingState.Compute_All(true);
-
- XXZ_Bethe_State SeedScanState;
- if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState;
- else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState);
- else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState);
- else ABACUSerror("Unknown whichDSF in Scan_Heis.");
-
- // Now the scan itself
- General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
- defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
-
- }
-
- // Scanning on a previously-defined AveragingState
- Scan_Info Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, string defaultScanStatename,
- int iKmin, int iKmax,
- int Max_Secs, DP target_sumrule, bool refine,
- int paralevel, Vect<int> rank, Vect<int> nr_processors)
- {
- // General state scanning for Heisenberg chains
-
- // PRECONDITIONS:
- // - the Ix2 of AveragingState are properly set.
-
- // Prepare the AveragingState:
- AveragingState.Compute_All(true);
-
- XXX_Bethe_State SeedScanState;
- if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState;
- else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState);
- else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState);
- else ABACUSerror("Unknown whichDSF in Scan_Heis.");
-
- // Now the scan itself
- return General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
- defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
-
- }
-
- // Scanning on a previously-defined AveragingState
- Scan_Info Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, string defaultScanStatename,
- int iKmin, int iKmax,
- int Max_Secs, DP target_sumrule, bool refine,
- int paralevel, Vect<int> rank, Vect<int> nr_processors)
- {
- // General state scanning for Heisenberg chains
-
- // PRECONDITIONS:
- // - the Ix2 of AveragingState are properly set.
-
- // Prepare the AveragingState:
- AveragingState.Compute_All(true);
-
- XXZ_gpd_Bethe_State SeedScanState;
- if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState;
- else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState);
- else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState);
- else ABACUSerror("Unknown whichDSF in Scan_Heis.");
-
- // Now the scan itself
- return General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
- defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
-
- }
-
-
- Scan_Info Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
- int Max_Secs, DP target_sumrule, bool refine,
- int paralevel, Vect<int> rank, Vect<int> nr_processors)
- {
- // This function scans the Hilbert space of the Heisenberg spin-1/2 chain
- // for the function identified by whichDSF.
-
- // whichDSF == 'Z': canonical partition function
- // whichDSF == 'm': S^{-+}
- // whichDSF == 'z': S^{zz}
- // whichDSF == 'p': S^{+-}
- // whichDSF == 'a': < S^z_j S^z_{j+1} S^z_l S^z_{l+1} > for RIXS
- // whichDSF == 'b': < S^z_j S^-_{j+1} S^-_l S^z_{l+1} > + (m <-> z) for RIXS
- // whichDSF == 'c': < S^-_j S^-_{j+1} S^-_l S^-_{l+1} > for RIXS
-
- Heis_Chain BD1(1.0, Delta, 0.0, N);
-
- Vect_INT Nrapidities_groundstate(0, BD1.Nstrings);
-
- Nrapidities_groundstate[0] = M;
-
- Heis_Base baseconfig_groundstate(BD1, Nrapidities_groundstate);
-
-
- if ((Delta > 0.0) && (Delta < 1.0)) {
-
- XXZ_Bethe_State GroundState(BD1, baseconfig_groundstate);
- GroundState.Compute_All(true);
- // The ground state is now fully defined.
-
- XXZ_Bethe_State SeedScanState;
- if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = GroundState;
- else if (whichDSF == 'm') SeedScanState = XXZ_Bethe_State(GroundState.chain, M - 1);
- else if (whichDSF == 'p') SeedScanState = XXZ_Bethe_State(GroundState.chain, M + 1);
- else ABACUSerror("Unknown whichDSF in Scan_Heis.");
-
- // Now the scan itself
- General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "",
- Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
-
- }
-
- else if (Delta == 1.0) {
-
- XXX_Bethe_State GroundState(BD1, baseconfig_groundstate);
- GroundState.Compute_All(true);
- // The ground state is now fully defined.
-
- XXX_Bethe_State SeedScanState;
- if (whichDSF == 'Z' || whichDSF == 'z' || whichDSF == 'a' || whichDSF == 'q') SeedScanState = GroundState;
- else if (whichDSF == 'm') SeedScanState = XXX_Bethe_State(GroundState.chain, M - 1);
- else if (whichDSF == 'p') SeedScanState = XXX_Bethe_State(GroundState.chain, M + 1);
- else if (whichDSF == 'c') SeedScanState = XXX_Bethe_State(GroundState.chain, M - 2);
- else ABACUSerror("Unknown whichDSF in Scan_Heis.");
-
- // Now the scan itself
- General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "",
- Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
- }
-
- else if (Delta > 1.0) {
-
- XXZ_gpd_Bethe_State GroundState(BD1, baseconfig_groundstate);
- GroundState.Compute_All(true);
- // The ground state is now fully defined.
-
- XXZ_gpd_Bethe_State SeedScanState;
- if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = GroundState;
- else if (whichDSF == 'm') SeedScanState = XXZ_gpd_Bethe_State(GroundState.chain, M - 1);
- else if (whichDSF == 'p') SeedScanState = XXZ_gpd_Bethe_State(GroundState.chain, M + 1);
- else ABACUSerror("Unknown whichDSF in Scan_Heis.");
-
- // Now the scan itself
- return General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "",
- Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
- }
-
- else ABACUSerror("Delta out of range in Heis_Structure_Factor");
-
- return Scan_Info();
- }
-
- Scan_Info Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
- int Max_Secs, DP target_sumrule, bool refine)
- {
- int paralevel = 0;
- Vect<int> rank(0,1);
- Vect<int> nr_processors(0,1);
-
- return Scan_Heis (whichDSF, Delta, N, M, iKmin, iKmax,
- Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
- }
-
-
- } // namespace ABACUS
|