123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436 |
- /**********************************************************
-
- This software is part of J.-S. Caux's ABACUS library.
-
- Copyright (c) J.-S. Caux.
-
- -----------------------------------------------------------
-
- File: ABACUS_Heis.h
-
- Purpose: Declares Heisenberg chain classes and functions.
-
- ***********************************************************/
-
- #ifndef ABACUS_HEIS_H
- #define ABACUS_HEIS_H
-
- #include "ABACUS.h"
-
- namespace ABACUS {
-
- // First, some global constants...
-
- const long long int ID_UPPER_LIMIT = 10000000LL; // max size we can define without seg fault
- const int INTERVALS_SIZE = 100000; // size of Scan_Intervals arrays
- const int NBASESMAX = 1000; // max number of bases kept
-
- const DP ITER_REQ_PREC = 100.0 * MACHINE_EPS_SQ;
-
- // Cutoffs on particle numbers
- const int MAXSTRINGS = 20; // maximal number of particle types we allow in bases
-
- const int NEXC_MAX_HEIS = 16; // max nr of excitations (string binding/unbinding, particle-hole) considered
-
-
- //***********************************************************************
-
- class Heis_Chain {
-
- public:
- DP J;
- DP Delta;
- DP anis; // acos(Delta) if Delta < 1.0, 0 if Delta == 1.0, acosh(Delta) if Delta > 1.0
- DP hz;
- int Nsites;
- int Nstrings; // how many possible strings. The following two arrays have Nstrings nonzero elements.
- int* Str_L; // vector (length M) containing the allowed string lengths. Elements that are 0 have no meaning.
- int* par; // vector (length M) containing the parities of the strings. Elements that are 0 have no meaning.
- // Parities are all +1 except for gapless XXZ subcases
- DP* si_n_anis_over_2; // for optimization: sin for XXZ, sinh for XXZ_gpd
- DP* co_n_anis_over_2; // for optimization
- DP* ta_n_anis_over_2; // for optimization
- DP prec; // precision required for computations, always put to ITER_REQ_PREC
-
- public:
- Heis_Chain ();
- Heis_Chain (DP JJ, DP DD, DP hh, int NN); // contructor: simply initializes
- Heis_Chain (const Heis_Chain& RefChain); // copy constructor;
- Heis_Chain& operator= (const Heis_Chain& RefChain);
- bool operator== (const Heis_Chain& RefChain);
- bool operator!= (const Heis_Chain& RefChain);
- ~Heis_Chain(); // destructor
-
- };
-
-
- //****************************************************************************
-
- // Objects in class Heis_Base are a checked vector containing the number of rapidities
- // of allowable types for a given state
-
- class Heis_Base {
-
- public:
- int Mdown; // total number of down spins
- Vect<int> Nrap; // Nrap[i] contains the number of rapidities of type i, i = 0, Nstrings - 1.
- int Nraptot; // total number of strings in this state
- Vect<DP> Ix2_infty; // Ix2_infty[i] contains the max of BAE function for the (half-)integer I[i], i = 0, Nstrings - 1.
- Vect<int> Ix2_min;
- Vect<int> Ix2_max; // Ix2_max[i] contains the integer part of 2*I_infty, with correct parity for base.
- double dimH; // dimension of sub Hilbert space associated to this base; use double to avoid max int problems.
- std::string baselabel; // base label
-
- public:
- Heis_Base ();
- Heis_Base (const Heis_Base& RefBase); // copy constructor
- Heis_Base (const Heis_Chain& RefChain, int M); // constructs configuration with all Mdown in one-string of +1 parity
- Heis_Base (const Heis_Chain& RefChain, const Vect<int>& Nrapidities); // sets to Nrapidities vector, and checks consistency
- Heis_Base (const Heis_Chain& RefChain, std::string baselabel_ref);
- inline int& operator[] (const int i);
- inline const int& operator[] (const int i) const;
- Heis_Base& operator= (const Heis_Base& RefBase);
- bool operator== (const Heis_Base& RefBase);
- bool operator!= (const Heis_Base& RefBase);
-
- void Compute_Ix2_limits(const Heis_Chain& RefChain); // computes the Ix2_infty and Ix2_max
-
- };
-
- inline int& Heis_Base::operator[] (const int i)
- {
- return Nrap[i];
- }
-
- inline const int& Heis_Base::operator[] (const int i) const
- {
- return Nrap[i];
- }
-
-
- //****************************************************************************
-
- // Objects in class Lambda carry all rapidities of a state
-
- class Lambda {
-
- private:
- int Nstrings;
- Vect<int> Nrap;
- int Nraptot;
- DP** lambda;
-
- public:
- Lambda ();
- Lambda (const Heis_Chain& RefChain, int M); // constructor, puts all lambda's to zero
- Lambda (const Heis_Chain& RefChain, const Heis_Base& base); // constructor, putting I's to lowest-energy config
- // consistent with Heis_Base configuration for chain RefChain
- Lambda& operator= (const Lambda& RefConfig);
- inline DP* operator[] (const int i);
- inline const DP* operator[] (const int i) const;
- ~Lambda();
-
- };
-
- inline DP* Lambda::operator[] (const int i)
- {
- return lambda[i];
- }
-
- inline const DP* Lambda::operator[] (const int i) const
- {
- return lambda[i];
- }
-
-
- //****************************************************************************
-
- // Objects in class Heis_Bethe_State carry all information about an eigenstate
-
- // Derived classes include XXZ_Bethe_State, XXX_Bethe_State, XXZ_gpd_Bethe_State
- // These contain subclass-specific functions and data.
-
- class Heis_Bethe_State {
-
- public:
- Heis_Chain chain;
- Heis_Base base;
- Vect<Vect<int> > Ix2;
- Lambda lambda;
- Lambda deviation; // string deviations
- Lambda BE; // Bethe equation for relevant rapidity, in the form BE = theta - (1/N)\sum ... - \pi I/N = 0
- DP diffsq; // sum of squares of rapidity differences in last iteration
- int conv; // convergence status
- DP dev; // sum of absolute values of string deviations
- int iter; // number of iterations necessary for convergence
- int iter_Newton; // number of iterations necessary for convergence (Newton method)
- DP E; // total energy
- int iK; // K = 2.0*PI * iK/Nsites
- DP K; // total momentum
- DP lnnorm; // ln of norm of reduced Gaudin matrix
- std::string label;
-
- public:
- Heis_Bethe_State ();
- Heis_Bethe_State (const Heis_Bethe_State& RefState); // copy constructor
- Heis_Bethe_State (const Heis_Chain& RefChain, int M); // constructor to ground-state configuration
- Heis_Bethe_State (const Heis_Chain& RefChain, const Heis_Base& base); // constructor to lowest-energy config with base
- virtual ~Heis_Bethe_State () {};
-
- public:
- int Charge () { return(base.Mdown); };
- void Set_to_Label (std::string label_ref, const Vect<Vect<int> >& OriginIx2);
- void Set_Label_from_Ix2 (const Vect<Vect<int> >& OriginIx2);
- bool Check_Symmetry (); // checks whether the I's are symmetrically distributed
- void Compute_diffsq (); // \sum BE[j][alpha]^2
- void Find_Rapidities (bool reset_rapidities); // Finds the rapidities
- void Find_Rapidities_Twisted (bool reset_rapidities, DP twist); // Finds the rapidities with twist added to RHS of logBE
- void Solve_BAE_bisect (int j, int alpha, DP req_prec, int itermax);
- void Iterate_BAE (DP iter_factor); // Finds new set of lambda[j][alpha] from previous one by simple iteration
- void Solve_BAE_straight_iter (DP straight_prec, int max_iter_interp, DP iter_factor);
- void Solve_BAE_extrap (DP extrap_prec, int max_iter_extrap, DP iter_factor);
- void Iterate_BAE_Newton (); // Finds new set of lambda[j][alpha] from previous one by a Newton step
- void Solve_BAE_Newton (DP Newton_prec, int max_iter_Newton);
- void Solve_BAE_with_silk_gloves (DP silk_prec, int max_iter_silk, DP iter_factor);
- void Compute_lnnorm ();
- void Compute_Momentum ();
- void Compute_All (bool reset_rapidities); // solves BAE, computes E, K and lnnorm
-
- inline bool Set_to_Inner_Skeleton (int iKneeded, const Vect<Vect<int> >& OriginStateIx2)
- {
- Ix2[0][0] = Ix2[0][1] - 2;
- Ix2[0][base.Nrap[0] - 1] = Ix2[0][base.Nrap[0] - 2] + 2;
- (*this).Compute_Momentum();
- if (base.Nrap[0] == 0) return(false);
- if (iKneeded >= iK) Ix2[0][base.Nrap[0]-1] += 2*(iKneeded - iK);
- else Ix2[0][0] += 2*(iKneeded - iK);
- if (Ix2[0][0] < base.Ix2_min[0] || Ix2[0][base.Nrap[0]-1] > base.Ix2_max[0]) return(false);
- (*this).Set_Label_from_Ix2 (OriginStateIx2);
- return(true);
- }
- void Set_to_Outer_Skeleton (const Vect<Vect<int> >& OriginStateIx2) {
- Ix2[0][0] = base.Ix2_min[0] - 4;
- Ix2[0][base.Nrap[0]-1] = base.Ix2_max[0] + 4;
- (*this).Set_Label_from_Ix2 (OriginStateIx2);
- };
-
- void Set_to_Closest_Matching_Ix2_fixed_Base (const Heis_Bethe_State& StateToMatch); // defined in Heis.cc
-
-
- // Virtual functions, all defined in the derived classes
-
- public:
- virtual void Set_Free_lambdas() { ABACUSerror("Heis_Bethe_State::..."); } // sets the rapidities to solutions of BAEs without scattering terms
- virtual bool Check_Admissibility(char option) { ABACUSerror("Heis_Bethe_State::..."); return(false); }
- // verifies that we don't have a symmetrical Ix2 config with a Ix2 == 0 for a string of even length >= 2.
- virtual void Compute_BE (int j, int alpha) { ABACUSerror("Heis_Bethe_State::..."); }
- virtual void Compute_BE () { ABACUSerror("Heis_Bethe_State::..."); }
- virtual DP Iterate_BAE(int i, int alpha) { ABACUSerror("Heis_Bethe_State::..."); return(0.0);}
- virtual bool Check_Rapidities() { ABACUSerror("Heis_Bethe_State::..."); return(false); }
- virtual DP String_delta () { ABACUSerror("Heis_Bethe_State::..."); return(0.0); }
- virtual void Compute_Energy () { ABACUSerror("Heis_Bethe_State::..."); }
- virtual void Build_Reduced_Gaudin_Matrix (SQMat<std::complex<DP> >& Gaudin_Red) { ABACUSerror("Heis_Bethe_State::..."); }
- };
-
- inline bool Is_Inner_Skeleton (Heis_Bethe_State& State) {
- return (State.base.Nrap[0] >= 2 && (State.Ix2[0][0] == State.Ix2[0][1] - 2 || State.Ix2[0][State.base.Nrap[0]-1] == State.Ix2[0][State.base.Nrap[0]-2] + 2));
- };
- inline bool Is_Outer_Skeleton (Heis_Bethe_State& State) {
- return (State.Ix2[0][0] == State.base.Ix2_min[0] - 4 && State.Ix2[0][State.base.Nrap[0]-1] == State.base.Ix2_max[0] + 4);
- };
-
- inline bool Force_Descent (char whichDSF, Heis_Bethe_State& ScanState, Heis_Bethe_State& RefState, int desc_type_required, int iKmod, DP Chem_Pot)
- {
- bool force_descent = false;
-
- // Force descent for all DSFs if we're at K = 0 or PI and not conserving momentum upon descent:
- if (desc_type_required > 8 && (2*(ScanState.iK - RefState.iK) % iKmod == 0)) force_descent = true; // type_req > 8 means that we don't conserve momentum
-
- return(force_descent);
- }
-
-
- std::ostream& operator<< (std::ostream& s, const Heis_Bethe_State& state);
-
-
- //****************************************************************************
-
- // Objects in class XXZ_Bethe_State carry all extra information pertaining to XXZ gapless
-
- class XXZ_Bethe_State : public Heis_Bethe_State {
-
- public:
- Lambda sinhlambda;
- Lambda coshlambda;
- Lambda tanhlambda;
-
- public:
- XXZ_Bethe_State ();
- XXZ_Bethe_State (const XXZ_Bethe_State& RefState); // copy constructor
- XXZ_Bethe_State (const Heis_Chain& RefChain, int M); // constructor to ground-state configuration
- XXZ_Bethe_State (const Heis_Chain& RefChain, const Heis_Base& base); // constructor to lowest-energy config with base
-
- public:
- XXZ_Bethe_State& operator= (const XXZ_Bethe_State& RefState);
-
- public:
- void Set_Free_lambdas(); // sets the rapidities to solutions of BAEs without scattering terms
- void Compute_sinhlambda();
- void Compute_coshlambda();
- void Compute_tanhlambda();
- bool Check_Admissibility(char option); // verifies that we don't have a symmetrical Ix2 config with a Ix2 == 0 for a string of even length >= 2.
- void Compute_BE (int j, int alpha);
- void Compute_BE ();
- DP Iterate_BAE(int i, int j);
- bool Check_Rapidities(); // checks that all rapidities are not nan
- DP String_delta ();
- void Compute_Energy ();
- void Build_Reduced_Gaudin_Matrix (SQMat<std::complex<DP> >& Gaudin_Red);
-
- // XXZ specific functions:
- public:
-
- };
-
- XXZ_Bethe_State Add_Particle_at_Center (const XXZ_Bethe_State& RefState);
- XXZ_Bethe_State Remove_Particle_at_Center (const XXZ_Bethe_State& RefState);
-
- // Defined in XXZ_Bethe_State.cc
- inline DP fbar_XXZ (DP lambda, int par, DP tannzetaover2);
- DP Theta_XXZ (DP lambda, int nj, int nk, int parj, int park, DP* tannzetaover2);
- DP hbar_XXZ (DP lambda, int n, int par, DP* si_n_anis_over_2);
- DP ddlambda_Theta_XXZ (DP lambda, int nj, int nk, int parj, int park, DP* si_n_anis_over_2);
-
-
-
- //****************************************************************************
-
- // Objects in class XXX_Bethe_State carry all extra information pertaining to XXX antiferromagnet
-
- class XXX_Bethe_State : public Heis_Bethe_State {
-
- public:
- XXX_Bethe_State ();
- XXX_Bethe_State (const XXX_Bethe_State& RefState); // copy constructor
- XXX_Bethe_State (const Heis_Chain& RefChain, int M); // constructor to ground-state configuration
- XXX_Bethe_State (const Heis_Chain& RefChain, const Heis_Base& base); // constructor to lowest-energy config with base
-
- public:
- XXX_Bethe_State& operator= (const XXX_Bethe_State& RefState);
-
- public:
- void Set_Free_lambdas(); // sets the rapidities to solutions of BAEs without scattering terms
- bool Check_Admissibility(char option); // verifies that we don't have a symmetrical Ix2 config with a Ix2 == 0 for a string of even length >= 2.
- void Compute_BE (int j, int alpha);
- void Compute_BE ();
- DP Iterate_BAE(int i, int j);
- bool Check_Rapidities(); // checks that all rapidities are not nan
- DP String_delta ();
- void Compute_Energy ();
- void Build_Reduced_Gaudin_Matrix (SQMat<std::complex<DP> >& Gaudin_Red);
-
- // XXX specific functions
- public:
- bool Check_Finite_rap ();
-
- };
-
- XXX_Bethe_State Add_Particle_at_Center (const XXX_Bethe_State& RefState);
- XXX_Bethe_State Remove_Particle_at_Center (const XXX_Bethe_State& RefState);
-
- //****************************************************************************
-
- // Objects in class XXZ_gpd_Bethe_State carry all extra information pertaining to XXZ gapped antiferromagnets
-
- class XXZ_gpd_Bethe_State : public Heis_Bethe_State {
-
- public:
- Lambda sinlambda;
- Lambda coslambda;
- Lambda tanlambda;
-
- public:
- XXZ_gpd_Bethe_State ();
- XXZ_gpd_Bethe_State (const XXZ_gpd_Bethe_State& RefState); // copy constructor
- XXZ_gpd_Bethe_State (const Heis_Chain& RefChain, int M); // constructor to ground-state configuration
- XXZ_gpd_Bethe_State (const Heis_Chain& RefChain, const Heis_Base& base); // constructor to lowest-energy config with base
-
- public:
- XXZ_gpd_Bethe_State& operator= (const XXZ_gpd_Bethe_State& RefState);
-
- public:
- void Set_Free_lambdas(); // sets the rapidities to solutions of BAEs without scattering terms
- void Compute_sinlambda();
- void Compute_coslambda();
- void Compute_tanlambda();
- int Weight(); // weight function for contributions cutoff
- bool Check_Admissibility(char option); // verifies that we don't have a symmetrical Ix2 config with a Ix2 == 0 for a string of even length >= 2.
- void Compute_BE (int j, int alpha);
- void Compute_BE ();
- DP Iterate_BAE(int i, int j);
- void Iterate_BAE_Newton();
- bool Check_Rapidities(); // checks that all rapidities are not nan and are in interval ]-PI/2, PI/2]
- DP String_delta ();
- void Compute_Energy ();
- void Build_Reduced_Gaudin_Matrix (SQMat<std::complex<DP> >& Gaudin_Red);
-
- // XXZ_gpd specific functions
- public:
-
- };
-
- XXZ_gpd_Bethe_State Add_Particle_at_Center (const XXZ_gpd_Bethe_State& RefState);
- XXZ_gpd_Bethe_State Remove_Particle_at_Center (const XXZ_gpd_Bethe_State& RefState);
-
-
- //***********************************************
-
- // Function declarations
-
- // in M_vs_H.cc
- DP Ezero (DP Delta, int N, int M);
- DP H_vs_M (DP Delta, int N, int M);
- DP HZmin (DP Delta, int N, int M, Vect_DP& Ezero_ref);
- int M_vs_H (DP Delta, int N, DP HZ);
-
- DP X_avg (char xyorz, DP Delta, int N, int M);
-
- DP Chemical_Potential (const Heis_Bethe_State& RefState);
- DP Particle_Hole_Excitation_Cost (char whichDSF, Heis_Bethe_State& AveragingState);
-
- //DP Sumrule_Factor (char whichDSF, Heis_Bethe_State& RefState, DP Chem_Pot, bool fixed_iK, int iKneeded);
- DP Sumrule_Factor (char whichDSF, Heis_Bethe_State& RefState, DP Chem_Pot, int iKmin, int iKmax);
- void Evaluate_F_Sumrule (std::string prefix, char whichDSF, const Heis_Bethe_State& RefState, DP Chem_Pot, int iKmin, int iKmax);
-
- std::complex<DP> ln_Sz_ME (XXZ_Bethe_State& A, XXZ_Bethe_State& B);
- std::complex<DP> ln_Smin_ME (XXZ_Bethe_State& A, XXZ_Bethe_State& B);
-
- std::complex<DP> ln_Sz_ME (XXX_Bethe_State& A, XXX_Bethe_State& B);
- std::complex<DP> ln_Smin_ME (XXX_Bethe_State& A, XXX_Bethe_State& B);
-
- // From Antoine Klauser:
- std::complex<DP> ln_Szz_ME (XXX_Bethe_State& A, XXX_Bethe_State& B);
- std::complex<DP> ln_Szm_p_Smz_ME (XXX_Bethe_State& A, XXX_Bethe_State& B);
- std::complex<DP> ln_Smm_ME (XXX_Bethe_State& A, XXX_Bethe_State& B);
-
- std::complex<DP> ln_Sz_ME (XXZ_gpd_Bethe_State& A, XXZ_gpd_Bethe_State& B);
- std::complex<DP> ln_Smin_ME (XXZ_gpd_Bethe_State& A, XXZ_gpd_Bethe_State& B);
-
-
- DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_Bethe_State& LeftState,
- XXZ_Bethe_State& RefState, DP Chem_Pot, std::stringstream& DAT_outfile);
- DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXX_Bethe_State& LeftState,
- XXX_Bethe_State& RefState, DP Chem_Pot, std::stringstream& DAT_outfile);
- DP Compute_Matrix_Element_Contrib (char whichDSF, int iKmin, int iKmax, XXZ_gpd_Bethe_State& LeftState,
- XXZ_gpd_Bethe_State& RefState, DP Chem_Pot, std::stringstream& DAT_outfile);
-
- // For geometric quench:
- std::complex<DP> ln_Overlap (XXX_Bethe_State& A, XXX_Bethe_State& B);
-
- void Scan_Heis_Geometric_Quench (DP Delta, int N_1, int M, long long int base_id_1, long long int type_id_1, long long int id_1,
- int N_2, int iKmin, int iKmax, int Max_Secs, bool refine);
-
-
- } // namespace ABACUS
-
- #endif
|