HILA
Loading...
Searching...
No Matches
energy_and_topo_charge_clover.h
Go to the documentation of this file.
1/** @file energy_and_topo_charge_clover.h */
2
3#ifndef ENERGY_AND_TOPO_CHARGE_CLOVER_H_
4#define ENERGY_AND_TOPO_CHARGE_CLOVER_H_
5
6#include "hila.h"
7
8
9template <typename group, typename atype = hila::arithmetic_type<group>>
10void measure_topo_charge_and_energy_clover(const GaugeField<group> &U, atype &qtopo_out,
11 atype &energy_out) {
12 // measure topological charge and field strength energy of the gauge field, using the
13 // clover matrices as components of the field strength tensor
14
15 Reduction<double> qtopo = 0;
16 Reduction<double> energy = 0;
17 qtopo.allreduce(false).delayed(true);
18 energy.allreduce(false).delayed(true);
19
20#if NDIM == 4
21 Field<group> F[6];
22 // F[0]: F[0][1], F[1]: F[0][2], F[2]: F[0][3],
23 // F[3]: F[1][2], F[4]: F[1][3], F[5]: F[2][3]
24 /*
25 int k=0;
26 foralldir(dir1) foralldir(dir2) if(dir1<dir2) {
27 onsites(ALL) {
28 // clover operator as in eq. (2.9) of Nuclear Physics B259 (1985) 572-596
29 F[k][X]=0.25*(U[dir1][X]*U[dir2][X+dir1]*(U[dir2][X]*U[dir1][X+dir2]).dagger()
30 +U[dir2][X]*(U[dir2][X-dir1]*U[dir1][X-dir1+dir2]).dagger()*U[dir1][X-dir1]
31 +(U[dir2][X-dir1-dir2]*U[dir1][X-dir1]).dagger()*U[dir1][X-dir1-dir2]*U[dir2][X-dir2]
32 +U[dir2][X-dir2].dagger()*U[dir1][X-dir2]*U[dir2][X+dir1-dir2]*U[dir1][X].dagger());
33 // Lie-algebra projection:
34 //F[k][X]=0.5*(F[k][X]-F[k][X].dagger());
35 F[k][X]-=F[k][X].dagger();
36 F[k][X]*=0.5;
37 F[k][X]-=trace(F[k][X])/group::size();
38 }
39 ++k;
40 }
41 */
42
43 Field<group> tF0, tF1;
44
45 int k = 0;
46 foralldir(dir1) foralldir(dir2) if (dir1 < dir2) {
47 U[dir2].start_gather(dir1, ALL);
48 U[dir1].start_gather(dir2, ALL);
49
50 onsites(ALL) {
51 // dir1-dir2-plaquette that starts and ends at X; corresponds to F[dir1][dir2]
52 // at center location X+dir1/2+dir2/2 of plaquette:
53 tF0[X] = U[dir1][X] * U[dir2][X + dir1] * (U[dir2][X] * U[dir1][X + dir2]).dagger();
54 // project to Lie-algebra (anti-hermitian trace-free)
55 tF0[X] -= tF0[X].dagger();
56 tF0[X] *= 0.5;
57 tF0[X] -= trace(tF0[X]) / group::size();
58 // parallel transport to X+dir1
59 tF1[X] = U[dir1][X].dagger() * tF0[X] * U[dir1][X];
60 }
61
62 tF1.start_gather(-dir1, ALL);
63 onsites(ALL) {
64 tF0[X] += tF1[X - dir1];
65 }
66
67 U[dir2].start_gather(-dir2, ALL);
68 tF0.start_gather(-dir2, ALL);
69 onsites(ALL) {
70 // get F[dir1][dir2] at X from average of the (parallel transported) F[dir1][dir2] from
71 // the centers of all dir1-dir2-plaquettes that touch X :
72 F[k][X] =
73 (tF0[X] + U[dir2][X - dir2].dagger() * tF0[X - dir2] * U[dir2][X - dir2]) * 0.25;
74 }
75 ++k;
76 }
77
78 onsites(ALL) {
79 qtopo += real(mul_trace(F[0][X], F[5][X]));
80 qtopo += -real(mul_trace(F[1][X], F[4][X]));
81 qtopo += real(mul_trace(F[2][X], F[3][X]));
82
83 energy += F[0][X].squarenorm();
84 energy += F[1][X].squarenorm();
85 energy += F[2][X].squarenorm();
86 energy += F[3][X].squarenorm();
87 energy += F[4][X].squarenorm();
88 energy += F[5][X].squarenorm();
89 }
90#endif
91 qtopo_out = (atype)qtopo.value() / (4.0 * M_PI * M_PI);
92 energy_out = (atype)energy.value();
93}
94
95#endif
Array< n, m, hila::arithmetic_type< T > > real(const Array< n, m, T > &arg)
Return real part of Array.
Definition array.h:688
The field class implements the standard methods for accessing Fields. Hilapp replaces the parity acce...
Definition field.h:61
Field< A > dagger() const
Returns dagger or Hermitian conjugate of Field depending on how it is defined for Field type T.
Definition field.h:1123
double squarenorm() const
Squarenorm.
Definition field.h:1083
dir_mask_t start_gather(Direction d, Parity p=ALL) const
Definition field_comm.h:262
Gauge field class.
Definition gaugefield.h:22
Special reduction class: enables delayed and non-blocking reductions, which are not possible with the...
Definition reduction.h:14
const T value()
Return value of the reduction variable. Wait for the comms if needed.
Definition reduction.h:157
Reduction & allreduce(bool b=true)
allreduce(bool) turns allreduce on or off. By default on.
Definition reduction.h:130
Reduction & delayed(bool b=true)
deferred(bool) turns deferred on or off. By default turns on.
Definition reduction.h:148
Complex< T > dagger(const Complex< T > &val)
Return dagger of Complex number.
Definition cmplx.h:1223
#define foralldir(d)
Macro to loop over (all) Direction(s)
Definition coordinates.h:78
constexpr Parity ALL
bit pattern: 011
auto mul_trace(const Mtype1 &a, const Mtype2 &b)
Returns trace of product of two matrices.
Definition matrix.h:2309