From a77dc0fcf7ca21f041a7e084425dc49dbf112fea Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Thu, 28 Jul 2022 13:55:34 +0200
Subject: [PATCH 01/17] sparse domain improved memory footprint

---
 lib/gpt/core/domain/sparse.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/lib/gpt/core/domain/sparse.py b/lib/gpt/core/domain/sparse.py
index 7022afd7..fae4e919 100644
--- a/lib/gpt/core/domain/sparse.py
+++ b/lib/gpt/core/domain/sparse.py
@@ -28,7 +28,7 @@ def get_cache(c, l):
 
 
 class sparse_kernel:
-    def __init__(self, grid, local_coordinates):
+    def __init__(self, grid, local_coordinates, divisible_by):
         assert grid.cb.n == 1
         self.grid = grid
         self.local_coordinates = local_coordinates
@@ -37,9 +37,9 @@ def __init__(self, grid, local_coordinates):
         n = len(local_coordinates)
         N = self.grid.Nprocessors
         l = np.zeros(N, dtype=np.uint64)
-        l[self.grid.processor] = 2 ** int(np.ceil(np.log(n) / np.log(2)))
+        l[self.grid.processor] = 2 ** int(np.ceil(np.log(n / divisible_by) / np.log(2)))
         l = grid.globalsum(l)
-        self.L = [int(np.max(l)) * self.grid.mpi[0]] + self.grid.mpi[1:]
+        self.L = [int(np.max(l) * divisible_by) * self.grid.mpi[0]] + self.grid.mpi[1:]
 
         cb_simd_only_first_dimension = gpt.general(1, [0] * grid.nd, [1] + [0] * (grid.nd - 1))
 
@@ -130,12 +130,12 @@ def coordinate_lattices(self, mark_empty=0):
 
 
 class sparse:
-    def __init__(self, grid, local_coordinates):
+    def __init__(self, grid, local_coordinates, divisible_by=1):
         self.local_coordinates = local_coordinates
         self.grid = grid
 
         # kernel to avoid circular references through captures below
-        kernel = sparse_kernel(grid, local_coordinates)
+        kernel = sparse_kernel(grid, local_coordinates, divisible_by)
 
         def _project(dst, src):
             for d, s in zip(dst, src):

From 9a3ea848254c60dbc37001cb67d744109ce135ee Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Tue, 23 Aug 2022 14:53:25 +0200
Subject: [PATCH 02/17] extended functionality of divisible_by in sparse domain

---
 lib/gpt/core/domain/sparse.py | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/lib/gpt/core/domain/sparse.py b/lib/gpt/core/domain/sparse.py
index fae4e919..34d8374b 100644
--- a/lib/gpt/core/domain/sparse.py
+++ b/lib/gpt/core/domain/sparse.py
@@ -37,9 +37,12 @@ def __init__(self, grid, local_coordinates, divisible_by):
         n = len(local_coordinates)
         N = self.grid.Nprocessors
         l = np.zeros(N, dtype=np.uint64)
-        l[self.grid.processor] = 2 ** int(np.ceil(np.log(n / divisible_by) / np.log(2)))
+        m = np.prod(divisible_by)
+        l[self.grid.processor] = 2 ** int(np.ceil(np.log(n / m) / np.log(2)))
         l = grid.globalsum(l)
-        self.L = [int(np.max(l) * divisible_by) * self.grid.mpi[0]] + self.grid.mpi[1:]
+        self.L = [int(np.max(l) * divisible_by[0]) * self.grid.mpi[0]] + [
+            self.grid.mpi[mu] * divisible_by[mu] for mu in range(1, self.grid.nd)
+        ]
 
         cb_simd_only_first_dimension = gpt.general(1, [0] * grid.nd, [1] + [0] * (grid.nd - 1))
 
@@ -130,10 +133,15 @@ def coordinate_lattices(self, mark_empty=0):
 
 
 class sparse:
-    def __init__(self, grid, local_coordinates, divisible_by=1):
+    def __init__(self, grid, local_coordinates, divisible_by=None):
         self.local_coordinates = local_coordinates
         self.grid = grid
 
+        if type(divisible_by) is int:
+            divisible_by = [divisible_by] + [1] * (grid.nd - 1)
+        if divisible_by is None:
+            divisible_by = [1] * grid.nd
+
         # kernel to avoid circular references through captures below
         kernel = sparse_kernel(grid, local_coordinates, divisible_by)
 

From 7138082945a1221ba48b6dbe149c1ae6f87d4b02 Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Tue, 29 Aug 2023 15:54:39 +0200
Subject: [PATCH 03/17] fixing merge with lehner master

---
 lib/gpt/core/domain/sparse.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/lib/gpt/core/domain/sparse.py b/lib/gpt/core/domain/sparse.py
index 44bc666f..a66a4719 100644
--- a/lib/gpt/core/domain/sparse.py
+++ b/lib/gpt/core/domain/sparse.py
@@ -145,11 +145,6 @@ class sparse:
     def __init__(self, grid, unmasked_local_coordinates, dimensions_divisible_by=None, mask=None):
         self.grid = grid
 
-        if type(divisible_by) is int:
-            divisible_by = [divisible_by] + [1] * (grid.nd - 1)
-        if divisible_by is None:
-            divisible_by = [1] * grid.nd
-
         # kernel to avoid circular references through captures below
         kernel = sparse_kernel(grid, unmasked_local_coordinates, dimensions_divisible_by, mask)
 

From b481e170b1b6200cceb9977b2c7f95152f589722 Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Fri, 1 Sep 2023 14:36:21 +0200
Subject: [PATCH 04/17] added support for CP^N-1 model

---
 applications/hmc/cpn.py               | 120 ++++++++++++++++++++++++++
 lib/gpt/qcd/scalar/action/__init__.py |   1 +
 lib/gpt/qcd/scalar/action/cpn.py      | 118 +++++++++++++++++++++++++
 3 files changed, 239 insertions(+)
 create mode 100755 applications/hmc/cpn.py
 create mode 100644 lib/gpt/qcd/scalar/action/cpn.py

diff --git a/applications/hmc/cpn.py b/applications/hmc/cpn.py
new file mode 100755
index 00000000..cbeab866
--- /dev/null
+++ b/applications/hmc/cpn.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python3
+#
+# Authors: Christoph Lehner, Mattia Bruno, Gabriele Morandi 2023
+#
+# HMC for 2D CP^N-1 theory
+#
+import gpt as g
+import sys, os
+import numpy
+
+grid = g.grid([42, 42], g.double)
+rng = g.random("hmc-cpn-model")
+
+# action for conj. momenta:
+g.message(f"Lattice = {grid.fdimensions}")
+g.message("Actions:")
+a0 = g.qcd.scalar.action.mass_term()
+g.message(f" - {a0.__name__}")
+
+# CP^N-1 action
+beta = 0.70
+N = 10
+a1 = g.qcd.scalar.action.cpn(N, beta)
+g.message(f" - {a1.__name__}")
+
+# fields
+z = g.vcomplex(grid, N)
+a1.draw(z, rng)
+l = [g.u1(grid) for _ in range(grid.nd)] 
+rng.element(l)
+fields = [z] + l
+
+# conjugate momenta
+mom_z = g.group.cartesian(z)
+mom_l = g.group.cartesian(l)
+moms = [mom_z] + mom_l
+
+def hamiltonian():
+    return a0(mom_z) + a0(mom_l) + a1(fields)
+
+
+# molecular dynamics
+sympl = g.algorithms.integrator.symplectic
+
+ip_z = sympl.update_p(mom_z, lambda: a1.gradient(fields, z))
+ip_l = sympl.update_p(mom_l, lambda: a1.gradient(fields, l))
+
+class constrained_iq(sympl.symplectic_base):
+    def __init__(self, fields, moms):
+        z = fields[0]
+        l = fields[1:]
+        mom_z = moms[0]
+        mom_l = moms[1:]
+
+        iq_l = sympl.update_q(l, lambda: a0.gradient(mom_l, mom_l))
+
+        def inner(eps):
+            a1.constrained_leap_frog(eps, z, mom_z)
+            iq_l(eps)
+        super().__init__(1, [], [inner], None, "constrained iq")
+        
+iq = constrained_iq(fields, moms)
+
+# integrator
+mdint = sympl.leap_frog(50, [ip_z, ip_l], iq)
+g.message(f"Integration scheme:\n{mdint}")
+
+# metropolis
+metro = g.algorithms.markov.metropolis(rng)
+
+# MD units
+tau = 1.0
+g.message(f"tau = {tau} MD units")
+
+
+def hmc(tau, moms):
+    mom_z = moms[0]
+    mom_l = moms[1:]
+    
+    rng.normal_element(mom_l)
+    a1.draw(mom_z, rng, z)
+    
+    accrej = metro(fields)
+    h0 = hamiltonian()
+    mdint(tau)
+    h1 = hamiltonian()
+    return [accrej(h1, h0), h1 - h0]
+
+
+# thermalization
+for i in range(1, 21):
+    h = []
+    for _ in range(10):
+        h += [hmc(tau, moms)]
+    h = numpy.array(h)
+    g.message(f"{i*5} % of thermalization completed")
+    g.message(
+        f"Action = {a1(fields)}, Acceptance = {numpy.mean(h[:,0]):.2f}, |dH| = {numpy.mean(numpy.abs(h[:,1])):.4e}"
+    )
+
+# measure action
+def measure():
+    return [a1(fields) / (N * beta * grid.fsites * grid.nd)]
+
+
+# production
+history = []
+data = []
+for i in range(100):
+    for k in range(10):
+        history += [hmc(tau, moms)]
+    data += [measure()]
+    g.message(f"Trajectory {i}")
+
+history = numpy.array(history)
+g.message(f"Acceptance rate = {numpy.mean(history[:,0]):.2f}")
+g.message(f"<|dH|> = {numpy.mean(numpy.abs(history[:,1])):.4e}")
+
+data = numpy.array(data)
+g.message(f"Energy density   = {numpy.mean(data[:,0])}")
diff --git a/lib/gpt/qcd/scalar/action/__init__.py b/lib/gpt/qcd/scalar/action/__init__.py
index 81709bcc..d203f1a7 100644
--- a/lib/gpt/qcd/scalar/action/__init__.py
+++ b/lib/gpt/qcd/scalar/action/__init__.py
@@ -18,4 +18,5 @@
 #    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 #
 from gpt.qcd.scalar.action.phi4 import phi4
+from gpt.qcd.scalar.action.cpn import cpn
 from gpt.qcd.scalar.action.mass_term import mass_term
diff --git a/lib/gpt/qcd/scalar/action/cpn.py b/lib/gpt/qcd/scalar/action/cpn.py
new file mode 100644
index 00000000..749c1221
--- /dev/null
+++ b/lib/gpt/qcd/scalar/action/cpn.py
@@ -0,0 +1,118 @@
+#
+#    GPT - Grid Python Toolkit
+#    Copyright (C) 2023  Christoph Lehner (christoph.lehner@ur.de, https://github.com/lehner/gpt)
+#                  2023  Mattia Bruno, Gabriele Morandi
+#
+#    This program is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation; either version 2 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with this program; if not, write to the Free Software Foundation, Inc.,
+#    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+import gpt as g
+from gpt.core.group import differentiable_functional
+import numpy
+
+# CP^{N-1} model
+
+# S[z,l] = -N * beta * sum_n,mu (z_{n+mu}^dag * z_n * l_{n,mu} + z_n^dag * z_{n+mu} * l_{n,mu}^dag)
+#        = -2 * N * beta * sum_n,mu [Re(z_{n+mu}^dag * z_n * l_{n,mu}) - 1]
+class cpn(differentiable_functional):
+    def __init__(self, N, b):
+        self.N = N
+        self.beta = b
+        self.__name__ = f"cpn(N = {self.N}, beta = {self.beta})"
+
+    # z = fields[0], l = fields[1:]
+    def split(self, fields):
+        return fields[0], fields[1:]
+
+    def __call__(self, fields):
+        z, l = self.split(fields)
+
+        J = g.lattice(z)
+        J[:] = 0.0
+        for mu in range(z.grid.nd):
+            J += g.cshift(z, mu, +1) * g.adj(l[mu])
+            
+        action = -2 * self.N * self.beta * (g.inner_product(J, z).real - z.grid.fsites * z.grid.nd)
+        return action
+
+
+    @differentiable_functional.multi_field_gradient
+    def gradient(self, fields, dfields):
+        def gradient_l(z, l, mu):
+            frc  = g.lattice(l[0])
+            frc @= 2 * self.beta * self.N * g.component.imag(g.trace(z * g.adj(g.cshift(z, mu, +1))) * l[mu])
+            frc.otype = l[0].otype.cartesian()
+            return frc
+         
+        def gradient_z(z, l):
+            J = g.lattice(z)
+            J[:] = 0.0
+            for mu in range(z.grid.nd):
+                J += g.cshift(z, mu, +1) * g.adj(l[mu]) 
+                J += g.cshift(z * l[mu], mu, -1)
+
+            frc = g.lattice(z)
+            frc @= - 2 * self.beta * self.N * J
+
+            frc -= g.trace(frc * g.adj(z)) * z
+            return frc
+        
+        z, l = self.split(fields)
+        frc = []
+        for df in g.core.util.to_list(dfields):
+            k = fields.index(df)
+            if k == 0:
+                frc.append(gradient_z(z, l))
+            else:
+                frc.append(gradient_l(z,l,k-1))
+        return frc
+
+    # https://arxiv.org/abs/1102.1852
+    def constrained_leap_frog(self, eps, z, mom_z):
+        # TO DO: replace with g.adj(v1) * v2
+        def dot(v1, v2):
+            return g.trace(v2 * g.adj(v1))
+
+        n  = g.real(z.grid)
+        n @= g.component.sqrt(g.component.real(dot(mom_z, mom_z))) 
+
+        # z'      =  cos(alpha) z + (1/|pi|) sin(alpha) mom_z
+        # mom_z'  = -|pi| sin(alpha) z + cos(alpha) mom_z
+        # alpha = eps |pi|
+        _z = g.lattice(z)
+        _z @= z
+        
+        cos  = g.real(z.grid) 
+        cos @= g.component.cos(eps * n)
+
+        sin  = g.real(z.grid)
+        sin @= g.component.sin(eps * n)
+        
+        z @= cos * _z + g(g.component.inv(n) * sin) * mom_z
+        mom_z @= - g(n * sin) * _z + cos * mom_z
+        del _z, cos, sin, n
+        
+    # https://arxiv.org/abs/1102.1852
+    def draw(self, field, rng, constraint=None):
+        if constraint is None:
+            z = field
+            rng.element(z)
+            n = g.component.real(g.trace(z * g.adj(z)))
+            z @= z * g.component.inv(g.component.sqrt(n))
+        else:
+            mom_z = field
+            z = constraint
+            rng.normal_element(mom_z)
+            # TO DO change to z * g(g.adj(z) * mom_z)
+            mom_z @= mom_z - g(z * g.adj(z)) * mom_z
\ No newline at end of file

From eacc2ff4daf3322ceba5b59d7859d8294fc028b9 Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Fri, 1 Sep 2023 14:43:18 +0200
Subject: [PATCH 05/17] introduced step_size verbosity for symplectic
 integrators

---
 applications/hmc/cpn.py                     | 1 +
 lib/gpt/algorithms/integrator/symplectic.py | 5 ++++-
 lib/gpt/default.py                          | 2 +-
 3 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/applications/hmc/cpn.py b/applications/hmc/cpn.py
index cbeab866..3ec72254 100755
--- a/applications/hmc/cpn.py
+++ b/applications/hmc/cpn.py
@@ -8,6 +8,7 @@
 import sys, os
 import numpy
 
+g.default.set_verbose("step_size", False)
 grid = g.grid([42, 42], g.double)
 rng = g.random("hmc-cpn-model")
 
diff --git a/lib/gpt/algorithms/integrator/symplectic.py b/lib/gpt/algorithms/integrator/symplectic.py
index 6516eeb3..bcbb6c98 100644
--- a/lib/gpt/algorithms/integrator/symplectic.py
+++ b/lib/gpt/algorithms/integrator/symplectic.py
@@ -81,8 +81,11 @@ def __mul__(self, f):
         return step(self.funcs, [c * f for c in self.c], self.n)
 
     def __call__(self, eps):
+        verbose = gpt.default.is_verbose("step_size")
+        
         for i in range(self.nf):
-            gpt.message(f"call eps = {eps}, {i} / {self.nf}")
+            if verbose:
+                gpt.message(f"call eps = {eps}, {i} / {self.nf}")
             self.funcs[i](self.c[i] * eps**self.n)
 
 
diff --git a/lib/gpt/default.py b/lib/gpt/default.py
index 80427ec4..00717c8e 100644
--- a/lib/gpt/default.py
+++ b/lib/gpt/default.py
@@ -76,7 +76,7 @@ def get_ivec(tag, default, ndim):
     "io,bicgstab,cg,defect_correcting,cagcr,fgcr,fgmres,mr,irl,repository,arnoldi,power_iteration,"
     + "checkpointer,modes,random,split,coarse_grid,gradient_descent,adam,non_linear_cg,"
     + "coarsen,qis_map,metropolis,su2_heat_bath,u1_heat_bath,fom,chronological,minimal_residual_extrapolation,"
-    + "subspace_minimal_residual"
+    + "subspace_minimal_residual,step_size"
 )
 verbose_additional = "eval,merge,orthogonalize,copy_plan"
 verbose = set()

From ea2e23913175be2f54f4bbc9d99f8c661b59f6fa Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Sun, 3 Sep 2023 15:49:37 +0200
Subject: [PATCH 06/17] towards inner product iVSinglet4-10-30 [missing expr
 template]

---
 lib/cgpt/lib/eval/mul_vlat_vlat.h | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/lib/cgpt/lib/eval/mul_vlat_vlat.h b/lib/cgpt/lib/eval/mul_vlat_vlat.h
index 022b88f9..62ff122a 100644
--- a/lib/cgpt/lib/eval/mul_vlat_vlat.h
+++ b/lib/cgpt/lib/eval/mul_vlat_vlat.h
@@ -59,9 +59,9 @@ void eval_mul_vlat_vlat(std::vector<cgpt_Lattice_base*> & dst_vl,
 
   // VV -> S/M
   if (lhs_singlet_rank == 1 && rhs_singlet_rank == 1) {
+    ASSERT(lhs_singlet_dim == rhs_singlet_dim);
     if (lhs_unary == 0 && rhs_unary == (BIT_TRANS|BIT_CONJ)) {
       // outer product -> M
-      ASSERT(lhs_singlet_dim == rhs_singlet_dim);
       dst_vl.resize(lhs_singlet_dim*rhs_singlet_dim, 0);
       for (int i=0;i<lhs_singlet_dim;i++) {
 	for (int j=0;j<rhs_singlet_dim;j++) {
@@ -71,9 +71,8 @@ void eval_mul_vlat_vlat(std::vector<cgpt_Lattice_base*> & dst_vl,
       }
       return;
     } else if (lhs_unary == (BIT_TRANS|BIT_CONJ) && rhs_unary == 0) {
-      ERR("Not implemented");
+      //ERR("Not implemented");
       // inner product -> S
-      /*ASSERT(lhs_singlet_dim == rhs_singlet_dim);
       dst_vl.resize(1, 0);
       bool _ac = ac;
       for (int i=0;i<lhs_singlet_dim;i++) {
@@ -82,7 +81,7 @@ void eval_mul_vlat_vlat(std::vector<cgpt_Lattice_base*> & dst_vl,
 	  _ac = true;
 	}
       }
-      return;*/
+      return;
     } else {
       ERR("Invalid combination of two vectors");
     }

From bbfd183f19f220d228e91d441df705372ba6e366 Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Sun, 3 Sep 2023 15:55:45 +0200
Subject: [PATCH 07/17] style

---
 applications/hmc/cpn.py                     | 12 ++++---
 lib/gpt/algorithms/integrator/symplectic.py |  2 +-
 lib/gpt/qcd/scalar/action/cpn.py            | 40 +++++++++++----------
 3 files changed, 31 insertions(+), 23 deletions(-)

diff --git a/applications/hmc/cpn.py b/applications/hmc/cpn.py
index 3ec72254..59e6599f 100755
--- a/applications/hmc/cpn.py
+++ b/applications/hmc/cpn.py
@@ -27,7 +27,7 @@
 # fields
 z = g.vcomplex(grid, N)
 a1.draw(z, rng)
-l = [g.u1(grid) for _ in range(grid.nd)] 
+l = [g.u1(grid) for _ in range(grid.nd)]
 rng.element(l)
 fields = [z] + l
 
@@ -36,6 +36,7 @@
 mom_l = g.group.cartesian(l)
 moms = [mom_z] + mom_l
 
+
 def hamiltonian():
     return a0(mom_z) + a0(mom_l) + a1(fields)
 
@@ -46,6 +47,7 @@ def hamiltonian():
 ip_z = sympl.update_p(mom_z, lambda: a1.gradient(fields, z))
 ip_l = sympl.update_p(mom_l, lambda: a1.gradient(fields, l))
 
+
 class constrained_iq(sympl.symplectic_base):
     def __init__(self, fields, moms):
         z = fields[0]
@@ -58,8 +60,10 @@ def __init__(self, fields, moms):
         def inner(eps):
             a1.constrained_leap_frog(eps, z, mom_z)
             iq_l(eps)
+
         super().__init__(1, [], [inner], None, "constrained iq")
-        
+
+
 iq = constrained_iq(fields, moms)
 
 # integrator
@@ -77,10 +81,10 @@ def inner(eps):
 def hmc(tau, moms):
     mom_z = moms[0]
     mom_l = moms[1:]
-    
+
     rng.normal_element(mom_l)
     a1.draw(mom_z, rng, z)
-    
+
     accrej = metro(fields)
     h0 = hamiltonian()
     mdint(tau)
diff --git a/lib/gpt/algorithms/integrator/symplectic.py b/lib/gpt/algorithms/integrator/symplectic.py
index bcbb6c98..b2e262f5 100644
--- a/lib/gpt/algorithms/integrator/symplectic.py
+++ b/lib/gpt/algorithms/integrator/symplectic.py
@@ -82,7 +82,7 @@ def __mul__(self, f):
 
     def __call__(self, eps):
         verbose = gpt.default.is_verbose("step_size")
-        
+
         for i in range(self.nf):
             if verbose:
                 gpt.message(f"call eps = {eps}, {i} / {self.nf}")
diff --git a/lib/gpt/qcd/scalar/action/cpn.py b/lib/gpt/qcd/scalar/action/cpn.py
index 749c1221..3dd3f01e 100644
--- a/lib/gpt/qcd/scalar/action/cpn.py
+++ b/lib/gpt/qcd/scalar/action/cpn.py
@@ -42,32 +42,36 @@ def __call__(self, fields):
         J[:] = 0.0
         for mu in range(z.grid.nd):
             J += g.cshift(z, mu, +1) * g.adj(l[mu])
-            
+
         action = -2 * self.N * self.beta * (g.inner_product(J, z).real - z.grid.fsites * z.grid.nd)
         return action
 
-
     @differentiable_functional.multi_field_gradient
     def gradient(self, fields, dfields):
         def gradient_l(z, l, mu):
-            frc  = g.lattice(l[0])
-            frc @= 2 * self.beta * self.N * g.component.imag(g.trace(z * g.adj(g.cshift(z, mu, +1))) * l[mu])
+            frc = g.lattice(l[0])
+            frc @= (
+                2
+                * self.beta
+                * self.N
+                * g.component.imag(g.trace(z * g.adj(g.cshift(z, mu, +1))) * l[mu])
+            )
             frc.otype = l[0].otype.cartesian()
             return frc
-         
+
         def gradient_z(z, l):
             J = g.lattice(z)
             J[:] = 0.0
             for mu in range(z.grid.nd):
-                J += g.cshift(z, mu, +1) * g.adj(l[mu]) 
+                J += g.cshift(z, mu, +1) * g.adj(l[mu])
                 J += g.cshift(z * l[mu], mu, -1)
 
             frc = g.lattice(z)
-            frc @= - 2 * self.beta * self.N * J
+            frc @= -2 * self.beta * self.N * J
 
             frc -= g.trace(frc * g.adj(z)) * z
             return frc
-        
+
         z, l = self.split(fields)
         frc = []
         for df in g.core.util.to_list(dfields):
@@ -75,7 +79,7 @@ def gradient_z(z, l):
             if k == 0:
                 frc.append(gradient_z(z, l))
             else:
-                frc.append(gradient_l(z,l,k-1))
+                frc.append(gradient_l(z, l, k - 1))
         return frc
 
     # https://arxiv.org/abs/1102.1852
@@ -84,25 +88,25 @@ def constrained_leap_frog(self, eps, z, mom_z):
         def dot(v1, v2):
             return g.trace(v2 * g.adj(v1))
 
-        n  = g.real(z.grid)
-        n @= g.component.sqrt(g.component.real(dot(mom_z, mom_z))) 
+        n = g.real(z.grid)
+        n @= g.component.sqrt(g.component.real(dot(mom_z, mom_z)))
 
         # z'      =  cos(alpha) z + (1/|pi|) sin(alpha) mom_z
         # mom_z'  = -|pi| sin(alpha) z + cos(alpha) mom_z
         # alpha = eps |pi|
         _z = g.lattice(z)
         _z @= z
-        
-        cos  = g.real(z.grid) 
+
+        cos = g.real(z.grid)
         cos @= g.component.cos(eps * n)
 
-        sin  = g.real(z.grid)
+        sin = g.real(z.grid)
         sin @= g.component.sin(eps * n)
-        
+
         z @= cos * _z + g(g.component.inv(n) * sin) * mom_z
-        mom_z @= - g(n * sin) * _z + cos * mom_z
+        mom_z @= -g(n * sin) * _z + cos * mom_z
         del _z, cos, sin, n
-        
+
     # https://arxiv.org/abs/1102.1852
     def draw(self, field, rng, constraint=None):
         if constraint is None:
@@ -115,4 +119,4 @@ def draw(self, field, rng, constraint=None):
             z = constraint
             rng.normal_element(mom_z)
             # TO DO change to z * g(g.adj(z) * mom_z)
-            mom_z @= mom_z - g(z * g.adj(z)) * mom_z
\ No newline at end of file
+            mom_z @= mom_z - g(z * g.adj(z)) * mom_z

From c947436538a8a705e947908b1ec2cf0a86333c62 Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Fri, 20 Oct 2023 17:57:00 +0200
Subject: [PATCH 08/17] added support site inner product vcomplex

---
 lib/cgpt/lib/expression/unary.h               |  2 +-
 .../expression_mul_double_iVSinglet10.cc      |  2 +-
 .../expression_mul_double_iVSinglet30.cc      |  2 +-
 .../expression_mul_double_iVSinglet4.cc       |  2 +-
 .../expression_mul_single_iVSinglet10.cc      |  2 +-
 .../expression_mul_single_iVSinglet30.cc      |  2 +-
 .../expression_mul_single_iVSinglet4.cc       |  2 +-
 .../iVSinglet/expression_mul.template         |  2 +-
 tests/core/vcomplex.py                        | 21 +++++++++++++++++++
 9 files changed, 29 insertions(+), 8 deletions(-)
 create mode 100644 tests/core/vcomplex.py

diff --git a/lib/cgpt/lib/expression/unary.h b/lib/cgpt/lib/expression/unary.h
index a06d9e60..c80bc7cf 100644
--- a/lib/cgpt/lib/expression/unary.h
+++ b/lib/cgpt/lib/expression/unary.h
@@ -78,7 +78,7 @@ cgpt_Lattice_base* lattice_lat(cgpt_Lattice_base* dst, bool ac, const A& lat, Co
 template<typename A>
 cgpt_Lattice_base* lattice_unary_lat(cgpt_Lattice_base* dst, bool ac, const A& la,int unary_expr,ComplexD coef) {
   if (unary_expr == 0) {
-    return lattice_lat(dst, ac, la, coef);
+    return lattice_lat(dst, ac, closure(ToSinglet(la)), coef);
   } else if (unary_expr == (BIT_SPINTRACE|BIT_COLORTRACE)) {
     return lattice_expr(dst, ac, coef*ToSinglet(trace(la)));
   } else if (unary_expr == BIT_SPINTRACE) {
diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc
index 5b98f13e..d139d1e7 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexD vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet10);
-  //  _INNER_PRODUCT_(iVSinglet10);
+  _INNER_PRODUCT_(iVSinglet10);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc
index 90bc07f7..492c1d37 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexD vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet30);
-  //  _INNER_PRODUCT_(iVSinglet30);
+  _INNER_PRODUCT_(iVSinglet30);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc
index 760f1028..7f632587 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexD vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet4);
-  //  _INNER_PRODUCT_(iVSinglet4);
+  _INNER_PRODUCT_(iVSinglet4);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc
index ff1b993d..f2be2da1 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexF vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet10);
-  //  _INNER_PRODUCT_(iVSinglet10);
+  _INNER_PRODUCT_(iVSinglet10);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc
index 8dd56db1..9391e0c7 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexF vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet30);
-  //  _INNER_PRODUCT_(iVSinglet30);
+  _INNER_PRODUCT_(iVSinglet30);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc
index 7abd5db9..904e7fae 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexF vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet4);
-  //  _INNER_PRODUCT_(iVSinglet4);
+  _INNER_PRODUCT_(iVSinglet4);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template b/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template
index 092ee1b8..d7eb6f27 100644
--- a/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template
+++ b/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template
@@ -40,6 +40,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef {precision_vector} vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet{tensor_arg_1});
-  //  _INNER_PRODUCT_(iVSinglet{tensor_arg_1});
+  _INNER_PRODUCT_(iVSinglet{tensor_arg_1});
   ERR("Not implemented");
 }
diff --git a/tests/core/vcomplex.py b/tests/core/vcomplex.py
new file mode 100644
index 00000000..d1b5097f
--- /dev/null
+++ b/tests/core/vcomplex.py
@@ -0,0 +1,21 @@
+#!/usr/bin/env python3
+#
+# Authors: Christoph Lehner 2020, Mattia Bruno 2023
+#
+# Desc.: Illustrate core concepts and features
+#
+import gpt as g
+
+L = [8, 4, 4, 4]
+grid = g.grid(L, g.double)
+rng = g.random("vcomplex")
+
+def test_local_inner_product(v):
+    i0 = g(g.adj(v) * v)
+    i1 = g(g.trace(v * g.adj(v)))
+    assert g.norm2(i0 - i1) < 1e-15
+
+for N in [4, 10, 30]:
+    v1 = g.vcomplex(grid, N)
+    test_local_inner_product(v1)
+    del v1

From 67b6ed12a4aa941c9b680e5e1815fd80df9b286d Mon Sep 17 00:00:00 2001
From: Christoph Lehner <christoph@lhnr.de>
Date: Wed, 31 May 2023 08:12:30 +0300
Subject: [PATCH 09/17] timings

---
 lib/gpt/qcd/gauge/stencil/staple.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/gpt/qcd/gauge/stencil/staple.py b/lib/gpt/qcd/gauge/stencil/staple.py
index 2397efe2..3d14c228 100644
--- a/lib/gpt/qcd/gauge/stencil/staple.py
+++ b/lib/gpt/qcd/gauge/stencil/staple.py
@@ -35,8 +35,8 @@ def staple_sum(U, rho, mu=None, cache=default_staple_cache):
 
     Ntarget = len(target_mu)
 
+    t = g.timer("staples")
     if tag not in cache:
-
         code = []
         nwr = [0] * Ntarget
         for idx in range(Ntarget):

From 6920f168b4c77bd08d1de61c9e946db6248f5100 Mon Sep 17 00:00:00 2001
From: Christoph Lehner <christoph@lhnr.de>
Date: Wed, 31 May 2023 16:03:50 +0200
Subject: [PATCH 10/17] new parallel_transport_matrix

---
 lib/gpt/qcd/gauge/stencil/staple.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/lib/gpt/qcd/gauge/stencil/staple.py b/lib/gpt/qcd/gauge/stencil/staple.py
index 3d14c228..1269da6b 100644
--- a/lib/gpt/qcd/gauge/stencil/staple.py
+++ b/lib/gpt/qcd/gauge/stencil/staple.py
@@ -35,7 +35,6 @@ def staple_sum(U, rho, mu=None, cache=default_staple_cache):
 
     Ntarget = len(target_mu)
 
-    t = g.timer("staples")
     if tag not in cache:
         code = []
         nwr = [0] * Ntarget

From 58397426075e3b9ef615c26e27cb192ac766408e Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Fri, 1 Sep 2023 14:36:21 +0200
Subject: [PATCH 11/17] added support for CP^N-1 model

---
 applications/hmc/cpn.py               | 120 ++++++++++++++++++++++++++
 lib/gpt/qcd/scalar/action/__init__.py |   1 +
 lib/gpt/qcd/scalar/action/cpn.py      | 118 +++++++++++++++++++++++++
 3 files changed, 239 insertions(+)
 create mode 100755 applications/hmc/cpn.py
 create mode 100644 lib/gpt/qcd/scalar/action/cpn.py

diff --git a/applications/hmc/cpn.py b/applications/hmc/cpn.py
new file mode 100755
index 00000000..cbeab866
--- /dev/null
+++ b/applications/hmc/cpn.py
@@ -0,0 +1,120 @@
+#!/usr/bin/env python3
+#
+# Authors: Christoph Lehner, Mattia Bruno, Gabriele Morandi 2023
+#
+# HMC for 2D CP^N-1 theory
+#
+import gpt as g
+import sys, os
+import numpy
+
+grid = g.grid([42, 42], g.double)
+rng = g.random("hmc-cpn-model")
+
+# action for conj. momenta:
+g.message(f"Lattice = {grid.fdimensions}")
+g.message("Actions:")
+a0 = g.qcd.scalar.action.mass_term()
+g.message(f" - {a0.__name__}")
+
+# CP^N-1 action
+beta = 0.70
+N = 10
+a1 = g.qcd.scalar.action.cpn(N, beta)
+g.message(f" - {a1.__name__}")
+
+# fields
+z = g.vcomplex(grid, N)
+a1.draw(z, rng)
+l = [g.u1(grid) for _ in range(grid.nd)] 
+rng.element(l)
+fields = [z] + l
+
+# conjugate momenta
+mom_z = g.group.cartesian(z)
+mom_l = g.group.cartesian(l)
+moms = [mom_z] + mom_l
+
+def hamiltonian():
+    return a0(mom_z) + a0(mom_l) + a1(fields)
+
+
+# molecular dynamics
+sympl = g.algorithms.integrator.symplectic
+
+ip_z = sympl.update_p(mom_z, lambda: a1.gradient(fields, z))
+ip_l = sympl.update_p(mom_l, lambda: a1.gradient(fields, l))
+
+class constrained_iq(sympl.symplectic_base):
+    def __init__(self, fields, moms):
+        z = fields[0]
+        l = fields[1:]
+        mom_z = moms[0]
+        mom_l = moms[1:]
+
+        iq_l = sympl.update_q(l, lambda: a0.gradient(mom_l, mom_l))
+
+        def inner(eps):
+            a1.constrained_leap_frog(eps, z, mom_z)
+            iq_l(eps)
+        super().__init__(1, [], [inner], None, "constrained iq")
+        
+iq = constrained_iq(fields, moms)
+
+# integrator
+mdint = sympl.leap_frog(50, [ip_z, ip_l], iq)
+g.message(f"Integration scheme:\n{mdint}")
+
+# metropolis
+metro = g.algorithms.markov.metropolis(rng)
+
+# MD units
+tau = 1.0
+g.message(f"tau = {tau} MD units")
+
+
+def hmc(tau, moms):
+    mom_z = moms[0]
+    mom_l = moms[1:]
+    
+    rng.normal_element(mom_l)
+    a1.draw(mom_z, rng, z)
+    
+    accrej = metro(fields)
+    h0 = hamiltonian()
+    mdint(tau)
+    h1 = hamiltonian()
+    return [accrej(h1, h0), h1 - h0]
+
+
+# thermalization
+for i in range(1, 21):
+    h = []
+    for _ in range(10):
+        h += [hmc(tau, moms)]
+    h = numpy.array(h)
+    g.message(f"{i*5} % of thermalization completed")
+    g.message(
+        f"Action = {a1(fields)}, Acceptance = {numpy.mean(h[:,0]):.2f}, |dH| = {numpy.mean(numpy.abs(h[:,1])):.4e}"
+    )
+
+# measure action
+def measure():
+    return [a1(fields) / (N * beta * grid.fsites * grid.nd)]
+
+
+# production
+history = []
+data = []
+for i in range(100):
+    for k in range(10):
+        history += [hmc(tau, moms)]
+    data += [measure()]
+    g.message(f"Trajectory {i}")
+
+history = numpy.array(history)
+g.message(f"Acceptance rate = {numpy.mean(history[:,0]):.2f}")
+g.message(f"<|dH|> = {numpy.mean(numpy.abs(history[:,1])):.4e}")
+
+data = numpy.array(data)
+g.message(f"Energy density   = {numpy.mean(data[:,0])}")
diff --git a/lib/gpt/qcd/scalar/action/__init__.py b/lib/gpt/qcd/scalar/action/__init__.py
index 81709bcc..d203f1a7 100644
--- a/lib/gpt/qcd/scalar/action/__init__.py
+++ b/lib/gpt/qcd/scalar/action/__init__.py
@@ -18,4 +18,5 @@
 #    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 #
 from gpt.qcd.scalar.action.phi4 import phi4
+from gpt.qcd.scalar.action.cpn import cpn
 from gpt.qcd.scalar.action.mass_term import mass_term
diff --git a/lib/gpt/qcd/scalar/action/cpn.py b/lib/gpt/qcd/scalar/action/cpn.py
new file mode 100644
index 00000000..749c1221
--- /dev/null
+++ b/lib/gpt/qcd/scalar/action/cpn.py
@@ -0,0 +1,118 @@
+#
+#    GPT - Grid Python Toolkit
+#    Copyright (C) 2023  Christoph Lehner (christoph.lehner@ur.de, https://github.com/lehner/gpt)
+#                  2023  Mattia Bruno, Gabriele Morandi
+#
+#    This program is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation; either version 2 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with this program; if not, write to the Free Software Foundation, Inc.,
+#    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+import gpt as g
+from gpt.core.group import differentiable_functional
+import numpy
+
+# CP^{N-1} model
+
+# S[z,l] = -N * beta * sum_n,mu (z_{n+mu}^dag * z_n * l_{n,mu} + z_n^dag * z_{n+mu} * l_{n,mu}^dag)
+#        = -2 * N * beta * sum_n,mu [Re(z_{n+mu}^dag * z_n * l_{n,mu}) - 1]
+class cpn(differentiable_functional):
+    def __init__(self, N, b):
+        self.N = N
+        self.beta = b
+        self.__name__ = f"cpn(N = {self.N}, beta = {self.beta})"
+
+    # z = fields[0], l = fields[1:]
+    def split(self, fields):
+        return fields[0], fields[1:]
+
+    def __call__(self, fields):
+        z, l = self.split(fields)
+
+        J = g.lattice(z)
+        J[:] = 0.0
+        for mu in range(z.grid.nd):
+            J += g.cshift(z, mu, +1) * g.adj(l[mu])
+            
+        action = -2 * self.N * self.beta * (g.inner_product(J, z).real - z.grid.fsites * z.grid.nd)
+        return action
+
+
+    @differentiable_functional.multi_field_gradient
+    def gradient(self, fields, dfields):
+        def gradient_l(z, l, mu):
+            frc  = g.lattice(l[0])
+            frc @= 2 * self.beta * self.N * g.component.imag(g.trace(z * g.adj(g.cshift(z, mu, +1))) * l[mu])
+            frc.otype = l[0].otype.cartesian()
+            return frc
+         
+        def gradient_z(z, l):
+            J = g.lattice(z)
+            J[:] = 0.0
+            for mu in range(z.grid.nd):
+                J += g.cshift(z, mu, +1) * g.adj(l[mu]) 
+                J += g.cshift(z * l[mu], mu, -1)
+
+            frc = g.lattice(z)
+            frc @= - 2 * self.beta * self.N * J
+
+            frc -= g.trace(frc * g.adj(z)) * z
+            return frc
+        
+        z, l = self.split(fields)
+        frc = []
+        for df in g.core.util.to_list(dfields):
+            k = fields.index(df)
+            if k == 0:
+                frc.append(gradient_z(z, l))
+            else:
+                frc.append(gradient_l(z,l,k-1))
+        return frc
+
+    # https://arxiv.org/abs/1102.1852
+    def constrained_leap_frog(self, eps, z, mom_z):
+        # TO DO: replace with g.adj(v1) * v2
+        def dot(v1, v2):
+            return g.trace(v2 * g.adj(v1))
+
+        n  = g.real(z.grid)
+        n @= g.component.sqrt(g.component.real(dot(mom_z, mom_z))) 
+
+        # z'      =  cos(alpha) z + (1/|pi|) sin(alpha) mom_z
+        # mom_z'  = -|pi| sin(alpha) z + cos(alpha) mom_z
+        # alpha = eps |pi|
+        _z = g.lattice(z)
+        _z @= z
+        
+        cos  = g.real(z.grid) 
+        cos @= g.component.cos(eps * n)
+
+        sin  = g.real(z.grid)
+        sin @= g.component.sin(eps * n)
+        
+        z @= cos * _z + g(g.component.inv(n) * sin) * mom_z
+        mom_z @= - g(n * sin) * _z + cos * mom_z
+        del _z, cos, sin, n
+        
+    # https://arxiv.org/abs/1102.1852
+    def draw(self, field, rng, constraint=None):
+        if constraint is None:
+            z = field
+            rng.element(z)
+            n = g.component.real(g.trace(z * g.adj(z)))
+            z @= z * g.component.inv(g.component.sqrt(n))
+        else:
+            mom_z = field
+            z = constraint
+            rng.normal_element(mom_z)
+            # TO DO change to z * g(g.adj(z) * mom_z)
+            mom_z @= mom_z - g(z * g.adj(z)) * mom_z
\ No newline at end of file

From fd3aed5ee0a67c02a1d9c32905e7e2f5760771f7 Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Fri, 1 Sep 2023 14:43:18 +0200
Subject: [PATCH 12/17] introduced step_size verbosity for symplectic
 integrators

---
 applications/hmc/cpn.py                     | 1 +
 lib/gpt/algorithms/integrator/symplectic.py | 5 ++++-
 lib/gpt/default.py                          | 2 +-
 3 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/applications/hmc/cpn.py b/applications/hmc/cpn.py
index cbeab866..3ec72254 100755
--- a/applications/hmc/cpn.py
+++ b/applications/hmc/cpn.py
@@ -8,6 +8,7 @@
 import sys, os
 import numpy
 
+g.default.set_verbose("step_size", False)
 grid = g.grid([42, 42], g.double)
 rng = g.random("hmc-cpn-model")
 
diff --git a/lib/gpt/algorithms/integrator/symplectic.py b/lib/gpt/algorithms/integrator/symplectic.py
index 3241219b..bcbb6c98 100644
--- a/lib/gpt/algorithms/integrator/symplectic.py
+++ b/lib/gpt/algorithms/integrator/symplectic.py
@@ -81,8 +81,11 @@ def __mul__(self, f):
         return step(self.funcs, [c * f for c in self.c], self.n)
 
     def __call__(self, eps):
+        verbose = gpt.default.is_verbose("step_size")
+        
         for i in range(self.nf):
-            # gpt.message(f"call eps = {eps}, {i} / {self.nf}")
+            if verbose:
+                gpt.message(f"call eps = {eps}, {i} / {self.nf}")
             self.funcs[i](self.c[i] * eps**self.n)
 
 
diff --git a/lib/gpt/default.py b/lib/gpt/default.py
index 80427ec4..00717c8e 100644
--- a/lib/gpt/default.py
+++ b/lib/gpt/default.py
@@ -76,7 +76,7 @@ def get_ivec(tag, default, ndim):
     "io,bicgstab,cg,defect_correcting,cagcr,fgcr,fgmres,mr,irl,repository,arnoldi,power_iteration,"
     + "checkpointer,modes,random,split,coarse_grid,gradient_descent,adam,non_linear_cg,"
     + "coarsen,qis_map,metropolis,su2_heat_bath,u1_heat_bath,fom,chronological,minimal_residual_extrapolation,"
-    + "subspace_minimal_residual"
+    + "subspace_minimal_residual,step_size"
 )
 verbose_additional = "eval,merge,orthogonalize,copy_plan"
 verbose = set()

From eac6d4581b35980797065ab9b8ab937d7b3de6a1 Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Sun, 3 Sep 2023 15:49:37 +0200
Subject: [PATCH 13/17] towards inner product iVSinglet4-10-30 [missing expr
 template]

---
 lib/cgpt/lib/eval/mul_vlat_vlat.h | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/lib/cgpt/lib/eval/mul_vlat_vlat.h b/lib/cgpt/lib/eval/mul_vlat_vlat.h
index 022b88f9..62ff122a 100644
--- a/lib/cgpt/lib/eval/mul_vlat_vlat.h
+++ b/lib/cgpt/lib/eval/mul_vlat_vlat.h
@@ -59,9 +59,9 @@ void eval_mul_vlat_vlat(std::vector<cgpt_Lattice_base*> & dst_vl,
 
   // VV -> S/M
   if (lhs_singlet_rank == 1 && rhs_singlet_rank == 1) {
+    ASSERT(lhs_singlet_dim == rhs_singlet_dim);
     if (lhs_unary == 0 && rhs_unary == (BIT_TRANS|BIT_CONJ)) {
       // outer product -> M
-      ASSERT(lhs_singlet_dim == rhs_singlet_dim);
       dst_vl.resize(lhs_singlet_dim*rhs_singlet_dim, 0);
       for (int i=0;i<lhs_singlet_dim;i++) {
 	for (int j=0;j<rhs_singlet_dim;j++) {
@@ -71,9 +71,8 @@ void eval_mul_vlat_vlat(std::vector<cgpt_Lattice_base*> & dst_vl,
       }
       return;
     } else if (lhs_unary == (BIT_TRANS|BIT_CONJ) && rhs_unary == 0) {
-      ERR("Not implemented");
+      //ERR("Not implemented");
       // inner product -> S
-      /*ASSERT(lhs_singlet_dim == rhs_singlet_dim);
       dst_vl.resize(1, 0);
       bool _ac = ac;
       for (int i=0;i<lhs_singlet_dim;i++) {
@@ -82,7 +81,7 @@ void eval_mul_vlat_vlat(std::vector<cgpt_Lattice_base*> & dst_vl,
 	  _ac = true;
 	}
       }
-      return;*/
+      return;
     } else {
       ERR("Invalid combination of two vectors");
     }

From e000a8a71fa5d2496e903fb1a0d8b2ce0d6449df Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Sun, 3 Sep 2023 15:55:45 +0200
Subject: [PATCH 14/17] style

---
 applications/hmc/cpn.py                     | 12 ++++---
 lib/gpt/algorithms/integrator/symplectic.py |  2 +-
 lib/gpt/qcd/scalar/action/cpn.py            | 40 +++++++++++----------
 3 files changed, 31 insertions(+), 23 deletions(-)

diff --git a/applications/hmc/cpn.py b/applications/hmc/cpn.py
index 3ec72254..59e6599f 100755
--- a/applications/hmc/cpn.py
+++ b/applications/hmc/cpn.py
@@ -27,7 +27,7 @@
 # fields
 z = g.vcomplex(grid, N)
 a1.draw(z, rng)
-l = [g.u1(grid) for _ in range(grid.nd)] 
+l = [g.u1(grid) for _ in range(grid.nd)]
 rng.element(l)
 fields = [z] + l
 
@@ -36,6 +36,7 @@
 mom_l = g.group.cartesian(l)
 moms = [mom_z] + mom_l
 
+
 def hamiltonian():
     return a0(mom_z) + a0(mom_l) + a1(fields)
 
@@ -46,6 +47,7 @@ def hamiltonian():
 ip_z = sympl.update_p(mom_z, lambda: a1.gradient(fields, z))
 ip_l = sympl.update_p(mom_l, lambda: a1.gradient(fields, l))
 
+
 class constrained_iq(sympl.symplectic_base):
     def __init__(self, fields, moms):
         z = fields[0]
@@ -58,8 +60,10 @@ def __init__(self, fields, moms):
         def inner(eps):
             a1.constrained_leap_frog(eps, z, mom_z)
             iq_l(eps)
+
         super().__init__(1, [], [inner], None, "constrained iq")
-        
+
+
 iq = constrained_iq(fields, moms)
 
 # integrator
@@ -77,10 +81,10 @@ def inner(eps):
 def hmc(tau, moms):
     mom_z = moms[0]
     mom_l = moms[1:]
-    
+
     rng.normal_element(mom_l)
     a1.draw(mom_z, rng, z)
-    
+
     accrej = metro(fields)
     h0 = hamiltonian()
     mdint(tau)
diff --git a/lib/gpt/algorithms/integrator/symplectic.py b/lib/gpt/algorithms/integrator/symplectic.py
index bcbb6c98..b2e262f5 100644
--- a/lib/gpt/algorithms/integrator/symplectic.py
+++ b/lib/gpt/algorithms/integrator/symplectic.py
@@ -82,7 +82,7 @@ def __mul__(self, f):
 
     def __call__(self, eps):
         verbose = gpt.default.is_verbose("step_size")
-        
+
         for i in range(self.nf):
             if verbose:
                 gpt.message(f"call eps = {eps}, {i} / {self.nf}")
diff --git a/lib/gpt/qcd/scalar/action/cpn.py b/lib/gpt/qcd/scalar/action/cpn.py
index 749c1221..3dd3f01e 100644
--- a/lib/gpt/qcd/scalar/action/cpn.py
+++ b/lib/gpt/qcd/scalar/action/cpn.py
@@ -42,32 +42,36 @@ def __call__(self, fields):
         J[:] = 0.0
         for mu in range(z.grid.nd):
             J += g.cshift(z, mu, +1) * g.adj(l[mu])
-            
+
         action = -2 * self.N * self.beta * (g.inner_product(J, z).real - z.grid.fsites * z.grid.nd)
         return action
 
-
     @differentiable_functional.multi_field_gradient
     def gradient(self, fields, dfields):
         def gradient_l(z, l, mu):
-            frc  = g.lattice(l[0])
-            frc @= 2 * self.beta * self.N * g.component.imag(g.trace(z * g.adj(g.cshift(z, mu, +1))) * l[mu])
+            frc = g.lattice(l[0])
+            frc @= (
+                2
+                * self.beta
+                * self.N
+                * g.component.imag(g.trace(z * g.adj(g.cshift(z, mu, +1))) * l[mu])
+            )
             frc.otype = l[0].otype.cartesian()
             return frc
-         
+
         def gradient_z(z, l):
             J = g.lattice(z)
             J[:] = 0.0
             for mu in range(z.grid.nd):
-                J += g.cshift(z, mu, +1) * g.adj(l[mu]) 
+                J += g.cshift(z, mu, +1) * g.adj(l[mu])
                 J += g.cshift(z * l[mu], mu, -1)
 
             frc = g.lattice(z)
-            frc @= - 2 * self.beta * self.N * J
+            frc @= -2 * self.beta * self.N * J
 
             frc -= g.trace(frc * g.adj(z)) * z
             return frc
-        
+
         z, l = self.split(fields)
         frc = []
         for df in g.core.util.to_list(dfields):
@@ -75,7 +79,7 @@ def gradient_z(z, l):
             if k == 0:
                 frc.append(gradient_z(z, l))
             else:
-                frc.append(gradient_l(z,l,k-1))
+                frc.append(gradient_l(z, l, k - 1))
         return frc
 
     # https://arxiv.org/abs/1102.1852
@@ -84,25 +88,25 @@ def constrained_leap_frog(self, eps, z, mom_z):
         def dot(v1, v2):
             return g.trace(v2 * g.adj(v1))
 
-        n  = g.real(z.grid)
-        n @= g.component.sqrt(g.component.real(dot(mom_z, mom_z))) 
+        n = g.real(z.grid)
+        n @= g.component.sqrt(g.component.real(dot(mom_z, mom_z)))
 
         # z'      =  cos(alpha) z + (1/|pi|) sin(alpha) mom_z
         # mom_z'  = -|pi| sin(alpha) z + cos(alpha) mom_z
         # alpha = eps |pi|
         _z = g.lattice(z)
         _z @= z
-        
-        cos  = g.real(z.grid) 
+
+        cos = g.real(z.grid)
         cos @= g.component.cos(eps * n)
 
-        sin  = g.real(z.grid)
+        sin = g.real(z.grid)
         sin @= g.component.sin(eps * n)
-        
+
         z @= cos * _z + g(g.component.inv(n) * sin) * mom_z
-        mom_z @= - g(n * sin) * _z + cos * mom_z
+        mom_z @= -g(n * sin) * _z + cos * mom_z
         del _z, cos, sin, n
-        
+
     # https://arxiv.org/abs/1102.1852
     def draw(self, field, rng, constraint=None):
         if constraint is None:
@@ -115,4 +119,4 @@ def draw(self, field, rng, constraint=None):
             z = constraint
             rng.normal_element(mom_z)
             # TO DO change to z * g(g.adj(z) * mom_z)
-            mom_z @= mom_z - g(z * g.adj(z)) * mom_z
\ No newline at end of file
+            mom_z @= mom_z - g(z * g.adj(z)) * mom_z

From 936490d92b3542c314d07138519a284533ecfaee Mon Sep 17 00:00:00 2001
From: Mattia Bruno <38449948+mbruno46@users.noreply.github.com>
Date: Fri, 20 Oct 2023 17:57:00 +0200
Subject: [PATCH 15/17] added support site inner product vcomplex

---
 lib/cgpt/lib/expression/unary.h               |  2 +-
 .../expression_mul_double_iVSinglet10.cc      |  2 +-
 .../expression_mul_double_iVSinglet30.cc      |  2 +-
 .../expression_mul_double_iVSinglet4.cc       |  2 +-
 .../expression_mul_single_iVSinglet10.cc      |  2 +-
 .../expression_mul_single_iVSinglet30.cc      |  2 +-
 .../expression_mul_single_iVSinglet4.cc       |  2 +-
 .../iVSinglet/expression_mul.template         |  2 +-
 tests/core/vcomplex.py                        | 21 +++++++++++++++++++
 9 files changed, 29 insertions(+), 8 deletions(-)
 create mode 100644 tests/core/vcomplex.py

diff --git a/lib/cgpt/lib/expression/unary.h b/lib/cgpt/lib/expression/unary.h
index a06d9e60..c80bc7cf 100644
--- a/lib/cgpt/lib/expression/unary.h
+++ b/lib/cgpt/lib/expression/unary.h
@@ -78,7 +78,7 @@ cgpt_Lattice_base* lattice_lat(cgpt_Lattice_base* dst, bool ac, const A& lat, Co
 template<typename A>
 cgpt_Lattice_base* lattice_unary_lat(cgpt_Lattice_base* dst, bool ac, const A& la,int unary_expr,ComplexD coef) {
   if (unary_expr == 0) {
-    return lattice_lat(dst, ac, la, coef);
+    return lattice_lat(dst, ac, closure(ToSinglet(la)), coef);
   } else if (unary_expr == (BIT_SPINTRACE|BIT_COLORTRACE)) {
     return lattice_expr(dst, ac, coef*ToSinglet(trace(la)));
   } else if (unary_expr == BIT_SPINTRACE) {
diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc
index 5b98f13e..d139d1e7 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet10.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexD vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet10);
-  //  _INNER_PRODUCT_(iVSinglet10);
+  _INNER_PRODUCT_(iVSinglet10);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc
index 90bc07f7..492c1d37 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet30.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexD vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet30);
-  //  _INNER_PRODUCT_(iVSinglet30);
+  _INNER_PRODUCT_(iVSinglet30);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc
index 760f1028..7f632587 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_double_iVSinglet4.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexD vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet4);
-  //  _INNER_PRODUCT_(iVSinglet4);
+  _INNER_PRODUCT_(iVSinglet4);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc
index ff1b993d..f2be2da1 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet10.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexF vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet10);
-  //  _INNER_PRODUCT_(iVSinglet10);
+  _INNER_PRODUCT_(iVSinglet10);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc
index 8dd56db1..9391e0c7 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet30.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexF vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet30);
-  //  _INNER_PRODUCT_(iVSinglet30);
+  _INNER_PRODUCT_(iVSinglet30);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc
index 7abd5db9..904e7fae 100644
--- a/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc
+++ b/lib/cgpt/lib/instantiate/expression_mul_single_iVSinglet4.cc
@@ -41,6 +41,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef vComplexF vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet4);
-  //  _INNER_PRODUCT_(iVSinglet4);
+  _INNER_PRODUCT_(iVSinglet4);
   ERR("Not implemented");
 }
diff --git a/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template b/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template
index 092ee1b8..d7eb6f27 100644
--- a/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template
+++ b/lib/cgpt/lib/instantiate/templates/iVSinglet/expression_mul.template
@@ -40,6 +40,6 @@ cgpt_Lattice_base* cgpt_lattice_mul(cgpt_Lattice_base* dst, bool ac, int unary_a
   typedef {precision_vector} vtype;
   _COMPATIBLE_MSR_(iSinglet);
   _OUTER_PRODUCT_(iVSinglet{tensor_arg_1});
-  //  _INNER_PRODUCT_(iVSinglet{tensor_arg_1});
+  _INNER_PRODUCT_(iVSinglet{tensor_arg_1});
   ERR("Not implemented");
 }
diff --git a/tests/core/vcomplex.py b/tests/core/vcomplex.py
new file mode 100644
index 00000000..d1b5097f
--- /dev/null
+++ b/tests/core/vcomplex.py
@@ -0,0 +1,21 @@
+#!/usr/bin/env python3
+#
+# Authors: Christoph Lehner 2020, Mattia Bruno 2023
+#
+# Desc.: Illustrate core concepts and features
+#
+import gpt as g
+
+L = [8, 4, 4, 4]
+grid = g.grid(L, g.double)
+rng = g.random("vcomplex")
+
+def test_local_inner_product(v):
+    i0 = g(g.adj(v) * v)
+    i1 = g(g.trace(v * g.adj(v)))
+    assert g.norm2(i0 - i1) < 1e-15
+
+for N in [4, 10, 30]:
+    v1 = g.vcomplex(grid, N)
+    test_local_inner_product(v1)
+    del v1

From 504ce22c25346fab08f233572829c40dd821b15a Mon Sep 17 00:00:00 2001
From: Christoph Lehner <christoph@lhnr.de>
Date: Wed, 31 May 2023 08:12:30 +0300
Subject: [PATCH 16/17] timings

---
 lib/gpt/qcd/gauge/stencil/staple.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/lib/gpt/qcd/gauge/stencil/staple.py b/lib/gpt/qcd/gauge/stencil/staple.py
index 64b043dd..41820bd5 100644
--- a/lib/gpt/qcd/gauge/stencil/staple.py
+++ b/lib/gpt/qcd/gauge/stencil/staple.py
@@ -34,6 +34,7 @@ def staple_sum(U, rho, mu=None, cache=default_staple_cache):
 
     Ntarget = len(target_mu)
 
+    t = g.timer("staples")
     if tag not in cache:
         code = []
         nwr = [0] * Ntarget

From a9b32480c1dbe31c370e70201e50c47df45a3497 Mon Sep 17 00:00:00 2001
From: Christoph Lehner <christoph@lhnr.de>
Date: Wed, 31 May 2023 16:03:50 +0200
Subject: [PATCH 17/17] new parallel_transport_matrix

---
 lib/gpt/qcd/gauge/stencil/staple.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/lib/gpt/qcd/gauge/stencil/staple.py b/lib/gpt/qcd/gauge/stencil/staple.py
index 41820bd5..64b043dd 100644
--- a/lib/gpt/qcd/gauge/stencil/staple.py
+++ b/lib/gpt/qcd/gauge/stencil/staple.py
@@ -34,7 +34,6 @@ def staple_sum(U, rho, mu=None, cache=default_staple_cache):
 
     Ntarget = len(target_mu)
 
-    t = g.timer("staples")
     if tag not in cache:
         code = []
         nwr = [0] * Ntarget