4#include "datatypes/sun.h"
5#include "hmc/gauge_field.h"
9void exp_and_derivative(sun &Q, sun &m0, sun &lambda, sun &eQ,
int exp_steps) {
22 for (
int k = 2; k <= exp_steps; k++) {
23 n = n * 1.0 / ((double)k);
24 m1 = m0 * qn + Q * m1;
29 lambda = lambda + m1 * n;
35template <
typename matrix,
typename forcetype>
42 matrix
U1, U2, U3, U4, L, L2;
44 U2 = basegauge2[X + dir1];
45 U3 = basegauge1[X + dir2];
48 L2 = Lambda[X + dir2];
51 result2[X] += (L * U2 * U3.dagger()).
dagger();
52 stapleder2[X] = U3.
dagger() * U4.dagger() * L;
53 stapleder3[X] = (U4.dagger() * L * U2).
dagger();
56 stapleder2[X] = stapleder2[X] + L2.
dagger() * U4.dagger() *
U1;
57 result1[X] += U2 * L2.
dagger() * U4.dagger();
58 result2[X] += L2 * U2.
dagger() *
U1.dagger();
63 result2[X] += stapleder2[X - dir1];
64 result1[X] += stapleder3[X - dir2];
78 using gauge_type = sun;
79 using fund_type = sun;
80 using basetype = hila::arithmetic_type<sun>;
81 static constexpr int Nf = sun::size;
82 static constexpr int N = sun::size;
93 : base_field(f), c(coeff) {
98 : base_field(f), c(coeff), smear_steps(nsteps) {
103 : base_field(f), c(coeff), smear_steps(nsteps), exp_steps(expsteps) {
108 : base_field(r.base_field), c(r.c), smear_steps(r.smear_steps),
109 exp_steps(r.exp_steps) {
117 for (
int step = 0; step < smear_steps - 1; step++) {
121 staples[smear_steps - 1] =
new Field<sun>[NDIM];
122 smeared_fields[smear_steps - 1] = &(this->
gauge[0]);
126 for (
int step = 0; step < smear_steps - 1; step++) {
127 delete[] staples[step];
128 delete[] smeared_fields[step];
130 delete[] staples[smear_steps - 1];
132 free(smeared_fields);
138 previous = &base_field.
gauge[0];
140 for (
int step = 0; step < smear_steps; step++) {
143 staples[step][dir] = calc_staples(previous, dir);
146 Q = -c * previous[dir][X] * staples[step][dir][X];
147 project_antihermitean(Q);
148 Q = Q.exp(exp_steps);
149 smeared_fields[step][dir][X] = previous[dir][X] * Q;
152 previous = smeared_fields[step];
170 foralldir(dir) { storage1[dir] = force[dir]; }
177 for (
int step = smear_steps - 1; step >= 0; step--) {
181 basegauge = &base_field.
gauge[0];
183 basegauge = smeared_fields[step - 1];
188 result[dir][
ALL] = 0;
190 element<sun> m0, m1, qn, eQ, Q;
191 Q = -c * basegauge[dir][X] * staples[step][dir][X];
192 project_antihermitean(Q);
194 m0 = previous[dir][X] * basegauge[dir][X];
195 exp_and_derivative(Q, m0, Lambda[dir][X], eQ, exp_steps);
197 project_antihermitean(Lambda[dir][X]);
200 result[dir][X] = eQ * previous[dir][X];
203 result[dir][X] -= c * staples[step][dir][X] * Lambda[dir][X];
206 Lambda[dir][X] = -c * Lambda[dir][X] * basegauge[dir][X];
212 staple_dir_derivative(basegauge[dir1], basegauge[dir2], Lambda[dir1],
213 result[dir1], result[dir2], dir1, dir2);
220 basegauge = &smeared_fields[step][0];
241template <
typename sun>
struct HEX_smeared_field :
public gauge_field_base<sun> {
242 using gauge_type = sun;
243 using fund_type = sun;
244 using basetype = hila::arithmetic_type<sun>;
245 static constexpr int Nf = sun::size;
246 static constexpr int N = sun::size;
265 : base_field(f), exp_steps(nsteps) {
269 : base_field(f), c1(_c1), c2(_c2), c3(_c3) {
274 : base_field(f), c1(_c1), c2(_c2), c3(_c3), exp_steps(nsteps) {
281 previous = &base_field.
gauge[0];
288 staples3[nu][mu] = calc_staples(base_field.
gauge, base_field.
gauge, mu, nu);
291 Q = -c3 * base_field.
gauge[mu][X] * staples3[nu][mu][X];
292 project_antihermitean(Q);
293 Q = Q.exp(exp_steps);
294 level3[nu][mu][X] = base_field.
gauge[mu][X] * Q;
301 staples2[nu][mu][
ALL] = 0;
302 foralldir(rho)
if (rho != mu)
if (rho != nu) {
304 foralldir(e)
if (e != mu && e != nu && e != rho) eta = e;
305 Field<sun> stp = calc_staples(level3[eta], level3[eta], mu, rho);
306 staples2[nu][mu][
ALL] = staples2[nu][mu][X] + stp[X];
310 Q = -c2 * base_field.
gauge[mu][X] * staples2[nu][mu][X];
311 project_antihermitean(Q);
312 Q = Q.exp(exp_steps);
313 level2[nu][mu][X] = base_field.
gauge[mu][X] * Q;
320 staples1[mu][
ALL] = 0;
322 Field<sun> stp = calc_staples(level2[nu], level2[mu], mu, nu);
323 staples1[mu][
ALL] = staples1[mu][X] + stp[X];
327 Q = -c1 * base_field.
gauge[mu][X] * staples1[mu][X];
328 project_antihermitean(Q);
329 Q = Q.exp(exp_steps);
330 this->
gauge[mu][X] = base_field.
gauge[mu][X] * Q;
364 element<sun> m0, m1, qn, eQ, Q;
365 Q = -c1 * base_field.
gauge[mu][X] * staples1[mu][X];
366 project_antihermitean(Q);
368 m0 = force[mu][X] * base_field.
gauge[mu][X];
369 exp_and_derivative(Q, m0, lambda1[mu][X], eQ, exp_steps);
370 project_antihermitean(lambda1[mu][X]);
373 result[mu][X] = eQ * force[mu][X];
376 result[mu][X] -= c1 * staples1[mu][X] * lambda1[mu][X];
379 lambda1[mu][X] = -c1 * lambda1[mu][X] * base_field.
gauge[mu][X];
385 staple_dir_derivative(level2[nu][mu], level2[mu][nu], lambda1[mu],
386 result1[nu][mu], result1[mu][nu], mu, nu);
392 element<sun> m0, m1, qn, eQ, Q;
393 Q = -c2 * base_field.
gauge[mu][X] * staples2[nu][mu][X];
394 project_antihermitean(Q);
396 m0 = result1[nu][mu][X] * base_field.
gauge[mu][X];
397 exp_and_derivative(Q, m0, lambda2[nu][mu][X], eQ, exp_steps);
398 project_antihermitean(lambda2[nu][mu][X]);
401 result[mu][X] += eQ * result1[nu][mu][X];
404 result[mu][X] -= c2 * staples2[nu][mu][X] * lambda2[nu][mu][X];
407 lambda2[nu][mu][X] = -c2 * lambda2[nu][mu][X] * base_field.
gauge[mu][X];
413 foralldir(rho)
if (rho != mu)
if (rho != nu) {
415 foralldir(e)
if (e != mu && e != nu && e != rho) eta = e;
416 staple_dir_derivative(level3[eta][mu], level3[eta][rho], lambda2[nu][mu],
417 result2[eta][mu], result2[eta][rho], mu, rho);
424 element<sun> m0, m1, qn, eQ, Q;
425 Q = -c3 * base_field.
gauge[mu][X] * staples3[nu][mu][X];
426 project_antihermitean(Q);
428 m0 = result2[nu][mu][X] * base_field.
gauge[mu][X];
429 exp_and_derivative(Q, m0, lambda3[nu][mu][X], eQ, exp_steps);
430 project_antihermitean(lambda3[nu][mu][X]);
433 result[mu][X] += eQ * result2[nu][mu][X];
436 result[mu][X] -= c3 * staples3[nu][mu][X] * lambda3[nu][mu][X];
439 lambda3[nu][mu][X] = -c3 * lambda3[nu][mu][X] * base_field.
gauge[mu][X];
445 staple_dir_derivative(base_field.
gauge[mu], base_field.
gauge[nu],
446 lambda3[nu][mu], result[mu], result[nu], mu, nu);
The field class implements the standard methods for accessing Fields. Hilapp replaces the parity acce...
Field< A > dagger() const
Returns dagger or Hermitian conjugate of Field depending on how it is defined for Field type T.
void check_alloc()
Allocate Field if it is not already allocated.
Matrix class which defines matrix operations.
virtual void restore_backup()
virtual void refresh()
Recalculate represented or smeared field.
virtual void draw_momentum()
Draw gaussian random momentum.
virtual void random()
Draw a random gauge field.
Field< sun > gauge[4]
A matrix field for each Direction.
virtual void set_unity()
Set the field to unity.
virtual void backup()
Make a copy of the gauge field for HMC.
virtual void add_momentum(Field< SquareMatrix< N, Complex< basetype > > > *force)
Update the momentum by given force.
virtual void zero_momentum()
Set the momentum to zero.
void restore_backup()
Restore the previous backup.
void zero_momentum()
Set the momentum to zero.
void set_unity()
Set the gauge field to unity.
void random()
Draw a random gauge field.
void backup()
Make a copy of fields updated in a trajectory.
Field< gauge_type > & get_gauge(int dir)
Return a reference to the gauge field.
void draw_momentum()
Gaussian random momentum for each element.
void add_momentum(Field< SquareMatrix< N, Complex< basetype > > > *force)
Field< gauge_type > & get_momentum(int dir)
Return a reference to the momentum field.
void random()
Draw a random gauge field.
void draw_momentum()
Draw gaussian random momentum.
void set_unity()
Set the field to unity.
void zero_momentum()
Set the momentum to zero.
void add_momentum(Field< SquareMatrix< N, Complex< basetype > > > *force)
Update the momentum by given force.
void backup()
Make a copy of the gauge field for HMC.
void refresh()
Recalculate represented or smeared field.
Complex< T > dagger(const Complex< T > &val)
Return dagger of Complex number.
#define foralldir(d)
Macro to loop over (all) Direction(s)
Direction
Enumerator for direction that assigns integer to direction to be interpreted as unit vector.
constexpr Parity ALL
bit pattern: 011