From ccc67c46c2f535e9d6cd352eeb43b8cd348f35b7 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Wed, 15 Feb 2023 17:08:47 +0100 Subject: [PATCH 01/13] DO NOT MERGE implementing support for SU(3) fields --- include/adjoint_su2.hh | 68 +++++++++++++++++++ include/adjoint_su3.hh | 79 ++++++++++++++++++++++ include/adjoint_u1.hh | 47 +++++++++++++ include/adjointfield.hh | 98 +++++---------------------- include/su3.hh | 146 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 358 insertions(+), 80 deletions(-) create mode 100644 include/adjoint_su2.hh create mode 100644 include/adjoint_su3.hh create mode 100644 include/adjoint_u1.hh create mode 100644 include/su3.hh diff --git a/include/adjoint_su2.hh b/include/adjoint_su2.hh new file mode 100644 index 0000000..92d1e23 --- /dev/null +++ b/include/adjoint_su2.hh @@ -0,0 +1,68 @@ +/** + * @file adjoint_su2.hh + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief class and routines for an su(2) matrix in the adjoint representation + * @version 0.1 + * @date 2023-02-15 + * + * @copyright Copyright (c) 2023 + * + */ + +#pragma once + +#include "geometry.hh" +#include "su2.hh" + +#include +#include +#include +#include +#include + + +template class adjointsu2 { +public: + adjointsu2(Float _a, Float _b, Float _c) : a(_a), b(_b), c(_c) {} + adjointsu2() : a(0.), b(0.), c(0.) {} + void flipsign() { + a = -a; + b = -b; + c = -c; + } + Float geta() const { return a; } + Float getb() const { return b; } + Float getc() const { return c; } + void seta(Float _a) { a = _a; } + void setb(Float _a) { b = _a; } + void setc(Float _a) { c = _a; } + void setzero() { a = b = c = 0.; } + adjointsu2 round(size_t n) const { + Float dn = n; + return adjointsu2(std::round(a * dn) / dn, std::round(b * dn) / dn, + std::round(c * dn) / dn); + } + void operator=(const adjointsu2 &A) { + a = A.geta(); + b = A.getb(); + c = A.getc(); + } + void operator+=(const adjointsu2 &A) { + a += A.geta(); + b += A.getb(); + c += A.getc(); + } + void operator-=(const adjointsu2 &A) { + a -= A.geta(); + b -= A.getb(); + c -= A.getc(); + } + +private: + Float a, b, c; +}; + +template inline adjointsu2 get_deriv(su2 &A) { + const Complex a = A.geta(), b = A.getb(); + return adjointsu2(2. * std::imag(b), 2. * std::real(b), 2. * std::imag(a)); +} diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh new file mode 100644 index 0000000..705b0c9 --- /dev/null +++ b/include/adjoint_su3.hh @@ -0,0 +1,79 @@ +/** + * @file adjoint_su3.hh + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief class and routines for an su(3) matrix in the adjoint representation + * @version 0.1 + * @date 2023-02-15 + * + * @copyright Copyright (c) 2023 + * + */ + +#pragma once + +#include "geometry.hh" +#include "su3.hh" + +#include +#include +#include +#include +#include + +/** + * @brief element of the su(3) algebra + * + * the element is parametrized by 8 real numbers (number of generators for SU(3)) + * + * @tparam Float + */ +template class adjointsu3 { +public: + adjointsu3(const std::array &_arr) : arr(_arr) {} + adjointsu3() { this->setzero(); } + void flipsign() { + for (size_t i = 0; i < 8; i++) { + arr[i] = -arr[i]; + } + } + std::array get_arr() const { return arr; } + void set_arr(const std::array &_arr) { arr = _arr; } + void setzero() { + const std::array _arr{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + arr = _arr; + } + adjointsu3 round(const size_t &n) const { + Float dn = n; + std::array arr2; + for (size_t i = 0; i < 8; i++) { + arr2[i] = std::round(arr[i] * dn) / dn; + } + return adjointsu3(arr2); + } + + void operator=(const adjointsu3 &A) { arr = A.get_arr(); } + void operator+=(const adjointsu3 &A) { + const std::array arr_A = A.get_arr(); + for (size_t i = 0; i < 8; i++) { + arr[i] += arr_A[i]; + } + } + void operator-=(const adjointsu3 &A) { + const std::array arr_A = A.get_arr(); + for (size_t i = 0; i < 8; i++) { + arr[i] -= arr_A[i]; + } + } + +private: + std::array arr; +}; + +template inline adjointsu3 get_deriv(const su3 &A) { + const std::array, 9> arr_A = A.get_arr(); + const std::array arr; + for (size_t i = 0; i < 8; i++) { + arr[i] = 2.0 * std::imag(arr_A[i+1]); + } + return adjointsu3(arr); +} diff --git a/include/adjoint_u1.hh b/include/adjoint_u1.hh new file mode 100644 index 0000000..8789a05 --- /dev/null +++ b/include/adjoint_u1.hh @@ -0,0 +1,47 @@ +/** + * @file adjoint_u1.hh + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief class and routines for an u(1) matrix in the adjoint representation + * (constant*Identity) + * @version 0.1 + * @date 2023-02-15 + * + * @copyright Copyright (c) 2023 + * + */ + +#pragma once + +#include "geometry.hh" +#include "u1.hh" + +#include +#include +#include +#include +#include + +template class adjointu1 { +public: + adjointu1(Float _a) : a(_a) {} + adjointu1() : a(0.) {} + void flipsign() { a = -a; } + Float geta() const { return a; } + void seta(Float _a) { a = _a; } + void setzero() { a = 0.; } + + adjointu1 round(size_t n) const { + Float dn = n; + return adjointu1(std::round(a * dn) / dn); + } + void operator=(const adjointu1 &A) { a = A.geta(); } + void operator+=(const adjointu1 &A) { a += A.geta(); } + void operator-=(const adjointu1 &A) { a -= A.geta(); } + +private: + Float a; +}; + +template inline adjointu1 get_deriv(Complex &A) { + return adjointu1(std::imag(A)); +} diff --git a/include/adjointfield.hh b/include/adjointfield.hh index 2c5c717..1de3f82 100644 --- a/include/adjointfield.hh +++ b/include/adjointfield.hh @@ -1,88 +1,21 @@ #pragma once #include "geometry.hh" -#include "su2.hh" -#include "u1.hh" +#include #include #include #include #include -#include - -template class adjointsu2 { -public: - adjointsu2(Float _a, Float _b, Float _c) : a(_a), b(_b), c(_c) {} - adjointsu2() : a(0.), b(0.), c(0.) {} - void flipsign() { - a = -a; - b = -b; - c = -c; - } - Float geta() const { return a; } - Float getb() const { return b; } - Float getc() const { return c; } - void seta(Float _a) { a = _a; } - void setb(Float _a) { b = _a; } - void setc(Float _a) { c = _a; } - void setzero() { a = b = c = 0.; } - adjointsu2 round(size_t n) const { - Float dn = n; - return adjointsu2(std::round(a * dn) / dn, std::round(b * dn) / dn, - std::round(c * dn) / dn); - } - void operator=(const adjointsu2 &A) { - a = A.geta(); - b = A.getb(); - c = A.getc(); - } - void operator+=(const adjointsu2 &A) { - a += A.geta(); - b += A.getb(); - c += A.getc(); - } - void operator-=(const adjointsu2 &A) { - a -= A.geta(); - b -= A.getb(); - c -= A.getc(); - } - -private: - Float a, b, c; -}; - -template inline adjointsu2 get_deriv(su2 &A) { - const Complex a = A.geta(), b = A.getb(); - return adjointsu2(2. * std::imag(b), 2. * std::real(b), 2. * std::imag(a)); -} - -template class adjointu1 { -public: - adjointu1(Float _a) : a(_a) {} - adjointu1() : a(0.) {} - void flipsign() { a = -a; } - Float geta() const { return a; } - void seta(Float _a) { a = _a; } - void setzero() { a = 0.; } - - adjointu1 round(size_t n) const { - Float dn = n; - return adjointu1(std::round(a * dn) / dn); - } - void operator=(const adjointu1 &A) { a = A.geta(); } - void operator+=(const adjointu1 &A) { a += A.geta(); } - void operator-=(const adjointu1 &A) { a -= A.geta(); } - -private: - Float a; -}; -template inline adjointu1 get_deriv(Complex &A) { - return adjointu1(std::imag(A)); -} +#include "adjoint_su2.hh" +#include "adjoint_su3.hh" +#include "adjoint_u1.hh" // The following class will be used to deliver the // adjoint type depending on the gauge group -template struct adjoint_type { typedef Group type; }; +template struct adjoint_type { + typedef Group type; +}; template struct adjoint_type { typedef adjointsu2 type; @@ -92,11 +25,16 @@ template struct adjoint_type { typedef adjointu1 type; }; +/** + * @brief field of adjoint matrices + * + * @tparam Float + * @tparam Group + */ template class adjointfield { public: using value_type = typename adjoint_type::type; - template - using nd_max_arr = typename spacetime_lattice::nd_max_arr; + template using nd_max_arr = typename spacetime_lattice::nd_max_arr; adjointfield(const size_t Lx, const size_t Ly, @@ -164,13 +102,13 @@ public: return data[getIndex(coords[0], coords[1], coords[2], coords[3], mu)]; } - template - value_type &operator()(const nd_max_arr &x, const size_t& mu) { + template + value_type &operator()(const nd_max_arr &x, const size_t &mu) { return data[getIndex(x[0], x[1], x[2], x[3], mu)]; } - template - const value_type &operator()(const nd_max_arr &x, const size_t& mu) const { + template + const value_type &operator()(const nd_max_arr &x, const size_t &mu) const { return data[getIndex(x[0], x[1], x[2], x[3], mu)]; } diff --git a/include/su3.hh b/include/su3.hh new file mode 100644 index 0000000..3623518 --- /dev/null +++ b/include/su3.hh @@ -0,0 +1,146 @@ +/** + * @file su3.hh + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief class and routines for an SU(3) matrix in the fundamental representation + * @version 0.1 + * @date 2023-02-15 + * + * @copyright Copyright (c) 2023 + * + */ + +#pragma once + +#include +// #include +#include + +#include "accum_type.hh" +#include "dagger.hh" +#include "traceless_antiherm.hh" + +using Complex = std::complex; + +/** + * @brief representation of an SU(3) matrix as U = a_0*1 + \sum_{i=1}^{8} a_i * + \lambda_i, where lambda_i are the Gell-Mann matrices normalized such that + tr(lambda_i*lambda_j)=2*delta_ij + * + */ +class _su3 { +public: + const size_t N_c = 3; + explicit _su3() { + const std::array _arr{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + arr = _arr; + } + explicit _su3(const std::array &_arr) : arr(_arr) {} + _su3(const _su3 &U) : arr(U.get_arr()) {} + + // friend inline _su3 operator+(const _su3 &U1, const _su3 &U2); + // friend inline _su3 operator-(const _su3 &U1, const _su3 &U2); + friend inline _su3 operator*(const _su3 &U1, const _su3 &U2); + friend inline _su3 operator*(const Complex &U1, const _su3 &U2); + friend inline _su3 operator*(const _su3 &U1, const Complex &U2); + // _su3 &operator*=(const _su3 &U1) { + // // CHANGE + // // 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); + // // return *this; + // } + _su3 round(size_t n) const { + double dn = n; + std::array arr2; + for (size_t i = 0; i < 9; i++) { + arr2[i] = Complex(std::round(std::real(arr[i]) * dn) / dn, + std::round(std::imag(arr[i]) * dn) / dn); + } + return _su3(arr2); + } + + inline std::array get_arr() const { return arr; } + inline void operator=(const _su3 &U) { arr = U.get_arr(); } + // inline void operator+=(const _su3 &U) { + // a += U.a; + // b += U.b; + // } + void set(const std::array &_arr) { arr = _arr; } + inline _su3 dagger() const { + std::array arr2; + for (size_t i = 0; i < 9; i++) { + arr2[i] = std::conj(arr[i]); + } + return _su3(arr2); + } + inline Complex trace() { return double(N_c) * arr[0]; } + inline double retrace() { return std::real(this->trace()); } + // 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)); + // a /= r; + // b /= r; + // } + +private: + std::array arr; // {a_0, a_1, ..., a_8} +}; + +inline Complex trace(_su3 const &U) { + return double(U.N_c) * U.get_arr()[0]; +} + +inline double retrace(_su3 const &U) { + return std::real(trace(U)); +} + +template <> inline _su3 dagger(const _su3 &u) { + return u.dagger(); +} + +template <> inline _su3 traceless_antiherm(const _su3 &x) { + std::array arr = x.get_arr(); + // make it anti-hermitian (loop can start from 1) + for (size_t i = 1; i < 9; i++) { + arr[i] = (arr[i] - std::conj(arr[i])); + } + arr[0] = 0.0; // make it traceless (trace comes from the 0-th element only) + return _su3(arr); +} + +// inline _su3 operator*(const _su3 &U1, const _su3 &U2) { +// _su3 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 _su3 operator+(const _su3 &U1, const _su3 &U2) { +// _su3 res; +// res.a = U1.a + U2.a; +// res.b = U1.b + U2.b; +// return (res); +// } + +// inline _su3 operator-(const _su3 &U1, const _su3 &U2) { +// _su3 res; +// res.a = U1.a - U2.a; +// res.b = U1.b - U2.b; +// return (res); +// } + +// inline _su3 operator*(const Complex &U1, const _su3 &U2) { +// _su3 res; +// res.a = U2.a * U1; +// res.b = U2.b * U1; +// return (res); +// } + +// inline _su3 operator*(const _su3 &U1, const Complex &U2) { +// _su3 res; +// res.a = U1.a * U2; +// res.b = U1.b * U2; +// return (res); +// } + +using su3 = _su3; From 279dad4cee855de74d23e4467cca53746a2562a8 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Thu, 16 Feb 2023 20:55:50 +0100 Subject: [PATCH 02/13] DO NOT MERGE implementing eq. 4.26 of Gattringer&Lang --- include/su2.hh | 4 +- include/su3.hh | 252 ++++++++++++++++++++++++++++++++----------------- 2 files changed, 169 insertions(+), 87 deletions(-) diff --git a/include/su2.hh b/include/su2.hh index c340ea6..88028ca 100644 --- a/include/su2.hh +++ b/include/su2.hh @@ -81,8 +81,8 @@ inline Complex trace(_su2 const &U) { 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; + _su2 udag(std::conj(a), -b); + return udag; } template <> inline _su2 traceless_antiherm(const _su2 &x) { diff --git a/include/su3.hh b/include/su3.hh index 3623518..fd7cca2 100644 --- a/include/su3.hh +++ b/include/su3.hh @@ -12,7 +12,6 @@ #pragma once #include -// #include #include #include "accum_type.hh" @@ -22,125 +21,208 @@ using Complex = std::complex; /** - * @brief representation of an SU(3) matrix as U = a_0*1 + \sum_{i=1}^{8} a_i * - \lambda_i, where lambda_i are the Gell-Mann matrices normalized such that - tr(lambda_i*lambda_j)=2*delta_ij + * @brief representation of an SU(2) matrix in the fundamental representation with the + * convention of eq. (4.26) of Gattringer&Lang + * https://link.springer.com/book/10.1007/978-3-642-01850-3 * */ class _su3 { public: const size_t N_c = 3; explicit _su3() { - const std::array _arr{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - arr = _arr; + // default: identity matrix + u = std::array({1.0, 0.0, 0.0}); + v = std::array({0.0, 1.0, 0.0}); } - explicit _su3(const std::array &_arr) : arr(_arr) {} - _su3(const _su3 &U) : arr(U.get_arr()) {} + explicit _su3(const std::array &_u, const std::array &_v) + : u(_u), v(_v) {} + _su3(const _su3 &U) : u(U.get_u()), v(U.get_v()) {} - // friend inline _su3 operator+(const _su3 &U1, const _su3 &U2); - // friend inline _su3 operator-(const _su3 &U1, const _su3 &U2); + friend inline _su3 operator+(const _su3 &U1, const _su3 &U2); + friend inline _su3 operator-(const _su3 &U1, const _su3 &U2); friend inline _su3 operator*(const _su3 &U1, const _su3 &U2); friend inline _su3 operator*(const Complex &U1, const _su3 &U2); friend inline _su3 operator*(const _su3 &U1, const Complex &U2); - // _su3 &operator*=(const _su3 &U1) { - // // CHANGE - // // 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); - // // return *this; - // } + _su3 &operator*=(const _su3 &U2) { + const std::array u1 = (*this).u; + const std::array v1 = (*this).v; + const std::array u2 = U2.u; + const std::array v2 = U2.v; + // vector components computed using sympy + (*this).u = {u1[0] * u2[0] + u1[1] * v2[0] + + u1[2] * (std::conj(u2[1]) * std::conj(v2[2]) - + std::conj(u2[2]) * std::conj(v2[1])), + u1[0] * u2[1] + u1[1] * v2[1] + + u1[2] * (-std::conj(u2[0]) * std::conj(v2[2]) + + std::conj(u2[2]) * std::conj(v2[0])), + u1[0] * u2[2] + u1[1] * v2[2] + + u1[2] * (std::conj(u2[0]) * std::conj(v2[1]) - + std::conj(u2[1]) * std::conj(v2[0]))}; + (*this).v = {u2[0] * v1[0] + v1[1] * v2[0] + + v1[2] * (std::conj(u2[1]) * std::conj(v2[2]) - + std::conj(u2[2]) * std::conj(v2[1])), + u2[1] * v1[0] + v1[1] * v2[1] + + v1[2] * (-std::conj(u2[0]) * std::conj(v2[2]) + + std::conj(u2[2]) * std::conj(v2[0])), + u2[2] * v1[0] + v1[1] * v2[2] + + v1[2] * (std::conj(u2[0]) * std::conj(v2[1]) - + std::conj(u2[1]) * std::conj(v2[0]))}; + return *this; + } _su3 round(size_t n) const { double dn = n; - std::array arr2; - for (size_t i = 0; i < 9; i++) { - arr2[i] = Complex(std::round(std::real(arr[i]) * dn) / dn, - std::round(std::imag(arr[i]) * dn) / dn); + std::array u2, v2; + for (size_t i = 0; i < N_c; i++) { + u2[i] = Complex(std::round(std::real(u[i]) * dn) / dn, + std::round(std::imag(u[i]) * dn) / dn); + v2[i] = Complex(std::round(std::real(v[i]) * dn) / dn, + std::round(std::imag(v[i]) * dn) / dn); } - return _su3(arr2); + return _su3(u2, v2); } - inline std::array get_arr() const { return arr; } - inline void operator=(const _su3 &U) { arr = U.get_arr(); } - // inline void operator+=(const _su3 &U) { - // a += U.a; - // b += U.b; - // } - void set(const std::array &_arr) { arr = _arr; } + inline std::array get_u() const { return u; } + inline std::array get_v() const { return v; } + inline void operator=(const _su3 &U) { + u = U.get_u(); + v = U.get_v(); + } + inline void operator+=(const _su3 &U) { + for (size_t i = 0; i < N_c; i++) { + u[i] += U.u[i]; + v[i] += U.v[i]; + } + } + inline void operator-=(const _su3 &U) { + for (size_t i = 0; i < N_c; i++) { + u[i] -= U.u[i]; + v[i] -= U.v[i]; + } + } + void set(const std::array &_u, const std::array &_v) { + u = _u; + v = _v; + } inline _su3 dagger() const { - std::array arr2; - for (size_t i = 0; i < 9; i++) { - arr2[i] = std::conj(arr[i]); + // vectors computed using sympy + const std::array u2 = {std::conj(u[0]), std::conj(v[0]), + u[1] * v[2] - u[2] * v[1]}; + const std::array v2 = {std::conj(u[1]), std::conj(v[1]), + -u[0] * v[2] + u[2] * v[0]}; + + return _su3(u2, v2); + } + inline Complex trace() const { + // expression computed using sympy + const Complex res = + u[0] + v[1] + std::conj(u[0]) * std::conj(v[1]) - std::conj(u[1]) * std::conj(v[0]); + return res; + } + inline double retrace() const { return std::real(this->trace()); } + Complex det() const { + // expression computed using sympy + Complex res = u[0] * v[1] * std::conj(u[0]) * std::conj(v[1]) - + u[0] * v[1] * std::conj(u[1]) * std::conj(v[0]) + + u[0] * v[2] * std::conj(u[0]) * std::conj(v[2]) - + u[0] * v[2] * std::conj(u[2]) * std::conj(v[0]) - + u[1] * v[0] * std::conj(u[0]) * std::conj(v[1]) + + u[1] * v[0] * std::conj(u[1]) * std::conj(v[0]) + + u[1] * v[2] * std::conj(u[1]) * std::conj(v[2]) - + u[1] * v[2] * std::conj(u[2]) * std::conj(v[1]) - + u[2] * v[0] * std::conj(u[0]) * std::conj(v[2]) + + u[2] * v[0] * std::conj(u[2]) * std::conj(v[0]) - + u[2] * v[1] * std::conj(u[1]) * std::conj(v[2]) + + u[2] * v[1] * std::conj(u[2]) * std::conj(v[1]); + + return res; + } + /** + * @brief formula (4.27) of Gattringer&Lang + * https://link.springer.com/book/10.1007/978-3-642-01850-3 + * + */ + void restoreSU() { + double n_u = 0.0; // norms of 'u' and 'v2' + for (size_t i = 0; i < N_c; i++) { + n_u += std::pow(std::abs((*this).u[i]), 2.0); + } + n_u = std::sqrt(n_u); + for (size_t i = 0; i < N_c; i++) { + (*this).u[i] /= n_u; } - return _su3(arr2); + const Complex v_ustar = + v[0] * std::conj(u[0]) + v[1] * std::conj(u[1]) + v[2] * std::conj(u[2]); + std::array v2 = (*this).v; + double n_v2 = 0.0; + for (size_t i = 0; i < N_c; i++) { + v2[i] -= v_ustar * u[i]; + } + for (size_t i = 0; i < N_c; i++) { + n_v2 += std::pow(std::abs(v2[i]), 2.0); + } + n_v2 = std::sqrt(n_v2); + for (size_t i = 0; i < N_c; i++) { + v2[i] /= n_v2; + } + (*this).v = v2; + return; } - inline Complex trace() { return double(N_c) * arr[0]; } - inline double retrace() { return std::real(this->trace()); } - // 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)); - // a /= r; - // b /= r; - // } private: - std::array arr; // {a_0, a_1, ..., a_8} + std::array u, v; }; inline Complex trace(_su3 const &U) { - return double(U.N_c) * U.get_arr()[0]; + return U.trace(); } inline double retrace(_su3 const &U) { - return std::real(trace(U)); + return U.retrace(); } template <> inline _su3 dagger(const _su3 &u) { return u.dagger(); } -template <> inline _su3 traceless_antiherm(const _su3 &x) { - std::array arr = x.get_arr(); - // make it anti-hermitian (loop can start from 1) - for (size_t i = 1; i < 9; i++) { - arr[i] = (arr[i] - std::conj(arr[i])); +inline _su3 operator*(const _su3 &U1, const _su3 &U2) { + _su3 U = U1; + U *= U2; + return U; +} + +inline _su3 operator+(const _su3 &U1, const _su3 &U2) { + _su3 U = U1; + U += U2; + return U; +} + +inline _su3 operator-(const _su3 &U1, const _su3 &U2) { + _su3 U = U1; + U -= U2; + return U; +} + +inline _su3 operator*(const Complex &z, const _su3 &U2) { + const size_t N_c = U2.N_c; + std::array u, v; + for (size_t i = 0; i < N_c; i++) { + u[i] *= z; + v[i] *= z; } - arr[0] = 0.0; // make it traceless (trace comes from the 0-th element only) - return _su3(arr); + return _su3(u, v); } -// inline _su3 operator*(const _su3 &U1, const _su3 &U2) { -// _su3 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 _su3 operator+(const _su3 &U1, const _su3 &U2) { -// _su3 res; -// res.a = U1.a + U2.a; -// res.b = U1.b + U2.b; -// return (res); -// } - -// inline _su3 operator-(const _su3 &U1, const _su3 &U2) { -// _su3 res; -// res.a = U1.a - U2.a; -// res.b = U1.b - U2.b; -// return (res); -// } - -// inline _su3 operator*(const Complex &U1, const _su3 &U2) { -// _su3 res; -// res.a = U2.a * U1; -// res.b = U2.b * U1; -// return (res); -// } - -// inline _su3 operator*(const _su3 &U1, const Complex &U2) { -// _su3 res; -// res.a = U1.a * U2; -// res.b = U1.b * U2; -// return (res); -// } +inline _su3 operator*(const _su3 &U1, const Complex &z) { + return z * U1; +} + +template <> inline _su3 traceless_antiherm(const _su3 &x0) { + const _su3 Id({1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}); // default constructor: identity operator + _su3 x = x0; + x = x - (trace(x) / double(x.N_c)) * Id; + x = 0.5 * (x - x.dagger()); + return x; +} using su3 = _su3; From 741b8ceceb7ce06a22b133c230961c9df2b1f24b Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Thu, 16 Feb 2023 21:09:19 +0100 Subject: [PATCH 03/13] implemented the adjointsu3 class --- include/adjoint_su3.hh | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh index 705b0c9..1a89f55 100644 --- a/include/adjoint_su3.hh +++ b/include/adjoint_su3.hh @@ -23,7 +23,9 @@ /** * @brief element of the su(3) algebra * - * the element is parametrized by 8 real numbers (number of generators for SU(3)) + * any su(3) matrix can be expressed as the anti-hermitian traceless part of an SU(3) + * matrix therefore, the object is parametrized by 8 real numbers (number of generators + * for SU(3)) * * @tparam Float */ @@ -69,11 +71,21 @@ private: std::array arr; }; +/** + * @brief returns (-i)*(B - B.dagger()), where B = A - A - (tr(a)/N_c) * 1, parametrized + * as an element of su(3) + * + * @tparam Float + * @param A + * @return adjointsu3 + */ template inline adjointsu3 get_deriv(const su3 &A) { - const std::array, 9> arr_A = A.get_arr(); - const std::array arr; - for (size_t i = 0; i < 8; i++) { - arr[i] = 2.0 * std::imag(arr_A[i+1]); - } + const su3 A_thh = traceless_antiherm(A); + const std::array u = A_thh.get_u(); + const std::array v = A_thh.get_v(); + const std::array, 8> arr = { + 2.0 * std::imag(u[0]), 2.0 * std::imag(u[1]), 2.0 * std::imag(u[2]), + 2.0 * std::imag(v[0]), + 2.0 * std::imag(v[1])}; // u*v=0 --> v[2] depends on v[0] and v[1] return adjointsu3(arr); } From d760027343153847900a2f66f6eeba42ee1b8ac7 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Thu, 16 Feb 2023 22:09:37 +0100 Subject: [PATCH 04/13] debugging --- CMakeLists.txt | 6 +++-- include/adjoint_su3.hh | 2 +- include/adjointfield.hh | 25 +++++++++++++++--- include/random_element.hh | 54 ++++++++++++++++++++++++++------------- include/staggered.hpp | 28 ++++++++++++++++---- include/su3.hh | 27 +++++++++++++++----- include/u1.hh | 2 ++ main-su3.cpp | 19 ++++++++++++++ 8 files changed, 127 insertions(+), 36 deletions(-) create mode 100644 main-su3.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ab37fcd..1b8e8cf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -88,6 +88,7 @@ endif() set(${CMAKE_BINARY_DIR} ${CMAKE_INSTALL_PREFIX}/bin/) add_executable(u1-main main-u1.cpp) add_executable(su2-main main-su2.cpp) +add_executable(su3-main main-su3.cpp) add_executable(su2-kramers kramers.cc) add_executable(test-groups test.cc) @@ -118,7 +119,8 @@ 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 ${test_progs_fullname}) +foreach(target u1-main su2-main su3-main su2-kramers test-groups scaling try ${test_progs_fullname}) + target_link_libraries(${target} su2 ${YAML_CPP_LIBRARIES}) target_link_libraries(${target} su2 ${YAML_CPP_LIBRARIES}) if(Boost_FOUND) @@ -133,7 +135,7 @@ endforeach(target) install(TARGETS ${PROJECT_NAME} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} ) -install(TARGETS u1-main su2-main) +install(TARGETS u1-main su2-main su3-main) install(TARGETS su2-kramers test-groups scaling try) # tests diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh index 1a89f55..4c64f2f 100644 --- a/include/adjoint_su3.hh +++ b/include/adjoint_su3.hh @@ -83,7 +83,7 @@ template inline adjointsu3 get_deriv(const su3 const su3 A_thh = traceless_antiherm(A); const std::array u = A_thh.get_u(); const std::array v = A_thh.get_v(); - const std::array, 8> arr = { + const std::array arr = { 2.0 * std::imag(u[0]), 2.0 * std::imag(u[1]), 2.0 * std::imag(u[2]), 2.0 * std::imag(v[0]), 2.0 * std::imag(v[1])}; // u*v=0 --> v[2] depends on v[0] and v[1] diff --git a/include/adjointfield.hh b/include/adjointfield.hh index 1de3f82..bfa6e6e 100644 --- a/include/adjointfield.hh +++ b/include/adjointfield.hh @@ -17,12 +17,16 @@ template struct adjoint_type { typedef Group type; }; +template struct adjoint_type { + typedef adjointu1 type; +}; + template struct adjoint_type { typedef adjointsu2 type; }; -template struct adjoint_type { - typedef adjointu1 type; +template struct adjoint_type { + typedef adjointsu3 type; }; /** @@ -175,6 +179,15 @@ Float operator*(const adjointfield &A, const adjointfield +void initnormal(URNG &engine, adjointfield &A) { + std::normal_distribution normal(0., 1.); + for (size_t i = 0; i < A.getSize(); i++) { + A[i].seta(Float(normal(engine))); + } + return; +} + template void initnormal(URNG &engine, adjointfield &A) { std::normal_distribution normal(0., 1.); @@ -187,10 +200,14 @@ void initnormal(URNG &engine, adjointfield &A) { } template -void initnormal(URNG &engine, adjointfield &A) { +void initnormal(URNG &engine, adjointfield &A) { std::normal_distribution normal(0., 1.); for (size_t i = 0; i < A.getSize(); i++) { - A[i].seta(Float(normal(engine))); + std::array arr; + for (size_t i = 0; i < 8; i++) { + arr[i] = Float(normal(engine)); + } + A[i].set_arr(arr); } return; } diff --git a/include/random_element.hh b/include/random_element.hh index 7a9182b..4204a0f 100644 --- a/include/random_element.hh +++ b/include/random_element.hh @@ -1,34 +1,52 @@ #pragma once -#include"su2.hh" -#include"u1.hh" -#include +#include "su2.hh" +#include "su3.hh" +#include "u1.hh" -constexpr double pi() { return std::atan(1)*4; } +#include -template void random_element(T &U, URNG &engine, - const double delta = 1.) { - +constexpr double pi() { + return std::atan(1) * 4; +} + +template +void random_element(T &U, URNG &engine, const double delta = 1.); + +template void random_element(_u1 &U, URNG &engine, const double delta = 1.) { + std::uniform_real_distribution dist(-pi() * delta, pi() * delta); + + U = _u1(dist(engine)); + return; +} + +template +void random_element(_su2 &U, URNG &engine, const double delta = 1.) { std::uniform_real_distribution dist1(-1., 1.); - std::uniform_real_distribution dist2(0., 2*pi()); - std::uniform_real_distribution dist3(0., delta*2*pi()); + std::uniform_real_distribution dist2(0., 2 * pi()); + std::uniform_real_distribution dist3(0., delta * 2 * pi()); const double alpha = dist3(engine); const double u = dist1(engine); const double theta = dist2(engine); - - const double r = sqrt(1-u*u); + + const double r = sqrt(1 - u * u); const double salpha = sin(alpha); - - U = T(Complex(cos(alpha), salpha*u), - salpha*r*Complex(sin(theta), cos(theta))); + + U = _su2(Complex(cos(alpha), salpha * u), salpha * r * Complex(sin(theta), cos(theta))); return; } -template void random_element(_u1 &U, URNG &engine, - const double delta = 1.) { +template +void random_element(_su3 &U, URNG &engine, const double delta = 1.) { + std::uniform_real_distribution dist(-delta, delta); + std::array u, v; + for (size_t i = 0; i < U.N_c; i++) { + u[i] = dist(engine); + v[i] = dist(engine); + } - std::uniform_real_distribution dist(-pi()*delta, pi()*delta); + U = _su3(u, v); + U.restoreSU(); - U = _u1(dist(engine)); return; } diff --git a/include/staggered.hpp b/include/staggered.hpp index f6da441..dee60f9 100644 --- a/include/staggered.hpp +++ b/include/staggered.hpp @@ -23,6 +23,7 @@ #include "hamiltonian_field.hh" #include "monomial.hh" #include "su2.hh" +#include "su3.hh" #include "u1.hh" #include "solver_type.hh" // generic type of solver @@ -281,7 +282,7 @@ namespace staggered { * * @tparam Float : double * @tparam Type : std::complex - * @tparam Group : su1 or su2 + * @tparam Group : u1 or su2 * @param M : wrapper for D*D^{\dagger} matrix. * Contains the information only about the configuration U and the quark mass * @param psi : spinor to which we apply D*D^{\dagger} @@ -303,9 +304,8 @@ namespace staggered { // returns the result of D*psi, where D is the Dirac operator template - spinor_lat apply_D(const gaugeconfig<_u1> &U, - const Float &m, - const spinor_lat &psi) { + spinor_lat + apply_D(const gaugeconfig<_u1> &U, const Float &m, const spinor_lat &psi) { const size_t Lt = U.getLt(), Lx = U.getLx(), Ly = U.getLy(), Lz = U.getLz(); const size_t nd = U.getndims(); const nd_max_arr dims = psi.get_dims(); // vector of dimensions @@ -313,7 +313,7 @@ namespace staggered { const int N = psi.size(); spinor_lat phi(dims); -//#pragma omp target teams distribute parallel for //collapse(4) +// #pragma omp target teams distribute parallel for //collapse(4) #pragma omp parallel for for (int x0 = 0; x0 < Lt; x0++) { for (int x1 = 0; x1 < Lx; x1++) { @@ -349,6 +349,15 @@ namespace staggered { return psi; } + template + spinor_lat apply_D(const gaugeconfig<_su3> &U, + const Float &m, + const spinor_lat &psi) { + spacetime_lattice::fatal_error("Staggered fermions not supported for SU(3)", + __func__); + return psi; + } + // returns the reslut of D^{\dagger}*psi, where D is the Dirac operator template spinor_lat apply_Ddag(const gaugeconfig<_u1> &U, @@ -397,4 +406,13 @@ namespace staggered { return psi; } + template + spinor_lat apply_Ddag(const gaugeconfig<_su3> &U, + const Float &m, + const spinor_lat &psi) { + spacetime_lattice::fatal_error("Staggered fermions not supported for SU(3)", + __func__); + return psi; + } + } // namespace staggered diff --git a/include/su3.hh b/include/su3.hh index fd7cca2..3c958ee 100644 --- a/include/su3.hh +++ b/include/su3.hh @@ -11,6 +11,7 @@ #pragma once +#include #include #include @@ -202,18 +203,32 @@ inline _su3 operator-(const _su3 &U1, const _su3 &U2) { return U; } -inline _su3 operator*(const Complex &z, const _su3 &U2) { - const size_t N_c = U2.N_c; +inline _su3 operator*(const double &z, const _su3 &U) { + const size_t N_c = U.N_c; std::array u, v; for (size_t i = 0; i < N_c; i++) { - u[i] *= z; - v[i] *= z; + u[i] *= double(z); + v[i] *= double(z); } return _su3(u, v); } -inline _su3 operator*(const _su3 &U1, const Complex &z) { - return z * U1; +inline _su3 operator*(const _su3 &U, const double &z) { + return z * U; +} + +inline _su3 operator*(const Complex &z, const _su3 &U) { + const size_t N_c = U.N_c; + std::array u, v; + for (size_t i = 0; i < N_c; i++) { + u[i] *= Complex(z); + v[i] *= Complex(z); + } + return _su3(u, v); +} + +inline _su3 operator*(const _su3 &U, const Complex &z) { + return z * U; } template <> inline _su3 traceless_antiherm(const _su3 &x0) { diff --git a/include/u1.hh b/include/u1.hh index 98ea897..62cd835 100644 --- a/include/u1.hh +++ b/include/u1.hh @@ -105,3 +105,5 @@ inline void operator+=(Complex &U1, const _u1 &U2) { inline void operator*=(Complex &U1, const _u1 &U2) { U1 *= std::exp(U2.geta() * Complex(0., 1.)); } + +using u1 = _u1; diff --git a/main-su3.cpp b/main-su3.cpp new file mode 100644 index 0000000..8e950f4 --- /dev/null +++ b/main-su3.cpp @@ -0,0 +1,19 @@ +/** + * @file main-su3.hpp + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief main programm running any simulation of this library for the SU(3) theory + * @version 0.1 + * @date 2022-10-03 + * + * @copyright Copyright (c) 2022 + * + */ + +#include "run_program.hpp" + +int main(int argc, char *argv[]) { + typedef _su3 Group; + run_program(argc, argv); + + return (0); +} From b7548d0e496732464cdfc1704484e258b354f483 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Fri, 17 Feb 2023 18:14:17 +0100 Subject: [PATCH 05/13] DO NOT MERGE 1st version that compiles successfully. E and dH in the hmc give nan, still debugging. --- exp_gauge.cc | 104 +++++++++++++++++++++++++++++++---- include/adjoint_su2.hh | 7 ++- include/adjoint_su3.hh | 15 ++++- include/adjoint_u1.hh | 14 +++++ include/adjointfield.hh | 28 ++++++---- include/detDDdag_monomial.hh | 24 ++++++++ include/exp_gauge.hh | 6 +- include/glueballs.hpp | 2 +- include/operators.hpp | 5 +- include/su2.hh | 4 ++ include/su3.hh | 4 ++ include/u1.hh | 24 ++++---- include/wilsonloop.hh | 6 +- 13 files changed, 198 insertions(+), 45 deletions(-) diff --git a/exp_gauge.cc b/exp_gauge.cc index 3708758..4d8bbeb 100644 --- a/exp_gauge.cc +++ b/exp_gauge.cc @@ -1,20 +1,100 @@ -#include"exp_gauge.hh" -#include"su2.hh" -#include"u1.hh" -#include"adjointfield.hh" -#include -#include - -_su2 exp(adjointsu2 const & x) { +#include "exp_gauge.hh" +#include "adjointfield.hh" +#include "su2.hh" +#include "u1.hh" +#include +#include +#include + +_su2 exp(adjointsu2 const &x) { double a = x.geta(), b = x.getb(), c = x.getc(); // normalise the vector - const double alpha = sqrt(a*a+b*b+c*c); - std::vector n = {a/alpha, b/alpha, c/alpha}; + const double alpha = sqrt(a * a + b * b + c * c); + std::vector n = {a / alpha, b / alpha, c / alpha}; const double salpha = sin(alpha); - _su2 res(Complex(cos(alpha), salpha*n[2]), - salpha*Complex(n[1], n[0])); + _su2 res(Complex(cos(alpha), salpha * n[2]), salpha * Complex(n[1], n[0])); return res; } +// eq. (4) of https://arxiv.org/pdf/1508.00868v2.pdf +// the Gell-Mann decomposition coefficients have been computed with sympy +_su3 exp(const adjointsu3 &x) { + std::array arr = x.get_arr(); + double theta = 0.0; + for (size_t i = 0; i < 8; i++) { + theta += std::pow(arr[i], 2.0); + } + theta = std::sqrt(theta); + + // normalised vector + std::array n = arr; + for (size_t i = 0; i < 8; i++) { + n[i] /= theta; + } + + // determinant of H + const double detH = + 2.0 * sqrt(3.0) * std::pow(n[0], 2.0) * n[7] / 3.0 + 2.0 * n[0] * n[3] * n[5] + + 2.0 * n[0] * n[4] * n[6] + 2.0 * sqrt(3.0) * std::pow(n[1], 2.0) * n[7] / 3.0 - + 2.0 * n[1] * n[3] * n[6] + 2.0 * n[1] * n[4] * n[5] + + 2.0 * sqrt(3.0) * std::pow(n[2], 2.0) * n[7] / 3.0 + n[2] * std::pow(n[3], 2.0) + + n[2] * std::pow(n[4], 2.0) - n[2] * std::pow(n[5], 2.0) - n[2] * std::pow(n[6], 2.0) - + sqrt(3.0) * std::pow(n[3], 2.0) * n[7] / 3.0 - + sqrt(3.0) * std::pow(n[4], 2.0) * n[7] / 3.0 - + sqrt(3.0) * std::pow(n[5], 2.0) * n[7] / 3.0 - + sqrt(3.0) * std::pow(n[6], 2.0) * n[7] / 3.0 - + 2.0 * sqrt(3.0) * std::pow(n[7], 3.0 / 9.0); + const double phi = (1.0 / 3.0) * (acos((3.0 / 2.0) * sqrt(3.0) * detH) - M_PI / 2.0); + + const Complex I(0.0, 1.0); // imaginary unit + + // 1st and 2nd row of H + const std::array H_1 = {n[2] + sqrt(3.0) * n[7] / 3.0, n[0] - I * n[1], + n[3] - I * n[4]}; + const std::array H_2 = {n[0] + I * n[1], -n[2] + sqrt(3.0) * n[7] / 3.0, + n[5] - I * n[6]}; + + // 1st and 2nd row of H^2 + const std::array H2_1 = { + (n[0] - I * n[1]) * (n[0] + I * n[1]) + std::pow(n[2] + sqrt(3.0) * n[7] / 3.0, 2.0) + + (n[3] - I * n[4]) * (n[3] + I * n[4]), + (n[0] - I * n[1]) * (-n[2] + sqrt(3.0) * n[7] / 3.0) + + (n[0] - I * n[1]) * (n[2] + sqrt(3.0) * n[7] / 3.0) + + (n[3] - I * n[4]) * (n[5] + I * n[6]), + -2.0 * sqrt(3.0) * n[7] * (n[3] - I * n[4]) / 3.0 + + (n[0] - I * n[1]) * (n[5] - I * n[6]) + + (n[2] + sqrt(3.0) * n[7] / 3.0) * (n[3] - I * n[4])}; + const std::array H2_2 = { + (n[0] + I * n[1]) * (-n[2] + sqrt(3.0) * n[7] / 3.0) + + (n[0] + I * n[1]) * (n[2] + sqrt(3.0) * n[7] / 3.0) + + (n[3] + I * n[4]) * (n[5] - I * n[6]), + (n[0] - I * n[1]) * (n[0] + I * n[1]) + + std::pow(-n[2] + sqrt(3.0) * n[7] / 3.0, 2.0) + + (n[5] - I * n[6]) * (n[5] + I * n[6]), + -2.0 * sqrt(3.0) * n[7] * (n[5] - I * n[6]) / 3.0 + + (n[0] + I * n[1]) * (n[3] - I * n[4]) + + (-n[2] + sqrt(3.0) * n[7] / 3.0) * (n[5] - I * n[6])}; + std::array u, v; + for (size_t k = 0; k <= 2; k++) { + const double pk = phi + (2.0 / 3.0) * M_PI * k; + const double arg_num_k = (2.0 / sqrt(3.0)) * theta * sin(pk); + const Complex num_k = std::exp(I * arg_num_k); + const double den_k = 1 - 2.0 * cos(2.0 * pk); + const Complex fact_k = num_k / den_k; + for (size_t i = 0; i < 3; i++) { + const Complex u_ik = H2_1[i] + (2.0 / sqrt(3.0)) * H_1[i] * sin(pk) - + (1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)); + u[i] += u_ik * fact_k; + + const Complex v_ik = H2_2[i] + (2.0 / sqrt(3.0)) * H_2[i] * sin(pk) - + (1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)); + + v[i] += v_ik; + } + } + + _su3 res(u, v); + return res; +} diff --git a/include/adjoint_su2.hh b/include/adjoint_su2.hh index 92d1e23..8abf345 100644 --- a/include/adjoint_su2.hh +++ b/include/adjoint_su2.hh @@ -20,7 +20,6 @@ #include #include - template class adjointsu2 { public: adjointsu2(Float _a, Float _b, Float _c) : a(_a), b(_b), c(_c) {} @@ -62,7 +61,13 @@ private: Float a, b, c; }; + template inline adjointsu2 get_deriv(su2 &A) { const Complex a = A.geta(), b = A.getb(); return adjointsu2(2. * std::imag(b), 2. * std::real(b), 2. * std::imag(a)); } + +template +inline adjointsu2 operator*(const Float &x, const adjointsu2 &A) { + return adjointsu2(x * A.geta(), x * A.getb(), x * A.getc()); +} diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh index 4c64f2f..6dc1c07 100644 --- a/include/adjoint_su3.hh +++ b/include/adjoint_su3.hh @@ -72,7 +72,7 @@ private: }; /** - * @brief returns (-i)*(B - B.dagger()), where B = A - A - (tr(a)/N_c) * 1, parametrized + * @brief returns (-i)*(B - B.dagger()), where B = A - (tr(a)/N_c) * 1, parametrized * as an element of su(3) * * @tparam Float @@ -81,11 +81,20 @@ private: */ template inline adjointsu3 get_deriv(const su3 &A) { const su3 A_thh = traceless_antiherm(A); - const std::array u = A_thh.get_u(); - const std::array v = A_thh.get_v(); + const std::array u = A_thh.get_u(); + const std::array v = A_thh.get_v(); const std::array arr = { 2.0 * std::imag(u[0]), 2.0 * std::imag(u[1]), 2.0 * std::imag(u[2]), 2.0 * std::imag(v[0]), 2.0 * std::imag(v[1])}; // u*v=0 --> v[2] depends on v[0] and v[1] return adjointsu3(arr); } + +template +inline adjointsu3 operator*(const Float &x, const adjointsu3 &A) { + std::array arr = A.get_arr(); + for (size_t i = 0; i < 8; i++) { + arr[i] *= x; + } + return adjointsu3(arr); +} diff --git a/include/adjoint_u1.hh b/include/adjoint_u1.hh index 8789a05..90bc155 100644 --- a/include/adjoint_u1.hh +++ b/include/adjoint_u1.hh @@ -21,6 +21,15 @@ #include #include +/** + * @brief element of the adjount representation of U(1), i.e. just a number proportional + * to the Identity. + * + * Note: We use the convention such that tr(t_i*t_j) = delta_{ij} (no factor 1/2), so the + * generator of the algebra is the Identity, not divided by sqrt(2) + * + * @tparam Float + */ template class adjointu1 { public: adjointu1(Float _a) : a(_a) {} @@ -45,3 +54,8 @@ private: template inline adjointu1 get_deriv(Complex &A) { return adjointu1(std::imag(A)); } + +template +inline adjointu1 operator*(const Float &x, const adjointu1 &A) { + return adjointu1(x * A.geta()); +} diff --git a/include/adjointfield.hh b/include/adjointfield.hh index bfa6e6e..6e46de3 100644 --- a/include/adjointfield.hh +++ b/include/adjointfield.hh @@ -139,16 +139,6 @@ private: } }; -template -inline adjointsu2 operator*(const Float &x, const adjointsu2 &A) { - return adjointsu2(x * A.geta(), x * A.getb(), x * A.getc()); -} - -template -inline adjointu1 operator*(const Float &x, const adjointu1 &A) { - return adjointu1(x * A.geta()); -} - template adjointfield operator*(const Float &x, const adjointfield &A) { adjointfield res(A.getLx(), A.getLy(), A.getLz(), A.getLt()); @@ -158,6 +148,16 @@ adjointfield operator*(const Float &x, const adjointfield +Float operator*(const adjointfield &A, const adjointfield &B) { + Float res = 0.; + assert(A.getSize() == B.getSize()); + for (size_t i = 0; i < A.getSize(); i++) { + res += A[i].geta() * B[i].geta(); + } + return res; +} + template Float operator*(const adjointfield &A, const adjointfield &B) { Float res = 0.; @@ -170,11 +170,15 @@ Float operator*(const adjointfield &A, const adjointfield -Float operator*(const adjointfield &A, const adjointfield &B) { +Float operator*(const adjointfield &A, const adjointfield &B) { Float res = 0.; assert(A.getSize() == B.getSize()); for (size_t i = 0; i < A.getSize(); i++) { - res += A[i].geta() * B[i].geta(); + const std::array arr_A = A[i].get_arr(); + const std::array arr_B = B[i].get_arr(); + for (size_t k = 0; k < 8; k++) { + res += arr_A[i] * arr_B[i]; + } } return res; } diff --git a/include/detDDdag_monomial.hh b/include/detDDdag_monomial.hh index 040d854..70a23af 100644 --- a/include/detDDdag_monomial.hh +++ b/include/detDDdag_monomial.hh @@ -191,4 +191,28 @@ namespace staggered { } }; + template + class detDDdag_monomial : public monomial { + public: + detDDdag_monomial(unsigned int _timescale, + const Float &m0_val, + const std::string &solver, + const Float &tolerance, + const size_t &seed, + const size_t &verb) + : monomial::monomial(_timescale) { + spacetime_lattice::fatal_error("SU(3) not supported yet.", __func__); + } + + void heatbath(const hamiltonian_field &h) override { return; } + + void accept(const hamiltonian_field &h) override { return; } + + void derivative(adjointfield &deriv, + const hamiltonian_field &h, + const Float fac = 1.) const override { + return; + } + }; + } // namespace staggered diff --git a/include/exp_gauge.hh b/include/exp_gauge.hh index 68231f9..cb1891d 100644 --- a/include/exp_gauge.hh +++ b/include/exp_gauge.hh @@ -1,11 +1,13 @@ #pragma once #include "su2.hh" +#include "su3.hh" #include "u1.hh" #include "adjointfield.hh" -_su2 exp(adjointsu2 const & x); - inline _u1 exp(adjointu1 const & x) { return _u1(x.geta()); } + +_su2 exp(adjointsu2 const & x); +_su3 exp(adjointsu3 const & x); diff --git a/include/glueballs.hpp b/include/glueballs.hpp index d1bb247..7d33372 100644 --- a/include/glueballs.hpp +++ b/include/glueballs.hpp @@ -31,7 +31,7 @@ namespace glueballs { const bool &Px) { typedef typename accum_type::type accum; - accum L = U(x, mu) * (U(x, mu).dagger()); // "1", independently of the group + Group L = U(x, mu) * (U(x, mu).dagger()); // "1", independently of the group for (size_t s = 0; s < a; s++) { L *= operators::parity(Px, 0, U, x, mu); x[mu] += 1; diff --git a/include/operators.hpp b/include/operators.hpp index 87956c7..bf4e0a2 100644 --- a/include/operators.hpp +++ b/include/operators.hpp @@ -12,8 +12,8 @@ #include #include -#include "gradient_flow.hh" #include "gaugeconfig.hh" +#include "gradient_flow.hh" #include "links.hpp" #include "parameters.hh" #include "propagator.hpp" @@ -111,7 +111,8 @@ namespace operators { abort(); } std::vector xrun = x; - accum L(1., 0.); + Group L; + L.set_to_identity(); // L = 1.0 for (size_t i = 0; i < lengths.size(); i++) { if (sign[i]) { for (size_t j = 0; j < lengths[i]; j++) { diff --git a/include/su2.hh b/include/su2.hh index 88028ca..117881b 100644 --- a/include/su2.hh +++ b/include/su2.hh @@ -55,6 +55,10 @@ public: a = _a; b = _b; } + void set_to_identity() { + a = 1.0; + b = 0.0; + } 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)); } diff --git a/include/su3.hh b/include/su3.hh index 3c958ee..2af648c 100644 --- a/include/su3.hh +++ b/include/su3.hh @@ -104,6 +104,10 @@ public: u = _u; v = _v; } + void set_to_identity() { + u = {1.0, 0.0, 0.0}; + v = {0.0, 1.0, 0.0}; + } inline _su3 dagger() const { // vectors computed using sympy const std::array u2 = {std::conj(u[0]), std::conj(v[0]), diff --git a/include/u1.hh b/include/u1.hh index 62cd835..5146f2c 100644 --- a/include/u1.hh +++ b/include/u1.hh @@ -35,7 +35,8 @@ public: double geta() const { return (a); } void operator=(const _u1 &U) { a = U.a; } - void set(const double _a) { a = _a; } + void set(const double& _a) { a = _a; } + void set_to_identity() { a = 0; } _u1 dagger() const { return (_u1(-a)); } double retrace() const { return (std::cos(a)); } Complex det() const { return (std::exp(a * Complex(0., 1.))); } @@ -65,21 +66,24 @@ template <> inline Complex dagger(const Complex &x) { return std::conj(x); } +template <> inline _u1 dagger(const _u1 &x) { + return x.dagger(); +} + + 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); +_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); inline _u1 operator*(const _u1 &U1, const _u1 &U2) { - // _u1 res; - // res.a = U1.a + U2.a; return _u1(U1.a + U2.a); } diff --git a/include/wilsonloop.hh b/include/wilsonloop.hh index 47ce1d6..eb3656a 100644 --- a/include/wilsonloop.hh +++ b/include/wilsonloop.hh @@ -63,7 +63,8 @@ double planar_wilsonloop_dir(const gaugeconfig &U, for (x[2] = 0; x[2] < U.getLy(); x[2]++) { for (x[3] = 0; x[3] < U.getLz(); x[3]++) { std::vector xrun = x; - accum L(1., 0.); + Group L; + L.set_to_identity(); // L = 1.0 for (size_t _t = 0; _t < t; _t++) { L *= U(xrun, nu); xrun[nu] += 1; @@ -113,7 +114,8 @@ double wilsonloop_non_planar(const gaugeconfig &U, std::vector r) for (size_t x2 = 0; x2 < U.getLy(); x2++) { for (size_t x3 = 0; x3 < U.getLz(); x3++) { std::vector xrun = {x0, x1, x2, x3}; - accum L(1., 0.); + Group L; + L.set_to_identity(); // L = 1.0 // needed if vector with directions contains more than 4 entries/if another // order than t-x-y-z is wanted size_t directionloop; From ee97d5553dc4296cd421a45cfeeedf773199e79c Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Fri, 17 Feb 2023 19:37:16 +0100 Subject: [PATCH 06/13] fact_k missing for v vector --- exp_gauge.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/exp_gauge.cc b/exp_gauge.cc index 4d8bbeb..c79fcf2 100644 --- a/exp_gauge.cc +++ b/exp_gauge.cc @@ -76,7 +76,7 @@ _su3 exp(const adjointsu3 &x) { (n[0] + I * n[1]) * (n[3] - I * n[4]) + (-n[2] + sqrt(3.0) * n[7] / 3.0) * (n[5] - I * n[6])}; - std::array u, v; + std::array u = {0.0, 0.0, 0.0}, v = {0.0, 0.0, 0.0}; for (size_t k = 0; k <= 2; k++) { const double pk = phi + (2.0 / 3.0) * M_PI * k; const double arg_num_k = (2.0 / sqrt(3.0)) * theta * sin(pk); @@ -91,7 +91,7 @@ _su3 exp(const adjointsu3 &x) { const Complex v_ik = H2_2[i] + (2.0 / sqrt(3.0)) * H_2[i] * sin(pk) - (1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)); - v[i] += v_ik; + v[i] += v_ik * fact_k; } } From afe37dfc6bea96637c2d342deaf4baa41d8c644f Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Fri, 17 Feb 2023 20:02:39 +0100 Subject: [PATCH 07/13] to do: check line 78 of exp_gauge.cc --- exp_gauge.cc | 53 ++++++++++++++++++++++++++-------------------------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/exp_gauge.cc b/exp_gauge.cc index c79fcf2..50e8f02 100644 --- a/exp_gauge.cc +++ b/exp_gauge.cc @@ -20,10 +20,12 @@ _su2 exp(adjointsu2 const &x) { // eq. (4) of https://arxiv.org/pdf/1508.00868v2.pdf // the Gell-Mann decomposition coefficients have been computed with sympy _su3 exp(const adjointsu3 &x) { + const double sqrt_of_3 = sqrt(3.0); // avoid computing it many times + std::array arr = x.get_arr(); double theta = 0.0; for (size_t i = 0; i < 8; i++) { - theta += std::pow(arr[i], 2.0); + theta += (arr[i] * arr[i]); } theta = std::sqrt(theta); @@ -35,60 +37,59 @@ _su3 exp(const adjointsu3 &x) { // determinant of H const double detH = - 2.0 * sqrt(3.0) * std::pow(n[0], 2.0) * n[7] / 3.0 + 2.0 * n[0] * n[3] * n[5] + - 2.0 * n[0] * n[4] * n[6] + 2.0 * sqrt(3.0) * std::pow(n[1], 2.0) * n[7] / 3.0 - + 2.0 * sqrt_of_3 * (n[0] * n[0]) * n[7] / 3.0 + 2.0 * n[0] * n[3] * n[5] + + 2.0 * n[0] * n[4] * n[6] + 2.0 * sqrt_of_3 * (n[1] * n[1]) * n[7] / 3.0 - 2.0 * n[1] * n[3] * n[6] + 2.0 * n[1] * n[4] * n[5] + - 2.0 * sqrt(3.0) * std::pow(n[2], 2.0) * n[7] / 3.0 + n[2] * std::pow(n[3], 2.0) + - n[2] * std::pow(n[4], 2.0) - n[2] * std::pow(n[5], 2.0) - n[2] * std::pow(n[6], 2.0) - - sqrt(3.0) * std::pow(n[3], 2.0) * n[7] / 3.0 - - sqrt(3.0) * std::pow(n[4], 2.0) * n[7] / 3.0 - - sqrt(3.0) * std::pow(n[5], 2.0) * n[7] / 3.0 - - sqrt(3.0) * std::pow(n[6], 2.0) * n[7] / 3.0 - - 2.0 * sqrt(3.0) * std::pow(n[7], 3.0 / 9.0); - const double phi = (1.0 / 3.0) * (acos((3.0 / 2.0) * sqrt(3.0) * detH) - M_PI / 2.0); + 2.0 * sqrt_of_3 * (n[2] * n[2]) * n[7] / 3.0 + n[2] * (n[3] * n[3]) + + n[2] * (n[4] * n[4]) - n[2] * (n[5] * n[5]) - n[2] * (n[6] * n[6]) - + sqrt_of_3 * (n[3] * n[3]) * n[7] / 3.0 - sqrt_of_3 * (n[4] * n[4]) * n[7] / 3.0 - + sqrt_of_3 * (n[5] * n[5]) * n[7] / 3.0 - sqrt_of_3 * (n[6] * n[6]) * n[7] / 3.0 - + 2.0 * sqrt_of_3 * (n[7] * n[7]*n[7]) / 9.0; + const double phi = (1.0 / 3.0) * (acos((3.0 / 2.0) * sqrt_of_3 * detH) - M_PI / 2.0); const Complex I(0.0, 1.0); // imaginary unit // 1st and 2nd row of H - const std::array H_1 = {n[2] + sqrt(3.0) * n[7] / 3.0, n[0] - I * n[1], + const std::array H_1 = {n[2] + sqrt_of_3 * n[7] / 3.0, n[0] - I * n[1], n[3] - I * n[4]}; - const std::array H_2 = {n[0] + I * n[1], -n[2] + sqrt(3.0) * n[7] / 3.0, + const std::array H_2 = {n[0] + I * n[1], -n[2] + sqrt_of_3 * n[7] / 3.0, n[5] - I * n[6]}; // 1st and 2nd row of H^2 const std::array H2_1 = { - (n[0] - I * n[1]) * (n[0] + I * n[1]) + std::pow(n[2] + sqrt(3.0) * n[7] / 3.0, 2.0) + + (n[0] - I * n[1]) * (n[0] + I * n[1]) + + (n[2] + sqrt_of_3 * n[7] / 3.0) * (n[2] + sqrt_of_3 * n[7] / 3.0) + (n[3] - I * n[4]) * (n[3] + I * n[4]), - (n[0] - I * n[1]) * (-n[2] + sqrt(3.0) * n[7] / 3.0) + - (n[0] - I * n[1]) * (n[2] + sqrt(3.0) * n[7] / 3.0) + + (n[0] - I * n[1]) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + + (n[0] - I * n[1]) * (n[2] + sqrt_of_3 * n[7] / 3.0) + (n[3] - I * n[4]) * (n[5] + I * n[6]), - -2.0 * sqrt(3.0) * n[7] * (n[3] - I * n[4]) / 3.0 + + -2.0 * sqrt_of_3 * n[7] * (n[3] - I * n[4]) / 3.0 + (n[0] - I * n[1]) * (n[5] - I * n[6]) + - (n[2] + sqrt(3.0) * n[7] / 3.0) * (n[3] - I * n[4])}; + (n[2] + sqrt_of_3 * n[7] / 3.0) * (n[3] - I * n[4])}; const std::array H2_2 = { - (n[0] + I * n[1]) * (-n[2] + sqrt(3.0) * n[7] / 3.0) + - (n[0] + I * n[1]) * (n[2] + sqrt(3.0) * n[7] / 3.0) + + (n[0] + I * n[1]) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + + (n[0] + I * n[1]) * (n[2] + sqrt_of_3 * n[7] / 3.0) + (n[3] + I * n[4]) * (n[5] - I * n[6]), (n[0] - I * n[1]) * (n[0] + I * n[1]) + - std::pow(-n[2] + sqrt(3.0) * n[7] / 3.0, 2.0) + + (-n[2] + sqrt_of_3 * n[7] / 3.0) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + (n[5] - I * n[6]) * (n[5] + I * n[6]), - -2.0 * sqrt(3.0) * n[7] * (n[5] - I * n[6]) / 3.0 + + -2.0 * sqrt_of_3 * n[7] * (n[5] - I * n[6]) / 3.0 + (n[0] + I * n[1]) * (n[3] - I * n[4]) + - (-n[2] + sqrt(3.0) * n[7] / 3.0) * (n[5] - I * n[6])}; + (-n[2] + sqrt_of_3 * n[7] / 3.0) * (n[5] - I * n[6])}; std::array u = {0.0, 0.0, 0.0}, v = {0.0, 0.0, 0.0}; for (size_t k = 0; k <= 2; k++) { const double pk = phi + (2.0 / 3.0) * M_PI * k; - const double arg_num_k = (2.0 / sqrt(3.0)) * theta * sin(pk); + const double arg_num_k = (2.0 / sqrt_of_3) * theta * sin(pk); const Complex num_k = std::exp(I * arg_num_k); const double den_k = 1 - 2.0 * cos(2.0 * pk); const Complex fact_k = num_k / den_k; for (size_t i = 0; i < 3; i++) { - const Complex u_ik = H2_1[i] + (2.0 / sqrt(3.0)) * H_1[i] * sin(pk) - + const Complex u_ik = H2_1[i] + (2.0 / sqrt_of_3) * H_1[i] * sin(pk) - (1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)); u[i] += u_ik * fact_k; - const Complex v_ik = H2_2[i] + (2.0 / sqrt(3.0)) * H_2[i] * sin(pk) - + const Complex v_ik = H2_2[i] + (2.0 / sqrt_of_3) * H_2[i] * sin(pk) - (1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)); v[i] += v_ik * fact_k; From 2dd2f44de31a695cc57957fa7d67c6f489aeba22 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Sat, 18 Feb 2023 13:19:50 +0100 Subject: [PATCH 08/13] still debugging, get_deriv not implemented yet --- exp_gauge.cc | 12 +++++++++++- include/adjoint_su3.hh | 7 +++---- include/su3.hh | 12 ++++++++++++ 3 files changed, 26 insertions(+), 5 deletions(-) diff --git a/exp_gauge.cc b/exp_gauge.cc index 50e8f02..b20decf 100644 --- a/exp_gauge.cc +++ b/exp_gauge.cc @@ -44,7 +44,7 @@ _su3 exp(const adjointsu3 &x) { n[2] * (n[4] * n[4]) - n[2] * (n[5] * n[5]) - n[2] * (n[6] * n[6]) - sqrt_of_3 * (n[3] * n[3]) * n[7] / 3.0 - sqrt_of_3 * (n[4] * n[4]) * n[7] / 3.0 - sqrt_of_3 * (n[5] * n[5]) * n[7] / 3.0 - sqrt_of_3 * (n[6] * n[6]) * n[7] / 3.0 - - 2.0 * sqrt_of_3 * (n[7] * n[7]*n[7]) / 9.0; + 2.0 * sqrt_of_3 * (n[7] * n[7] * n[7]) / 9.0; const double phi = (1.0 / 3.0) * (acos((3.0 / 2.0) * sqrt_of_3 * detH) - M_PI / 2.0); const Complex I(0.0, 1.0); // imaginary unit @@ -96,6 +96,16 @@ _su3 exp(const adjointsu3 &x) { } } + // Complex dp = 0.0; + // for (size_t i = 0; i < 3; i++) { + // dp += std::conj(u[i]) * v[i]; + // } + // std::cout << "Perpendicular? " << dp << "\n"; + // // std::abort(); + _su3 res(u, v); +// res.restoreSU(); + _su3 Uc = res*(res.dagger()); + Uc.print(); return res; } diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh index 6dc1c07..af02a7a 100644 --- a/include/adjoint_su3.hh +++ b/include/adjoint_su3.hh @@ -80,13 +80,12 @@ private: * @return adjointsu3 */ template inline adjointsu3 get_deriv(const su3 &A) { + std::cout << "ERROR: please implement " << __func__ << " for SU(3) \n"; + std::abort(); const su3 A_thh = traceless_antiherm(A); const std::array u = A_thh.get_u(); const std::array v = A_thh.get_v(); - const std::array arr = { - 2.0 * std::imag(u[0]), 2.0 * std::imag(u[1]), 2.0 * std::imag(u[2]), - 2.0 * std::imag(v[0]), - 2.0 * std::imag(v[1])}; // u*v=0 --> v[2] depends on v[0] and v[1] + const std::array arr = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // = ??? return adjointsu3(arr); } diff --git a/include/su3.hh b/include/su3.hh index 2af648c..c65bf4c 100644 --- a/include/su3.hh +++ b/include/su3.hh @@ -173,6 +173,18 @@ public: return; } + void print() const { + std::array w = { + u[1] * v[2] - v[1] * u[2], -u[0] * v[2] + v[0] * u[2], u[0] * v[1] - v[0] * u[1]}; + w = {std::conj(w[0]), std::conj(w[1]), std::conj(w[2])}; + + std::cout << "----------------------------------\n"; + std::cout << u[0] << " " << u[1] << " " << u[2] << "\n"; + std::cout << v[0] << " " << v[1] << " " << v[2] << "\n"; + std::cout << w[0] << " " << w[1] << " " << w[2] << "\n"; + std::cout << "----------------------------------\n"; + } + private: std::array u, v; }; From 4f50c3cca9d40841781dcda276af455fda5d9ba8 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Tue, 21 Feb 2023 14:53:44 +0100 Subject: [PATCH 09/13] DO NOT MERGE - exponentiation of su(3) works. The formulas were found using sympy and I made a numerical test of the unitarity of the matrix. - get_deriv() still implemented incorrectly --- exp_gauge.cc | 28 +++++++++------------------- include/adjoint_su3.hh | 4 ++-- include/adjointfield.hh | 9 +++++++-- 3 files changed, 18 insertions(+), 23 deletions(-) diff --git a/exp_gauge.cc b/exp_gauge.cc index b20decf..203b37b 100644 --- a/exp_gauge.cc +++ b/exp_gauge.cc @@ -57,13 +57,12 @@ _su3 exp(const adjointsu3 &x) { // 1st and 2nd row of H^2 const std::array H2_1 = { - (n[0] - I * n[1]) * (n[0] + I * n[1]) + - (n[2] + sqrt_of_3 * n[7] / 3.0) * (n[2] + sqrt_of_3 * n[7] / 3.0) + + (n[0] - I * n[1]) * (n[0] + I * n[1]) + std::pow(n[2] + sqrt_of_3 * n[7] / 3.0, 2.0) + (n[3] - I * n[4]) * (n[3] + I * n[4]), (n[0] - I * n[1]) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + (n[0] - I * n[1]) * (n[2] + sqrt_of_3 * n[7] / 3.0) + (n[3] - I * n[4]) * (n[5] + I * n[6]), - -2.0 * sqrt_of_3 * n[7] * (n[3] - I * n[4]) / 3.0 + + -2 * sqrt_of_3 * n[7] * (n[3] - I * n[4]) / 3.0 + (n[0] - I * n[1]) * (n[5] - I * n[6]) + (n[2] + sqrt_of_3 * n[7] / 3.0) * (n[3] - I * n[4])}; const std::array H2_2 = { @@ -71,9 +70,9 @@ _su3 exp(const adjointsu3 &x) { (n[0] + I * n[1]) * (n[2] + sqrt_of_3 * n[7] / 3.0) + (n[3] + I * n[4]) * (n[5] - I * n[6]), (n[0] - I * n[1]) * (n[0] + I * n[1]) + - (-n[2] + sqrt_of_3 * n[7] / 3.0) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + + std::pow(-n[2] + sqrt_of_3 * n[7] / 3.0, 2.0) + (n[5] - I * n[6]) * (n[5] + I * n[6]), - -2.0 * sqrt_of_3 * n[7] * (n[5] - I * n[6]) / 3.0 + + -2 * sqrt_of_3 * n[7] * (n[5] - I * n[6]) / 3.0 + (n[0] + I * n[1]) * (n[3] - I * n[4]) + (-n[2] + sqrt_of_3 * n[7] / 3.0) * (n[5] - I * n[6])}; @@ -85,27 +84,18 @@ _su3 exp(const adjointsu3 &x) { const double den_k = 1 - 2.0 * cos(2.0 * pk); const Complex fact_k = num_k / den_k; for (size_t i = 0; i < 3; i++) { - const Complex u_ik = H2_1[i] + (2.0 / sqrt_of_3) * H_1[i] * sin(pk) - - (1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)); + const Complex u_ik = H2_1[i] + (2.0 / sqrt_of_3) * H_1[i] * sin(pk); u[i] += u_ik * fact_k; - const Complex v_ik = H2_2[i] + (2.0 / sqrt_of_3) * H_2[i] * sin(pk) - - (1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)); + const Complex v_ik = H2_2[i] + (2.0 / sqrt_of_3) * H_2[i] * sin(pk); v[i] += v_ik * fact_k; } + // contibution from the identity + u[0] += -(1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)) * fact_k; + v[1] += -(1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)) * fact_k; } - // Complex dp = 0.0; - // for (size_t i = 0; i < 3; i++) { - // dp += std::conj(u[i]) * v[i]; - // } - // std::cout << "Perpendicular? " << dp << "\n"; - // // std::abort(); - _su3 res(u, v); -// res.restoreSU(); - _su3 Uc = res*(res.dagger()); - Uc.print(); return res; } diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh index af02a7a..03f5b69 100644 --- a/include/adjoint_su3.hh +++ b/include/adjoint_su3.hh @@ -80,8 +80,8 @@ private: * @return adjointsu3 */ template inline adjointsu3 get_deriv(const su3 &A) { - std::cout << "ERROR: please implement " << __func__ << " for SU(3) \n"; - std::abort(); +// std::cout << "ERROR: please implement " << __func__ << " for SU(3) \n"; +// std::abort(); const su3 A_thh = traceless_antiherm(A); const std::array u = A_thh.get_u(); const std::array v = A_thh.get_v(); diff --git a/include/adjointfield.hh b/include/adjointfield.hh index 6e46de3..e84323e 100644 --- a/include/adjointfield.hh +++ b/include/adjointfield.hh @@ -46,7 +46,7 @@ public: const size_t Lt, const size_t ndims = 4) : Lx(Lx), Ly(Ly), Lz(Lz), Lt(Lt), volume(Lx * Ly * Lz * Lt), ndims(ndims) { - data.resize(volume * 4); + data.resize(volume * ndims); } adjointfield(const adjointfield &U) : Lx(U.getLx()), @@ -120,6 +120,11 @@ public: const value_type &operator[](size_t const index) const { return data[index]; } + std::vector get_data() const { + std::vector data2 = data; + return data2; + } + private: size_t Lx, Ly, Lz, Lt, volume, ndims; @@ -177,7 +182,7 @@ Float operator*(const adjointfield &A, const adjointfield arr_A = A[i].get_arr(); const std::array arr_B = B[i].get_arr(); for (size_t k = 0; k < 8; k++) { - res += arr_A[i] * arr_B[i]; + res += arr_A[k] * arr_B[k]; } } return res; From 7a63e5acc1eb2f3be06059ab12bb232f20848e24 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Tue, 21 Feb 2023 15:03:27 +0100 Subject: [PATCH 10/13] better file structure + added some comments about the implementation --- exp_gauge.cc | 114 ++++++--------------------------------------------- exp_su2.cc | 40 ++++++++++++++++++ exp_su3.cc | 110 +++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 163 insertions(+), 101 deletions(-) create mode 100644 exp_su2.cc create mode 100644 exp_su3.cc diff --git a/exp_gauge.cc b/exp_gauge.cc index 203b37b..6739160 100644 --- a/exp_gauge.cc +++ b/exp_gauge.cc @@ -1,101 +1,13 @@ -#include "exp_gauge.hh" -#include "adjointfield.hh" -#include "su2.hh" -#include "u1.hh" -#include -#include -#include - -_su2 exp(adjointsu2 const &x) { - double a = x.geta(), b = x.getb(), c = x.getc(); - // normalise the vector - const double alpha = sqrt(a * a + b * b + c * c); - std::vector n = {a / alpha, b / alpha, c / alpha}; - - const double salpha = sin(alpha); - _su2 res(Complex(cos(alpha), salpha * n[2]), salpha * Complex(n[1], n[0])); - return res; -} - -// eq. (4) of https://arxiv.org/pdf/1508.00868v2.pdf -// the Gell-Mann decomposition coefficients have been computed with sympy -_su3 exp(const adjointsu3 &x) { - const double sqrt_of_3 = sqrt(3.0); // avoid computing it many times - - std::array arr = x.get_arr(); - double theta = 0.0; - for (size_t i = 0; i < 8; i++) { - theta += (arr[i] * arr[i]); - } - theta = std::sqrt(theta); - - // normalised vector - std::array n = arr; - for (size_t i = 0; i < 8; i++) { - n[i] /= theta; - } - - // determinant of H - const double detH = - 2.0 * sqrt_of_3 * (n[0] * n[0]) * n[7] / 3.0 + 2.0 * n[0] * n[3] * n[5] + - 2.0 * n[0] * n[4] * n[6] + 2.0 * sqrt_of_3 * (n[1] * n[1]) * n[7] / 3.0 - - 2.0 * n[1] * n[3] * n[6] + 2.0 * n[1] * n[4] * n[5] + - 2.0 * sqrt_of_3 * (n[2] * n[2]) * n[7] / 3.0 + n[2] * (n[3] * n[3]) + - n[2] * (n[4] * n[4]) - n[2] * (n[5] * n[5]) - n[2] * (n[6] * n[6]) - - sqrt_of_3 * (n[3] * n[3]) * n[7] / 3.0 - sqrt_of_3 * (n[4] * n[4]) * n[7] / 3.0 - - sqrt_of_3 * (n[5] * n[5]) * n[7] / 3.0 - sqrt_of_3 * (n[6] * n[6]) * n[7] / 3.0 - - 2.0 * sqrt_of_3 * (n[7] * n[7] * n[7]) / 9.0; - const double phi = (1.0 / 3.0) * (acos((3.0 / 2.0) * sqrt_of_3 * detH) - M_PI / 2.0); - - const Complex I(0.0, 1.0); // imaginary unit - - // 1st and 2nd row of H - const std::array H_1 = {n[2] + sqrt_of_3 * n[7] / 3.0, n[0] - I * n[1], - n[3] - I * n[4]}; - const std::array H_2 = {n[0] + I * n[1], -n[2] + sqrt_of_3 * n[7] / 3.0, - n[5] - I * n[6]}; - - // 1st and 2nd row of H^2 - const std::array H2_1 = { - (n[0] - I * n[1]) * (n[0] + I * n[1]) + std::pow(n[2] + sqrt_of_3 * n[7] / 3.0, 2.0) + - (n[3] - I * n[4]) * (n[3] + I * n[4]), - (n[0] - I * n[1]) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + - (n[0] - I * n[1]) * (n[2] + sqrt_of_3 * n[7] / 3.0) + - (n[3] - I * n[4]) * (n[5] + I * n[6]), - -2 * sqrt_of_3 * n[7] * (n[3] - I * n[4]) / 3.0 + - (n[0] - I * n[1]) * (n[5] - I * n[6]) + - (n[2] + sqrt_of_3 * n[7] / 3.0) * (n[3] - I * n[4])}; - const std::array H2_2 = { - (n[0] + I * n[1]) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + - (n[0] + I * n[1]) * (n[2] + sqrt_of_3 * n[7] / 3.0) + - (n[3] + I * n[4]) * (n[5] - I * n[6]), - (n[0] - I * n[1]) * (n[0] + I * n[1]) + - std::pow(-n[2] + sqrt_of_3 * n[7] / 3.0, 2.0) + - (n[5] - I * n[6]) * (n[5] + I * n[6]), - -2 * sqrt_of_3 * n[7] * (n[5] - I * n[6]) / 3.0 + - (n[0] + I * n[1]) * (n[3] - I * n[4]) + - (-n[2] + sqrt_of_3 * n[7] / 3.0) * (n[5] - I * n[6])}; - - std::array u = {0.0, 0.0, 0.0}, v = {0.0, 0.0, 0.0}; - for (size_t k = 0; k <= 2; k++) { - const double pk = phi + (2.0 / 3.0) * M_PI * k; - const double arg_num_k = (2.0 / sqrt_of_3) * theta * sin(pk); - const Complex num_k = std::exp(I * arg_num_k); - const double den_k = 1 - 2.0 * cos(2.0 * pk); - const Complex fact_k = num_k / den_k; - for (size_t i = 0; i < 3; i++) { - const Complex u_ik = H2_1[i] + (2.0 / sqrt_of_3) * H_1[i] * sin(pk); - u[i] += u_ik * fact_k; - - const Complex v_ik = H2_2[i] + (2.0 / sqrt_of_3) * H_2[i] * sin(pk); - - v[i] += v_ik * fact_k; - } - // contibution from the identity - u[0] += -(1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)) * fact_k; - v[1] += -(1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)) * fact_k; - } - - _su3 res(u, v); - return res; -} +/** + * @file exp_gauge.cc + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief exponentiation of su(2) and su(3) matrices (Lie algebra) + * @version 0.1 + * @date 2023-02-21 + * + * @copyright Copyright (c) 2023 + * + */ + +#include "exp_su2.cc" +#include "exp_su3.cc" diff --git a/exp_su2.cc b/exp_su2.cc new file mode 100644 index 0000000..b5fa109 --- /dev/null +++ b/exp_su2.cc @@ -0,0 +1,40 @@ +/** + * @file exp_su2.cc + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief exponentiation of an element of the su(2) algebra: A=alpha_k*sigma_k --> + * exp(i*alpha_k*sigma_k) + * @version 0.1 + * @date 2023-02-21 + * + * @copyright Copyright (c) 2023 + * + */ + +#include +#include +#include + +#include "adjointfield.hh" +#include "exp_gauge.hh" +#include "su2.hh" + +/** + * @brief exp(i*alpha_k*sigma_k), where the sigma_k are the Pauli matrices and alpha_k are + * real numbers. + * + * The code uses the explicit formula following by (4.23) of Gattringer&Lang + * https://link.springer.com/book/10.1007/978-3-642-01850-3 + * + * @param x + * @return _su2 + */ +_su2 exp(adjointsu2 const &x) { + double a = x.geta(), b = x.getb(), c = x.getc(); + // normalise the vector + const double alpha = sqrt(a * a + b * b + c * c); + std::vector n = {a / alpha, b / alpha, c / alpha}; + + const double salpha = sin(alpha); + _su2 res(Complex(cos(alpha), salpha * n[2]), salpha * Complex(n[1], n[0])); + return res; +} diff --git a/exp_su3.cc b/exp_su3.cc new file mode 100644 index 0000000..377bc7e --- /dev/null +++ b/exp_su3.cc @@ -0,0 +1,110 @@ +/** + * @file exp_su3.cc + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief exponentiation of an element of the su(3) algebra: A=a_k*sigma_k --> + * exp(i*a_k*lambda_k) + * @version 0.1 + * @date 2023-02-21 + * + * @copyright Copyright (c) 2023 + * + */ + +#include "adjointfield.hh" +#include "exp_gauge.hh" +#include "su3.hh" + +#include +#include +#include + +/** + * @brief exp(i*a_k*lambda_k), where the lambda_k are the Gell-Mann matrices and a_k are + * real numbers. + * + * eq. (4) of https://arxiv.org/pdf/1508.00868v2.pdf + * the Gell-Mann decomposition coefficients have been computed with sympy + * + * @param x + * @return _su2 + */ +_su3 exp(const adjointsu3 &x) { + const double sqrt_of_3 = sqrt(3.0); // avoid computing it many times + + std::array arr = x.get_arr(); + double theta = 0.0; + for (size_t i = 0; i < 8; i++) { + theta += (arr[i] * arr[i]); + } + theta = std::sqrt(theta); + + // normalised vector + std::array n = arr; + for (size_t i = 0; i < 8; i++) { + n[i] /= theta; + } + + // determinant of H + const double detH = + 2.0 * sqrt_of_3 * (n[0] * n[0]) * n[7] / 3.0 + 2.0 * n[0] * n[3] * n[5] + + 2.0 * n[0] * n[4] * n[6] + 2.0 * sqrt_of_3 * (n[1] * n[1]) * n[7] / 3.0 - + 2.0 * n[1] * n[3] * n[6] + 2.0 * n[1] * n[4] * n[5] + + 2.0 * sqrt_of_3 * (n[2] * n[2]) * n[7] / 3.0 + n[2] * (n[3] * n[3]) + + n[2] * (n[4] * n[4]) - n[2] * (n[5] * n[5]) - n[2] * (n[6] * n[6]) - + sqrt_of_3 * (n[3] * n[3]) * n[7] / 3.0 - sqrt_of_3 * (n[4] * n[4]) * n[7] / 3.0 - + sqrt_of_3 * (n[5] * n[5]) * n[7] / 3.0 - sqrt_of_3 * (n[6] * n[6]) * n[7] / 3.0 - + 2.0 * sqrt_of_3 * (n[7] * n[7] * n[7]) / 9.0; + const double phi = (1.0 / 3.0) * (acos((3.0 / 2.0) * sqrt_of_3 * detH) - M_PI / 2.0); + + const Complex I(0.0, 1.0); // imaginary unit + + // 1st and 2nd row of H + const std::array H_1 = {n[2] + sqrt_of_3 * n[7] / 3.0, n[0] - I * n[1], + n[3] - I * n[4]}; + const std::array H_2 = {n[0] + I * n[1], -n[2] + sqrt_of_3 * n[7] / 3.0, + n[5] - I * n[6]}; + + // 1st and 2nd row of H^2 + const std::array H2_1 = { + (n[0] - I * n[1]) * (n[0] + I * n[1]) + std::pow(n[2] + sqrt_of_3 * n[7] / 3.0, 2.0) + + (n[3] - I * n[4]) * (n[3] + I * n[4]), + (n[0] - I * n[1]) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + + (n[0] - I * n[1]) * (n[2] + sqrt_of_3 * n[7] / 3.0) + + (n[3] - I * n[4]) * (n[5] + I * n[6]), + -2 * sqrt_of_3 * n[7] * (n[3] - I * n[4]) / 3.0 + + (n[0] - I * n[1]) * (n[5] - I * n[6]) + + (n[2] + sqrt_of_3 * n[7] / 3.0) * (n[3] - I * n[4])}; + const std::array H2_2 = { + (n[0] + I * n[1]) * (-n[2] + sqrt_of_3 * n[7] / 3.0) + + (n[0] + I * n[1]) * (n[2] + sqrt_of_3 * n[7] / 3.0) + + (n[3] + I * n[4]) * (n[5] - I * n[6]), + (n[0] - I * n[1]) * (n[0] + I * n[1]) + + std::pow(-n[2] + sqrt_of_3 * n[7] / 3.0, 2.0) + + (n[5] - I * n[6]) * (n[5] + I * n[6]), + -2 * sqrt_of_3 * n[7] * (n[5] - I * n[6]) / 3.0 + + (n[0] + I * n[1]) * (n[3] - I * n[4]) + + (-n[2] + sqrt_of_3 * n[7] / 3.0) * (n[5] - I * n[6])}; + + std::array u = {0.0, 0.0, 0.0}, v = {0.0, 0.0, 0.0}; + for (size_t k = 0; k <= 2; k++) { + const double pk = phi + (2.0 / 3.0) * M_PI * k; + const double arg_num_k = (2.0 / sqrt_of_3) * theta * sin(pk); + const Complex num_k = std::exp(I * arg_num_k); + const double den_k = 1 - 2.0 * cos(2.0 * pk); + const Complex fact_k = num_k / den_k; + for (size_t i = 0; i < 3; i++) { + const Complex u_ik = H2_1[i] + (2.0 / sqrt_of_3) * H_1[i] * sin(pk); + u[i] += u_ik * fact_k; + + const Complex v_ik = H2_2[i] + (2.0 / sqrt_of_3) * H_2[i] * sin(pk); + + v[i] += v_ik * fact_k; + } + // contibution from the identity + u[0] += -(1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)) * fact_k; + v[1] += -(1.0 / 3.0) * (1.0 + 2.0 * cos(2.0 * pk)) * fact_k; + } + + _su3 res(u, v); + return res; +} From 40add1fe1ec52d91d4722d642ca27900a76185a2 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Tue, 21 Feb 2023 17:25:27 +0100 Subject: [PATCH 11/13] DO NOT MERGE still understanding how to implement get_deriv of SU(3) --- include/adjoint_su3.hh | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh index 03f5b69..9135822 100644 --- a/include/adjoint_su3.hh +++ b/include/adjoint_su3.hh @@ -79,13 +79,21 @@ private: * @param A * @return adjointsu3 */ -template inline adjointsu3 get_deriv(const su3 &A) { -// std::cout << "ERROR: please implement " << __func__ << " for SU(3) \n"; -// std::abort(); - const su3 A_thh = traceless_antiherm(A); - const std::array u = A_thh.get_u(); - const std::array v = A_thh.get_v(); - const std::array arr = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // = ??? +template inline adjointsu3 get_deriv(const su3 &U) { + const std::array u = U.get_u(); + const std::array v = U.get_v(); + + std::array arr; + +arr[0] = std::imag(u[1]) + std::imag(v[0]) ; +arr[1] = std::real(u[1]) - std::real(v[0]) ; +arr[2] = std::imag(u[0]) - std::imag(v[1]) ; +arr[3] = -std::real(u[1])*std::imag(v[2]) + std::real(u[2])*std::imag(v[1]) + std::real(v[1])*std::imag(u[2]) - std::real(v[2])*std::imag(u[1]) + std::imag(u[2]) ; +arr[4] = -std::real(u[1])*std::real(v[2]) + std::real(u[2])*std::real(v[1]) + std::real(u[2]) + std::imag(u[1])*std::imag(v[2]) - std::imag(u[2])*std::imag(v[1]) ; +arr[5] = std::real(u[0])*std::imag(v[2]) - std::real(u[2])*std::imag(v[0]) - std::real(v[0])*std::imag(u[2]) + std::real(v[2])*std::imag(u[0]) + std::imag(v[2]) ; +arr[6] = std::real(u[0])*std::real(v[2]) - std::real(u[2])*std::real(v[0]) + std::real(v[2]) - std::imag(u[0])*std::imag(v[2]) + std::imag(u[2])*std::imag(v[0]) ; +arr[7] = 2*sqrt(3.0)*std::real(u[0])*std::imag(v[1])/3.0 - 2*sqrt(3.0)*std::real(u[1])*std::imag(v[0])/3.0 - 2*sqrt(3.0)*std::real(v[0])*std::imag(u[1])/3.0 + 2*sqrt(3.0)*std::real(v[1])*std::imag(u[0])/3.0 + sqrt(3.0)*std::imag(u[0])/3.0 + sqrt(3.0)*std::imag(v[1])/3.0 ; + return adjointsu3(arr); } From a544b9beba6d52dfb55cbd6030341bb098187232 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Fri, 24 Feb 2023 19:13:02 +0100 Subject: [PATCH 12/13] clearer definition and use of accum types --- include/accum_type.hh | 22 +++- include/adjoint_su2.hh | 4 +- include/adjoint_su3.hh | 33 +++--- include/adjoint_u1.hh | 2 +- include/gaugeconfig.hh | 1 + include/gaugemonomial.hh | 2 + include/get_staples.hh | 4 +- include/integrator.hh | 2 + include/md_update.hh | 52 +++++---- include/obc_gaugemonomial.hh | 4 +- include/random_element.hh | 21 +++- include/smearape.hh | 8 +- include/su2.hh | 41 +++++++ include/su3.hh | 81 +++----------- include/su3_accum.hh | 204 +++++++++++++++++++++++++++++++++++ include/u1.hh | 44 +++++--- 16 files changed, 392 insertions(+), 133 deletions(-) create mode 100644 include/su3_accum.hh diff --git a/include/accum_type.hh b/include/accum_type.hh index 5c6384c..c546c27 100644 --- a/include/accum_type.hh +++ b/include/accum_type.hh @@ -1,5 +1,23 @@ #pragma once -template struct accum_type { - typedef T type; +/** + * @brief generic accumulation type for elements of the group + * + * accumulations of group elements (e.g. sums of plaquettes) don't belong to the group + * anymore. This the "type" attribute of this class defines a templated class representing + * these quantities. + * + * By default, the same class used for the elements of the group is used. It is the case + * of SU(2). However it becomes more expensive with U(1) and other SU(N) in general. + * + * The specialization of this templated struct have to be written + * explicitly for these other groups. + * See the implementation of U(1) and SU(3) accum types for an example. + * + * @tparam Group + */ +template struct accum_type { + typedef Group type; }; + + diff --git a/include/adjoint_su2.hh b/include/adjoint_su2.hh index 8abf345..0e16eda 100644 --- a/include/adjoint_su2.hh +++ b/include/adjoint_su2.hh @@ -61,8 +61,8 @@ private: Float a, b, c; }; - -template inline adjointsu2 get_deriv(su2 &A) { +template +inline adjointsu2 get_deriv(const su2_accum &A) { const Complex a = A.geta(), b = A.getb(); return adjointsu2(2. * std::imag(b), 2. * std::real(b), 2. * std::imag(a)); } diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh index 9135822..b2f0578 100644 --- a/include/adjoint_su3.hh +++ b/include/adjoint_su3.hh @@ -13,6 +13,7 @@ #include "geometry.hh" #include "su3.hh" +#include "su3_accum.hh" #include #include @@ -72,27 +73,33 @@ private: }; /** - * @brief returns (-i)*(B - B.dagger()), where B = A - (tr(a)/N_c) * 1, parametrized - * as an element of su(3) + * @brief returns A = B - (tr(B)/N_c) * 1_{3x3}, where B = U - U^{\dagger}. + * The matrix A is returned as an element of the algebra su(3), i.e. giving the + * coefficients in front of the Gell-Mall matrices * * @tparam Float * @param A * @return adjointsu3 */ -template inline adjointsu3 get_deriv(const su3 &U) { - const std::array u = U.get_u(); - const std::array v = U.get_v(); +template inline adjointsu3 get_deriv(const su3_accum &U) { + su3_accum A = U - U.dagger(); + _su3 Id; // by default is the identity + A = A - (trace(A) / double(U.N_c)) * Id; + const std::array u = A.get_u(); + const std::array v = A.get_v(); + const std::array w = A.get_w(); std::array arr; -arr[0] = std::imag(u[1]) + std::imag(v[0]) ; -arr[1] = std::real(u[1]) - std::real(v[0]) ; -arr[2] = std::imag(u[0]) - std::imag(v[1]) ; -arr[3] = -std::real(u[1])*std::imag(v[2]) + std::real(u[2])*std::imag(v[1]) + std::real(v[1])*std::imag(u[2]) - std::real(v[2])*std::imag(u[1]) + std::imag(u[2]) ; -arr[4] = -std::real(u[1])*std::real(v[2]) + std::real(u[2])*std::real(v[1]) + std::real(u[2]) + std::imag(u[1])*std::imag(v[2]) - std::imag(u[2])*std::imag(v[1]) ; -arr[5] = std::real(u[0])*std::imag(v[2]) - std::real(u[2])*std::imag(v[0]) - std::real(v[0])*std::imag(u[2]) + std::real(v[2])*std::imag(u[0]) + std::imag(v[2]) ; -arr[6] = std::real(u[0])*std::real(v[2]) - std::real(u[2])*std::real(v[0]) + std::real(v[2]) - std::imag(u[0])*std::imag(v[2]) + std::imag(u[2])*std::imag(v[0]) ; -arr[7] = 2*sqrt(3.0)*std::real(u[0])*std::imag(v[1])/3.0 - 2*sqrt(3.0)*std::real(u[1])*std::imag(v[0])/3.0 - 2*sqrt(3.0)*std::real(v[0])*std::imag(u[1])/3.0 + 2*sqrt(3.0)*std::real(v[1])*std::imag(u[0])/3.0 + sqrt(3.0)*std::imag(u[0])/3.0 + sqrt(3.0)*std::imag(v[1])/3.0 ; + arr[0] = std::imag(u[1]) + std::imag(v[0]); + arr[1] = std::real(u[1]) - std::real(v[0]); + arr[2] = std::imag(u[0]) - std::imag(v[1]); + arr[3] = std::imag(u[2]) + std::imag(w[0]); + arr[4] = std::real(u[2]) - std::real(w[0]); + arr[5] = std::imag(v[2]) + std::imag(w[1]); + arr[6] = std::real(v[2]) - std::real(w[1]); + arr[7] = sqrt(3.0) * std::imag(u[0]) / 3.0 + sqrt(3.0) * std::imag(v[1]) / 3.0 - + 2.0 * sqrt(3.0) * std::imag(w[2]) / 3.0; return adjointsu3(arr); } diff --git a/include/adjoint_u1.hh b/include/adjoint_u1.hh index 90bc155..2b67d20 100644 --- a/include/adjoint_u1.hh +++ b/include/adjoint_u1.hh @@ -51,7 +51,7 @@ private: Float a; }; -template inline adjointu1 get_deriv(Complex &A) { +template inline adjointu1 get_deriv(const u1_accum &A) { return adjointu1(std::imag(A)); } diff --git a/include/gaugeconfig.hh b/include/gaugeconfig.hh index bada64d..2dc92a5 100644 --- a/include/gaugeconfig.hh +++ b/include/gaugeconfig.hh @@ -14,6 +14,7 @@ #include "geometry.hh" #include "random_element.hh" +#include "su3.hh" #include "su2.hh" #include "u1.hh" diff --git a/include/gaugemonomial.hh b/include/gaugemonomial.hh index 5d307cb..db912c7 100644 --- a/include/gaugemonomial.hh +++ b/include/gaugemonomial.hh @@ -19,6 +19,8 @@ #include "hamiltonian_field.hh" #include "monomial.hh" #include "su2.hh" +#include "su3.hh" +#include "su3_accum.hh" #include "u1.hh" #include diff --git a/include/get_staples.hh b/include/get_staples.hh index aea2440..6b8157f 100644 --- a/include/get_staples.hh +++ b/include/get_staples.hh @@ -94,9 +94,9 @@ RetType get_staples_down(gaugeconfig &U, const Arr &x, const size_t &mu) * @param anisotropic boolean flag: true when considering anisotropy * @param spatial_only boolean flag: true when summing over spatial staples only */ -template +template void get_staples_MCMC_step(T &K, - gaugeconfig &U, + gaugeconfig &U, Arr const x, const size_t mu, const double xi = 1.0, diff --git a/include/integrator.hh b/include/integrator.hh index 1d4c297..577c0e9 100644 --- a/include/integrator.hh +++ b/include/integrator.hh @@ -46,6 +46,8 @@ public: update_momenta(monomial_list, deriv, h, dtau/2.); // restore SU if(restore) h.U->restoreSU(); + + } }; diff --git a/include/md_update.hh b/include/md_update.hh index e4737d5..76f88dd 100644 --- a/include/md_update.hh +++ b/include/md_update.hh @@ -1,26 +1,26 @@ #pragma once -#include"gaugeconfig.hh" -#include"adjointfield.hh" -#include"hamiltonian_field.hh" -#include"monomial.hh" -#include"md_params.hh" -#include"integrator.hh" -#include -#include -#include -#include -#include +#include "adjointfield.hh" +#include "gaugeconfig.hh" +#include "hamiltonian_field.hh" +#include "integrator.hh" +#include "md_params.hh" +#include "monomial.hh" +#include +#include +#include +#include using std::vector; - -template void md_update(gaugeconfig &U, - URNG &engine, - md_params ¶ms, - std::list*> &monomial_list, - integrator &md_integ) { - adjointfield momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), U.getndims()); +template +void md_update(gaugeconfig &U, + URNG &engine, + md_params ¶ms, + std::list *> &monomial_list, + integrator &md_integ) { + adjointfield momenta(U.getLx(), U.getLy(), U.getLz(), U.getLt(), + U.getndims()); // generate standard normal distributed random momenta // normal distribution checked! initnormal(engine, momenta); @@ -31,7 +31,7 @@ template void md_update(gaugeconfigheatbath(h); + (*it)->heatbath(h); } // keep a copy of original gauge field @@ -45,28 +45,28 @@ template void md_update(gaugeconfigaccept(h); + (*it)->accept(h); delta_H += (*it)->getDeltaH(); } params.setdeltaH(delta_H); // accept/reject step, if needed params.setaccept(true); - if(delta_H > 0) { - if(uniform(engine) > exp(-delta_H)) { + if (delta_H > 0) { + if (uniform(engine) > exp(-delta_H)) { params.setaccept(false); } } // if wanted, perform a reversibility violation test. - if(params.getrevtest()) { + if (params.getrevtest()) { delta_H = 0.; gaugeconfig U_save(U); h.momenta->flipsign(); md_integ.integrate(monomial_list, h, params); for (auto it = monomial_list.begin(); it != monomial_list.end(); it++) { - (*it)->accept(h); + (*it)->accept(h); delta_H += (*it)->getDeltaH(); } params.setdeltadeltaH(delta_H); @@ -74,10 +74,8 @@ template void md_update(gaugeconfig(S); + deriv(x, mu) += num_fact_i * get_deriv(S); xpmu[mu]--; // x } diff --git a/include/random_element.hh b/include/random_element.hh index 4204a0f..83d87cb 100644 --- a/include/random_element.hh +++ b/include/random_element.hh @@ -36,15 +36,30 @@ void random_element(_su2 &U, URNG &engine, const double delta = 1.) { return; } +/** + * @brief initialized the configuration to some random SU(3) matrix + * + * In the (u,v) representation of U (see eq. 4.26 of Gattringer&Lang), there's no unique + * way of defining "v", because the orthogonal space to a vector "u" is 2-dimensional. + * Here we adopt our custom "u_2"-convention, according to which: + * v = (u_2, -u_1-u_3, u_2)^{\dagger} + * At the end we normalize the vectors using the restoreSU() method. + * + * @tparam URNG + * @param U + * @param engine + * @param delta + */ template void random_element(_su3 &U, URNG &engine, const double delta = 1.) { + std::uniform_real_distribution dist(-delta, delta); - std::array u, v; + std::array u; for (size_t i = 0; i < U.N_c; i++) { u[i] = dist(engine); - v[i] = dist(engine); } - + std::array v = {std::conj(u[1]), -std::conj(u[0]) - std::conj(u[2]), + std::conj(u[1])}; U = _su3(u, v); U.restoreSU(); diff --git a/include/smearape.hh b/include/smearape.hh index 430b2d1..e33e53e 100644 --- a/include/smearape.hh +++ b/include/smearape.hh @@ -5,6 +5,9 @@ #include "random_element.hh" #include "su2.hh" #include "u1.hh" +#include "su3.hh" +#include "su3_accum.hh" + #include #include #include @@ -114,9 +117,10 @@ void APEsmearing(gaugeconfig &U, const double &alpha, const bool spatial= // K is intialized to (0,0) even if not explicitly specified accum K; get_staples_APE(K, Uold, x, i, spatial); - const Group Uprime(alpha*Uold(x, i) + beta*K); + K = alpha*accum(Uold(x, i)) + beta*K; + const Group Uprime = accum_to_Group(K); U(x, i) = Uprime; - U(x, i).restoreSU(); //only necessary for su2 + U(x, i).restoreSU(); } } } diff --git a/include/su2.hh b/include/su2.hh index 117881b..cbcfdc5 100644 --- a/include/su2.hh +++ b/include/su2.hh @@ -1,3 +1,14 @@ +/** + * @file su2.hh + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief class representing an SU(2) matrix in the fundamental representation + * @version 0.1 + * @date 2023-02-24 + * + * @copyright Copyright (c) 2023 + * + */ + #pragma once #include @@ -68,6 +79,13 @@ public: b /= r; } + void print(){ + std::cout << "--------------------\n"; + std::cout << a << " " << b << "\n"; + std::cout << -std::conj(b) << " " << std::conj(a) << "\n"; + std::cout << "--------------------\n"; + } + private: Complex a, b; }; @@ -128,4 +146,27 @@ inline _su2 operator*(const _su2 &U1, const Complex &U2) { return (res); } + +using su2_accum = _su2; + +inline _su2 accum_to_Group(const su2_accum &x) { + _su2 U(x.geta(), x.getb()); + U.restoreSU(); + return U; +} + +/** + * @brief accumulation type for SU(2) matrices + * + * Incidentally, for SU(2) linear combinations of products its matrices can be still + * parametrized by 2 numbers only. Therefore we use the same class for SU(2) elements and + * their accumulations. + * + * @tparam + */ + +template <> struct accum_type<_su2> { + typedef _su2 type; +}; + using su2 = _su2; diff --git a/include/su3.hh b/include/su3.hh index c65bf4c..4b9aa36 100644 --- a/include/su3.hh +++ b/include/su3.hh @@ -22,7 +22,7 @@ using Complex = std::complex; /** - * @brief representation of an SU(2) matrix in the fundamental representation with the + * @brief representation of an SU(3) matrix in the fundamental representation with the * convention of eq. (4.26) of Gattringer&Lang * https://link.springer.com/book/10.1007/978-3-642-01850-3 * @@ -41,8 +41,8 @@ public: friend inline _su3 operator+(const _su3 &U1, const _su3 &U2); friend inline _su3 operator-(const _su3 &U1, const _su3 &U2); - friend inline _su3 operator*(const _su3 &U1, const _su3 &U2); - friend inline _su3 operator*(const Complex &U1, const _su3 &U2); + // friend inline _su3 operator*(const _su3 &U1, const _su3 &U2); + // friend inline _su3 operator*(const Complex &U1, const _su3 &U2); friend inline _su3 operator*(const _su3 &U1, const Complex &U2); _su3 &operator*=(const _su3 &U2) { const std::array u1 = (*this).u; @@ -84,22 +84,18 @@ public: inline std::array get_u() const { return u; } inline std::array get_v() const { return v; } + inline std::array get_w() const { + std::array w = {u[1] * v[2] - v[1] * u[2], -u[0] * v[2] + v[0] * u[2], + u[0] * v[1] - v[0] * u[1]}; + w = {std::conj(w[0]), std::conj(w[1]), std::conj(w[2])}; + return w; + } + inline void operator=(const _su3 &U) { u = U.get_u(); v = U.get_v(); } - inline void operator+=(const _su3 &U) { - for (size_t i = 0; i < N_c; i++) { - u[i] += U.u[i]; - v[i] += U.v[i]; - } - } - inline void operator-=(const _su3 &U) { - for (size_t i = 0; i < N_c; i++) { - u[i] -= U.u[i]; - v[i] -= U.v[i]; - } - } + void set(const std::array &_u, const std::array &_v) { u = _u; v = _v; @@ -119,8 +115,7 @@ public: } inline Complex trace() const { // expression computed using sympy - const Complex res = - u[0] + v[1] + std::conj(u[0]) * std::conj(v[1]) - std::conj(u[1]) * std::conj(v[0]); + const Complex res = u[0] + v[1] + std::conj(u[0] * v[1] - u[1] * v[0]); return res; } inline double retrace() const { return std::real(this->trace()); } @@ -174,8 +169,7 @@ public: } void print() const { - std::array w = { - u[1] * v[2] - v[1] * u[2], -u[0] * v[2] + v[0] * u[2], u[0] * v[1] - v[0] * u[1]}; + std::array w = this->get_w(); w = {std::conj(w[0]), std::conj(w[1]), std::conj(w[2])}; std::cout << "----------------------------------\n"; @@ -207,53 +201,4 @@ inline _su3 operator*(const _su3 &U1, const _su3 &U2) { return U; } -inline _su3 operator+(const _su3 &U1, const _su3 &U2) { - _su3 U = U1; - U += U2; - return U; -} - -inline _su3 operator-(const _su3 &U1, const _su3 &U2) { - _su3 U = U1; - U -= U2; - return U; -} - -inline _su3 operator*(const double &z, const _su3 &U) { - const size_t N_c = U.N_c; - std::array u, v; - for (size_t i = 0; i < N_c; i++) { - u[i] *= double(z); - v[i] *= double(z); - } - return _su3(u, v); -} - -inline _su3 operator*(const _su3 &U, const double &z) { - return z * U; -} - -inline _su3 operator*(const Complex &z, const _su3 &U) { - const size_t N_c = U.N_c; - std::array u, v; - for (size_t i = 0; i < N_c; i++) { - u[i] *= Complex(z); - v[i] *= Complex(z); - } - return _su3(u, v); -} - -inline _su3 operator*(const _su3 &U, const Complex &z) { - return z * U; -} - -template <> inline _su3 traceless_antiherm(const _su3 &x0) { - const _su3 Id({1.0, 0.0, 0.0}, - {0.0, 1.0, 0.0}); // default constructor: identity operator - _su3 x = x0; - x = x - (trace(x) / double(x.N_c)) * Id; - x = 0.5 * (x - x.dagger()); - return x; -} - using su3 = _su3; diff --git a/include/su3_accum.hh b/include/su3_accum.hh new file mode 100644 index 0000000..0ba95a5 --- /dev/null +++ b/include/su3_accum.hh @@ -0,0 +1,204 @@ +/** + * @file su3_accum.hh + * @author Simone Romiti (simone.romiti@uni-bonn.de) + * @brief accumulation type for SU(3) matrices + * @version 0.1 + * @date 2023-02-24 + * + * @copyright Copyright (c) 2023 + * + */ + +#pragma once + +#include +#include +#include + +#include "accum_type.hh" +#include "dagger.hh" +#include "su3.hh" +#include "traceless_antiherm.hh" + +/** + * @brief class describing a 3x3 matrix with rows (u, v, w). Its used to accumulate sums + * and differences of plaquettes and/or staples (which is not unitary anymore) + * + */ +class _su3_accum { +public: + const size_t N_c = 3; + explicit _su3_accum() { + this->set_to_zero(); // default: zero matrix + } + explicit _su3_accum(const std::array &_u, + const std::array &_v, + const std::array &_w) + : u(_u), v(_v), w(_w) {} + _su3_accum(const _su3_accum &U) : u(U.get_u()), v(U.get_v()), w(U.get_w()) {} + _su3_accum(const _su3 &U) : u(U.get_u()), v(U.get_v()), w(U.get_w()) {} + + friend inline _su3_accum operator+(const _su3_accum &U1, const _su3_accum &U2); + friend inline _su3_accum operator-(const _su3_accum &U1, const _su3_accum &U2); + friend inline _su3 operator*(const _su3 &U1, const _su3 &U2); + + template _su3_accum &operator*=(const matr_33 &U2) { + const std::array u1 = (*this).u, v1 = (*this).v, w1 = (*this).w; + const std::array u2 = U2.get_u(), v2 = U2.get_v(), w2 = U2.get_w(); + for (size_t i = 0; i < this->N_c; i++) { + (*this).u[i] = u1[0] * u2[i] + u1[1] * v2[i] + u1[2] * w2[i]; + (*this).v[i] = v1[0] * u2[i] + v1[1] * v2[i] + v1[2] * w2[i]; + (*this).w[i] = w1[0] * u2[i] + w1[1] * v2[i] + w1[2] * w2[i]; + } + return *this; + } + + inline std::array get_u() const { return u; } + inline std::array get_v() const { return v; } + inline std::array get_w() const { return w; } + + inline void operator=(const _su3 &U) { + u = U.get_u(); + v = U.get_v(); + w = U.get_w(); + } + + inline void operator=(const _su3_accum &U) { + u = U.get_u(); + v = U.get_v(); + w = U.get_w(); + } + + inline void operator+=(const _su3_accum &U) { + for (size_t i = 0; i < N_c; i++) { + u[i] += U.u[i]; + v[i] += U.v[i]; + w[i] += U.w[i]; + } + } + inline void operator-=(const _su3_accum &U) { + for (size_t i = 0; i < N_c; i++) { + u[i] -= U.u[i]; + v[i] -= U.v[i]; + w[i] -= U.w[i]; + } + } + void set(const std::array &_u, + const std::array &_v, + const std::array &_w) { + u = _u; + v = _v; + w = _w; + } + void set_to_zero() { + u = {0.0, 0.0, 0.0}; + v = u; + w = u; + } + inline _su3_accum dagger() const { + // vectors computed using sympy + const std::array u2 = {std::conj(u[0]), std::conj(v[0]), std::conj(w[0])}, + v2 = {std::conj(u[1]), std::conj(v[1]), std::conj(w[1])}, + w2 = {std::conj(u[2]), std::conj(v[2]), std::conj(w[2])}; + return _su3_accum(u2, v2, w2); + } + inline Complex trace() const { return (u[0] + v[1] + w[2]); } + inline double retrace() const { return std::real(this->trace()); } + + void print() const { + std::cout << "----------------------------------\n"; + std::cout << u[0] << " " << u[1] << " " << u[2] << "\n"; + std::cout << v[0] << " " << v[1] << " " << v[2] << "\n"; + std::cout << w[0] << " " << w[1] << " " << w[2] << "\n"; + std::cout << "----------------------------------\n"; + } + +private: + std::array u, v, w; +}; + +inline _su3_accum operator+(const _su3_accum &U1, const _su3_accum &U2) { + _su3_accum U = U1; + U += U2; + return U; +} + +inline _su3_accum operator-(const _su3_accum &U1, const _su3_accum &U2) { + _su3_accum U = U1; + U -= U2; + return U; +} + +inline Complex trace(_su3_accum const &U) { + return U.trace(); +} + +inline double retrace(_su3_accum const &U) { + return U.retrace(); +} + +template <> inline _su3_accum dagger(const _su3_accum &u) { + return u.dagger(); +} + +inline _su3_accum operator*(const _su3_accum &U1, const _su3_accum &U2) { + _su3_accum U = U1; + U *= U2; + return U; +} + +inline _su3_accum operator*(const double &z, const _su3_accum &U) { + std::array u = U.get_u(), v = U.get_v(), w = U.get_w(); + for (size_t i = 0; i < U.N_c; i++) { + u[i] *= double(z); + v[i] *= double(z); + w[i] *= double(z); + } + return _su3_accum(u, v, w); +} + +inline _su3_accum operator*(const _su3_accum &U, const double &z) { + return z * U; +} + +inline _su3_accum operator*(const Complex &z, const _su3_accum &U) { + std::array u = U.get_u(), v = U.get_v(), w = U.get_w(); + for (size_t i = 0; i < U.N_c; i++) { + u[i] *= Complex(z); + v[i] *= Complex(z); + w[i] *= Complex(z); + } + return _su3_accum(u, v, w); +} + +inline _su3_accum operator*(const _su3_accum &U, const Complex &z) { + return z * U; +} + +template <> inline _su3_accum traceless_antiherm(const _su3_accum &x0) { + _su3 Id; // by default is the identity + _su3_accum x = x0; + x = x - (trace(x) / double(x.N_c)) * Id; + x = 0.5 * (x - x.dagger()); + return x; +} + +using su3_accum = _su3_accum; + +inline su3 accum_to_Group(const su3_accum &x) { + su3 U(x.get_u(), x.get_v()); + U.restoreSU(); + return U; +} + +/** + * @brief accumulation type for SU(3) + * + * when summing SU(3) matrices the result is not unitary anymore. This type serves to + * define the sum of the plaquettes appearing in the action and monomial forces (covariant + * derivatives) + * + */ +template <> struct accum_type { + typedef su3_accum type; +}; diff --git a/include/u1.hh b/include/u1.hh index 5146f2c..d035bf9 100644 --- a/include/u1.hh +++ b/include/u1.hh @@ -35,13 +35,19 @@ public: double geta() const { return (a); } void operator=(const _u1 &U) { a = U.a; } - void set(const double& _a) { a = _a; } + void set(const double &_a) { a = _a; } void set_to_identity() { a = 0; } _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() {} + void print() { + std::cout << "---\n"; + std::cout << std::exp(a * std::complex(0.0, 1.0)) << "\n"; + std::cout << "---\n"; + } + private: double a; }; @@ -58,10 +64,6 @@ inline Complex trace(const Complex c) { return (c); // for U(1) the trace operator acts trivially } -template <> struct accum_type<_u1> { - typedef Complex type; -}; - template <> inline Complex dagger(const Complex &x) { return std::conj(x); } @@ -70,18 +72,17 @@ template <> inline _u1 dagger(const _u1 &x) { return x.dagger(); } - 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) { return _u1(U1.a + U2.a); @@ -110,4 +111,23 @@ inline void operator*=(Complex &U1, const _u1 &U2) { U1 *= std::exp(U2.geta() * Complex(0., 1.)); } +using u1_accum = Complex; + +inline _u1 accum_to_Group(const u1_accum &x) { + _u1 U(x); + U.restoreSU(); + return U; +} + +/** + * @brief accumulation type for U(1) + * + * sums of phases (U(1) elements) is just a complex numbers. + * No need to implement a separate class for it, we use std::complex. + * + */ +template <> struct accum_type<_u1> { + typedef Complex type; +}; + using u1 = _u1; From 1ff398d43603fb6099958bb288f50255601d2211 Mon Sep 17 00:00:00 2001 From: simone-romiti Date: Fri, 24 Feb 2023 20:04:12 +0100 Subject: [PATCH 13/13] fixed get_deriv finally get acceptance rate = 1 for SU(3) run. Check the plaquette, is low in value, maybe the normalization is wrong. --- include/adjoint_su3.hh | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/include/adjoint_su3.hh b/include/adjoint_su3.hh index b2f0578..b4469bc 100644 --- a/include/adjoint_su3.hh +++ b/include/adjoint_su3.hh @@ -75,19 +75,16 @@ private: /** * @brief returns A = B - (tr(B)/N_c) * 1_{3x3}, where B = U - U^{\dagger}. * The matrix A is returned as an element of the algebra su(3), i.e. giving the - * coefficients in front of the Gell-Mall matrices + * coefficients in front of the Gell-Mall matrices. The expression for them has been computed with sympy. * * @tparam Float * @param A * @return adjointsu3 */ template inline adjointsu3 get_deriv(const su3_accum &U) { - su3_accum A = U - U.dagger(); - _su3 Id; // by default is the identity - A = A - (trace(A) / double(U.N_c)) * Id; - const std::array u = A.get_u(); - const std::array v = A.get_v(); - const std::array w = A.get_w(); + const std::array u = U.get_u(); + const std::array v = U.get_v(); + const std::array w = U.get_w(); std::array arr; @@ -101,6 +98,8 @@ template inline adjointsu3 get_deriv(const su3_acc arr[7] = sqrt(3.0) * std::imag(u[0]) / 3.0 + sqrt(3.0) * std::imag(v[1]) / 3.0 - 2.0 * sqrt(3.0) * std::imag(w[2]) / 3.0; + // std::abort(); + return adjointsu3(arr); }