HILA
Loading...
Searching...
No Matches
clover_action.h
Go to the documentation of this file.
1/** @file clover_action.h */
2
3#ifndef CLOVER_ACTION_H_
4#define CLOVER_ACTION_H_
5
6#include "hila.h"
8
9// functions for clover gauge action
10
11template <typename group, typename atype = hila::arithmetic_type<group>>
12void get_clover_leaves(const GaugeField<group> &U, Direction d1, Direction d2,
13 Field<group>(out_only &C)[4], out_only Field<group> &Csum) {
14 // assignes the clover leaf matrices in counter-clockwise order to C[0],C[1],C[2],C[3] and sets
15 // Csum to the full Clover matrix, i.e. to the anti-hermitian tracless part of
16 // (C[0]+C[1]+C[2]+C[3])/4
17
18 U[d2].start_gather(d1, ALL);
19 U[d1].start_gather(d2, ALL);
20 Field<group> tF;
21
22 onsites(ALL) {
23 // d1-d2-plaquette that starts and ends at X
24 tF[X] = U[d2][X + d1] * (U[d2][X] * U[d1][X + d2]).dagger();
25 tF[X] *= 0.25;
26 C[0][X] = U[d1][X] * tF[X];
27
28 // parallel transport C[0][X] to X+d1 (use C[2][X] as temp. storage)
29 C[2][X] = tF[X] * U[d1][X];
30 }
31
32 C[2].start_gather(-d1, ALL);
33
34 onsites(ALL) {
35 // parallel transport C[0][X] to X+d2 (use Csum as temp. storage)
36 Csum[X] = U[d2][X].dagger() * C[0][X] * U[d2][X];
37 }
38
39 Csum.start_gather(-d2, ALL);
40
41 onsites(ALL) {
42 C[1][X] = C[2][X - d1];
43
44 // parallel transport C[1][X] to X+d2 (use tF as temp. storage)
45 tF[X] = U[d2][X].dagger() * C[1][X] * U[d2][X];
46 }
47
48 tF.start_gather(-d2, ALL);
49
50 onsites(ALL) {
51 C[3][X] = Csum[X - d2];
52 }
53
54 onsites(ALL) {
55 C[2][X] = tF[X - d2];
56
57 Csum[X] = C[0][X];
58 Csum[X] += C[1][X];
59 Csum[X] += C[2][X];
60 Csum[X] += C[3][X];
61
62 // Lie-algebra project Csum[X]
63 Csum[X] -= Csum[X].dagger();
64 Csum[X] *= 0.5;
65 Csum[X] -= trace(Csum[X]) / group::size();
66 }
67}
68
69template <typename group>
70double measure_s_clover(const GaugeField<group> &U) {
71 // measure the clover action for dir1<dir2
72 // (just to have same normalization as with plaquette action)
73 Reduction<double> stot = 0;
74 stot.allreduce(false).delayed(true);
75 Field<group> C[4];
76 Field<group> Csum;
77 foralldir(dir1) foralldir(dir2) if (dir1 < dir2) {
78 get_clover_leaves(U, dir1, dir2, C, Csum);
79 onsites(ALL) {
80 stot += 0.5 * Csum[X].squarenorm();
81 }
82 }
83 return stot.value() / group::size();
84}
85
86template <typename group, typename atype = hila::arithmetic_type<group>>
87void get_force_clover_add(const GaugeField<group> &U, VectorField<Algebra<group>> &K,
88 atype eps = 1.0) {
89 // compute gauge force for clover action and add result to K
90 Field<group> C[4];
91 Field<group> Csum;
92
93 Field<group> tC;
94 foralldir(dir1) foralldir(dir2) if (dir1 < dir2) {
95 // get the 4 clover leaf matrices C[][X] (ordered counter-clockwise in (dir1,dir2)-plane)
96 // and their anti-hermitian traceless sum Csum[X]
97 get_clover_leaves(U, dir1, dir2, C, Csum);
98
99 // define for each clover leaf matrix the correspondin path of gauge links
100 std::vector<Direction> paths[4] = {{dir1, dir2, -dir1, -dir2},
101 {dir2, -dir1, -dir2, dir1},
102 {-dir1, -dir2, dir1, dir2},
103 {-dir2, dir1, dir2, -dir1}};
104
105 // the force on the link (Y,dir) from the clover term with foot point X is:
106 // F[dir][Y] = d_{Y,dir}(0.5*Tr(Csum*Csum))
107 // = 0.5*Tr((d_{Y,dir}Csum)*Csum) + 0.5*Tr(Csum*(d_{Y,dir}Csum))
108 // = Tr((d_{Y,dir}Csum)*Csum)
109 // = 0.5*\sum_k (Tr((d_{Y,dir}C[k][X])*Csum[X]) -
110 // Tr((d_{Y,dir}C[k][X].dagger())*Csum[X])) = \sum_k
111 // Tr((d_{Y,dir}C[k][X])*Csum[X])
112 for (int k = 0; k < 4; ++k) {
113 // multiply the clover leave matrix C[k][X] by Csum[X] --> yields (weighted) sum
114 // of Wilson loops with foot point X and whose first 4 links are always those of C[k][X]
115 onsites(ALL) {
116 tC[X] = C[k][X] * Csum[X];
117 }
118 // compute the gauge force from the first 4 links of all the Wilson loops that are
119 // summed in tC[k] (is valid to do so since first 4 links are the same for all Wilson
120 // loops and projection on Lie-algebra is linear)
121 get_wloop_force_from_wl_add(U, paths[k], tC, eps, K);
122 }
123 }
124}
125
126template <typename group, typename atype = hila::arithmetic_type<group>>
127void get_force_clover(const GaugeField<group> &U, out_only VectorField<Algebra<group>> &K,
128 atype eps = 1.0) {
129 // determine gauge force for clover action and store result in K
130 foralldir(d1) {
131 K[d1][ALL] = 0;
132 }
133 get_force_clover_add(U, K, eps);
134}
135
136
137#endif
Definition su2.h:7
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
Direction
Enumerator for direction that assigns integer to direction to be interpreted as unit vector.
Definition coordinates.h:34
constexpr Parity ALL
bit pattern: 011