HILA
Loading...
Searching...
No Matches
improved_action.h
Go to the documentation of this file.
1/** @file improved_action.h */
2
3#ifndef IMPROVED_ACTION_H_
4#define IMPROVED_ACTION_H_
5
6#include "hila.h"
8
9/*
10* parameters for for different improved actions:
11// DBW2 action parameters:
12//ftype c12=-1.4088; // rectangle weight
13//ftype c11=1.0-8.0*c12; // plaquette weight
14
15// Iwasaki action parameters:
16//ftype c12=-0.331; // rectangle weight
17//ftype c11=1.0-8.0*c12; // plaquette weight
18
19// LW action parameters:
20//ftype c12=-1.0/12.0; // rectangle weight
21//ftype c11=1.0-8.0*c12; // plaquette weight
22
23// Wilson plaquette action parameters:
24//ftype c12=0;
25//ftype c11=1.0;
26*/
27
28
29// functions for improved action -S_{impr} = \p.beta/N*( c11*\sum_{plaquettes P} ReTr(P) +
30// c12*\sum_{1x2-rectangles R} ReTr(R) )
31/*
32template <typename group,typename atype=hila::arithmetic_type<group>>
33void get_force_impr_add(const GaugeField<group>& U,VectorField<Algebra<group>>& K,atype c11,atype
34c12) {
35 // compute the force for the improved action -S_{impr}=\beta/N*(c11*ReTr(plaq)+c12*ReTr(rect))
36 // and write it to K
37
38 // plaquette part:
39 foralldir(dir1) foralldir(dir2) if(dir1<dir2) {
40 Direction path[4]={dir1,dir2,-dir1,-dir2};
41 get_wloop_force_add(U,path,c11,K);
42 }
43
44 if(c12!=0) {
45 // rectangle part:
46 foralldir(dir1) foralldir(dir2) if(dir1!=dir2) {
47 Direction path[6]={dir1,dir2,dir2,-dir1,-dir2,-dir2};
48 get_wloop_force_add(U,path,c12,K);
49 }
50 }
51}
52*/
53
54template <typename group, typename atype = hila::arithmetic_type<group>>
55double measure_s_impr(const GaugeField<group> &U, atype c11, atype c12) {
56 // measure the improved action for dir1<dir2
57 Reduction<double> stot = 0;
58 stot.allreduce(false).delayed(true);
59 Field<group> tP;
60 foralldir(dir1) foralldir(dir2) if (dir1 < dir2) {
61 U[dir2].start_gather(dir1, ALL);
62 U[dir1].start_gather(dir2, ALL);
63 onsites(ALL) {
64 // plaquettes:
65 tP[X] = U[dir1][X] * U[dir2][X + dir1] * (U[dir2][X] * U[dir1][X + dir2]).dagger();
66
67 // plaquette part:
68 stot += c11 * (1.0 - real(trace(tP[X])) / group::size());
69 }
70
71 tP.start_gather(dir1, ALL);
72 tP.start_gather(dir2, ALL);
73
74 onsites(ALL) {
75 // 2x1-rectangle part:
76 stot +=
77 c12 * (1.0 - real(trace(U[dir1][X] * tP[X + dir1] * U[dir1][X].dagger() * tP[X])) /
78 group::size());
79 // 1x2-rectangle part:
80 stot +=
81 c12 * (1.0 - real(trace(tP[X] * U[dir2][X] * tP[X + dir2] * U[dir2][X].dagger())) /
82 group::size());
83 }
84 }
85 return stot.value();
86}
87
88template <typename group, typename atype = hila::arithmetic_type<group>>
89void get_force_impr_add(const GaugeField<group> &U, VectorField<Algebra<group>> &K, atype c11,
90 atype c12) {
91 // compute the force for the improved action -S_{impr}=\beta/N*(c11*ReTr(plaq)+c12*ReTr(rect))
92 // in an even faster way and add the results to K
93 Field<group> ustap;
94 Field<group> lstap;
95 Field<group> tstap;
96 bool first;
97 foralldir(dir1) {
98 first = true;
99 foralldir(dir2) if (dir1 != dir2) {
100 U[dir2].start_gather(dir1, ALL);
101 U[dir1].start_gather(dir2, ALL);
102
103 // get upper (dir1,dir2) and lower (dir1,-dir2) staples
104 onsites(ALL) {
105 lstap[X] = (U[dir1][X] * U[dir2][X + dir1]).dagger() * U[dir2][X];
106 ustap[X] = U[dir2][X + dir1] * (U[dir2][X] * U[dir1][X + dir2]).dagger();
107 }
108
109 lstap.start_gather(-dir2, ALL);
110
111 // sum the staples for dir1-link
112 if (first) {
113 onsites(ALL) {
114 tstap[X] = ustap[X];
115 tstap[X] += lstap[X - dir2];
116 }
117 first = false;
118 } else {
119 onsites(ALL) {
120 tstap[X] += ustap[X];
121 tstap[X] += lstap[X - dir2];
122 }
123 }
124
125 // compute rectangle contribution to force (if any) and add it to K
126 if (c12 != 0) {
127 // compose rectangle and store it in ustap
128 onsites(ALL) ustap[X] = lstap[X - dir2].dagger() * ustap[X];
129
130 // corresponding path
131 std::vector<Direction> path = {-dir2, dir1, dir2, dir2, -dir1, -dir2};
132
133 // compute rectangle force and add it to K
134 get_wloop_force_from_wl_add(U, path, ustap, c12, K);
135 }
136 }
137 // compute plaquette contribution to force and add it to K
138 onsites(ALL) {
139 K[dir1][X] -= (U[dir1][X] * tstap[X]).project_to_algebra_scaled(c11);
140 }
141 }
142}
143
144template <typename group, typename atype = hila::arithmetic_type<group>>
145void get_force_impr(const GaugeField<group> &U, out_only VectorField<Algebra<group>> &K, atype c11,
146 atype c12) {
147 // compute the force for the improved action -S_{impr}=\beta/N*(c11*ReTr(plaq)+c12*ReTr(rect))
148 // in an even faster way and write the results to K
149 foralldir(d1) {
150 K[d1][ALL] = 0;
151 }
152 get_force_impr_add(U, K, c11, c12);
153}
154
155
156#endif
Array< n, m, hila::arithmetic_type< T > > real(const Array< n, m, T > &arg)
Return real part of Array.
Definition array.h:688
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
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