HILA
Loading...
Searching...
No Matches
sun_heatbath.h
Go to the documentation of this file.
1#ifndef HILA_SUN_HEATBATH_H
2#define HILA_SUN_HEATBATH_H
3
4/** @file sun_heatbath.h */
5
6/* Kennedy-Pendleton quasi heat bath on SU(2) subgroups */
7
8/* This does the heatbath with
9 * beta ( 1 - 1/Ncol Tr U S^+ ) S is the sum of staples.
10 * -> -beta/Ncol Tr u (US+) = -beta/Ncol u_ab (US+)_ba
11 * here u is a SU(2) embedded in SU(N); u = su2 * 1.
12 */
13
14#include "sun_matrix.h"
15#include "su2.h"
16
17/**
18 * @brief \f$ SU(N) \f$ heatbath
19 * @details Kennedy-Pendleton quasi heat bath on \f$ SU(2)\f$ subgroups
20 * @tparam T Group element type such as Real or Complex
21 * @tparam N Number of colors
22 * @param U \f$ SU(N) \f$ link to perform heatbath on
23 * @param staple Staple to compute heatbath with
24 * @param beta
25 * @return double
26 */
27template <typename T, int N>
28double suN_heatbath(SU<N, T> &U, const SU<N, T> &staple, double beta) {
29 // K-P quasi-heat bath by SU(2) subgroups
30
31 double utry, uhit;
32 double xr1, xr2, xr3, xr4;
33
34 double r, r2, rho, z;
35 double al, d, xl, xd;
36 int k, test;
37 double b3;
38 SU<N, T> action;
39 SU<2, T> h2x2;
40 SU2<T> v, a, h;
41 double pi2 = M_PI * 2.0;
42
43 b3 = beta / N;
44
45 utry = uhit = 0.0;
46
47 /* now for the qhb updating */
48
49 action = U * staple.dagger();
50
51 /* SU(2) subgroups */
52 for (int ina = 0; ina < N - 1; ina++)
53 for (int inb = ina + 1; inb < N; inb++) {
54
55 /* decompose the action into SU(2) subgroups using
56 * Pauli matrix expansion
57 * The SU(2) hit matrix is represented as
58 * a0 + i * Sum j (sigma j * aj)
59 * Let us actually store dagger, v = action(a,b).dagger
60 */
61 v.d = action.e(ina, ina).re + action.e(inb, inb).re;
62 v.c = -(action.e(ina, ina).im - action.e(inb, inb).im);
63 v.a = -(action.e(ina, inb).im + action.e(inb, ina).im);
64 v.b = -(action.e(ina, inb).re - action.e(inb, ina).re);
65
66 z = sqrt(v.det());
67
68 v /= z;
69
70 /* end norm check--trial SU(2) matrix is a0 + i a(j)sigma(j)*/
71 /* test */
72 /* now begin qhb */
73 /* get four random numbers */
74
75 xr1 = log(1.0 - hila::random());
76 xr2 = log(1.0 - hila::random());
77 xr3 = hila::random();
78 xr4 = hila::random();
79
80 xr3 = cos(pi2 * xr3);
81
82
83 // generate a0 component of suN matrix
84 //
85 // first consider generating an su(2) matrix h
86 // according to exp(bg/Ncol * re tr(h*s))
87 // rewrite re tr(h*s) as re tr(h*v)z where v is
88 // an su(2) matrix and z is a real normalization constant
89 // let v = z*v. (z is 2*xi in k-p notation)
90 // v is represented in the form v(0) + i*sig*v (sig are pauli)
91 // v(0) and vector v are real
92 //
93 // let a = h*v and now generate a
94 // rewrite beta/3 * re tr(h*v) * z as al*a0
95 // a0 has prob(a0) = n0 * sqrt(1 - a0**2) * exp(al * a0)
96
97
98 al = b3 * z;
99
100 // let a0 = 1 - del**2
101 // get d = del**2
102 // such that prob2(del) = n1 * del**2 * exp(-al*del**2)
103
104
105 d = -(xr2 + xr1 * xr3 * xr3) / al;
106
107 /* monte carlo prob1(del) = n2 * sqrt(1 - 0.5*del**2)
108 * then prob(a0) = n3 * prob1(a0)*prob2(a0)
109 */
110
111 if ((1.00 - 0.5 * d) <= xr4 * xr4) {
112 if (al > 2.0) { /* k-p algorithm */
113 test = 0;
114 for (k = 0; k < 40 && !test; k++) {
115 /* get four random numbers */
116 xr1 = log(1.0 - hila::random());
117 xr2 = log(1.0 - hila::random());
118 xr3 = hila::random();
119 xr4 = hila::random();
120
121 xr3 = cos(pi2 * xr3);
122
123 d = -(xr2 + xr1 * xr3 * xr3) / al;
124 if ((1.00 - 0.5 * d) > xr4 * xr4)
125 test = 1;
126 }
127 utry += k;
128 uhit++;
129#ifndef __CUDA_ARCH__
130 // TODO: Make parallel-warning out of this
131 // hila::out0 << "site took 20 kp hits\n";
132#endif
133 } else { /* now al <= 2.0 */
134
135 /* creutz algorithm */
136 xl = exp((double)(-2.0 * al));
137 xd = 1.0 - xl;
138 test = 0;
139 double a0;
140 for (k = 0; k < 40 && test == 0; k++) {
141 /* get two random numbers */
142 xr1 = hila::random();
143 xr2 = hila::random();
144
145 r = xl + xd * xr1;
146 a0 = 1.00 + log((double)r) / al;
147 if ((1.0 - a0 * a0) > xr2 * xr2)
148 test = 1;
149 }
150 d = 1.0 - a0;
151 utry += k;
152 uhit++;
153#ifndef __CUDA_ARCH__
154 // hila::out0 << "site took 20 creutz hits\n";
155#endif
156 } /* endif al */
157
158 } else {
159 /* direct hit */
160 utry += 1.0;
161 uhit += 1.0;
162 }
163
164 /* generate full su(2) matrix and update link matrix*/
165
166 /* find a0 = 1 - d*/
167 a.d = 1.0 - d;
168 /* compute r */
169 r2 = 1.0 - a.d * a.d;
170 r2 = fabs(r2);
171 r = sqrt(r2);
172
173 /* compute a3 */
174 a.c = (2.0 * hila::random() - 1.0) * r;
175
176 /* compute a1 and a2 */
177 rho = r2 - a.c * a.c;
178 rho = sqrt(fabs(rho));
179
180 /*xr2 is a random number between 0 and 2*pi */
181 xr2 = pi2 * hila::random();
182 a.a = rho * cos((double)xr2);
183 a.b = rho * sin((double)xr2);
184
185 /* now do the updating. h = a*v^dagger, new u = h*u */
186 h = a * v; // v was daggerized above
187
188 /* Elements of SU(2) matrix */
189 h2x2 = h.convert_to_2x2_matrix();
190
191 /* update the link and 'action' */
192
193 U.mult_by_2x2_left(ina, inb, h2x2);
194 action.mult_by_2x2_left(ina, inb, h2x2);
195
196 } /* hits */
197
198 /* check_unitarity( U );
199 */
200 return (uhit / utry);
201
202} /* site */
203
204#endif
void mult_by_2x2_left(int p, int q, const Matrix_t< 2, 2, R, Mt > &M)
Multiply -matrix from the left by matrix defined by sub matrix.
Definition matrix.h:1595
Rtype dagger() const
Hermitian conjugate of matrix.
Definition matrix.h:1094
T e(const int i, const int j) const
Standard array indexing operation for matrices.
Definition matrix.h:279
Definition su2.h:14
Class for SU(N) matrix.
Definition sun_matrix.h:110
double random()
Real valued uniform random number generator.
Definition hila_gpu.cpp:121
double suN_heatbath(SU< N, T > &U, const SU< N, T > &staple, double beta)
heatbath
SU(N) Matrix definitions.