Browse Source

Restart work on code. Fix changes of last period.

master
J.-S. Caux 6 years ago
parent
commit
f035415b1b
5 changed files with 56 additions and 24 deletions
  1. 11
    0
      ABACUS.org
  2. 18
    11
      ABACUS_Usage_Example_Heis.cc
  3. 7
    0
      include/ABACUS_Heis.h
  4. 16
    9
      src/HEIS/Heis.cc
  5. 4
    4
      src/HEIS/XXZ_Bethe_State.cc

+ 11
- 0
ABACUS.org View File

@@ -32,6 +32,7 @@ Type your description here
32 32
 
33 33
 TIP: Search for the string BUGRISK in the codebase
34 34
 
35
+
35 36
 * Priority 					       :ABACUS:Dev:Priority:
36 37
   :PROPERTIES:
37 38
   :ARCHIVE:  %s_archive::* Priority
@@ -45,6 +46,16 @@ TIP: Search for the string BUGRISK in the codebase
45 46
   :CUSTOM_ID: ImplementationQueue
46 47
   :END:
47 48
 
49
+** CONCEPT Unittests for integration functions
50
+   - State "CONCEPT"    from ""           [2018-02-20 Tue 07:01]
51
+Select a number of standard functions with known definite integrals:
52
+- polynomial
53
+- rational
54
+- exponential
55
+- sinusoidal
56
+
57
+Write a unittest aiming to reproduce the exact result, and displaying the accuracy.
58
+
48 59
 
49 60
 ** CONCEPT Complex integration
50 61
    - State "CONCEPT"    from ""           [2018-02-10 Sat 06:28]

+ 18
- 11
ABACUS_Usage_Example_Heis.cc View File

@@ -6,7 +6,7 @@ Copyright (c) J.-S. Caux
6 6
 
7 7
 -----------------------------------------------------------
8 8
 
9
-File:  ABACUS+_Usage_Example_Heis.cc
9
+File:  ABACUS_Usage_Example_Heis.cc
10 10
 
11 11
 Purpose:  illustrates basic use of ABACUS for spin chains.
12 12
 
@@ -22,20 +22,21 @@ int main()
22 22
 {
23 23
   clock_t StartTime = clock();
24 24
 
25
-  // Refer to include/JSC_Heis.h for all class definitions
25
+  // Refer to include/ABACUS_Heis.h for all class definitions
26 26
   // and to src/HEIS/ for the actual implementations.
27 27
 
28 28
   // Set basic system parameters:
29 29
   DP Delta = 1.0; // anisotropy
30 30
   int N = 64; // chain length
31
-  int M = 32; // number of downturned spins as compared to all spins up; must be <= N/2 (on or above equator)
31
+  int M = 32; // number of downturned spins; must be <= N/2 (on or above equator)
32 32
 
33 33
   // Define the chain:  J, Delta, h, Nsites
34 34
   Heis_Chain chain(1.0, Delta, 0.0, N);
35 35
 
36 36
   // The Heis_Chain class so constructed contains information about all types of strings:
37 37
   cout << "Delta = " << Delta << endl << "Possible string lengths and parities: ";
38
-  for (int j = 0; j < chain.Nstrings; ++j) cout << "(" << chain.Str_L[j] << ", " << chain.par[j] << ")" << "\t";
38
+  for (int j = 0; j < chain.Nstrings; ++j)
39
+    cout << "(" << chain.Str_L[j] << ", " << chain.par[j] << ")" << "\t";
39 40
   cout << endl;
40 41
 
41 42
 
@@ -64,22 +65,26 @@ int main()
64 65
   // Now define some excited state.
65 66
 
66 67
   // First method:
67
-  // Start with defining a base using a base constructor with a vector of numbers of rapidities of each type as argument:
68
+  // Start with defining a base using a base constructor with a vector of
69
+  //  numbers of rapidities of each type as argument:
68 70
   // First define Nrapidities vector:
69 71
   Vect<int> Nrapidities(0, chain.Nstrings);
70
-  // 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).
72
+  // Assigning the numbers that follow, you have to ensure (yourself!) that
73
+  //  the \sum_j M_j n_j = M (where n_j is the length of type j).
71 74
   Nrapidities[0] = M-2; // number of one-strings
72 75
   Nrapidities[1] = 1; // one two-string
73 76
   // Define the base:
74 77
   Heis_Base ebase(chain, Nrapidities);
75 78
   // Once the base is defined, the limiting quantum numbers are automatically computed:
76 79
   cout << "ebase defined, data is (Mdown, Nrap, Nraptot, Ix2_infty, Ix2_min, Ix2_max, baselabel):"
77
-       << ebase.Mdown << endl << ebase.Nrap << endl << ebase.Nraptot << endl << ebase.Ix2_infty << endl
80
+       << ebase.Mdown << endl << ebase.Nrap << endl << ebase.Nraptot << endl
81
+       << ebase.Ix2_infty << endl
78 82
        << ebase.Ix2_min << endl << ebase.Ix2_max << endl << ebase.baselabel << endl;
79 83
 
80 84
   // An excited state can then be defined using this new base:
81 85
   XXX_Bethe_State estate(chain, ebase);
82
-  // Individual quantum numbers can be manipulated: this will NOT update the state label or verify range of Ix2
86
+  // Individual quantum numbers can be manipulated: this will NOT update the
87
+  //  state label or verify range of Ix2
83 88
   estate.Ix2[0][0] = M+1;
84 89
   estate.Compute_All(true);
85 90
   cout << endl << "estate: " << estate << endl;
@@ -91,11 +96,12 @@ int main()
91 96
 
92 97
 
93 98
   // Construct a new Bethe state
94
-  XXX_Bethe_State estateref (chain, ebase2); // this will contain the lowest-energy quantum numbers for this base
99
+  XXX_Bethe_State estateref (chain, ebase2); // defaults to the lowest-energy state for this base
95 100
   XXX_Bethe_State estate2 (chain, ebase2); // yet another state
96 101
 
97 102
   // Setting a state to a given label:
98
-  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
103
+  // The base of estate must coincide with the base in label. Label is relative to estateref.Ix2.
104
+  estate2.Set_to_Label ("31_1_nh", estateref.Ix2);
99 105
   estate2.Compute_All(true);
100 106
   cout << "estate2: " << estate2 << endl;
101 107
 
@@ -106,7 +112,8 @@ int main()
106 112
   cout << "Excitation momentum (mod 2 pi) = estate.K - gstate.K = " << estate.K - gstate.K << endl;
107 113
 
108 114
 
109
-  // Computing matrix elements: CAREFUL: magnetizations must be consistent with operator (error is flagged).
115
+  // Computing matrix elements:
116
+  // CAREFUL: magnetizations must be consistent with operator (error is flagged).
110 117
   cout << "Matrix elements: " << endl;
111 118
   // The logarithm of the matrix elements are computed as complex numbers:
112 119
   cout << "ln_Sz (gstate, estate) = " << ln_Sz_ME (gstate, estate) << endl;

+ 7
- 0
include/ABACUS_Heis.h View File

@@ -294,6 +294,13 @@ namespace ABACUS {
294 294
   XXZ_Bethe_State Add_Particle_at_Center (const XXZ_Bethe_State& RefState);
295 295
   XXZ_Bethe_State Remove_Particle_at_Center (const XXZ_Bethe_State& RefState);
296 296
 
297
+  // Defined in XXZ_Bethe_State.cc
298
+  inline DP fbar_XXZ (DP lambda, int par, DP tannzetaover2);
299
+  DP Theta_XXZ (DP lambda, int nj, int nk, int parj, int park, DP* tannzetaover2);
300
+  DP hbar_XXZ (DP lambda, int n, int par, DP* si_n_anis_over_2);
301
+  DP ddlambda_Theta_XXZ (DP lambda, int nj, int nk, int parj, int park, DP* si_n_anis_over_2);
302
+
303
+
297 304
 
298 305
   //****************************************************************************
299 306
 

+ 16
- 9
src/HEIS/Heis.cc View File

@@ -422,17 +422,20 @@ namespace ABACUS {
422 422
 
423 423
 	for (int k = 0; k < RefChain.Nstrings; ++k) {
424 424
 
425
-	  sum2 = 0.0;
425
+	  // sum2 = 0.0;
426 426
 
427
-	  sum2 += (RefChain.Str_L[j] == RefChain.Str_L[k]) ? 0.0 :
428
-	    2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
429
-			   - 0.5 * fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) * RefChain.anis));
430
-	  sum2 += 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
431
-				 - 0.5 * (RefChain.Str_L[j] + RefChain.Str_L[k]) * RefChain.anis));
427
+	  // sum2 += (RefChain.Str_L[j] == RefChain.Str_L[k]) ? 0.0 :
428
+	  //   2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
429
+	  // 		   - 0.5 * fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) * RefChain.anis));
430
+	  // sum2 += 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
431
+	  // 			 - 0.5 * (RefChain.Str_L[j] + RefChain.Str_L[k]) * RefChain.anis));
432 432
 
433
-	  for (int a = 1; a < ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]); ++a)
434
-	    sum2 += 2.0 * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
435
-					 - 0.5 * (fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) + 2.0*a) * RefChain.anis));
433
+	  // for (int a = 1; a < ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]); ++a)
434
+	  //   sum2 += 2.0 * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
435
+	  // 				 - 0.5 * (fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) + 2.0*a) * RefChain.anis));
436
+
437
+	  sum2 = Theta_XXZ(1.0, RefChain.Str_L[j], RefChain.Str_L[k], RefChain.par[j],
438
+			   RefChain.par[k], RefChain.ta_n_anis_over_2);
436 439
 
437 440
 	  sum1 += (Nrap[k] - ((j == k) ? 1 : 0)) * sum2;
438 441
 	}
@@ -459,6 +462,10 @@ namespace ABACUS {
459 462
 
460 463
 	if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;
461 464
 
465
+	// If Ix2_max equals Ix2_infty, we reduce it by 2:
466
+
467
+	if (Ix2_max[j] == int(Ix2_infty[j])) Ix2_max[j] -= 2;
468
+
462 469
 	while (Ix2_max[j] > RefChain.Nsites) {
463 470
 	  Ix2_max[j] -= 2;
464 471
 	}

+ 4
- 4
src/HEIS/XXZ_Bethe_State.cc View File

@@ -20,10 +20,10 @@ namespace ABACUS {
20 20
 
21 21
   // Function prototypes
22 22
 
23
-  inline DP fbar_XXZ (DP lambda, int par, DP tannzetaover2);
24
-  DP Theta_XXZ (DP lambda, int nj, int nk, int parj, int park, DP* tannzetaover2);
25
-  DP hbar_XXZ (DP lambda, int n, int par, DP* si_n_anis_over_2);
26
-  DP ddlambda_Theta_XXZ (DP lambda, int nj, int nk, int parj, int park, DP* si_n_anis_over_2);
23
+  // inline DP fbar_XXZ (DP lambda, int par, DP tannzetaover2);
24
+  // DP Theta_XXZ (DP lambda, int nj, int nk, int parj, int park, DP* tannzetaover2);
25
+  // DP hbar_XXZ (DP lambda, int n, int par, DP* si_n_anis_over_2);
26
+  // DP ddlambda_Theta_XXZ (DP lambda, int nj, int nk, int parj, int park, DP* si_n_anis_over_2);
27 27
 
28 28
 
29 29
   //***************************************************************************************************

Loading…
Cancel
Save