Fix plaquette action

src/action.hpp

@@ -125,59 +125,53 @@ inline Scalar S_H(const array::ArrayView<ConstScalar, Shape>& Qmat, const array:
     using array::Array;
     using array::Shape;
 
-    /* Neighbor plaquettes: (nth,dir) -> (di,dj)
+    /* Plaquettes: (dir) -> (di,dj)
      * where 'dir' is the direction index of the link,
      * starting from bottom and cycling counterclockwise.
      */
-    const int plaqc[4][4][2] = {
-        {{0,0}, {1,0}, {1,1}, {0,1}},
-        {{0,-1}, {1,-1}, {1,0}, {0,0}},
-        {{-1,0}, {0,0}, {0,1}, {-1,1}},
-        {{-1,-1}, {0,-1}, {0,0}, {0,1}},
-    };
+    const int plaq[4][2] = {{0,0}, {1,0}, {1,1}, {0,1}};
+
+    const size_t invdir[4] = {2, 3, 0, 1};
 
     Scalar S{0.0};
 
     const int nx = Qmat.dim(0);
     const int ny = Qmat.dim(1);
 
-    for (int i = 0; i < nx; ++i) {
-        for (int j = 0; j < ny; ++j) {
-
-            /* Sum over neighboring plaquettes */
-            for (int iplaq = 0; iplaq < 4; ++iplaq) {
-                const int i1 = (int)i + plaqc[iplaq][0][0];
-                const int j1 = (int)j + plaqc[iplaq][0][1];
-                const int i2 = (int)i + plaqc[iplaq][1][0];
-                const int j2 = (int)j + plaqc[iplaq][1][1];
-                const int i3 = (int)i + plaqc[iplaq][2][0];
-                const int j3 = (int)j + plaqc[iplaq][2][1]; 
-                const int i4 = (int)i + plaqc[iplaq][3][0];
-                const int j4 = (int)j + plaqc[iplaq][3][1];
-
-                if (i1 < 0 || i1 >= nx || j1 < 0 || j1 >= ny)
-                    continue;
-                if (i2 < 0 || i2 >= nx || j2 < 0 || j2 >= ny)
-                    continue;
-                if (i3 < 0 || i3 >= nx || j3 < 0 || j3 >= ny)
-                    continue;
-                if (i4 < 0 || i4 >= nx || j4 < 0 || j4 >= ny)
-                    continue;
+    for (int i = 0; i < nx-1; ++i) {
+        for (int j = 0; j < ny-1; ++j) {
+            /* Plaquette with lower left corner at (i,j) */
+            
+            /* Sum over cyclic permutations */
+            for (int iperm = 0; iperm < 4; ++iperm) {
+                const int k1 = (iperm+0)%4;
+                const int k2 = (iperm+1)%4;
+                const int k3 = (iperm+2)%4;
+                const int k4 = (iperm+3)%4;
+                
+                const int i1 = (int)i + plaq[k1][0];
+                const int j1 = (int)j + plaq[k1][1];
+                const int i2 = (int)i + plaq[k2][0];
+                const int j2 = (int)j + plaq[k2][1];
+                const int i3 = (int)i + plaq[k3][0];
+                const int j3 = (int)j + plaq[k3][1]; 
+                const int i4 = (int)i + plaq[k4][0];
+                const int j4 = (int)j + plaq[k4][1];
 
                 auto Q1 = Qmat.part(i1,j1).matrix();
                 auto Q2 = Qmat.part(i2,j2).matrix();
-                auto Q3 = Qmat.part(i3,j3).matrix();
+                //auto Q3 = Qmat.part(i3,j3).matrix();
                 auto Q4 = Qmat.part(i4,j4).matrix();
 
-                auto U21 = U.part(i1,j1,0).matrix();
-                auto U32 = U.part(i2,j2,1).matrix();
-                auto U43 = U.part(i3,j3,2).matrix();
-                auto U14 = U.part(i4,j4,3).matrix();
+                auto U21 = U.part(i1,j1,k1).matrix();
+                auto U32 = U.part(i2,j2,k2).matrix();
+                auto U43 = U.part(i3,j3,k3).matrix();
+                auto U14 = U.part(i4,j4,k4).matrix();
 
-                auto U12 = U.part(i2,j2,3).matrix();
-                auto U23 = U.part(i3,j3,2).matrix();
-                auto U34 = U.part(i4,j4,0).matrix();
-                auto U41 = U.part(i1,j1,1).matrix();
+                auto U12 = U.part(i2,j2,invdir[k1]).matrix();
+                //auto U23 = U.part(i3,j3,invdir[k2]).matrix();
+                //auto U34 = U.part(i4,j4,invdir[k3]).matrix();
+                auto U41 = U.part(i1,j1,invdir[k4]).matrix();
 
                 auto P = U14*U43*U32*U21;
 
@@ -188,7 +182,7 @@ inline Scalar S_H(const array::ArrayView<ConstScalar, Shape>& Qmat, const array:
                     0, 0, 1, 0,
                     0, 0, 0, 1;
 
-                S += ((P - I) * (Q1 - nambu_diag_mul(U12, Q2, U21) * nambu_diag_mul(U14, Q4, U41))).trace();
+                S += ((P - I) * (Q1 - nambu_diag_mul(U12, Q2, U21) * Q1 * nambu_diag_mul(U14, Q4, U41))).trace();
             }
         }
     }