32#include "gauge/stout_smear.h"
33#define STOUTSTEPS STOUTSMEAR
43int stout_nsteps = STOUTSTEPS;
60 std::string config_file;
68template <
typename group,
typename atype = hila::arithmetic_type<group>>
74 stout_smear(U, tU, stoutc, stout_nsteps);
83 double res = measure_s_bp(tU);
85 atype c12 = -1.0 / 12.0;
86 atype c11 = 1.0 - 8.0 * c12;
87 double res = measure_s_impr(tU, c11, c12);
90 atype c11 = 1.0 - 8.0 * c12;
91 double res = measure_s_impr(tU, c11, c12);
94 atype c11 = 1.0 - 8.0 * c12;
95 double res = measure_s_impr(tU, c11, c12);
97 double res = measure_s_log_plaq(tU);
99 double res = measure_s_wplaq(tU);
105template <
typename group,
typename atype = hila::arithmetic_type<group>>
108 atype eps = delta / group::size();
112 std::vector<GaugeField<group>> tUl(stout_nsteps + 1);
113 std::vector<VectorField<group>> tUKl(stout_nsteps);
114 std::vector<VectorField<group>> tstapl(stout_nsteps);
115 stout_smeark(U, tUl, tstapl, tUKl, stoutc);
131 get_force_bp_add(tU, tE, eps);
133 atype c12 = -1.0 / 12.0;
134 atype c11 = 1.0 - 8.0 * c12;
135 get_force_impr_add(tU, tE, eps * c11, eps * c12);
138 atype c11 = 1.0 - 8.0 * c12;
139 get_force_impr_add(tU, tE, eps * c11, eps * c12);
142 atype c11 = 1.0 - 8.0 * c12;
143 get_force_impr_add(tU, tE, eps * c11, eps * c12);
145 get_force_log_plaq_add(tU, tE, eps);
147 get_force_wplaq_add(tU, tE, eps);
153 stout_smeark_force(tUl, tstapl, tUKl, tE, KS, stoutc);
156 onsites(
ALL) E[d1][X] += KS[d1][X];
163template <
typename group,
typename atype = hila::arithmetic_type<group>>
167 onsites(
ALL) U[d][X] =
chexp(E[d][X].expand_scaled(delta)) * U[d][X];
171template <
typename group>
177 onsites(
ALL) e2 += E[d][X].squarenorm();
179 return e2.
value() / 2;
182template <
typename group>
184 const parameters &p,
double &plaq) {
186 plaq = p.beta * measure_s(U);
187 double e2 = measure_e2(E);
188 return plaq + e2 / 2;
191template <
typename group>
193 const parameters &p) {
196 return measure_action(U, E, p, plaq);
200template <
typename group,
typename atype = hila::arithmetic_type<group>>
205 update_U(U, E, p.dt / 2);
207 for (
int n = 0; n < p.trajlen - 1; ++n) {
208 update_E(U, E, p.beta * p.dt);
209 update_U(U, E, p.dt);
212 update_E(U, E, p.beta * p.dt);
213 update_U(U, E, p.dt / 2);
222template <
typename group>
226 static bool first =
true;
230 <<
"LMEAS: s-local plaq E^2 P.real P.imag\n";
233 auto slocal = measure_s(U) / (lattice.volume() * NDIM * (NDIM - 1) / 2);
234 auto plaq = measure_s_wplaq(U) / (lattice.volume() * NDIM * (NDIM - 1) / 2);
235 auto e2 = measure_e2(E) / (lattice.volume() * NDIM);
237 hila::out0 << string_format(
"MEAS % 0.6e % 0.6e % 0.6e % 0.6e % 0.6e", slocal, plaq, e2,
238 poly.real(), poly.imag())
246template <
typename group>
247void checkpoint(
const GaugeField<group> &U,
int trajectory,
const parameters &p) {
250 std::string config_file =
251 p.config_file +
"_" + std::to_string(
abs((trajectory + 1) / p.n_save) % 2);
257 outf.open(
"run_status", std::ios::out | std::ios::trunc);
258 outf <<
"trajectory " << trajectory + 1 <<
'\n';
259 outf <<
"seed " <<
static_cast<uint64_t
>(
hila::random() * (1UL << 61)) <<
'\n';
262 outf <<
"config name " << config_file <<
'\n';
265 std::stringstream msg;
267 hila::timestamp(msg.str().c_str());
270template <
typename group>
276 if (status.
open(
"run_status",
false,
false)) {
278 trajectory = status.
get(
"trajectory");
279 seed = status.
get(
"seed");
280 p.time_offset = status.
get(
"time");
282 std::string config_file = status.
get(
"config name");
285 U.config_read(config_file);
289 in.open(p.config_file, std::ios::in | std::ios::binary);
293 U.config_read(p.config_file);
306int main(
int argc,
char **argv) {
338 hila::out0 <<
"Using floating point epsilon: " << fp<ftype>::epsilon <<
"\n";
346 lsize = par.get(
"lattice size");
348 p.beta = par.get(
"beta");
350 p.dt = par.get(
"dt");
352 p.trajlen = par.get(
"trajectory length");
354 p.n_traj = par.get(
"number of trajectories");
356 p.n_therm = par.get(
"thermalization trajs");
358 p.gflow_freq = par.get(
"gflow freq");
360 p.gflow_max_l = par.get(
"gflow max lambda");
362 p.gflow_l_step = par.get(
"gflow lambda step");
364 p.gflow_a_accu = par.get(
"gflow abs. accuracy");
366 p.gflow_r_accu = par.get(
"gflow rel. accuracy");
368 long seed = par.get(
"random seed");
370 p.n_save = par.get(
"trajs/saved");
372 p.config_file = par.get(
"config name");
377 lattice.
setup(lsize);
388 Alg_gen<NCOLOR, ftype> genlist[NCOLOR * NCOLOR - 1];
391 for (
int i = 0; i < NCOLOR * NCOLOR - 1; ++i) {
392 hila::out0 << hila::prettyprint(genlist[i].to_matrix()) <<
"\n\n";
395 Matrix<NCOLOR * NCOLOR - 1, NCOLOR * NCOLOR - 1, ftype> omat;
398 chexp(2.0*genlist[NCOLOR+1].to_matrix()+1.5*genlist[1].to_matrix(), V, dV);
399 hila::out0 <<
"exp(A):\n" << hila::prettyprint(V) <<
"\n";
400 for (
int i = 0; i < NCOLOR; ++i) {
401 for (
int j = 0; j < NCOLOR; ++j) {
402 dV[i][j] = dV[i][j] * V.
dagger();
403 hila::out0 <<
"dexp(A)/dA[" << i <<
"][" << j <<
"]*exp(-A):\n"
404 << hila::prettyprint(dV[i][j]) <<
"\n";
407 project_to_algebra_bilinear(dV, omat, genlist);
408 hila::out0 <<
"dexp(A)^i/dA^j*exp(-A):\n" << hila::prettyprint(omat) <<
"\n";
411 Alg_gen<NCOLOR, ftype> genprodlist[NCOLOR * NCOLOR - 1][NCOLOR * NCOLOR - 1];
414 for (
int i = 0; i < NCOLOR * NCOLOR - 1; ++i) {
415 for (
int j = 0; j < NCOLOR * NCOLOR - 1; ++j) {
416 hila::out0 << hila::prettyprint(genprodlist[i][j].to_matrix()) <<
"\n\n";
419 project_to_algebra_bilinear(V, omat, genprodlist);
420 hila::out0 <<
"algebra_bilinear omat[][]=2*ReTr(t[].dagger()*V*t[]):\n"
421 << hila::prettyprint(omat) <<
"\n";
426 if (!restore_checkpoint(U, start_traj, p)) {
428 hila::out0 <<
"cold start: initial link variables set to identity\n";
436 auto orig_trajlen = p.trajlen;
440 double g_act_old, act_old, g_act_new, act_new;
441 g_act_old = p.beta * measure_s(U);
443 for (
int trajectory = start_traj; trajectory < p.n_traj; ++trajectory) {
444 if (trajectory < p.n_therm) {
448 if (trajectory < p.n_therm * 3.0 / 4.0) {
449 p.dt = orig_dt * (0.1 + 0.9 * 4.0 / 3.0 * trajectory / p.n_therm);
450 p.trajlen = orig_trajlen;
453 p.trajlen = orig_trajlen;
459 for (
int irej = 0; irej < nreject - 1; ++irej) {
463 hila::out0 <<
" thermalization step size(reduzed due to multiple reject) dt="
464 << std::setprecision(8) << p.dt <<
'\n';
466 hila::out0 <<
" thermalization step size dt=" << std::setprecision(8) << p.dt
469 }
else if (trajectory == p.n_therm) {
471 p.trajlen = orig_trajlen;
472 hila::out0 <<
" normal stepsize dt=" << std::setprecision(8) << p.dt <<
'\n';
475 update_timer.start();
482 act_old = g_act_old + measure_e2(E) / 2;
486 act_new = measure_action(U, E, p, g_act_new);
491 if (trajectory < p.n_therm) {
503 hila::out0 << std::setprecision(12) <<
"HMC " << trajectory <<
" S_TOT_start " << act_old
504 <<
" dS_TOT " << std::setprecision(6) << act_new - act_old
505 << std::setprecision(12);
507 hila::out0 <<
" REJECT" <<
" --> S_GAUGE " << g_act_old;
510 hila::out0 <<
" ACCEPT" <<
" --> S_GAUGE " << g_act_new;
511 g_act_old = g_act_new;
518 hila::out0 <<
"Measure_start " << trajectory <<
'\n';
520 measure_timer.start();
524 measure_timer.stop();
526 hila::out0 <<
"Measure_end " << trajectory <<
'\n';
528 if (trajectory >= p.n_therm) {
530 if (p.gflow_freq > 0 && trajectory % p.gflow_freq == 0) {
531 int gtrajectory = trajectory / p.gflow_freq;
532 if (p.gflow_l_step > 0) {
536 int nflow_steps = (int)(p.gflow_max_l / p.gflow_l_step);
539 hila::out0 <<
"Gflow_start " << gtrajectory <<
'\n';
543 ftype t_step = t_step0;
544 measure_gradient_flow_stuff(V, (ftype)0.0, t_step);
545 t_step = do_gradient_flow_adapt(V, (ftype)0.0, p.gflow_l_step, p.gflow_a_accu,
546 p.gflow_r_accu, t_step);
547 measure_gradient_flow_stuff(V, p.gflow_l_step, t_step);
550 for (
int i = 1; i < nflow_steps; ++i) {
553 do_gradient_flow_adapt(V, i * p.gflow_l_step, (i + 1) * p.gflow_l_step,
554 p.gflow_a_accu, p.gflow_r_accu, t_step);
556 measure_gradient_flow_stuff(V, (i + 1) * p.gflow_l_step, t_step);
561 hila::out0 <<
"Gflow_end " << gtrajectory <<
" time " << std::setprecision(3)
567 if (p.n_save > 0 && trajectory>=0 && (trajectory + 1) % p.n_save == 0) {
568 checkpoint(U, trajectory, p);
Array< n, m, T > exp(Array< n, m, T > a)
Exponential.
void reunitarize_gauge()
Reunitarize Gauge Field consisting of matrices.
void config_write(const std::string &filename) const
config_write writes the gauge field to file, with additional "verifying" header
static constexpr int size()
Returns size of Vector or square Matrix.
Rtype dagger() const
Hermitian conjugate of matrix.
Matrix class which defines matrix operations.
Special reduction class: enables delayed and non-blocking reductions, which are not possible with the...
const T value()
Return value of the reduction variable. Wait for the comms if needed.
Reduction & allreduce(bool b=true)
allreduce(bool) turns allreduce on or off. By default on.
Reduction & delayed(bool b=true)
deferred(bool) turns deferred on or off. By default turns on.
void setup(const CoordinateVector &siz)
General lattice setup.
T abs(const Complex< T > &a)
Return absolute value of Complex number.
#define foralldir(d)
Macro to loop over (all) Direction(s)
constexpr Parity ALL
bit pattern: 011
int chexp(const Matrix_t< n, m, T, MT > &mat, Matrix_t< n, m, T, MT > &omat, Mt(&pl)[n])
Calculate exp of a square matrix.
double random()
Real valued uniform random number generator.
int myrank()
rank of this node
std::ostream out0
This writes output only from main process (node 0)
void initialize(int argc, char **argv)
Read in command line arguments. Initialise default stream and MPI communication.
void seed_random(uint64_t seed, bool device_rng=true)
Seed random generators with 64-bit unsigned value. On MPI shuffles the seed so that different MPI ran...
T broadcast(T &var, int rank=0)
Broadcast the value of var to all MPI ranks from rank (default=0).
void finishrun()
Normal, controlled exit - all nodes must call this. Prints timing information and information about c...
Complex< double > measure_polyakov(const GaugeField< T > &U, Direction dir=Direction(NDIM - 1))
Measure Polyakov lines to direction dir.
void do_trajectory(GaugeField< group > &U, const parameters &p)
Evolve gauge field.
void measure_stuff(const GaugeField< group > &U, const parameters &p)
Measures Polyakov lines and Wilson action.