diff --git a/CMakeLists.txt b/CMakeLists.txt index 35b5f16..ab37fcd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -51,7 +51,6 @@ configure_file( set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -D_USE_OMP_") find_package(yaml-cpp) - include_directories(${YAML_CPP_INCLUDE_DIR} ${YAML_CPP_INCLUDE_DIR}/yaml-cpp/) find_package(xtl) @@ -60,6 +59,11 @@ configure_file( find_package(xtensor) include_directories(${xtensor_INCLUDE_DIRS} ${xtensor_INCLUDE_DIRS}/xtensor/) + + if(${HIGHFIVE_FORMAT}) + find_package(HighFive REQUIRED) + endif() + # We could just have two `add_executable` blocks listing almost all source # files. However, CMake would then compile the source files twice, once for # each executable. The rationale is that different flags can be used for @@ -90,9 +94,31 @@ add_executable(test-groups test.cc) add_executable(scaling scaling.cc) add_executable(try tryreduce.cc) +set(test_progs "links" "smearing" "yaml") + +if(${HIGHFIVE_FORMAT}) + set(test_progs ${test_progs} "h5-gen" "h5-dump") +endif() + +set(test_progs_fullname "") +foreach(target ${test_progs}) + add_executable(test-${target}.exe test/${target}.cpp) + list(APPEND test_progs_fullname test-${target}.exe) # now FOO="a;b" +endforeach(target) + + +target_link_libraries(test-yaml.exe ${YAML_CPP_LIBRARIES}) + +foreach(target ${test_progs_fullname}) + if(${HIGHFIVE_FORMAT}) + target_link_libraries(${target} HighFive) + endif() +endforeach(target) + + target_link_directories(${CMAKE_PROJECT_NAME} PUBLIC ${YAML_CPP_LIBRARY_DIR}) -foreach(target u1-main su2-main su2-kramers test-groups scaling try) +foreach(target u1-main su2-main su2-kramers test-groups scaling try ${test_progs_fullname}) target_link_libraries(${target} su2 ${YAML_CPP_LIBRARIES}) if(Boost_FOUND) @@ -110,4 +136,9 @@ install(TARGETS ${PROJECT_NAME} install(TARGETS u1-main su2-main) install(TARGETS su2-kramers test-groups scaling try) +# tests + +foreach(target ${test_progs}) + install(TARGETS test-${target}.exe LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) +endforeach(target) diff --git a/base_program.hpp b/base_program.hpp index f6a51a8..a568d3f 100644 --- a/base_program.hpp +++ b/base_program.hpp @@ -12,7 +12,7 @@ #pragma once #include "flat-energy_density.hh" -#include "flat-gauge_energy.hpp" +#include "gauge_energy.hh" #include "flat-sweep.hh" // flat spacetime #include "gaugeconfig.hh" #include "io.hh" @@ -21,7 +21,7 @@ #include "random_gauge_trafo.hh" // #include "rotating-energy_density.hpp" // rotating spacetime // #include "rotating-gauge_energy.hpp" // rotating spacetime -//#include "rotating-sweep.hpp" // rotating spacetime +// #include "rotating-sweep.hpp" // rotating spacetime #include "su2.hh" #include "u1.hh" #include "vectorfunctions.hh" @@ -423,10 +423,9 @@ template class base_program { if ((*this).omeas.verbosity > 0) { std::cout << "## online measuring: J^{PC} glueball correlators.\n"; } - if ((*this).omeas.glueball.correlator) { - omeasurements::meas_glueball_correlator(omeas.glueball.interpolator_type, - U, i, (*this).omeas); - } + + const std::string glb_interp = omeas.glueball.interpolator_type; + omeasurements::meas_glueball_correlator(glb_interp, U, i, (*this).omeas); } return; } diff --git a/doc/users/2d-materials.csl b/doc/2d-materials.csl similarity index 100% rename from doc/users/2d-materials.csl rename to doc/2d-materials.csl diff --git a/doc/users/Makefile b/doc/Makefile similarity index 65% rename from doc/users/Makefile rename to doc/Makefile index 3ed4250..724d76f 100644 --- a/doc/users/Makefile +++ b/doc/Makefile @@ -1,21 +1,18 @@ # Makefile # file paths and variables -STAGG = ./staggered_fermions.Rmd -u1 = ./u1.Rmd +main = ./main.Rmd define generate_html Rscript -e 'library(rmarkdown); rmarkdown::render("$(1)", "html_document")' + Rscript -e 'library(rmarkdown); rmarkdown::render("$(1)", "pdf_document")' endef -all: staggered u1 +all: main make clean -u1: $(u1) - $(call generate_html,$(u1)) - -staggered: $(STAGG) - $(call generate_html,$(STAGG)) +main: $(main) + $(call generate_html,$(main)) clean: diff --git a/doc/README.md b/doc/README.md new file mode 100644 index 0000000..973587a --- /dev/null +++ b/doc/README.md @@ -0,0 +1,4 @@ +# Users guide + +This directory contains the users guide to this project. +In order to produce the `html` and `pdf` output run `make all` from the present directory. diff --git a/doc/users/biblio.bib b/doc/biblio.bib similarity index 100% rename from doc/users/biblio.bib rename to doc/biblio.bib diff --git a/doc/users/u1.Rmd b/doc/howto.Rmd similarity index 76% rename from doc/users/u1.Rmd rename to doc/howto.Rmd index 47c332b..9d7d0f1 100644 --- a/doc/users/u1.Rmd +++ b/doc/howto.Rmd @@ -1,28 +1,7 @@ --- -title: u1 -author: Simone Romiti output: html_document -bibliography: biblio.bib -csl: 2d-materials.csl -link-citations: true --- -## Introduction and main features - -This repository contains a program called `main-u1.cpp`, which is capable of both generating gauge configurations and doing measurements on them for a $U(1)$ theory in $d=2,3,4$ dimensions. -At present the gauge configurations can be generated with the following Markov Chain Monte Carlo algorithms (see @gattringer2009quantum for the theoretical background): - -- Metropolis -- Hybrid Monte Carlo - -The measurements can be generated either online (during the MCMC) or offline (on previously generated configurations). At present the supported observables are: - -- Plaquette standard and clover-improved energy -- Glueball correlators -- "pion" correlator -- Wilson loops and static potential among $\bar{q} q$ pairs -- Wilson gradient flow - ## Program execution The `main-u1.cpp` program works using an input file `main.input` as follows: @@ -101,7 +80,4 @@ WRITE DOCUMENTATION HERE - `Wloop` - `potential` block. - `glueball` block. The user can specify the interpolators and the APE smearing parameters -- `gradient_flow` block - -## References - +- `gradient_flow` block \ No newline at end of file diff --git a/doc/main.Rmd b/doc/main.Rmd new file mode 100644 index 0000000..499b641 --- /dev/null +++ b/doc/main.Rmd @@ -0,0 +1,74 @@ +--- +title: Documentation `su2` project +author: Simone Romiti +output: + html_document: + toc: true + theme: united +bibliography: biblio.bib +csl: 2d-materials.csl +link-citations: true +--- + +The present repository contains routines for numerical simulations of $SU(2)$ and $U(1)$ lattice gauge theories in $d=2,3,4$ dimensions. + +# Library routines + +The gauge configurations can be generated with the following Markov Chain Monte Carlo algorithms (see @gattringer2009quantum for the theoretical background): + +- Metropolis +- Hybrid Monte Carlo + +The measurements can be generated either online (during the MCMC) or offline (on previously generated configurations). At present the supported observables are: + +- Plaquette standard and clover-improved energy +- Glueball correlators +- "pion" correlator +- Wilson loops and static potential among $\bar{q} q$ pairs +- Wilson gradient flow + +### Programs info + +- Programs: `main-u1.cpp` and `main-su2.cpp` +- Parameters parsing: `yaml` input file +- Input file example: see [./main.input](main.input) + Notes: + - Each run produces/updates a file named `nconf_counter.txt` in the directory storing the gauge configurations. This file contains a header line and the values below. An example: + ``` + heat i path_conf + 1 1140 ./confs/config_u1.4.4.1.16.b1.900000.1140 + ``` + - `restart` and `heat` cannot be passed simultaneously + - If `restart: true` is passed, the last configuration id and path are read from the `nconf_counter.txt` file. + - The program can do online and offline measurements + +##### `hadron` folder + +In the `hadron/` folder are found scripts for data manipulation such that online measurements are more easily readable with the `R` library `hadron`: https://github.com/HISKP-LQCD/hadron: + +- `glueball.py`: converts the data for the interpolators and correlators of the glueball. + 1. This is done in a separate step and not directly in the `C++` code on purpose. One may want to generate configurations and do offline measurements afterwards in parallel (multiple jobs, one for each bunch of measurements). + 2. The scripts takes care of deleting the files in the old format. For each choice of the parameters and interpolator type, only 4 files are saved in the corresponding folders: + + - `iconf.dat` : list of trajectory indices that have been put in `hadron` format + + - `pp.dat`, `pm.dat`, `mp.dat`, `mm.dat`: correlators + + 3. The storage space saved in this way is mostly negligible, but avoids potential filesystem slowdowns due to the access to large number of files. Besides, the analysis scripts won't have to cope with reading the configurations one by one. + 4. The script uses the same input file used for the runs + +### tests + +In the `test/` folder the user can find various simple programs implementing the routines of this repository. For further details please look at the `test/README.md` file. + +# How to use the library + +```{r child = './howto.Rmd'} +``` + +# Theoretical background + +Fermions are included as a degenerate doublet using the staggered discretization. + +```{r child = './staggered_fermions.Rmd'} +``` diff --git a/doc/users/main.input b/doc/main.input similarity index 100% rename from doc/users/main.input rename to doc/main.input diff --git a/doc/users/staggered_fermions.Rmd b/doc/staggered_fermions.Rmd similarity index 97% rename from doc/users/staggered_fermions.Rmd rename to doc/staggered_fermions.Rmd index 417a58c..28c566b 100644 --- a/doc/users/staggered_fermions.Rmd +++ b/doc/staggered_fermions.Rmd @@ -1,12 +1,14 @@ --- -title: Staggered fermions -author: Simone Romiti +#title: Staggered fermions +#author: Simone Romiti output: html_document -bibliography: biblio.bib -csl: 2d-materials.csl -link-citations: true +#bibliography: biblio.bib +#csl: 2d-materials.csl +#link-citations: true --- +## Staggered fermions + For staggered fermions the Dirac operator is: $$ diff --git a/doc/bibliography.bib b/doc/su2_partitioning/bibliography.bib similarity index 100% rename from doc/bibliography.bib rename to doc/su2_partitioning/bibliography.bib diff --git a/doc/charge-notes.xoj b/doc/su2_partitioning/charge-notes.xoj similarity index 100% rename from doc/charge-notes.xoj rename to doc/su2_partitioning/charge-notes.xoj diff --git a/doc/doc.Rmd b/doc/su2_partitioning/doc.Rmd similarity index 100% rename from doc/doc.Rmd rename to doc/su2_partitioning/doc.Rmd diff --git a/doc/ref.txt b/doc/su2_partitioning/ref.txt similarity index 100% rename from doc/ref.txt rename to doc/su2_partitioning/ref.txt diff --git a/doc/reference-values.txt b/doc/su2_partitioning/reference-values.txt similarity index 100% rename from doc/reference-values.txt rename to doc/su2_partitioning/reference-values.txt diff --git a/doc/su2strong.R b/doc/su2_partitioning/su2strong.R similarity index 100% rename from doc/su2strong.R rename to doc/su2_partitioning/su2strong.R diff --git a/doc/u1strong.R b/doc/su2_partitioning/u1strong.R similarity index 100% rename from doc/u1strong.R rename to doc/su2_partitioning/u1strong.R diff --git a/doc/users/README.md b/doc/users/README.md deleted file mode 100644 index 485436c..0000000 --- a/doc/users/README.md +++ /dev/null @@ -1,31 +0,0 @@ -# Users guide - -This directory contains the users guide to this project. -The latter is written in `Rmd` format. -In order to produce the `html` and `pdf` output run `make all` from the present directory. - -## Theoretical background - -* At the moment fermions are included as a degenerate doublet using the staggered discretization: see [./staggered_fermions.Rmd](./staggered_fermions.Rmd) - - -## Library usage - -The library contains code for both $SU(2)$ and $U(1)$ in $d=2,3,4$ dimensions. - -- Programs: `main-u1.cpp` and `main-su2.cpp` -- Parameters parsing: `yaml` input file -- Input file example: see [./main.input](main.input) - Notes: - - Each run produces/updates a file named `nconf_counter.txt` in the directory storing the gauge configurations. This file contains a header line and the values below. At the moment an example: - ``` - heat i path_conf - 1 1140 ./confs/config_u1.4.4.1.16.b1.900000.1140 - ``` - - `restart` and `heat` cannot be passed simultaneously - - If `restart: true` is passed, the last configuration id and path are read from the `nconf_counter.txt` file. - - The program can do online and offline measurements - -## tests - -In the `test/` folder the user can find various simple programs implementing the routines of this repository. For further details please look at the `test/README.md` file. diff --git a/hadron/README.md b/hadron/README.md new file mode 100644 index 0000000..c5f951e --- /dev/null +++ b/hadron/README.md @@ -0,0 +1,2 @@ +This folder contains scripts for formatting online (offline) measurements. +Tha data is formatted such that is more easily readable by the routines of the "R" library "hadron": https://github.com/HISKP-LQCD/hadron \ No newline at end of file diff --git a/hadron/glueball.py b/hadron/glueball.py new file mode 100644 index 0000000..fc6ff36 --- /dev/null +++ b/hadron/glueball.py @@ -0,0 +1,133 @@ + +import pandas as pd +import numpy as np +import os +import yaml +import argparse +import glob + + +parser = argparse.ArgumentParser( + prog='glueball', + description="""format glueball correlators and interpolator files in the "hadron" format""") + +parser.add_argument('-f', '--inputfile', help="Path to the same yaml input file used for the run") +args = parser.parse_args() + +# loading and parsing the input file for the needed info +nd = yaml.load(open(str(args.inputfile)), Loader=yaml.Loader) + +T = nd["geometry"]["T"] +nd_omeas = nd["omeas"] +resdir = nd_omeas["res_dir"] +nd_glueball = nd_omeas["glueball"] + +save_corr = nd_glueball["correlator"] +nd_interpolator = nd_glueball["interpolator"] + +do_APE_smearing = bool(nd_glueball["do_APE_smearing"]) +n_APE = [] +alpha_APE = 0.0 +if (do_APE_smearing): + nd_APE_smearing = nd_glueball["APE_smearing"] + n_APE = nd_APE_smearing["n"] + alpha = nd_APE_smearing["alpha"] +#### + +rmin = int(nd_interpolator["rmin"]) +rmax = int(nd_interpolator["rmax"]) +interp_type = nd_interpolator["type"] +interp_save = nd_interpolator["save"] + + +def hadronize(operator, name_obj, name_ij): + d1 = resdir+"/"+"glueball/" + d2 = d1 + operator + "/" + interp_type + "/" + d2 += "smearAPEn{n_APE}alpha{alpha}/".format(n_APE=n_APE, alpha=alpha)+"/" + d2 += name_ij+"/" + print(d2) + # list all files in d2 and extract configuration numbers from names + file_names = glob.glob(d2+"/"+name_obj+"_*") + file_index = [int(x.split(name_obj+"_")[1]) for x in file_names] + # + # load .dat file with indices of formatted configurations, if existing + iconf_file = d2 + "/iconfs.dat" + # already formatted configurations + file_index_old = [] + if os.path.exists(iconf_file): + print("Old configurations found") + file_index_old = pd.read_csv( + iconf_file, sep=" ", header=None)[0].tolist() + #### + n_conf_old = len(file_index_old) + # new matrices for a each J^{PC} + # dict({"pp": pd.DataFrame(), "pm": pd.DataFrame(), "mp": pd.DataFrame(), "mm":pd.DataFrame()}) + df_new = dict() + print("Reading old data if present") + for j_pc in ["pp", "pm", "mp", "mm"]: + # loop over all J^{PC} quantum numbers + df_new[j_pc] = pd.DataFrame(np.zeros(shape=(1, T+1)), columns=["i"]+[ + str(t) for t in range(T)], index=["remove"]) + df_new[j_pc] = df_new[j_pc].astype({"i": int}) + if n_conf_old > 0: + path_old = d2+"/"+j_pc+".dat" + if os.path.exists(iconf_file) and (not os.path.exists(path_old)): + print("Error! Something bad happened, please redo the online measurements:\n", + iconf_file, "\nexists, but\n", path_old, "\ndoesn't.") + raise ValueError + #### + df_old = pd.read_csv(path_old, sep=" ", dtype={"i": int}) + df_new[j_pc] = pd.concat([df_new[j_pc], df_old], axis=0) + # print(df_new[j_pc]) + #### + #### + print("loop over new and updated files") + for i_f in range(len(file_index)): + path_i = file_names[i_f] + print(i_f, path_i) + idx_i = file_index[i_f] + df_i = pd.read_csv(path_i, sep=" ") + for j_pc in ["pp", "pm", "mp", "mm"]: + # loop over all J^{PC} quantum numbers + arr1 = df_i[j_pc].to_numpy() + if idx_i in file_index_old: + # removing old configuration + df_new[j_pc] = df_new[j_pc][df_new[j_pc]["i"] != idx_i] + #### + df1 = pd.DataFrame(arr1).transpose() + df1.insert(0, "i", idx_i) + df1.columns = ["i"] + [str(t) for t in range(T)] + # update new array with the new data + df_new[j_pc] = pd.concat([df_new[j_pc], df1]) + #### + #### + #### + print("Cleaning up the dataframe format before exporting") + save_iconf = True + for j_pc in ["pp", "pm", "mp", "mm"]: + df_new[j_pc] = df_new[j_pc].drop("remove") + save_df_new = (df_new[j_pc].shape[0]>0) + save_iconf = save_iconf and save_df_new + # saving + if save_df_new: + df_new[j_pc] = df_new[j_pc].sort_values(by=['i']) + df_new[j_pc].to_csv(d2 + "/"+j_pc+".dat", sep=" ", index=False) + #### + print("Saving the new iconf file") + if save_iconf: + pd.DataFrame(sorted(file_index)).to_csv( + d2+"iconfs.dat", header=False, sep=" ", index=False) + #### + print("Removing data in the old format") + for p in file_names: + print("removing:", p) + os.remove(p) + #### +#### + + +for i in range(rmin, rmax+1): + hadronize("interpolator", "phi", str(i)) + for j in range(rmin, i+1): + hadronize("correlator", "C_glueball", str(i)+"_"+str(j)) + diff --git a/include/accum_type.hh b/include/accum_type.hh index 2990e5b..5c6384c 100644 --- a/include/accum_type.hh +++ b/include/accum_type.hh @@ -1,5 +1,5 @@ #pragma once -template struct accum_type { +template struct accum_type { typedef T type; }; diff --git a/include/dagger.hh b/include/dagger.hh new file mode 100644 index 0000000..3af0f9d --- /dev/null +++ b/include/dagger.hh @@ -0,0 +1,3 @@ +#pragma once + +template inline T dagger(const T& x); diff --git a/include/detDDdag_monomial.hh b/include/detDDdag_monomial.hh index ea21126..040d854 100644 --- a/include/detDDdag_monomial.hh +++ b/include/detDDdag_monomial.hh @@ -10,7 +10,7 @@ #pragma once #include "adjointfield.hh" -#include "flat-gauge_energy.hpp" +#include "gauge_energy.hpp" #include "gaugeconfig.hh" #include "get_staples.hh" #include "hamiltonian_field.hh" diff --git a/include/flat-energy_density.hh b/include/flat-energy_density.hh index fec480b..e370ac1 100644 --- a/include/flat-energy_density.hh +++ b/include/flat-energy_density.hh @@ -80,39 +80,39 @@ namespace flat_spacetime { accum G[4][4]; for (size_t mu = mu_start; mu < U.getndims() - 1; mu++) { for (size_t nu = mu + 1; nu < U.getndims(); nu++) { - x1[mu] += 1; - x2[nu] += 1; + x1[mu] += 1; // x + mu + x2[nu] += 1; // x + nu accum leaf = U(x, mu) * U(x1, nu) * U(x2, mu).dagger() * U(x, nu).dagger(); - x1[mu] -= 1; - x2[nu] -= 1; + x1[mu] -= 1; // x + x2[nu] -= 1; // x if (cloverdef) { - x1[mu] -= 1; - x1[nu] += 1; - x2[mu] -= 1; + x1[mu] -= 1; // x - mu + x1[nu] += 1; // x - mu + nu + x2[mu] -= 1; // x - mu leaf += U(x, nu) * U(x1, mu).dagger() * U(x2, nu).dagger() * U(x2, mu); - x1[mu] += 1; - x1[nu] -= 1; - x2[mu] += 1; + x1[mu] += 1; // x + nu + x1[nu] -= 1; // x + x2[mu] += 1; // x - x1[mu] -= 1; - x2[mu] -= 1; - x2[nu] -= 1; - x3[nu] -= 1; + x1[mu] -= 1; // x - mu + x2[mu] -= 1; // x - mu + x2[nu] -= 1; // x - mu -nu + x3[nu] -= 1; // x - nu leaf += U(x1, mu).dagger() * U(x2, nu).dagger() * U(x2, mu) * U(x3, nu); - x1[mu] += 1; - x2[mu] += 1; - x2[nu] += 1; - x3[nu] += 1; + x1[mu] += 1; // x + x2[mu] += 1; // x + nu + x2[nu] += 1; // x + x3[nu] += 1; // x - x1[nu] -= 1; - x2[nu] -= 1; - x2[mu] += 1; + x1[nu] -= 1; // x - nu + x2[nu] -= 1; // x - nu + x2[mu] += 1; // x + mu -nu leaf += U(x1, nu).dagger() * U(x1, mu) * U(x2, nu) * U(x, mu).dagger(); - x1[nu] += 1; - x2[nu] += 1; - x2[mu] -= 1; + x1[nu] += 1; // x + x2[nu] += 1; // x + mu + x2[mu] -= 1; // x } // traceless and anti-hermitian // here we include a factor 1/2 already diff --git a/include/flat-gauge_energy.hpp b/include/flat-gauge_energy.hpp deleted file mode 100644 index 0118cc3..0000000 --- a/include/flat-gauge_energy.hpp +++ /dev/null @@ -1,139 +0,0 @@ -/** - * @file flat_spacetime_gauge_energy.hpp - * @author Simone Romiti (simone.romiti@uni-bonn.de) - * @brief gauge energy in flat spacetime (euclidean metric) - * @version 0.1 - * @date 2022-05-20 - * - * @copyright Copyright (c) 2022 - * - */ -#pragma once - -#include "gaugeconfig.hh" -#ifdef _USE_OMP_ -#include -#endif - -namespace flat_spacetime { - - /** - * @brief real part of the trace of the Wilson plaquette - * - * \sum_mu \sum_nu 0 are calculated - * @param xi anisotropy in the action: S_G \supset (\beta/xi)*P_{0i} + (xi*\beta)*P_{ij} - * @return double - */ - template - double retr_sum_Wplaquettes(const gaugeconfig &U, - const double &xi = 1.0, - const bool &anisotropic = false, - const bool &spatial_only = false) { - double res = 0.; - size_t startmu = spatial_only; // 0 if spatial_only==false, 1 if spatial_only==true - - if (!anisotropic) { -#pragma omp parallel for reduction(+ : res) - for (size_t x0 = 0; x0 < U.getLt(); x0++) { - for (size_t x1 = 0; x1 < U.getLx(); x1++) { - for (size_t x2 = 0; x2 < U.getLy(); x2++) { - for (size_t x3 = 0; x3 < U.getLz(); x3++) { - std::vector x = {x0, x1, x2, x3}; - std::vector xplusmu = x; - std::vector xplusnu = x; - for (size_t mu = startmu; mu < U.getndims() - 1; mu++) { - for (size_t nu = mu + 1; nu < U.getndims(); nu++) { - xplusmu[mu] += 1; - xplusnu[nu] += 1; - res += retrace(U(x, mu) * U(xplusmu, nu) * U(xplusnu, mu).dagger() * - U(x, nu).dagger()); - xplusmu[mu] -= 1; - xplusnu[nu] -= 1; - } - } - } - } - } - } - } - // 2n option - anisotropic lattice present - if (anisotropic) { -#pragma omp parallel for reduction(+ : res) - for (size_t x0 = 0; x0 < U.getLt(); x0++) { - for (size_t x1 = 0; x1 < U.getLx(); x1++) { - for (size_t x2 = 0; x2 < U.getLy(); x2++) { - for (size_t x3 = 0; x3 < U.getLz(); x3++) { - std::vector x = {x0, x1, x2, x3}; - std::vector xplusmu = x; - std::vector xplusnu = x; - for (size_t mu = startmu; mu < U.getndims() - 1; mu++) { - for (size_t nu = mu + 1; nu < U.getndims(); nu++) { - xplusmu[mu] += 1; - xplusnu[nu] += 1; - double eta = xi; - if ((mu==0) ^ (nu==0)) { - // at least one direction is temporal (but not both) - eta = 1.0 / xi; - } - res += eta * retrace(U(x, mu) * U(xplusmu, nu) * - U(xplusnu, mu).dagger() * U(x, nu).dagger()); - xplusmu[mu] -= 1; - xplusnu[nu] -= 1; - } - } - } - } - } - } - } - return res; - } - - /** [DEPRECATED - USE retr_sum_Wplaquettes()] - * @brief Wilson plaquette gauge energy - * - * \sum_mu \sum_nu 0 are calculated - * - * @tparam T - * @param U - * @param spatial_only - * @return double - */ - template - double gauge_energy(const gaugeconfig &U, bool spatial_only = false) { - double res = 0.; - size_t startmu = spatial_only; // 0 if spatial_only==false, 1 if spatial_only==true - -#pragma omp parallel for reduction(+ : res) - for (size_t x0 = 0; x0 < U.getLt(); x0++) { - for (size_t x1 = 0; x1 < U.getLx(); x1++) { - for (size_t x2 = 0; x2 < U.getLy(); x2++) { - for (size_t x3 = 0; x3 < U.getLz(); x3++) { - std::vector x = {x0, x1, x2, x3}; - std::vector xplusmu = x; - std::vector xplusnu = x; - for (size_t mu = startmu; mu < U.getndims() - 1; mu++) { - for (size_t nu = mu + 1; nu < U.getndims(); nu++) { - xplusmu[mu] += 1; - xplusnu[nu] += 1; - res += retrace(U(x, mu) * U(xplusmu, nu) * U(xplusnu, mu).dagger() * - U(x, nu).dagger()); - xplusmu[mu] -= 1; - xplusnu[nu] -= 1; - } - } - } - } - } - } - return res; - } - -} // namespace flat_spacetime \ No newline at end of file diff --git a/include/gauge_energy.hh b/include/gauge_energy.hh index cb99a82..0118cc3 100644 --- a/include/gauge_energy.hh +++ b/include/gauge_energy.hh @@ -1,29 +1,139 @@ /** - * @file gauge_energy.hpp - * @author simone-romiti (simone.romiti@uni-bonn.de) - * @brief + * @file flat_spacetime_gauge_energy.hpp + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief gauge energy in flat spacetime (euclidean metric) * @version 0.1 * @date 2022-05-20 * * @copyright Copyright (c) 2022 * */ - #pragma once -#include "metric.hpp" -#include "flat-gauge_energy.hpp" -#include "rotating-gauge_energy.hpp" +#include "gaugeconfig.hh" +#ifdef _USE_OMP_ +#include +#endif + +namespace flat_spacetime { + + /** + * @brief real part of the trace of the Wilson plaquette + * + * \sum_mu \sum_nu 0 are calculated + * @param xi anisotropy in the action: S_G \supset (\beta/xi)*P_{0i} + (xi*\beta)*P_{ij} + * @return double + */ + template + double retr_sum_Wplaquettes(const gaugeconfig &U, + const double &xi = 1.0, + const bool &anisotropic = false, + const bool &spatial_only = false) { + double res = 0.; + size_t startmu = spatial_only; // 0 if spatial_only==false, 1 if spatial_only==true + + if (!anisotropic) { +#pragma omp parallel for reduction(+ : res) + for (size_t x0 = 0; x0 < U.getLt(); x0++) { + for (size_t x1 = 0; x1 < U.getLx(); x1++) { + for (size_t x2 = 0; x2 < U.getLy(); x2++) { + for (size_t x3 = 0; x3 < U.getLz(); x3++) { + std::vector x = {x0, x1, x2, x3}; + std::vector xplusmu = x; + std::vector xplusnu = x; + for (size_t mu = startmu; mu < U.getndims() - 1; mu++) { + for (size_t nu = mu + 1; nu < U.getndims(); nu++) { + xplusmu[mu] += 1; + xplusnu[nu] += 1; + res += retrace(U(x, mu) * U(xplusmu, nu) * U(xplusnu, mu).dagger() * + U(x, nu).dagger()); + xplusmu[mu] -= 1; + xplusnu[nu] -= 1; + } + } + } + } + } + } + } + // 2n option - anisotropic lattice present + if (anisotropic) { +#pragma omp parallel for reduction(+ : res) + for (size_t x0 = 0; x0 < U.getLt(); x0++) { + for (size_t x1 = 0; x1 < U.getLx(); x1++) { + for (size_t x2 = 0; x2 < U.getLy(); x2++) { + for (size_t x3 = 0; x3 < U.getLz(); x3++) { + std::vector x = {x0, x1, x2, x3}; + std::vector xplusmu = x; + std::vector xplusnu = x; + for (size_t mu = startmu; mu < U.getndims() - 1; mu++) { + for (size_t nu = mu + 1; nu < U.getndims(); nu++) { + xplusmu[mu] += 1; + xplusnu[nu] += 1; + double eta = xi; + if ((mu==0) ^ (nu==0)) { + // at least one direction is temporal (but not both) + eta = 1.0 / xi; + } + res += eta * retrace(U(x, mu) * U(xplusmu, nu) * + U(xplusnu, mu).dagger() * U(x, nu).dagger()); + xplusmu[mu] -= 1; + xplusnu[nu] -= 1; + } + } + } + } + } + } + } + return res; + } -template -double gauge_energy_g(const M& g_munu, const gaugeconfig &U, bool spatial_only = false); + /** [DEPRECATED - USE retr_sum_Wplaquettes()] + * @brief Wilson plaquette gauge energy + * + * \sum_mu \sum_nu 0 are calculated + * + * @tparam T + * @param U + * @param spatial_only + * @return double + */ + template + double gauge_energy(const gaugeconfig &U, bool spatial_only = false) { + double res = 0.; + size_t startmu = spatial_only; // 0 if spatial_only==false, 1 if spatial_only==true -template -double gauge_energy_g(const gaugeconfig &U, bool spatial_only = false) { - return flat_spacetime::gauge_energy(U, spatial_only); -} +#pragma omp parallel for reduction(+ : res) + for (size_t x0 = 0; x0 < U.getLt(); x0++) { + for (size_t x1 = 0; x1 < U.getLx(); x1++) { + for (size_t x2 = 0; x2 < U.getLy(); x2++) { + for (size_t x3 = 0; x3 < U.getLz(); x3++) { + std::vector x = {x0, x1, x2, x3}; + std::vector xplusmu = x; + std::vector xplusnu = x; + for (size_t mu = startmu; mu < U.getndims() - 1; mu++) { + for (size_t nu = mu + 1; nu < U.getndims(); nu++) { + xplusmu[mu] += 1; + xplusnu[nu] += 1; + res += retrace(U(x, mu) * U(xplusmu, nu) * U(xplusnu, mu).dagger() * + U(x, nu).dagger()); + xplusmu[mu] -= 1; + xplusnu[nu] -= 1; + } + } + } + } + } + } + return res; + } -template -double gauge_energy(const gaugeconfig &U, bool spatial_only = false) { - return gauge_energy_g(U, spatial_only); -} +} // namespace flat_spacetime \ No newline at end of file diff --git a/include/gauge_energy.hpp b/include/gauge_energy.hpp new file mode 100644 index 0000000..81a3b17 --- /dev/null +++ b/include/gauge_energy.hpp @@ -0,0 +1,35 @@ +/** + * @file gauge_energy.hpp + * @author simone-romiti (simone.romiti@uni-bonn.de) + * @brief + * @version 0.1 + * @date 2022-05-20 + * + * @copyright Copyright (c) 2022 + * + */ + +#pragma once +#include "metric.hpp" + +#include "gauge_energy.hh" +#include "rotating-gauge_energy.hpp" + +// template +// double gauge_energy_g(const M& g_munu, const gaugeconfig &U, bool spatial_only = false); + +// template +// double gauge_energy_g(const gaugeconfig &U, bool spatial_only = false) { +// return flat_spacetime::gauge_energy(U, spatial_only); +// } + +// template +// double gauge_energy(const gaugeconfig &U, bool spatial_only = false) { +// return gauge_energy_g(U, spatial_only); +// } + +// for the moment only the flat spacetime version is supported +template +double gauge_energy(const gaugeconfig &U, bool spatial_only = false) { + return flat_spacetime::gauge_energy(U, spatial_only); +} diff --git a/include/flat-gaugemonomial.hh b/include/gaugemonomial.hh similarity index 85% rename from include/flat-gaugemonomial.hh rename to include/gaugemonomial.hh index d1febf2..7dc234b 100644 --- a/include/flat-gaugemonomial.hh +++ b/include/gaugemonomial.hh @@ -12,7 +12,7 @@ #pragma once #include "adjointfield.hh" -#include "flat-gauge_energy.hpp" +#include "gauge_energy.hpp" #include "gaugeconfig.hh" #include "geometry.hh" #include "get_staples.hh" @@ -33,9 +33,9 @@ namespace flat_spacetime { * Returns: * (\beta/Nc)*\sum_{x} \sum_{\mu<\nu} \eta_{\mu \nu} Re(Tr(P_{\mu \nu}(x))) * with \eta_{0i}=1/\xi and \eta_{ij}=\xi - * + * * where beta = 2*N_c/g_0^2 - * + * * @tparam Float * @tparam Group * @param U gauge configuration @@ -59,7 +59,7 @@ namespace flat_spacetime { * * S_G = (\beta/Nc)*\sum_{x} \sum_{\mu<\nu} * \eta_{\mu \nu} Re(Tr(1 - P_{\mu\nu}(x))) - * + * * where beta = 2*N_c/g_0^2 * * @tparam Float @@ -102,6 +102,16 @@ namespace flat_spacetime { flat_spacetime::get_S_G_hmc(*h.U, (*this).xi, (*this).anisotropic); return; } + + /** + * @brief updating the derivative of the action with respect to the gauge field + * contribution in this monomial + * + * @param deriv derivative to be updated + * @param h hamiltonian field (Hybrid Monte Carlo evolution) + * @param fac multiplicative facator of the derivative (1.0 for the HMC but != 1 in + * other cases e.g. the gradient flow evolution) + */ void derivative(adjointfield &deriv, hamiltonian_field const &h, const Float fac = 1.) const override { @@ -115,13 +125,10 @@ namespace flat_spacetime { for (size_t mu = 0; mu < h.U->getndims(); mu++) { accum S; get_staples(S, *h.U, x, mu, (*this).xi, (*this).anisotropic); - S = (*h.U)(x, mu) * S; - // the antihermitian traceless part - // beta/N_c *(U*U^stap - (U*U^stap)^dagger) - // in get_deriv + S = (*h.U)(x, mu) * S; // U*A in eq. 8.40 in Gattringer&Lang - deriv(x, mu) += - fac * h.U->getBeta() / double(h.U->getNc()) * get_deriv(S); + const double num_fact_i = fac * h.U->getBeta() / double(h.U->getNc()); + deriv(x, mu) += num_fact_i * get_deriv(S); } } } diff --git a/include/glueballs.hpp b/include/glueballs.hpp index 063bfa8..d1bb247 100644 --- a/include/glueballs.hpp +++ b/include/glueballs.hpp @@ -49,6 +49,15 @@ namespace glueballs { L *= operators::parity(Px, 0, U, x, nu).dagger(); } + if (Px == true) { + // the loop operator L(x) is transformed by: + // 1. symmetrizing the loop (keeping the orientation) + // 2. daggering the result. + // What computed so far is L.dagger(), so we dagger it once more to get the correct + // result + L = dagger(L); + } + return trace(L); } @@ -68,7 +77,6 @@ namespace glueballs { const size_t &nu, const bool &Px) { typedef typename accum_type::type accum; - accum L = U(x, mu) * (U(x, mu).dagger()); // "1", independently of the group // start at x // go to x + (a+1)*\hat{mu} @@ -103,6 +111,15 @@ namespace glueballs { } // loop closed + if (Px == true) { + // the loop operator L(x) is transformed by: + // 1. symmetrizing the loop (keeping the orientation) + // 2. daggering the result. + // What computed so far is L.dagger(), so we dagger it once more to get the correct + // result + L = dagger(L); + } + return trace(L); } diff --git a/include/flat-gradient_flow.hh b/include/gradient_flow.hh similarity index 91% rename from include/flat-gradient_flow.hh rename to include/gradient_flow.hh index c0c125d..6b10f9f 100644 --- a/include/flat-gradient_flow.hh +++ b/include/gradient_flow.hh @@ -1,10 +1,10 @@ #pragma once -//#include"gradient_flow.hh" +// #include"gradient_flow.hh" #include "adjointfield.hh" #include "flat-energy_density.hh" -#include "flat-gauge_energy.hpp" -#include "flat-gaugemonomial.hh" +#include "gauge_energy.hpp" #include "gaugeconfig.hh" +#include "gaugemonomial.hh" #include "hamiltonian_field.hh" #include "monomial.hh" #include "su2.hh" @@ -32,17 +32,13 @@ void runge_kutta(hamiltonian_field &h, // before the three steps zero the derivative field zeroadjointfield(*(h.momenta)); + const double fact_f = 2.0 * h.Fact_Nc_force_Luscher * h.U->getNc() / h.U->getBeta(); for (int f = 0; f < 3; f++) { - // add to *(h.momenta) + // evolution of eq. C.2 of https://arxiv.org/pdf/1006.4518.pdf // we have to cancel beta/N_c from the derivative - // a factor two to obtain the correct normalisation of - // the Wilson plaquette action - // S_W = 1./g_0^2 \sum_x \sum_{p} Re Tr(1 - U(p)) - // where we sum over all oriented plaquettes - // we sum over unoriented plaquettes, so we have to multiply by 2 - // which is usually in beta - SW.derivative(*(h.momenta), h, 2. * h.U->getNc() * zfac[f] / h.U->getBeta()); - // The '-' comes from the action to be tr(1-U(p)) + + // updating the momenta *(h.momenta) + SW.derivative(*(h.momenta), h, fact_f * zfac[f]); // update the flowed gauge field Vt update_gauge(h, -eps * expfac[f]); } diff --git a/include/hamiltonian_field.hh b/include/hamiltonian_field.hh index aa2546f..300d758 100644 --- a/include/hamiltonian_field.hh +++ b/include/hamiltonian_field.hh @@ -1,12 +1,22 @@ #pragma once -#include"gaugeconfig.hh" -#include"adjointfield.hh" -#include - -template struct hamiltonian_field { - adjointfield * momenta; - gaugeconfig * U; - hamiltonian_field(adjointfield &momenta, gaugeconfig &U) : - momenta(&momenta), U(&U) {} +#include "adjointfield.hh" +#include "gaugeconfig.hh" +#include + +template struct hamiltonian_field { + adjointfield *momenta; + gaugeconfig *U; + + // Our normalization convention of the generators is different than + // https://arxiv.org/pdf/1006.4518.pdf this factor takes that into account when using + // this class for the gradiend flow evolution of the gauge fields. + double Fact_Nc_force_Luscher = 1.0 / 4.0; // default: non abelian case + + hamiltonian_field(adjointfield &momenta, gaugeconfig &U) + : momenta(&momenta), U(&U) { + if (U.getNc() == 1) { + Fact_Nc_force_Luscher = 1.0 / 2.0; + } + } }; diff --git a/include/monomial.hh b/include/monomial.hh index 2e489f7..ca6d98b 100644 --- a/include/monomial.hh +++ b/include/monomial.hh @@ -46,4 +46,4 @@ public: const Float fac) const {} }; -#include "flat-gaugemonomial.hh" +#include "gaugemonomial.hh" diff --git a/include/omeasurements.hpp b/include/omeasurements.hpp index 454ec2a..8a79fef 100644 --- a/include/omeasurements.hpp +++ b/include/omeasurements.hpp @@ -9,10 +9,11 @@ #include #include +#include #include #include -#include "flat-gradient_flow.hh" +#include "gradient_flow.hh" #include "glueballs.hpp" #include "io.hh" #include "links.hpp" @@ -22,7 +23,10 @@ #include "smearape.hh" #include "wilsonloop.hh" #include + #include +#include +#include namespace omeasurements { @@ -142,7 +146,7 @@ namespace omeasurements { */ template std::vector> - meas_glueball_interpolators_pp_mm(const std::string &type, + meas_glueball_interpolators(const std::string &type, const gaugeconfig &U, const size_t &i, const size_t &nAPEsmear, @@ -160,7 +164,7 @@ namespace omeasurements { if (S.glueball.doAPEsmear) { oss_dir << meas_details + "/"; } - fsys::create_directories(fsys::absolute(oss_dir.str())); + // fsys::create_directories(fsys::absolute(oss_dir.str())); if (S.glueball.doAPEsmear) { if (S.glueball.lengthy_file_name) { @@ -186,15 +190,7 @@ namespace omeasurements { for (size_t i1 = 0; i1 < nr; i1++) { const size_t r1 = i1 + rmin; - // directory path - const std::string dir_ij = - oss_dir.str() + std::to_string(r1) + "/"; - fsys::create_directories(fsys::absolute(dir_ij)); // creating directory - - const std::string path = - dir_ij + "phi" + oss_name.str(); // full path of output file - - phi[ii].resize({T_ext, 3}); // t, ++, -- + phi[ii].resize({T_ext, 5}); // t, ++, +-, -+, -- for (size_t t = 0; t < T_ext; t++) { phi[ii](t, 0) = t; @@ -203,13 +199,20 @@ namespace omeasurements { const std::complex Pm = glueballs::get_rest_trace_loop( type, t, U, r1, true, S.glueball.spatial_loops); phi[ii](t, 1) = (Pp + Pm).real() / 2.0; // PC=++ - phi[ii](t, 2) = (Pp + Pm).imag() / 2.0; // PC=-- + phi[ii](t, 3) = (Pp + Pm).imag() / 2.0; // PC=+- + phi[ii](t, 2) = (Pp - Pm).real() / 2.0; // PC=-+ + phi[ii](t, 4) = (Pp - Pm).imag() / 2.0; // PC=-- } if (save_interpolator) { + // directory path + const std::string dir_ij = oss_dir.str() + std::to_string(r1) + "/"; + fsys::create_directories(fsys::absolute(dir_ij)); + // full path of output file + const std::string path = dir_ij + "phi" + oss_name.str(); std::ofstream ofs(path, std::ios::out); ofs << std::scientific << std::setprecision(16); - io::xtensor_to_stream(ofs, phi[ii], " ", "t pp mm"); + io::xtensor_to_stream(ofs, phi[ii], " ", "t pp pm mp mm"); ofs.close(); } ++ii; @@ -228,7 +231,7 @@ namespace omeasurements { * @param S specific parameters */ template - void meas_glueball_correlator_pp_mm(const std::string &type, + void meas_glueball_correlator(const std::string &type, const gaugeconfig &U, const size_t &i, const size_t &nAPEsmear, @@ -243,7 +246,7 @@ namespace omeasurements { if (S.glueball.doAPEsmear) { oss_dir << meas_details + "/"; } - fsys::create_directories(fsys::absolute(oss_dir.str())); + // fsys::create_directories(fsys::absolute(oss_dir.str())); if (S.glueball.doAPEsmear) { if (S.glueball.lengthy_file_name) { @@ -262,8 +265,11 @@ namespace omeasurements { const size_t N_ops = rmax; // number of interpolatoing operators // phi_i(t)^{PC} - std::vector> phi = meas_glueball_interpolators_pp_mm( + std::vector> phi = meas_glueball_interpolators( type, U, i, nAPEsmear, S.glueball.save_interpolator, S); + if (!S.glueball.correlator) { + return; + } const size_t n_phi = phi.size(); @@ -277,14 +283,15 @@ namespace omeasurements { // full path of output file const std::string path = dir_ij + "C_glueball" + oss_name.str(); - const std::string header = "t pp mm"; + const std::string header = "t pp pm mp mm"; xt::xarray corr; - corr.resize({T_ext, 3}); // t, ++, -- + corr.resize({T_ext, 5}); // t, ++, +-, -+, -- for (size_t t = 0; t < T_ext; t++) { corr(t, 0) = t; size_t iPC = 1; // 1st column is time 't' - for (size_t PC = 0; PC <= 1; PC++) { + for (size_t P = 0; P <= 1; P++) { + for (size_t C = 0; C <= 1; C++) { corr(t, iPC) = 0.0; for (size_t tau = 0; tau < T_ext; tau++) { const size_t t1 = (t + tau) % T_ext; @@ -292,6 +299,7 @@ namespace omeasurements { } corr(t, iPC) /= double(T_ext); // average over all times ++iPC; + } } } @@ -319,7 +327,7 @@ namespace omeasurements { * @param S specific parameters */ template - void meas_glueball_interpolators_pp_mm(const std::string &type, + void meas_glueball_interpolators(const std::string &type, const gaugeconfig &U0, const size_t &i, const sparams &S) { @@ -334,7 +342,7 @@ namespace omeasurements { // other nsteps so that we reach the value of v_ns[is] spatial_APEsmearing(U, S.glueball.alphaAPEsmear); } - auto foo = from_smeared_field::meas_glueball_interpolators_pp_mm( + auto foo = from_smeared_field::meas_glueball_interpolators( type, U, i, v_ns[is], S.glueball.save_interpolator, S); } } @@ -365,7 +373,7 @@ namespace omeasurements { // other nsteps so that we reach the value of v_ns[is] spatial_APEsmearing(U, S.glueball.alphaAPEsmear); } - from_smeared_field::meas_glueball_correlator_pp_mm(type, U, i, v_ns[is], S); + from_smeared_field::meas_glueball_correlator(type, U, i, v_ns[is], S); } } diff --git a/include/operators.hpp b/include/operators.hpp index 3c37ac1..87956c7 100644 --- a/include/operators.hpp +++ b/include/operators.hpp @@ -12,7 +12,7 @@ #include #include -#include "flat-gradient_flow.hh" +#include "gradient_flow.hh" #include "gaugeconfig.hh" #include "links.hpp" #include "parameters.hh" diff --git a/include/parameters.hh b/include/parameters.hh index 35c2aff..1287d06 100644 --- a/include/parameters.hh +++ b/include/parameters.hh @@ -54,7 +54,7 @@ namespace global_parameters { // glueball interpolators are all loops of links of a certain size and shape bool interpolator = false; // true if measuring the glueball interpolators - std::string interpolator_type = "NONE"; // "squares", "rectangles", "fatL", etc. + std::string interpolator_type = "NONE"; // "rectangles" or "fatL" bool save_interpolator = false; // true if saving interpolators values bool spatial_loops = true; // true if using only spatial loops size_t rmin = 0, rmax = 0; // minum and maximum sizes of the interpolating loops diff --git a/include/staggered.hpp b/include/staggered.hpp index a8983f6..f6da441 100644 --- a/include/staggered.hpp +++ b/include/staggered.hpp @@ -16,7 +16,7 @@ #include #include "adjointfield.hh" -#include "flat-gauge_energy.hpp" +#include "gauge_energy.hpp" #include "gaugeconfig.hh" #include "geometry.hh" #include "get_staples.hh" diff --git a/include/su2.hh b/include/su2.hh index c48fcfd..c340ea6 100644 --- a/include/su2.hh +++ b/include/su2.hh @@ -1,43 +1,48 @@ -#pragma once +#pragma once -#include"accum_type.hh" -#include"traceless_antiherm.hh" -#include -//#include +#include +// #include +#include -using Complex = std::complex; +#include "accum_type.hh" +#include "dagger.hh" +#include "traceless_antiherm.hh" +using Complex = std::complex; +/** + * @brief representation of an SU(2) matrix in the fundamental representation with the + * convention of eq. (4.23) of Gattringer&Lang + * https://link.springer.com/book/10.1007/978-3-642-01850-3 + * + */ class _su2 { public: const size_t N_c = 2; explicit _su2() : a(0), b(0) {} explicit _su2(Complex a, Complex b) : a(a), b(b) {} - _su2(const _su2& U) : a(U.a), b(U.b) {} + _su2(const _su2 &U) : a(U.a), b(U.b) {} friend inline _su2 operator+(const _su2 &U1, const _su2 &U2); friend inline _su2 operator-(const _su2 &U1, const _su2 &U2); friend inline _su2 operator*(const _su2 &U1, const _su2 &U2); friend inline _su2 operator*(const Complex &U1, const _su2 &U2); friend inline _su2 operator*(const _su2 &U1, const Complex &U2); - _su2& operator*=(const _su2 &U1) { + _su2 &operator*=(const _su2 &U1) { Complex a = this->a; - this->a = a*U1.a - this->b*std::conj(U1.b); - this->b = a*U1.b + this->b*std::conj(U1.a); + this->a = a * U1.a - this->b * std::conj(U1.b); + this->b = a * U1.b + this->b * std::conj(U1.a); return *this; } _su2 round(size_t n) const { double dn = n; - return _su2(Complex(std::round(std::real(a)*dn), std::round(std::imag(a)*dn))/dn, - Complex(std::round(std::real(b)*dn), std::round(std::imag(b)*dn))/dn); + return _su2( + Complex(std::round(std::real(a) * dn), std::round(std::imag(a) * dn)) / dn, + Complex(std::round(std::real(b) * dn), std::round(std::imag(b) * dn)) / dn); } - inline Complex geta() const { - return(a); - } - inline Complex getb() const { - return(b); - } + inline Complex geta() const { return (a); } + inline Complex getb() const { return (b); } inline void operator=(const _su2 &U) { a = U.geta(); b = U.getb(); @@ -50,17 +55,11 @@ public: a = _a; b = _b; } - inline _su2 dagger() const { - return(_su2(std::conj(a), -b)); - } - inline double retrace() { - return(2.*std::real(a)); - } - Complex det() { - return(a*std::conj(a) + b*std::conj(b)); - } + inline _su2 dagger() const { return (_su2(std::conj(a), -b)); } + inline double retrace() { return (2. * std::real(a)); } + Complex det() { return (a * std::conj(a) + b * std::conj(b)); } void restoreSU() { - double r = sqrt(std::abs(a)*std::abs(a) + std::abs(b)*std::abs(b)); + double r = sqrt(std::abs(a) * std::abs(a) + std::abs(b) * std::abs(b)); a /= r; b /= r; } @@ -71,54 +70,58 @@ private: inline double retrace(_su2 const &U) { double a = std::real(U.geta()); - return(2*a); + return (2 * a); } inline Complex trace(_su2 const &U) { double a = std::real(U.geta()); - return(Complex(2*a, 0.)); + return (Complex(2 * a, 0.)); } +template <> inline _su2 dagger(const _su2 &u) { + const Complex a = u.geta(); + const Complex b = u.getb(); + _su2 ud(std::conj(a), -b); + return ud; +} -template<> inline _su2 traceless_antiherm(const _su2& x) { - return(_su2(0.5*(x.geta()-std::conj(x.geta())), x.getb())); +template <> inline _su2 traceless_antiherm(const _su2 &x) { + return (_su2(0.5 * (x.geta() - std::conj(x.geta())), x.getb())); } inline _su2 operator*(const _su2 &U1, const _su2 &U2) { _su2 res; - res.a = U1.a*U2.a - U1.b*std::conj(U2.b); - res.b = U1.a*U2.b + U1.b*std::conj(U2.a); - return(res); + res.a = U1.a * U2.a - U1.b * std::conj(U2.b); + res.b = U1.a * U2.b + U1.b * std::conj(U2.a); + return (res); } inline _su2 operator+(const _su2 &U1, const _su2 &U2) { _su2 res; res.a = U1.a + U2.a; res.b = U1.b + U2.b; - return(res); + return (res); } inline _su2 operator-(const _su2 &U1, const _su2 &U2) { _su2 res; res.a = U1.a - U2.a; res.b = U1.b - U2.b; - return(res); + return (res); } inline _su2 operator*(const Complex &U1, const _su2 &U2) { - _su2 res; - res.a = U2.a * U1; - res.b = U2.b * U1; - return(res); + _su2 res; + res.a = U2.a * U1; + res.b = U2.b * U1; + return (res); } inline _su2 operator*(const _su2 &U1, const Complex &U2) { - _su2 res; - res.a = U1.a * U2; - res.b = U1.b * U2; - return(res); + _su2 res; + res.a = U1.a * U2; + res.b = U1.b * U2; + return (res); } - using su2 = _su2; - diff --git a/include/u1.hh b/include/u1.hh index c6b4d89..98ea897 100644 --- a/include/u1.hh +++ b/include/u1.hh @@ -1,10 +1,12 @@ -#pragma once +#pragma once -#include"accum_type.hh" -#include"traceless_antiherm.hh" -#include -#include -#include +#include +#include +#include + +#include "accum_type.hh" +#include "dagger.hh" +#include "traceless_antiherm.hh" using Complex = std::complex; @@ -13,78 +15,67 @@ public: const size_t N_c = 1; explicit _u1() : a(0) {} explicit _u1(double _a) : a(_a) {} - _u1(const _u1& U) : a(U.a) {} + _u1(const _u1 &U) : a(U.a) {} _u1(Complex _a) : a(std::arg(_a)) {} friend Complex operator+(const _u1 &U1, const _u1 &U2); friend Complex operator-(const _u1 &U1, const _u1 &U2); friend _u1 operator*(const _u1 &U1, const _u1 &U2); - + // implicit conversion operator to complex - operator Complex() const { - return(std::exp(a*Complex(0., 1.))); - } - _u1& operator*=(const _u1 &U1) { + operator Complex() const { return (std::exp(a * Complex(0., 1.))); } + _u1 &operator*=(const _u1 &U1) { this->a += U1.a; return *this; } _u1 round(size_t n) const { double dn = n; - return _u1(double(std::round(a)*dn)/dn); + return _u1(double(std::round(a) * dn) / dn); } - double geta() const { - return(a); - } - void operator=(const _u1 &U) { - a = U.a; - } - void set(const double _a) { - a = _a; - } - _u1 dagger() const { - return(_u1(-a)); - } - double retrace() const { - return(std::cos(a)); - } - Complex det() const { - return(std::exp(a*Complex(0., 1.))); - } - void restoreSU() { - } + double geta() const { return (a); } + void operator=(const _u1 &U) { a = U.a; } + void set(const double _a) { a = _a; } + _u1 dagger() const { return (_u1(-a)); } + double retrace() const { return (std::cos(a)); } + Complex det() const { return (std::exp(a * Complex(0., 1.))); } + void restoreSU() {} private: double a; }; inline double retrace(_u1 const &U) { - return(std::cos(U.geta())); + return (std::cos(U.geta())); } inline double retrace(const Complex c) { - return(std::real(c)); + return (std::real(c)); } inline Complex trace(const Complex c) { - return(c); // for U(1) the trace operator acts trivially + return (c); // for U(1) the trace operator acts trivially } -template<> struct accum_type<_u1> { +template <> struct accum_type<_u1> { typedef Complex type; }; -template<> inline Complex traceless_antiherm(const Complex& x) { - return(Complex(0., std::imag(x))); +template <> inline Complex dagger(const Complex &x) { + return std::conj(x); +} + +template <> inline Complex traceless_antiherm(const Complex &x) { + return (Complex(0., std::imag(x))); } //_u1 operator*(const _u1 &U1, const _u1 &U2); -//Complex operator*(const _u1 &U1, const Complex &U2); -//Complex operator*(const Complex &U1, const _u1 &U2); -//Complex operator+(const _u1 &U1, const _u1 &U2); -//Complex operator-(const _u1 &U1, const _u1 &U2); -//void operator+=(Complex & U1, const _u1 & U2); -//void operator*=(Complex & U1, const _u1 & U2); +// Complex operator*(const _u1 &U1, const Complex &U2); +// Complex operator*(const Complex &U1, const _u1 &U2); +// Complex operator+(const _u1 &U1, const _u1 &U2); +// Complex operator-(const _u1 &U1, const _u1 &U2); +// void operator+=(Complex & U1, const _u1 & U2); +// void operator*=(Complex & U1, const _u1 & U2); inline _u1 operator*(const _u1 &U1, const _u1 &U2) { // _u1 res; @@ -93,27 +84,24 @@ inline _u1 operator*(const _u1 &U1, const _u1 &U2) { } inline Complex operator*(const _u1 &U1, const Complex &U2) { - return(std::exp(U1.geta()*Complex(0., 1.)) * U2); + return (std::exp(U1.geta() * Complex(0., 1.)) * U2); } inline Complex operator*(const Complex &U1, const _u1 &U2) { - return(U1 * std::exp(U2.geta()*Complex(0., 1.))); + return (U1 * std::exp(U2.geta() * Complex(0., 1.))); } - inline Complex operator+(const _u1 &U1, const _u1 &U2) { - return(std::exp(U1.a*Complex(0., 1.)) + - std::exp(U2.a*Complex(0., 1.))); + return (std::exp(U1.a * Complex(0., 1.)) + std::exp(U2.a * Complex(0., 1.))); } inline Complex operator-(const _u1 &U1, const _u1 &U2) { - return(std::exp(U1.a*Complex(0., 1.)) - - std::exp(U2.a*Complex(0., 1.))); + return (std::exp(U1.a * Complex(0., 1.)) - std::exp(U2.a * Complex(0., 1.))); } -inline void operator+=(Complex & U1, const _u1 & U2) { - U1 += std::exp(U2.geta()*Complex(0., 1.)); +inline void operator+=(Complex &U1, const _u1 &U2) { + U1 += std::exp(U2.geta() * Complex(0., 1.)); } -inline void operator*=(Complex & U1, const _u1 & U2) { - U1 *= std::exp(U2.geta()*Complex(0., 1.)); +inline void operator*=(Complex &U1, const _u1 &U2) { + U1 *= std::exp(U2.geta() * Complex(0., 1.)); } diff --git a/kramers.cc b/kramers.cc index e368d87..7692aa7 100644 --- a/kramers.cc +++ b/kramers.cc @@ -10,7 +10,7 @@ * */ -#include "flat-gauge_energy.hpp" +#include "gauge_energy.hh" #include "gaugeconfig.hh" #include "geometry.hh" #include "integrator.hh" diff --git a/parse_input_file.cc b/parse_input_file.cc index 02213ba..2bb3615 100644 --- a/parse_input_file.cc +++ b/parse_input_file.cc @@ -201,7 +201,9 @@ namespace input_file_parsing { std::abort(); } - mgparams.do_measure = true; + in.read_verb(mgparams.correlator, {"correlator"}); + mgparams.do_measure = + mgparams.interpolator || mgparams.correlator; // true if we save one of the two in.read_verb(mgparams.doAPEsmear, {"do_APE_smearing"}); if (mgparams.doAPEsmear) { in.read_sequence_verb(mgparams.vec_nAPEsmear, {"APE_smearing", "n"}); diff --git a/scaling.cc b/scaling.cc index 6c5e4e4..e026ee2 100644 --- a/scaling.cc +++ b/scaling.cc @@ -1,7 +1,7 @@ #include"su2.hh" #include"u1.hh" #include"gaugeconfig.hh" -#include"flat-gauge_energy.hpp" +#include"gauge_energy.hh" #include"random_gauge_trafo.hh" #include"flat-sweep.hh" #include"parse_commandline.hh" diff --git a/test.cc b/test.cc index b9cdfe0..1c0318a 100644 --- a/test.cc +++ b/test.cc @@ -2,7 +2,7 @@ #include"u1.hh" #include"random_element.hh" #include"gaugeconfig.hh" -#include"flat-gauge_energy.hpp" +#include"gauge_energy.hpp" #include"random_gauge_trafo.hh" #include"parse_commandline.hh" #include"flat-energy_density.hh" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt deleted file mode 100644 index aabe76f..0000000 --- a/test/CMakeLists.txt +++ /dev/null @@ -1,40 +0,0 @@ -cmake_minimum_required(VERSION 3.10) - -# set the project name -project(test VERSION 1.0) - -# specify the C++ standard -set(CMAKE_CXX_STANDARD 11) -set(CMAKE_CXX_STANDARD_REQUIRED True) - -# add the executables - -# conjugate gradient algorithms -include_directories(PUBLIC ../include/) - -add_executable(CG.exe CG.cpp) -add_executable(BiCGStab.exe BiCGStab.cpp) -add_executable(smearing.exe smearing.cpp) - -#yaml input file parsing - -set(yaml_src_dir ${YAML_SRC}) -set( - inc_dir - ${yaml_src_dir}/include/ - ${yaml_src_dir}/include/yaml-cpp/ - ) - -add_executable(yaml.exe yaml.cpp) -target_include_directories(yaml.exe PUBLIC ${inc_dir}) - -set( - link_dir - ${yaml_src_dir}/build/ - ${yaml_src_dir}/build/lib/ -) - - -target_link_directories(yaml.exe PUBLIC ${link_dir}) -target_link_libraries(yaml.exe libyaml-cpp.a) - diff --git a/test/Makefile b/test/Makefile deleted file mode 100644 index f2b54c5..0000000 --- a/test/Makefile +++ /dev/null @@ -1,18 +0,0 @@ -# Makefile - -all: CG_programs yaml links - -CG_programs: CG BiCGStab - -CG: - g++ CG.cpp -o CG.exe - - -BiCGStab: - g++ BiCGStab.cpp -o BiCGStab.exe - -links: - g++ links.cpp -o links.exe - -clean: - rm *.exe diff --git a/test/h5-dump.cpp b/test/h5-dump.cpp new file mode 100644 index 0000000..c72f1e7 --- /dev/null +++ b/test/h5-dump.cpp @@ -0,0 +1,38 @@ +// h5-dump.cpp +// reading .h5 file and dumping its content in the working directory + +#include +#include +#include +#include +#include +#include + +#include "boost/filesystem.hpp" +#include "boost/range/iterator_range.hpp" +#include + + +int main(int argc, char **argv) { + // check that argc == 3 + + H5Easy::File file(argv[1], H5Easy::File::ReadWrite); + std::vector path = H5Easy::load>(file, "path"); + + const int n = path.size(); + + for (std::string p : path) { + std::ofstream ofs(p); + + if (!boost::filesystem::exists(p)) { + std::vector strs; + boost::split(strs, p, boost::is_any_of("/")); + boost::filesystem::create_directory(strs[0]); + } + + const std::string result = H5Easy::load(file, p); + ofs << result; + } + + return 0; +} \ No newline at end of file diff --git a/test/h5-gen.cpp b/test/h5-gen.cpp new file mode 100644 index 0000000..a5bc87a --- /dev/null +++ b/test/h5-gen.cpp @@ -0,0 +1,57 @@ +// h5-gen.cpp +// generate .h5 file from directory path + +#include +#include +#include +#include +#include + +#include "boost/filesystem.hpp" +#include "boost/range/iterator_range.hpp" + +#include + +std::vector ls_recursive(const std::string &path) { + boost::filesystem::path dir = path; + boost::filesystem::recursive_directory_iterator it(dir), end; + + std::vector files; + for (auto &entry : boost::make_iterator_range(it, end)) { + if (boost::filesystem::is_regular(entry)) { + files.push_back(entry.path().native()); + } + } + + return files; +} + +int main(int argc, char **argv) { + // check that argc == 3 + + // list all files in directory argv[1] + + std::string h5_file = argv[2]; + std::vector path; // vector of paths + + if (boost::filesystem::exists(h5_file)) { + std::cerr << "Error. File already exists: " << h5_file << "\n"; + std::abort(); + } + + H5Easy::File file2(h5_file, H5Easy::File::Overwrite); + + const std::vector files = ls_recursive(argv[1]); + for (std::string f : files) { + path.push_back(f); + + std::string cnt; + boost::filesystem::load_string_file(f, cnt); + H5Easy::dump(file2, f, cnt); // dumping vector of paths + + } + + H5Easy::dump(file2, "path", path); // dumping vector of paths + + return 0; +} \ No newline at end of file diff --git a/test/init_test.sh b/test/init_test.sh deleted file mode 100755 index b5b2ff2..0000000 --- a/test/init_test.sh +++ /dev/null @@ -1,15 +0,0 @@ -# init_test.sh - -cd .. -d1=$(pwd -P) -cd test - -YAML_SRC_PATH=$d1/external/yaml-cpp/ - -rm -rf build # removing the previous build/ directory -mkdir build # creating a new one - -cp config.yaml build - -cd build -cmake ../ -D YAML_SRC=$YAML_SRC_PATH