/********************************************************** 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 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 Ix2_infty; // Ix2_infty[i] contains the max of BAE function for the (half-)integer I[i], i = 0, Nstrings - 1. Vect Ix2_min; Vect 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& 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 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 > 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 >& OriginIx2); void Set_Label_from_Ix2 (const Vect >& 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 >& 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 >& 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 >& 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 >& 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 >& 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 >& 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 ln_Sz_ME (XXZ_Bethe_State& A, XXZ_Bethe_State& B); std::complex ln_Smin_ME (XXZ_Bethe_State& A, XXZ_Bethe_State& B); std::complex ln_Sz_ME (XXX_Bethe_State& A, XXX_Bethe_State& B); std::complex ln_Smin_ME (XXX_Bethe_State& A, XXX_Bethe_State& B); // From Antoine Klauser: std::complex ln_Szz_ME (XXX_Bethe_State& A, XXX_Bethe_State& B); std::complex ln_Szm_p_Smz_ME (XXX_Bethe_State& A, XXX_Bethe_State& B); std::complex ln_Smm_ME (XXX_Bethe_State& A, XXX_Bethe_State& B); std::complex ln_Sz_ME (XXZ_gpd_Bethe_State& A, XXZ_gpd_Bethe_State& B); std::complex 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 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