123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202 |
- /**********************************************************
-
- This software is part of J.-S. Caux's ABACUS library.
-
- Copyright (c) J.-S. Caux.
-
- -----------------------------------------------------------
-
- File: ABACUS_LiebLin.h
-
- Purpose: Declares LiebLin gas-related classes and functions.
-
- ***********************************************************/
-
- #ifndef ABACUS_LIEBLIN_H
- #define ABACUS_LIEBLIN_H
-
- #include "ABACUS.h"
-
- namespace ABACUS {
-
- // First, some global constants...
-
- const DP ITER_REQ_PREC_LIEBLIN = 1.0e+4 * MACHINE_EPS_SQ;
-
- const int LIEBLIN_Ix2_MIN = -1000000; // Like a UV cutoff. Assumption: never reached in scanning.
- const int LIEBLIN_Ix2_MAX = -LIEBLIN_Ix2_MIN;
-
-
- //***********************************************************************
-
- /**
- Eigenstates of the Lieb-Liniger model in the repulsive \f$(c > 0)\f$ regime.
-
- For numerical convenience, rapidities are rescaled by the interaction parameter \f$c\f$.
- In equations throughout this documentation, \f$\tilde{\lambda} \equiv \lambda/c\f$.
- In the code, these rescaled rapidities are denoted `lambdaoc` (read: "lambda over c").
- */
-
- class LiebLin_Bethe_State {
-
- public:
- DP c_int; // interaction parameter
- DP L;
- DP cxL;
- int N;
- std::string label;
- Vect<int> Ix2_available; // quantum numbers which are allowable but not occupied
- Vect<int> index_first_hole_to_right;
- Vect<int> displacement;
- Vect<int> Ix2;
- Vect<DP> lambdaoc;
- Vect<DP> S; // scattering sum
- Vect<DP> dSdlambdaoc; // its derivative
- DP diffsq;
- DP prec;
- int conv;
- int iter_Newton;
- DP E;
- int iK;
- DP K;
- DP lnnorm;
-
- public:
- LiebLin_Bethe_State ();
- LiebLin_Bethe_State (DP c_int_ref, DP L_ref, int N_ref);
- LiebLin_Bethe_State& operator= (const LiebLin_Bethe_State& RefState);
-
- public:
- int Charge () { return(N); };
- void Set_to_Label (std::string label_ref, const Vect<int>& OriginStateIx2);
- void Set_to_Label (std::string label_ref); // assumes OriginState == GroundState
- void Set_Label_from_Ix2 (const Vect<int>& OriginStateIx2);
- void Set_Label_Internals_from_Ix2 (const Vect<int>& OriginStateIx2);
- bool Check_Admissibility(char whichDSF); // always returns true
- void Find_Rapidities (bool reset_rapidities);
- bool Check_Rapidities(); // checks that all rapidities are not nan
- DP String_delta(); // trivially returns 0; exists to mirror spin chain function
- bool Check_Symmetry (); // checks whether set of quantum numbers obeys { I } = { -I }
- void Compute_lnnorm ();
- void Compute_All (bool reset_rapidities); // solves BAE, computes E, K and lnnorm
- void Set_Free_lambdaocs();
- void Iterate_BAE(DP damping);
- void Iterate_BAE_S(DP damping);
- void Iterate_BAE_Newton (DP damping, Vect_DP& RHSBAE, Vect_DP& dlambdaoc, SQMat_DP& Gaudin, Vect_INT& indx);
- void Compute_Energy ();
- void Compute_Momentum ();
- DP Kernel (int a, int b);
- DP Kernel (DP lambdaoc_ref);
- void Build_Reduced_Gaudin_Matrix (SQMat<DP>& Gaudin_Red);
- void Build_Reduced_BEC_Quench_Gaudin_Matrix (SQMat<DP>& Gaudin_Red);
- void Annihilate_ph_pair (int ipart, int ihole, const Vect<int>& OriginStateIx2);
- void Parity_Flip (); // takes all lambdaoc to -lambdaoc
- inline bool Set_to_Inner_Skeleton (int iKneeded, const Vect<int>& OriginIx2)
- {
- if (N < 3) ABACUSerror("N<3 incompatible with fixed momentum scanning");
- Ix2[0] = Ix2[1] - 2;
- Ix2[N-1] = Ix2[N-2] + 2;
- (*this).Compute_Momentum();
- if (iKneeded >= iK) Ix2[N-1] += 2*(iKneeded - iK);
- else Ix2[0] += 2*(iKneeded - iK);
- if (Ix2[0] < LIEBLIN_Ix2_MIN || Ix2[N-1] > LIEBLIN_Ix2_MAX) return(false);
- (*this).Set_Label_from_Ix2 (OriginIx2);
- return(true);
- }
- void Set_to_Outer_Skeleton (const Vect<int>& OriginIx2)
- {
- Ix2[0] = LIEBLIN_Ix2_MIN + (N % 2) + 1;
- Ix2[N-1] = LIEBLIN_Ix2_MAX - (N % 2) - 1;
- (*this).Set_Label_from_Ix2 (OriginIx2);
- };
- };
-
- inline bool Is_Inner_Skeleton (LiebLin_Bethe_State& State) {
- return (State.N >= 2 && (State.Ix2[0] == State.Ix2[1] - 2
- || State.Ix2[State.N-1] == State.Ix2[State.N-2] + 2));
- };
- inline bool Is_Outer_Skeleton (LiebLin_Bethe_State& State) {
- return (State.N >= 2 && State.Ix2[0] == LIEBLIN_Ix2_MIN + (State.N % 2) + 1
- && State.Ix2[State.N-1] == LIEBLIN_Ix2_MAX - (State.N % 2) - 1);
- };
- inline bool Is_Outer_Skeleton (const LiebLin_Bethe_State& State) {
- return (State.N >= 2 && State.Ix2[0] == LIEBLIN_Ix2_MIN + (State.N % 2) + 1
- && State.Ix2[State.N-1] == LIEBLIN_Ix2_MAX - (State.N % 2) - 1);
- };
-
- inline bool Force_Descent (char whichDSF, LiebLin_Bethe_State& ScanState,
- LiebLin_Bethe_State& RefState,
- int desc_type_required, int iKmod, DP Chem_Pot)
- {
- bool forcedesc = false;
-
- // Force descent if we're computing density-density, we're at zero momentum
- // and we're descending with momentum preserved:
- if (whichDSF == 'd' && ScanState.iK == RefState.iK && desc_type_required > 8) forcedesc = true;
-
- // For BEC to c > 0 quench, g2(x=0): force first step
- else if (whichDSF == 'B' && ScanState.label == RefState.label) forcedesc = true;
- else if (whichDSF == 'C' && ScanState.label == RefState.label) forcedesc = true;
-
- return(forcedesc);
- }
-
- std::ostream& operator<< (std::ostream& s, const LiebLin_Bethe_State& state);
-
- // FUNCTION DECLARATIONS:
-
- DP Chemical_Potential (LiebLin_Bethe_State& RefState);
- DP Sumrule_Factor (char whichDSF, LiebLin_Bethe_State& RefState, DP Chem_Pot,
- int iKmin, int iKmax);
- void Evaluate_F_Sumrule (char whichDSF, const LiebLin_Bethe_State& RefState, DP Chem_Pot,
- int iKmin, int iKmax, const char* RAW_Cstr, const char* FSR_Cstr);
- void Evaluate_F_Sumrule (std::string prefix, char whichDSF, const LiebLin_Bethe_State& RefState,
- DP Chem_Pot, int iKmin, int iKmax);
- void Evaluate_F_Sumrule (char whichDSF, DP c_int, DP L, int N, DP kBT, int nstates_req,
- DP Chem_Pot, int iKmin, int iKmax, const char* FSR_Cstr);
-
- // in LiebLin_Utils.cc
- DP LiebLin_dE0_dc (DP c_int, DP L, int N);
- DP LiebLin_vs (DP c_int, DP L, int N);
- DP LiebLin_Dressed_Charge_N (DP c_int, DP L, int N);
- int Momentum_Right_Excitations (LiebLin_Bethe_State& ScanState);
- int Momentum_Left_Excitations (LiebLin_Bethe_State& ScanState);
- DP ln_Overlap_with_BEC (LiebLin_Bethe_State& lambda);
-
- DP Particle_Hole_Excitation_Cost (char whichDSF, LiebLin_Bethe_State& AveragingState);
-
- std::complex<DP> ln_Density_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate);
- std::complex<DP> ln_Psi_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate);
- std::complex<DP> ln_g2_ME (LiebLin_Bethe_State& mu, LiebLin_Bethe_State& lambda);
-
- DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax,
- LiebLin_Bethe_State& LeftState,
- LiebLin_Bethe_State& RefState, DP Chem_Pot,
- std::stringstream& DAT_outfile);
-
- DP LiebLin_Twisted_lnnorm (Vect<std::complex<DP> >& lambdaoc, double cxL);
- std::complex<DP> LiebLin_Twisted_ln_Overlap (DP expbeta, Vect<DP> lstate_lambdaoc,
- DP lstate_lnnorm,
- LiebLin_Bethe_State& rstate);
- std::complex<DP> LiebLin_Twisted_ln_Overlap (std::complex<DP> expbeta,
- Vect<std::complex<DP> > lstate_lambdaoc,
- DP lstate_lnnorm,
- LiebLin_Bethe_State& rstate);
- std::complex<DP> LiebLin_ln_Overlap (Vect<DP> lstate_lambdaoc, DP lstate_lnnorm,
- LiebLin_Bethe_State& rstate);
- std::complex<DP> LiebLin_ln_Overlap (Vect<std::complex<DP> > lstate_lambdaoc, DP lstate_lnnorm,
- LiebLin_Bethe_State& rstate);
-
- // In src/LiebLin_Tgt0.cc:
- DP Entropy (LiebLin_Bethe_State& RefState);
- DP Canonical_Free_Energy (LiebLin_Bethe_State& RefState, DP kBT);
- LiebLin_Bethe_State Canonical_Saddle_Point_State (DP c_int, DP L, int N, DP kBT);
- LiebLin_Bethe_State Add_Particle_at_Center (const LiebLin_Bethe_State& RefState);
- LiebLin_Bethe_State Remove_Particle_at_Center (const LiebLin_Bethe_State& RefState);
- DP rho_of_lambdaoc_1 (LiebLin_Bethe_State& RefState, DP lambdaoc, DP delta);
- DP rho_of_lambdaoc_2 (LiebLin_Bethe_State& RefState, DP lambdaoc, DP delta);
-
-
- } // namespace ABACUS
-
- #endif
|