$SU(3)$ support

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

exp_gauge.cc

@@ -1,20 +1,13 @@
-#include"exp_gauge.hh"
-#include"su2.hh"
-#include"u1.hh"
-#include"adjointfield.hh"
-#include<vector>
-#include<complex>
-
-_su2 exp(adjointsu2<double> 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<double> 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;
-}
-
-
+/**
+ * @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"

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 <cmath>
+#include <complex>
+#include <vector>
+
+#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<double> 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<double> 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;
+}

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 <cmath>
+#include <complex>
+#include <vector>
+
+/**
+ * @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<double> &x) {
+  const double sqrt_of_3 = sqrt(3.0); // avoid computing it many times
+
+  std::array<double, 8> 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<double, 8> 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<Complex, 3> H_1 = {n[2] + sqrt_of_3 * n[7] / 3.0, n[0] - I * n[1],
+                                      n[3] - I * n[4]};
+  const std::array<Complex, 3> 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<Complex, 3> 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<Complex, 3> 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<Complex, 3> 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;
+}

include/accum_type.hh

@@ -1,5 +1,23 @@
 #pragma once
 
-template <class T> 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 <class Group> struct accum_type {
+  typedef Group type;
 };
+
+

include/adjoint_su2.hh

@@ -0,0 +1,73 @@
+/**
+ * @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 <array>
+#include <cassert>
+#include <cmath>
+#include <random>
+#include <vector>
+
+template <typename Float> 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<Float> 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 <typename Float = double>
+inline adjointsu2<Float> get_deriv(const su2_accum &A) {
+  const Complex a = A.geta(), b = A.getb();
+  return adjointsu2<Float>(2. * std::imag(b), 2. * std::real(b), 2. * std::imag(a));
+}
+
+template <typename Float>
+inline adjointsu2<Float> operator*(const Float &x, const adjointsu2<Float> &A) {
+  return adjointsu2<Float>(x * A.geta(), x * A.getb(), x * A.getc());
+}