CP^N-1 model

applications/hmc/cpn.py

@@ -0,0 +1,125 @@
+#!/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
+
+g.default.set_verbose("step_size", False)
+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])}")

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");
     }

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) {

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");
 }

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");
 }

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");
 }

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");
 }

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");
 }

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");
 }

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");
 }

lib/gpt/core/domain/sparse.py

@@ -104,18 +104,30 @@ def slice(self, fields, ortho_dim):
         length = self.grid.gdimensions[ortho_dim]
         return gpt.indexed_sum(fields, self.coordinate_lattices()[ortho_dim], length)
 
-    def weight(self):
-        if self.weight_cache is not None:
-            return self.weight_cache
+    def global_coordinates(self):
+        # broadcast below uses grid.processor not gpt.rank()
+        rank_id = self.grid.processor
 
         # this function is mostly used for tests, so it is not performance critical
         for rank in range(self.grid.Nprocessors):
-            local_coordinates = np.copy(self.local_coordinates)
+            nlc = np.array([len(self.local_coordinates)])
+            self.grid.broadcast(rank, nlc)
+            if rank == rank_id:
+                local_coordinates = np.copy(self.local_coordinates)
+            else:
+                local_coordinates = np.zeros((nlc[0], self.grid.nd), dtype=self.local_coordinates.dtype)
             self.grid.broadcast(rank, local_coordinates)
             if rank == 0:
                 global_coordinates = local_coordinates
             else:
                 global_coordinates = np.concatenate((global_coordinates, local_coordinates))
+        return global_coordinates
+
+    def weight(self):
+        if self.weight_cache is not None:
+            return self.weight_cache
+
+        global_coordinates = self.global_coordinates()
 
         unique_coordinates, count = np.unique(global_coordinates, axis=0, return_counts=True)
         unique_coordinates = unique_coordinates.view(type(global_coordinates))
@@ -208,6 +220,27 @@ def converted(self, precision):
             mask=mask,
         )
 
+    # every rank must pass the same list of positions, otherwise code breaks
+    def restricted(self, positions):
+        # get emb_coor of wanted positions on each rank to build mask
+        ec = self.unique_embedded_coordinates(positions)
+
+        # creates full mask
+        mask = self.lattice(gpt.ot_real_additive_group())
+        mask[:] = 0
+        mask[ec] = 1
+
+        cl = self.coordinate_lattices(mark_empty=-1)
+        local_coordinates = np.hstack(tuple([x[:].real.astype(np.int32) for x in cl]))
+        mask_local = [m[0]==1 for m in mask[:]]
+
+        return gpt.core.domain.sparse(
+            self.kernel.grid,
+            local_coordinates,
+            dimensions_divisible_by=cl[0].grid.fdimensions,
+            mask = mask_local
+        )
+
     def weight(self):
         return self.kernel.weight()
 

lib/gpt/qcd/scalar/action/__init__.py

@@ -18,6 +18,7 @@
 #    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, fourier_mass_term, general_mass_term
 from gpt.qcd.scalar.action.coupling import coupling
 import gpt.qcd.scalar.action.hermitian_kernel