HILA
Loading...
Searching...
No Matches
log_action.h
Go to the documentation of this file.
1/** @file log_action.h */
2
3#ifndef LOG_ACTION_H_
4#define LOG_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_log_clover_mat(const GaugeField<group> &U, Direction dir1, Direction dir2,
13 Field<group> &Csum) {
14 // sets Csum=(C[0]+C[1]+C[2]+C[3])/4 where C[0],C[1],C[2],C[3] are the logs of the clover-leave
15 // matrices in counter-clockwise order
16
17 U[dir2].start_gather(dir1, ALL);
18 U[dir1].start_gather(dir2, ALL);
19
20 Field<group> tF0, tF1;
21
22 onsites(ALL) {
23 // log of dir1-dir2-plaquette that starts and ends at X; corresponds to F[dir1][dir2]
24 // at center location X+dir1/2+dir2/2 of plaquette:
25 tF0[X] = log((U[dir1][X] * U[dir2][X + dir1] * (U[dir2][X] * U[dir1][X + dir2]).dagger()))
26 .expand();
27 // parallel transport to X+dir1
28 tF1[X] = U[dir1][X].dagger() * tF0[X] * U[dir1][X];
29 }
30
31 tF1.start_gather(-dir1, ALL);
32 onsites(ALL) {
33 tF0[X] += tF1[X - dir1];
34 }
35
36 U[dir2].start_gather(-dir2, ALL);
37 tF0.start_gather(-dir2, ALL);
38 onsites(ALL) {
39 // get F[dir1][dir2] at X from average of the (parallel transported) F[dir1][dir2] from
40 // the centers of all dir1-dir2-plaquettes that touch X :
41 Csum[X] = (tF0[X] + U[dir2][X - dir2].dagger() * tF0[X - dir2] * U[dir2][X - dir2]) * 0.25;
42 }
43}
44
45template <typename group, typename atype = hila::arithmetic_type<group>>
46atype measure_s_log(const GaugeField<group> &U) {
47 // measure the log action for dir1<dir2
48 // (just to have same normalization as with plaquette action)
49 Reduction<double> stot = 0;
50 stot.allreduce(false).delayed(true);
51 Field<group> Csum;
52 foralldir(dir1) foralldir(dir2) if (dir1 < dir2) {
53 get_log_clover_mat(U, dir1, dir2, Csum);
54 onsites(ALL) {
55 stot += -0.5 * real(mul_trace(Csum[X], Csum[X]));
56 }
57 }
58 return (atype)stot.value() / group::size();
59}
60
61template <typename group, typename atype = hila::arithmetic_type<group>>
62void get_force_log_add(const GaugeField<group> &U, VectorField<Algebra<group>> &K,
63 atype eps = 1.0) {
64 // compute gauge force for clover action and add result to K
65 Field<group> Csum;
66
67 foralldir(dir1) foralldir(dir2) if (dir1 < dir2) {
68 // get the 4 clover leave matrices C[][X] (ordered counter-clockwise in (dir1,dir2)-plane)
69 // and their anti-hermitian traceless sum Csum[X]
70 get_log_clover_mat(U, dir1, dir2, Csum);
71
72 // define for each clover leave matrix the correspondin path of gauge links
73 std::vector<Direction> paths[4] = {{dir1, dir2, -dir1, -dir2},
74 {dir2, -dir1, -dir2, dir1},
75 {-dir1, -dir2, dir1, dir2},
76 {-dir2, dir1, dir2, -dir1}};
77
78 // the force on the link (Y,dir) from the clover term with foot point X is:
79 // F[dir][Y] = d_{Y,dir}(0.5*Tr(Csum*Csum))
80 // = 0.5*Tr((d_{Y,dir}Csum)*Csum) + 0.5*Tr(Csum*(d_{Y,dir}Csum))
81 // = Tr((d_{Y,dir}Csum)*Csum)
82 // = 0.25*\sum_k Tr((d_{Y,dir}log(C[k][X]))*Csum[X])
83 // with: d_{Y,dir}log(C[k][X]) = (d_{Y,dir}C[k][X])*C[k][X].dagger()
84 // where C[k][X] is the k-th clover leave matrix,
85 // Furthermore: (d^a_{Y,dir[0]}C[0][X])*C[0][X].dagger() = T^a
86 // (d^a_{Y,dir[1]}C[0][X])*C[0][X].dagger() = U[dir[0]][X].dagger() * T^a *
87 // U[dir[0]][X]
88 // ...
89 // ...
90 // We therefore need only Csum/4 to compute the log-action force:
91 Csum[ALL] *= 0.25;
92 for (int k = 0; k < 4; ++k) {
93 // compute the gauge force from the first 4 links of all the Wilson loops that are
94 // summed in tC[k] (is valid to do so since first 4 links are the same for all Wilson
95 // loops and projection on Lie-algebra is linear)
96 get_wloop_force_from_wl_add(U, paths[k], Csum, eps, K);
97 }
98 }
99}
100
101template <typename group, typename atype = hila::arithmetic_type<group>>
102void get_force_log(const GaugeField<group> &U, VectorField<Algebra<group>> &K, atype eps = 1.0) {
103 // determine gauge force for clover action and store result in K
104 foralldir(d1) {
105 K[d1][ALL] = 0;
106 }
107 get_force_log_add(U, K, eps);
108}
109
110
111#endif
Array< n, m, hila::arithmetic_type< T > > real(const Array< n, m, T > &arg)
Return real part of Array.
Definition array.h:688
Array< n, m, T > log(Array< n, m, T > a)
Logarithm.
Definition array.h:1069
Definition su2.h:7
The field class implements the standard methods for accessing Fields. Hilapp replaces the parity acce...
Definition field.h:61
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
auto mul_trace(const Mtype1 &a, const Mtype2 &b)
Returns trace of product of two matrices.
Definition matrix.h:2309