diff --git a/ABACUS.org b/ABACUS.org index 2cb6bdd..033b2b4 100644 --- a/ABACUS.org +++ b/ABACUS.org @@ -32,6 +32,7 @@ Type your description here TIP: Search for the string BUGRISK in the codebase + * Priority :ABACUS:Dev:Priority: :PROPERTIES: :ARCHIVE: %s_archive::* Priority @@ -45,6 +46,16 @@ TIP: Search for the string BUGRISK in the codebase :CUSTOM_ID: ImplementationQueue :END: +** CONCEPT Unittests for integration functions + - State "CONCEPT" from "" [2018-02-20 Tue 07:01] +Select a number of standard functions with known definite integrals: +- polynomial +- rational +- exponential +- sinusoidal + +Write a unittest aiming to reproduce the exact result, and displaying the accuracy. + ** CONCEPT Complex integration - State "CONCEPT" from "" [2018-02-10 Sat 06:28] diff --git a/ABACUS_Usage_Example_Heis.cc b/ABACUS_Usage_Example_Heis.cc index ef33418..9487667 100644 --- a/ABACUS_Usage_Example_Heis.cc +++ b/ABACUS_Usage_Example_Heis.cc @@ -6,7 +6,7 @@ Copyright (c) J.-S. Caux ----------------------------------------------------------- -File: ABACUS+_Usage_Example_Heis.cc +File: ABACUS_Usage_Example_Heis.cc Purpose: illustrates basic use of ABACUS for spin chains. @@ -22,20 +22,21 @@ int main() { clock_t StartTime = clock(); - // Refer to include/JSC_Heis.h for all class definitions + // Refer to include/ABACUS_Heis.h for all class definitions // and to src/HEIS/ for the actual implementations. // Set basic system parameters: DP Delta = 1.0; // anisotropy int N = 64; // chain length - int M = 32; // number of downturned spins as compared to all spins up; must be <= N/2 (on or above equator) + int M = 32; // number of downturned spins; must be <= N/2 (on or above equator) // Define the chain: J, Delta, h, Nsites Heis_Chain chain(1.0, Delta, 0.0, N); // The Heis_Chain class so constructed contains information about all types of strings: cout << "Delta = " << Delta << endl << "Possible string lengths and parities: "; - for (int j = 0; j < chain.Nstrings; ++j) cout << "(" << chain.Str_L[j] << ", " << chain.par[j] << ")" << "\t"; + for (int j = 0; j < chain.Nstrings; ++j) + cout << "(" << chain.Str_L[j] << ", " << chain.par[j] << ")" << "\t"; cout << endl; @@ -64,22 +65,26 @@ int main() // Now define some excited state. // First method: - // Start with defining a base using a base constructor with a vector of numbers of rapidities of each type as argument: + // Start with defining a base using a base constructor with a vector of + // numbers of rapidities of each type as argument: // First define Nrapidities vector: Vect Nrapidities(0, chain.Nstrings); - // Assigning the numbers that follow, you have to ensure (yourself!) the \sum_j M_j n_j = M (where n_j is the length of type j). + // Assigning the numbers that follow, you have to ensure (yourself!) that + // the \sum_j M_j n_j = M (where n_j is the length of type j). Nrapidities[0] = M-2; // number of one-strings Nrapidities[1] = 1; // one two-string // Define the base: Heis_Base ebase(chain, Nrapidities); // Once the base is defined, the limiting quantum numbers are automatically computed: cout << "ebase defined, data is (Mdown, Nrap, Nraptot, Ix2_infty, Ix2_min, Ix2_max, baselabel):" - << ebase.Mdown << endl << ebase.Nrap << endl << ebase.Nraptot << endl << ebase.Ix2_infty << endl + << ebase.Mdown << endl << ebase.Nrap << endl << ebase.Nraptot << endl + << ebase.Ix2_infty << endl << ebase.Ix2_min << endl << ebase.Ix2_max << endl << ebase.baselabel << endl; // An excited state can then be defined using this new base: XXX_Bethe_State estate(chain, ebase); - // Individual quantum numbers can be manipulated: this will NOT update the state label or verify range of Ix2 + // Individual quantum numbers can be manipulated: this will NOT update the + // state label or verify range of Ix2 estate.Ix2[0][0] = M+1; estate.Compute_All(true); cout << endl << "estate: " << estate << endl; @@ -91,11 +96,12 @@ int main() // Construct a new Bethe state - XXX_Bethe_State estateref (chain, ebase2); // this will contain the lowest-energy quantum numbers for this base + XXX_Bethe_State estateref (chain, ebase2); // defaults to the lowest-energy state for this base XXX_Bethe_State estate2 (chain, ebase2); // yet another state // Setting a state to a given label: - estate2.Set_to_Label ("31_1_nh", estateref.Ix2); // The base of estate must coincide with the base in the label. Label is relative to estateref.Ix2 + // The base of estate must coincide with the base in label. Label is relative to estateref.Ix2. + estate2.Set_to_Label ("31_1_nh", estateref.Ix2); estate2.Compute_All(true); cout << "estate2: " << estate2 << endl; @@ -106,7 +112,8 @@ int main() cout << "Excitation momentum (mod 2 pi) = estate.K - gstate.K = " << estate.K - gstate.K << endl; - // Computing matrix elements: CAREFUL: magnetizations must be consistent with operator (error is flagged). + // Computing matrix elements: + // CAREFUL: magnetizations must be consistent with operator (error is flagged). cout << "Matrix elements: " << endl; // The logarithm of the matrix elements are computed as complex numbers: cout << "ln_Sz (gstate, estate) = " << ln_Sz_ME (gstate, estate) << endl; diff --git a/include/ABACUS_Heis.h b/include/ABACUS_Heis.h index 8a7d187..6348aa7 100644 --- a/include/ABACUS_Heis.h +++ b/include/ABACUS_Heis.h @@ -294,6 +294,13 @@ namespace ABACUS { 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); + + //**************************************************************************** diff --git a/src/HEIS/Heis.cc b/src/HEIS/Heis.cc index dbef830..a6882ab 100644 --- a/src/HEIS/Heis.cc +++ b/src/HEIS/Heis.cc @@ -422,17 +422,20 @@ namespace ABACUS { for (int k = 0; k < RefChain.Nstrings; ++k) { - sum2 = 0.0; + // sum2 = 0.0; - sum2 += (RefChain.Str_L[j] == RefChain.Str_L[k]) ? 0.0 : - 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k]) - - 0.5 * fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) * RefChain.anis)); - sum2 += 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k]) - - 0.5 * (RefChain.Str_L[j] + RefChain.Str_L[k]) * RefChain.anis)); + // sum2 += (RefChain.Str_L[j] == RefChain.Str_L[k]) ? 0.0 : + // 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k]) + // - 0.5 * fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) * RefChain.anis)); + // sum2 += 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k]) + // - 0.5 * (RefChain.Str_L[j] + RefChain.Str_L[k]) * RefChain.anis)); - for (int a = 1; a < ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]); ++a) - sum2 += 2.0 * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k]) - - 0.5 * (fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) + 2.0*a) * RefChain.anis)); + // for (int a = 1; a < ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]); ++a) + // sum2 += 2.0 * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k]) + // - 0.5 * (fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) + 2.0*a) * RefChain.anis)); + + sum2 = Theta_XXZ(1.0, RefChain.Str_L[j], RefChain.Str_L[k], RefChain.par[j], + RefChain.par[k], RefChain.ta_n_anis_over_2); sum1 += (Nrap[k] - ((j == k) ? 1 : 0)) * sum2; } @@ -459,6 +462,10 @@ namespace ABACUS { if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1; + // If Ix2_max equals Ix2_infty, we reduce it by 2: + + if (Ix2_max[j] == int(Ix2_infty[j])) Ix2_max[j] -= 2; + while (Ix2_max[j] > RefChain.Nsites) { Ix2_max[j] -= 2; } diff --git a/src/HEIS/XXZ_Bethe_State.cc b/src/HEIS/XXZ_Bethe_State.cc index 55526a9..aeedf65 100644 --- a/src/HEIS/XXZ_Bethe_State.cc +++ b/src/HEIS/XXZ_Bethe_State.cc @@ -20,10 +20,10 @@ namespace ABACUS { // Function prototypes - 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); + // 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); //***************************************************************************************************