From 7356835026b58abc6687efd53725825212bc3536 Mon Sep 17 00:00:00 2001
From: JingJing Li <lijingjing416@gmail.com>
Date: Mon, 13 Apr 2026 16:35:34 +0200
Subject: [PATCH 1/8] rebase roman's app_context (direct file copy)

---
 ...-ptbc-hmc-rgmixedcg-tmcloverdetratio.input | 117 +++++
 src/bin/benchmark.c                           |  28 +-
 src/bin/deriv_mg_tune.c                       |   2 +-
 src/bin/hmc_tm.c                              |   6 +-
 src/bin/invert.c                              |   2 +-
 src/bin/offline_measurement.c                 |   2 +-
 src/bin/tests/hopping_test.c                  |   2 +-
 src/bin/tests/qphix_test_Dslash.c             |  16 +-
 src/bin/tests/test_eigenvalues.c              |   8 +-
 src/lib/CMakeLists.txt                        |  12 +
 src/lib/DDalphaAMG_interface.c                |   8 +-
 src/lib/alloc/alloc.c                         |  81 ++++
 src/lib/app_context/app_context.c             | 235 ++++++++++
 src/lib/compare_derivative.c                  |   4 +-
 src/lib/fatal_error.c                         |   2 +-
 src/lib/include/alloc.h                       |  49 +++
 src/lib/include/app.h                         | 108 +++++
 src/lib/include/mpi.h                         |  35 ++
 src/lib/init/init_parallel.c                  |   2 +-
 src/lib/io/deri_write_stdout.c                |   2 +-
 src/lib/io/dml.c                              |   2 +-
 src/lib/io/eospinor_read.c                    |   2 +-
 src/lib/io/eospinor_write.c                   |   8 +-
 src/lib/io/gauge_read_binary.c                |   2 +-
 src/lib/io/gauge_write.c                      |   2 +-
 src/lib/io/gauge_write_binary.c               |   2 +-
 src/lib/io/io_cm.c                            |   2 +-
 src/lib/io/spinor_read_binary.c               |   4 +-
 src/lib/io/spinor_write_binary.c              |   4 +-
 src/lib/io/spinor_write_stdout.c              |   2 +-
 src/lib/io/sw_write_stdout.c                  |   2 +-
 src/lib/io/utils_kill_with_error.c            |   2 +-
 src/lib/linalg/assign_mul_add_r_and_square.c  |   2 +-
 src/lib/linalg/diff_and_square_norm.c         |   2 +-
 src/lib/linalg/scalar_prod_body.inc           |   4 +-
 src/lib/linalg/scalar_prod_i.c                |   2 +-
 src/lib/linalg/scalar_prod_r.c                |   2 +-
 src/lib/linalg/scalar_prod_r_32.c             |   2 +-
 src/lib/linalg/square_and_minmax.c            |  32 +-
 src/lib/linalg/square_and_prod_r.c            |   4 +-
 src/lib/linalg/square_norm.c                  |   4 +-
 src/lib/linalg/square_norm_32.c               |   4 +-
 src/lib/meas/correlators.c                    |   4 +-
 ...easure_clover_field_strength_observables.c |   4 +-
 src/lib/meas/oriented_plaquettes.c            |   2 +-
 src/lib/meas/pion_norm.c                      |   2 +-
 src/lib/meas/polyakov_loop.c                  |   4 +-
 src/lib/measure_gauge_action.c                |   4 +-
 src/lib/measure_rectangles.c                  |   2 +-
 src/lib/monomial/moment_energy.c              |   2 +-
 src/lib/monomial/monitor_forces.c             |   4 +-
 src/lib/mpi_init.c                            |  12 +-
 src/lib/operator/clover_det.c                 |   4 +-
 src/lib/prepare_source.c                      |   8 +-
 src/lib/qphix/qphix_interface.cpp             |   4 +-
 src/lib/quda_interface.c                      |  14 +-
 src/lib/reweighting_factor.c                  |   2 +-
 src/lib/sighandler.c                          |   2 +-
 src/lib/solver/dfl_projector.c                |   6 +-
 src/lib/solver/gcr4complex_body.inc           |   6 +-
 src/lib/solver/little_project_eo_body.inc     |   2 +-
 src/lib/start.c                               |  12 +-
 src/lib/test/check_xchange.c                  | 400 +++++++++---------
 src/lib/test/measure_rectangles.debug.c       |   2 +-
 src/lib/update_tm.c                           |   4 +-
 src/lib/wrapper/lib_wrapper.c                 |   4 +-
 tests/test_buffers.c                          |   2 +-
 67 files changed, 981 insertions(+), 344 deletions(-)
 create mode 100644 doc/sample-input/sample-ptbc-hmc-rgmixedcg-tmcloverdetratio.input
 create mode 100644 src/lib/alloc/alloc.c
 create mode 100644 src/lib/app_context/app_context.c
 create mode 100644 src/lib/include/alloc.h
 create mode 100644 src/lib/include/app.h
 create mode 100644 src/lib/include/mpi.h

diff --git a/doc/sample-input/sample-ptbc-hmc-rgmixedcg-tmcloverdetratio.input b/doc/sample-input/sample-ptbc-hmc-rgmixedcg-tmcloverdetratio.input
new file mode 100644
index 000000000..89ff0511f
--- /dev/null
+++ b/doc/sample-input/sample-ptbc-hmc-rgmixedcg-tmcloverdetratio.input
@@ -0,0 +1,117 @@
+# this is a sample input file for cloverdet + cloverdetratio using
+# DDalphaAMG as a solver
+
+L=8
+T=8
+Measurements = 20
+Startcondition = hot
+2KappaMu = 0.01
+CSW = 1.00
+kappa = 0.138
+NSave = 500000
+ThetaT = 1.0
+UseEvenOdd = yes
+ReversibilityCheck = yes
+ReversibilityCheckIntervall = 4
+InitialStoreCounter = 0
+DebugLevel = 1
+
+BeginPTBC
+
+  BeginPTBCDefect 0
+    Extents = 2, 2, 2
+    Along = Y
+  EndPTBCDefect
+
+  BeginPTBCInstance 0
+    Coefficients = 1.0
+    Defects = 0
+  EndPTBCInstance
+
+  BeginPTBCInstance 1
+    Coefficients = 0.0
+    Defects = 0
+  EndPTBCInstance
+
+  BeginPTBCInstance 2
+    Coefficients = 0.75
+    Defects = 0
+  EndPTBCInstance
+
+  BeginPTBCInstance 3
+    Coefficients = 0.5
+    Defects = 0
+  EndPTBCInstance
+
+EndPTBC
+
+
+
+# since this is a test file, we employ the reproducible random numbers mode
+ReproduceRandomNumbers = yes
+Seed = 127782
+
+BeginMeasurement CORRELATORS
+  Frequency = 2
+EndMeasurement
+
+BeginMonomial GAUGE
+  Type = Wilson
+  beta = 5.60
+  Timescale = 0
+EndMonomial
+
+BeginMonomial CLOVERDET
+  Timescale = 1
+  2KappaMu = 0.01
+  CSW = 1.00
+  # nominator shift
+  rho = 0.1
+  kappa = 0.138
+  AcceptancePrecision =  1.e-20
+  ForcePrecision = 1.e-14
+  Name = cloverdet
+  solver = rgmixedcg
+  usesloppyprecision = single
+EndMonomial
+
+BeginMonomial CLOVERDETRATIO
+  Timescale = 2
+  2KappaMu = 0.01
+  # nominator shift
+  rho = 0.0
+  # denominator shift, should match CLOVERDET shift
+  rho2 = 0.1
+  CSW = 1.00
+  kappa = 0.138
+  AcceptancePrecision =  1.e-20
+  ForcePrecision = 1.e-16
+  Name = cloverdetratio
+  solver = rgmixedcg
+  usesloppyprecision = single
+EndMonomial
+
+BeginIntegrator 
+  Type0 = 2MNFG
+  Type1 = 2MNFG
+  Type2 = 2MNFG
+  IntegrationSteps0 = 1
+  IntegrationSteps1 = 1
+  IntegrationSteps2 = 4
+  tau = 1.00
+  Lambda0 = 0.16666667
+  Lambda1 = 0.16666667
+  Lambda2 = 0.16666667
+  NumberOfTimescales = 3
+EndIntegrator
+
+BeginOperator CLOVER
+  2KappaMu = 0.01
+  CSW = 1.00
+  kappa = 0.138
+  SolverPrecision = 1e-16
+  MaxSolverIterations = 1000
+  useevenodd = yes
+  solver = rgmixedcg
+  usesloppyprecision = single
+EndOperator
\ No newline at end of file
diff --git a/src/bin/benchmark.c b/src/bin/benchmark.c
index b2f4ee68c..f0b15ce24 100644
--- a/src/bin/benchmark.c
+++ b/src/bin/benchmark.c
@@ -102,7 +102,7 @@ int main(int argc, char *argv[]) {
 #else
   MPI_Init(&argc, &argv);
 #endif
-  MPI_Comm_rank(MPI_COMM_WORLD, &g_proc_id);
+  MPI_Comm_rank(app()->mpi.comm, &g_proc_id);
 
 #else
   g_proc_id = 0;
@@ -234,7 +234,7 @@ int main(int argc, char *argv[]) {
     antioptaway = 0.0;
     /* compute approximately how many applications we need to do to get a reliable measurement */
 #ifdef TM_USE_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 #endif
     t1 = gettime();
     for (j = 0; j < j_max; j++) {
@@ -248,14 +248,14 @@ int main(int argc, char *argv[]) {
     // division by g_nproc because we will average over processes
     j = (int)(ceil(j_max * 31.0 / dt / g_nproc));
 #ifdef TM_USE_MPI
-    MPI_Allreduce(&j, &j_max, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&j, &j_max, 1, MPI_INT, MPI_SUM, app()->mpi.comm);
 #else
     j_max = j;
 #endif
 
     /* perform the actual benchmark */
 #ifdef TM_USE_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 #endif
     t1 = gettime();
     antioptaway = 0.0;
@@ -268,14 +268,14 @@ int main(int argc, char *argv[]) {
     }
     dt = gettime() - t1;
 #ifdef TM_USE_MPI
-    MPI_Allreduce(&dt, &sdt, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&dt, &sdt, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 #else
     sdt = dt;
 #endif
 
     qdt = dt * dt;
 #ifdef TM_USE_MPI
-    MPI_Allreduce(&qdt, &sqdt, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&qdt, &sqdt, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 #else
     sqdt = qdt;
 #endif
@@ -321,9 +321,9 @@ int main(int argc, char *argv[]) {
     dt2 = t2 - t1;
     /* compute the bandwidth */
     dt = dts - dt2;
-    MPI_Allreduce(&dt, &sdt, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&dt, &sdt, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
     sdt = sdt / ((double)g_nproc);
-    MPI_Allreduce(&dt2, &dt, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&dt2, &dt, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
     dt = dt / ((double)g_nproc);
     dt = 1.0e6f * dt / ((double)(k_max * j_max * (VOLUME)));
     if (g_proc_id == 0) {
@@ -365,7 +365,7 @@ int main(int argc, char *argv[]) {
 
     /* estimate a reasonable number of applications to get a reliable measurement */
 #ifdef TM_USE_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 #endif
     t1 = gettime();
     for (j = 0; j < j_max; j++) {
@@ -379,14 +379,14 @@ int main(int argc, char *argv[]) {
     // division by g_nproc because we will average over processes using  MPI_SUM
     j = (int)(ceil(j_max * 31.0 / dt / g_nproc));
 #ifdef TM_USE_MPI
-    MPI_Allreduce(&j, &j_max, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&j, &j_max, 1, MPI_INT, MPI_SUM, app()->mpi.comm);
 #else
     j_max = j;
 #endif
 
     /* perform the actual measurement */
 #ifdef TM_USE_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 #endif
     t1 = gettime();
     for (j = 0; j < j_max; j++) {
@@ -398,13 +398,13 @@ int main(int argc, char *argv[]) {
     t2 = gettime();
     dt = t2 - t1;
 #ifdef TM_USE_MPI
-    MPI_Allreduce(&dt, &sdt, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&dt, &sdt, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 #else
     sdt = dt;
 #endif
     qdt = dt * dt;
 #ifdef TM_USE_MPI
-    MPI_Allreduce(&qdt, &sqdt, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&qdt, &sqdt, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 #else
     sqdt = qdt;
 #endif
@@ -451,7 +451,7 @@ int main(int argc, char *argv[]) {
   free_spinor_field();
   free_moment_field();
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
   MPI_Finalize();
 #endif
   return (0);
diff --git a/src/bin/deriv_mg_tune.c b/src/bin/deriv_mg_tune.c
index 0fd6ec529..d0e816309 100644
--- a/src/bin/deriv_mg_tune.c
+++ b/src/bin/deriv_mg_tune.c
@@ -357,7 +357,7 @@ int main(int argc, char *argv[]) {
   _endQuda();
 #endif
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
   MPI_Finalize();
 #endif
 
diff --git a/src/bin/hmc_tm.c b/src/bin/hmc_tm.c
index 017ab9cb3..6e3654c2c 100644
--- a/src/bin/hmc_tm.c
+++ b/src/bin/hmc_tm.c
@@ -479,7 +479,7 @@ int main(int argc, char *argv[]) {
 
           sleep(io_timeout);
 #ifdef TM_USE_MPI
-          MPI_Barrier(MPI_COMM_WORLD);
+          MPI_Barrier(app()->mpi.comm);
 #endif
         }
       /* Now move .conf.tmp into place */
@@ -528,7 +528,7 @@ int main(int argc, char *argv[]) {
     }
 
 #ifdef TM_USE_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 #endif
     if (ix == 0 && g_proc_id == 0) {
       countfile = fopen("history_hmc_tm", "a");
@@ -581,7 +581,7 @@ int main(int argc, char *argv[]) {
   _endQuda();
 #endif
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
   MPI_Finalize();
 #endif
 
diff --git a/src/bin/invert.c b/src/bin/invert.c
index b5040ba88..7e004188c 100644
--- a/src/bin/invert.c
+++ b/src/bin/invert.c
@@ -448,7 +448,7 @@ int main(int argc, char *argv[]) {
   _endQuda();
 #endif
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
   MPI_Finalize();
 #endif
   return (0);
diff --git a/src/bin/offline_measurement.c b/src/bin/offline_measurement.c
index 667dce3f9..f932019fe 100644
--- a/src/bin/offline_measurement.c
+++ b/src/bin/offline_measurement.c
@@ -294,7 +294,7 @@ int main(int argc, char *argv[]) {
   free(input_filename);
 
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
   MPI_Finalize();
 #endif
   return (0);
diff --git a/src/bin/tests/hopping_test.c b/src/bin/tests/hopping_test.c
index 0e5ff03e7..a09c0e46f 100644
--- a/src/bin/tests/hopping_test.c
+++ b/src/bin/tests/hopping_test.c
@@ -318,7 +318,7 @@ int main(int argc, char *argv[]) {
     }
 
 #ifdef TM_USE_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     MPI_Finalize();
 #endif
   }
diff --git a/src/bin/tests/qphix_test_Dslash.c b/src/bin/tests/qphix_test_Dslash.c
index 41e2602a4..b998fe52f 100644
--- a/src/bin/tests/qphix_test_Dslash.c
+++ b/src/bin/tests/qphix_test_Dslash.c
@@ -221,7 +221,7 @@ int main(int argc, char* argv[]) {
                    source_location, 12345 /* seed */);
 
 #ifdef TM_USE_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 #endif
 
     tm_t1 = gettime();
@@ -231,7 +231,7 @@ int main(int argc, char* argv[]) {
     tm_t2 = gettime();
 
 #ifdef TM_USE_MPI
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 #endif
     q_t1 = gettime();
     Mfull_qphix(qphix_out_cb_spinors[0], qphix_out_cb_spinors[1], op->sr0, op->sr1, op->type);
@@ -297,7 +297,7 @@ int main(int argc, char* argv[]) {
   free_spinor_field();
   free_moment_field();
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
   MPI_Finalize();
 #endif
   return (failed);
@@ -305,7 +305,7 @@ int main(int argc, char* argv[]) {
 
 double compare_spinors(spinor* s1, spinor* s2) {
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
 #endif
   int coords[4];
   int x, y, z, t, id = 0;
@@ -352,7 +352,7 @@ double compare_spinors(spinor* s1, spinor* s2) {
               }
             }
 #ifdef TM_USE_MPI
-            MPI_Barrier(MPI_COMM_WORLD);
+            MPI_Barrier(app()->mpi.comm);
 #endif
           }  // z
         }  // y
@@ -361,7 +361,7 @@ double compare_spinors(spinor* s1, spinor* s2) {
   }  // if( SourceInfo.type == SRC_TYPE_POINT )
 
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
 #endif
   if (g_proc_id == 0) {
     printf("\n");
@@ -376,7 +376,7 @@ double compare_spinors(spinor* s1, spinor* s2) {
   double squarenorm = diff_and_square_norm(s1, s2, VOLUME);
 
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
 #endif
   id = 0;
   for (int t_global = 0; t_global < g_nproc_t * T; t_global++) {
@@ -408,7 +408,7 @@ double compare_spinors(spinor* s1, spinor* s2) {
             }
           }
 #ifdef TM_USE_MPI
-          MPI_Barrier(MPI_COMM_WORLD);
+          MPI_Barrier(app()->mpi.comm);
 #endif
         }  // z
       }  // y
diff --git a/src/bin/tests/test_eigenvalues.c b/src/bin/tests/test_eigenvalues.c
index c52d29cf8..6d162117e 100644
--- a/src/bin/tests/test_eigenvalues.c
+++ b/src/bin/tests/test_eigenvalues.c
@@ -390,14 +390,14 @@ int main(int argc, char *argv[]) {
       }
       rlxd_get(rlxd_state);
 #ifdef TM_USE_MPI
-      MPI_Send(&rlxd_state[0], 105, MPI_INT, 1, 99, MPI_COMM_WORLD);
-      MPI_Recv(&rlxd_state[0], 105, MPI_INT, g_nproc - 1, 99, MPI_COMM_WORLD, &status);
+      MPI_Send(&rlxd_state[0], 105, MPI_INT, 1, 99, app()->mpi.comm);
+      MPI_Recv(&rlxd_state[0], 105, MPI_INT, g_nproc - 1, 99, app()->mpi.comm, &status);
       rlxd_reset(rlxd_state);
 #endif
     }
 #ifdef TM_USE_MPI
     else {
-      MPI_Recv(&rlxd_state[0], 105, MPI_INT, g_proc_id - 1, 99, MPI_COMM_WORLD, &status);
+      MPI_Recv(&rlxd_state[0], 105, MPI_INT, g_proc_id - 1, 99, app()->mpi.comm, &status);
       rlxd_reset(rlxd_state);
       /* hot */
       if (startoption == 1) {
@@ -408,7 +408,7 @@ int main(int argc, char *argv[]) {
         k = 0;
       }
       rlxd_get(rlxd_state);
-      MPI_Send(&rlxd_state[0], 105, MPI_INT, k, 99, MPI_COMM_WORLD);
+      MPI_Send(&rlxd_state[0], 105, MPI_INT, k, 99, app()->mpi.comm);
     }
 #endif
 
diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt
index 6974f4523..9c5322997 100644
--- a/src/lib/CMakeLists.txt
+++ b/src/lib/CMakeLists.txt
@@ -68,6 +68,16 @@ list(
   init/init_stout_smear_vars.c
   init/init_moment_field.c)
 
+list(
+  APPEND
+  APP_CONTEXT_SRC_C
+  app_context/app_context.c)
+
+list(
+  APPEND
+  ALLOC_SRC_C
+  alloc/alloc.c)
+
 list(
   APPEND
   SOLVER_SRC_C
@@ -368,6 +378,8 @@ list(
   ${LINALG_SRC_C}
   ${IO_SRC_C}
   ${INIT_SRC_C}
+  ${APP_CONTEXT_SRC_C}
+  ${ALLOC_SRC_C}
   ${SOLVER_SRC_C}
   ${TEST_SRC_C}
   ${MEAS_SRC_C}
diff --git a/src/lib/DDalphaAMG_interface.c b/src/lib/DDalphaAMG_interface.c
index a4b1d61cf..08b68ca07 100644
--- a/src/lib/DDalphaAMG_interface.c
+++ b/src/lib/DDalphaAMG_interface.c
@@ -347,7 +347,7 @@ static int MG_pre_solve(su3 **gf) {
     MG_init();
     mg_initialized = 1;
     if (g_proc_id == 0) printf("DDalphaAMG initialized\n");
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
   }
 
   if (mg_update_gauge == 1) {
@@ -1240,7 +1240,7 @@ int MG_solver(spinor *const phi_new, spinor *const phi_old, const double precisi
     if (g_proc_id == 0) printf("ERROR: solver didn't converge after two trials!! Aborting... \n");
     // TODO: handle abort
     DDalphaAMG_finalize();
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     MPI_Finalize();
     exit(1);
   }
@@ -1316,7 +1316,7 @@ int MG_solver_nd(spinor *const up_new, spinor *const dn_new, spinor *const up_ol
     if (g_proc_id == 0) printf("ERROR: solver didn't converge after two trials!! Aborting... \n");
     // TODO: handle abort
     DDalphaAMG_finalize();
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     MPI_Finalize();
     exit(1);
   }
@@ -1390,7 +1390,7 @@ int MG_mms_solver_nd(spinor **const up_new, spinor **const dn_new, spinor *const
     if (g_proc_id == 0) printf("ERROR: solver didn't converge after two trials!! Aborting... \n");
     // TODO: handle abort
     DDalphaAMG_finalize();
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     MPI_Finalize();
     exit(1);
   }
diff --git a/src/lib/alloc/alloc.c b/src/lib/alloc/alloc.c
new file mode 100644
index 000000000..a853d58d4
--- /dev/null
+++ b/src/lib/alloc/alloc.c
@@ -0,0 +1,81 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2026 Roman Gruber
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD 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 tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Allocation utils
+ *
+ * Author: Roman Gruber
+ *         roman.gruber@unibe.ch
+ *
+ *******************************************************************************/
+
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include <errno.h>
+
+
+/**
+ * @brief      Safe malloc implementation that checks malloc for NULL. Never use
+ *             this function instead use the macro safe_malloc()
+ *
+ * @param[in]  size  The allocation size in bytes
+ * @param[in]  file  __FILE__
+ * @param[in]  line  __LINE__
+ * @param[in]  func  __func__
+ *
+ * @return     Pointer returned by malloc()
+ */
+void *safe_malloc_impl(size_t size, const char *file, int line, const char *func)
+{
+    if (size <= 0) {
+        fprintf(stderr, "safe_malloc: zero-size allocation at %s:%d (%s)\n", file, line, func);
+        abort();
+    }
+
+    void *p = malloc(size);
+    if (p == NULL) {
+        fprintf(stderr, "safe_malloc: failed to allocate %zu bytes at %s:%d (%s): %s\n",
+                size, file, line, func, strerror(errno));
+        abort();
+    }
+    return p;
+}
+
+
+/**
+ * @brief      Identical to safe_malloc_impl above just that is call calloc
+ *             instead of malloc. Never use this function instead use the macro
+ *             safe_calloc()
+ */
+void *safe_calloc_impl(size_t size, const char *file, int line, const char *func)
+{
+   if (size <= 0) {
+      fprintf(stderr, "safe_calloc: zero-size allocation at %s:%d (%s)\n", file, line, func);
+      abort();
+   }
+
+   void *p = calloc(size, 1);
+   if (p == NULL) {
+      fprintf(stderr, "safe_calloc: failed to allocate %zu bytes at %s:%d (%s): %s\n",
+              size, file, line, func, strerror(errno));
+      abort();
+   }
+   return p;
+}
\ No newline at end of file
diff --git a/src/lib/app_context/app_context.c b/src/lib/app_context/app_context.c
new file mode 100644
index 000000000..0d2e514fa
--- /dev/null
+++ b/src/lib/app_context/app_context.c
@@ -0,0 +1,235 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2026 Roman Gruber
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD 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 tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * App context module
+ *
+ * Author: Roman Gruber
+ *         roman.gruber@unibe.ch
+ *
+ *******************************************************************************/
+
+#include <stdbool.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <mpi.h>
+#include "fatal_error.h"
+
+
+/**
+ * @brief      The rank topology struct
+ */
+typedef struct {
+    int number_of_nodes;    // total number of nodes in the job
+    int number_of_ranks;    // total number of processes in the job, i.e. size of the world communicator
+    int ranks_per_node;     // number of ranks per node
+    int node_index;         // index enumerating the node (unique per node)
+    int node_rank;          // rank number inside the node
+} RankTopology;
+
+
+/**
+ * @brief      Gracefully error with function name, file and line number along
+ *             the error message
+ *
+ * @param      test  The test
+ * @param      ...   Format parameters
+ */
+#define err(test, ...) err_impl(test, __func__, __FILE__, __LINE__, __VA_ARGS__)
+static void err_impl(const bool test, const char* func, const char* file, const int line, const char* format, ...)
+{
+    if (test) {
+        va_list args;
+        char message[1024];
+        va_start(args, format);
+        vsnprintf(message, 1024, format, args);
+        va_end(args);
+        char location[1024];
+        snprintf(location, 1024, "%s:%d %s", file, line, func);
+        fatal_error(message, location);
+    }
+}
+
+
+static void initialize(void);
+static AppContext app_instance = {
+    .mpi = {
+        .comm = MPI_COMM_WORLD, // default communicator
+        .world_comm = MPI_COMM_WORLD,
+    },
+    .ptbc = {
+        .instance_id = 0,
+        .n_instances = 1,
+        .n_defects = 0,
+        .active = false,
+        .initialize = initialize,
+        .instances = {{.active = false}},
+        .defects = {{.active = false}}
+    }
+};
+
+
+/**
+ * @brief      Return the global *immutable* application context struct. To be
+ *             used when reading parameters.
+ *
+ * @return     Global application context struct
+ */
+const AppContext* app(void)
+{
+    return &app_instance;
+}
+
+
+/**
+ * @brief      Return the global *mutable* application context struct. To be
+ *             used when initializing/setting parameters.
+ *
+ * @return     Global application context struct
+ */
+AppContext* appm(void)
+{
+    return &app_instance;
+}
+
+
+/**
+ * @brief      Return rank topology.
+ *
+ * @return     The topology.
+ */
+static RankTopology get_topology(void)
+{
+    int world_rank;
+    RankTopology topo;
+    MPI_Comm node_comm, leader_comm;
+
+    MPI_Comm_rank(app_instance.mpi.world_comm, &world_rank);
+    MPI_Comm_size(app_instance.mpi.world_comm, &topo.number_of_ranks);
+    MPI_Comm_split_type(app_instance.mpi.world_comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &node_comm);
+    MPI_Comm_rank(node_comm, &topo.node_rank);
+    bool is_leader = topo.node_rank == 0;
+    MPI_Comm_split(app_instance.mpi.world_comm, is_leader ? 0 : MPI_UNDEFINED, world_rank, &leader_comm);
+    if (is_leader) MPI_Comm_size(leader_comm, &topo.number_of_nodes);
+    if (is_leader) MPI_Comm_rank(leader_comm, &topo.node_index);
+    MPI_Bcast(&topo.number_of_nodes, 1, MPI_INT, 0, app_instance.mpi.world_comm);
+    MPI_Bcast(&topo.node_index, 1, MPI_INT, 0, node_comm);
+
+    topo.ranks_per_node = topo.number_of_ranks/topo.number_of_nodes;
+
+    return topo;
+}
+
+
+/**
+ * @brief      Initializes the difference application instances.
+ */
+static void initialize(void)
+{
+    printf("\033[0;31m[PTBC] Number of chains = %d\033[0m\n", app_instance.ptbc.n_instances);
+
+    int flag;
+    MPI_Initialized(&flag);
+    err(!flag, "Initialize has to be called *after* MPI_Init().");
+
+    MPI_Comm_rank(app_instance.mpi.world_comm, &app_instance.mpi.world_rank);
+
+    // do nothing in case of a single chain
+    if (app_instance.ptbc.n_instances == 1) return;
+
+    app_instance.ptbc.active = true;
+    RankTopology topo = get_topology();
+
+    int instance_size = topo.number_of_ranks / app_instance.ptbc.n_instances;
+
+    err(topo.number_of_ranks % app_instance.ptbc.n_instances != 0,
+        "PTBC_NCHAINS = %d must divide total number of ranks = %d",
+        app_instance.ptbc.n_instances, topo.number_of_ranks);
+    err(instance_size % topo.ranks_per_node != 0 && topo.ranks_per_node % instance_size != 0,
+        "The number of processes per node = %d and instance_size = %d: one must be divisible by the other",
+        topo.ranks_per_node, instance_size);
+
+    // We perform a topology-aware splitting of processes into instances using
+    // MPI_Comm_split. Processes within the same node should preferably be
+    // associated to the same instance. We have 3 cases:
+    //
+    // Case 1: If we have one instance per node, there is nothing special to
+    // consider.
+    //
+    // Case 2: If instances span multiple nodes, they should span over the
+    // minimal number of nodes possible. Instances only cover whole nodes. We
+    // have no notion of nodes being "close" to each other, alhtough we group
+    // nodes together with adjacent node indices. Node indices are inherited
+    // from rank numbers. If ranks with adjacent indices are "close", then nodes
+    // with adjacent node indices are "close".
+    //
+    // Case 3: If we have multiple instances per node, ranks in the same
+    // instance should have adjacent world rank numbers, i.e. they should be
+    // "close" to each other cache-wise. No instance covers more than one node.
+    int color;
+    int key = app_instance.mpi.world_rank; // order of the ranks is kept
+    if (instance_size == topo.ranks_per_node) { // case 1: one instance per node
+        color = topo.node_index;
+    } else if (instance_size % topo.ranks_per_node == 0) { // case 2: one instance spans multiple nodes
+        int nodes_per_instance = instance_size / topo.ranks_per_node;
+        int remainder = topo.node_index % nodes_per_instance;
+        color = (topo.node_index - remainder) / nodes_per_instance; // group nodes with adjacent node indices
+    } else if (topo.ranks_per_node % instance_size == 0) { // case 3: multiple instances per node
+        int instances_per_node = topo.ranks_per_node / instance_size;
+        int remainder = topo.node_rank % instance_size;
+        int per_node_instance_index = (topo.node_rank - remainder) / instance_size;
+        color = topo.node_index*instances_per_node + per_node_instance_index;
+    }
+
+    MPI_Comm_split(app_instance.mpi.world_comm, color, key, &app_instance.mpi.comm);
+
+    int n;
+    MPI_Comm_size(app_instance.mpi.comm, &n);
+    err(instance_size != n, "Rank topology is not uniform");
+
+    int instance_rank;
+    MPI_Comm_rank(app_instance.mpi.comm, &instance_rank);
+    app_instance.ptbc.instance_id = color;
+
+    printf("\033[0;31m[PTBC] world rank = %d/%d in instance_id = %d/%d, as instance_rank = %d/%d\033[0m\n",
+        app_instance.mpi.world_rank, topo.number_of_ranks,
+        app_instance.ptbc.instance_id, app_instance.ptbc.n_instances,
+        instance_rank, instance_size);
+
+
+    //err(true, "bailing out");
+    /*if (app_instance.ptbc.instance_id != 0) {
+        char logfile[1024];
+        snprintf(logfile, 1024, "logfile_%.2d.log", app_instance.ptbc.instance_id);
+        freopen(logfile, "w", stdout);
+    }*/
+
+    // Every instance just changes into a subdirectory "instance_xx". Relative
+    // paths work, absolute paths not.
+    struct stat st = {0};
+    char subdir[1024];
+    snprintf(subdir, 1024, "instance_%.2d", app_instance.ptbc.instance_id);
+    if (stat(subdir, &st) == -1)
+        mkdir(subdir, 0700);
+
+    chdir(subdir);
+}
\ No newline at end of file
diff --git a/src/lib/compare_derivative.c b/src/lib/compare_derivative.c
index 6ac39f222..b7802481f 100644
--- a/src/lib/compare_derivative.c
+++ b/src/lib/compare_derivative.c
@@ -64,8 +64,8 @@ void compare_derivative(monomial *mnl, su3adj **ext_lib, su3adj **native, const
 
   int red_n_diff = 0;
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
-  MPI_Reduce(&n_diff, &red_n_diff, 1, MPI_INT, MPI_MAX, 0, MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
+  MPI_Reduce(&n_diff, &red_n_diff, 1, MPI_INT, MPI_MAX, 0, app()->mpi.comm);
 #else
   red_n_diff = n_diff;
 #endif
diff --git a/src/lib/fatal_error.c b/src/lib/fatal_error.c
index b4f5d5be9..8e00dfbe7 100644
--- a/src/lib/fatal_error.c
+++ b/src/lib/fatal_error.c
@@ -46,7 +46,7 @@ void fatal_error(char const *error, char const *function) {
   }
 
 #ifdef TM_USE_MPI
-  MPI_Abort(MPI_COMM_WORLD, 1);
+  MPI_Abort(app()->mpi.comm, 1);
   MPI_Finalize();
 #endif
 
diff --git a/src/lib/include/alloc.h b/src/lib/include/alloc.h
new file mode 100644
index 000000000..8071e9a4c
--- /dev/null
+++ b/src/lib/include/alloc.h
@@ -0,0 +1,49 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2026 Roman Gruber
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD 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 tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Allocation utils
+ *
+ * Author: Roman Gruber
+ *         roman.gruber@unibe.ch
+ *
+ *******************************************************************************/
+
+#ifndef ALLOC_H
+#define ALLOC_H
+
+
+#include <stdbool.h>
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+void *safe_malloc_impl(size_t size, const char *file, int line, const char *func);
+void *safe_calloc_impl(size_t size, const char *file, int line, const char *func);
+#define safe_malloc(size) safe_malloc_impl((size), __FILE__, __LINE__, __func__)
+#define safe_calloc(size) safe_calloc_impl((size), __FILE__, __LINE__, __func__)
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
\ No newline at end of file
diff --git a/src/lib/include/app.h b/src/lib/include/app.h
new file mode 100644
index 000000000..169f48f22
--- /dev/null
+++ b/src/lib/include/app.h
@@ -0,0 +1,108 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2026 Roman Gruber
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD 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 tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Simple MPI header wrapper
+ *
+ * Author: Roman Gruber
+ *         roman.gruber@unibe.ch
+ *
+ *******************************************************************************/
+
+#ifndef APP_H
+#define APP_H
+
+
+#include <stdbool.h>
+
+
+#if defined(TM_USE_MPI)
+#include <mpi.h>
+#endif
+
+
+#ifndef MPI_VERSION
+typedef int MPI_Comm;
+#endif
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+#define MAX_N_DEFECTS 10
+#define MAX_N_INSTANCES 10
+
+
+typedef enum direction_t {
+   DIRECTION_T = 0,
+   DIRECTION_X = 1,
+   DIRECTION_Y = 2,
+   DIRECTION_Z = 3
+} direction_t;
+
+
+typedef struct {
+    MPI_Comm comm;          // MPI instance communicator
+    MPI_Comm world_comm;    // MPI world communicator
+    int world_rank;         // MPI world rank
+} MPIContext;
+
+
+typedef struct {
+    bool active;        // Whether the defect is active or not
+    int Ld[3];          // Extents of the defect
+    direction_t along;  // Along which dimension
+} PTBCDefect;
+
+
+typedef struct {
+    bool active;          // Whether the instance is active or not
+    int n_coeffs;         // Number of coefficients / defect this instance is associated to
+    PTBCDefect** defects; // List of defects where this instance is associated to
+    double* coefficients; // List of coefficients for the defects
+} PTBCInstance;
+
+
+typedef struct {
+    bool active;                                // Whether PTBC mode is active or not
+    int instance_id;                            // Instance ID
+    int n_instances;                            // Number of instances
+    int n_defects;                              // Number of defects
+    PTBCInstance instances[MAX_N_INSTANCES];    // List of all instances
+    PTBCDefect defects[MAX_N_DEFECTS];          // List of all defects
+    void (*initialize)(void);                   // PTBC algorithm initializer
+} PTBCContext;
+
+
+typedef struct {
+    MPIContext mpi;     // MPI context
+    PTBCContext ptbc;   // PTBC context
+} AppContext;
+
+
+const AppContext* app(void);
+AppContext* appm(void);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
\ No newline at end of file
diff --git a/src/lib/include/mpi.h b/src/lib/include/mpi.h
new file mode 100644
index 000000000..7dbf88df4
--- /dev/null
+++ b/src/lib/include/mpi.h
@@ -0,0 +1,35 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2026 Roman Gruber
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD 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 tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Simple MPI header wrapper
+ *
+ * Author: Roman Gruber
+ *         roman.gruber@unibe.ch
+ *
+ *******************************************************************************/
+
+#ifndef MPI_WRAPPER_H
+#define MPI_WRAPPER_H
+
+
+#include_next <mpi.h> // include *real* MPI header
+#include "app.h"
+
+
+#endif
\ No newline at end of file
diff --git a/src/lib/init/init_parallel.c b/src/lib/init/init_parallel.c
index 194404d7d..9dfdbb0c5 100644
--- a/src/lib/init/init_parallel.c
+++ b/src/lib/init/init_parallel.c
@@ -71,7 +71,7 @@ void init_parallel_and_read_input(int argc, char *argv[], const char input_filen
 #endif  // QPHIX_QMP_COMMS
 
 #if defined(TM_USE_MPI) || defined(QPHIX_QMP_COMMS)
-  MPI_Comm_rank(MPI_COMM_WORLD, &g_proc_id);
+  MPI_Comm_rank(app()->mpi.comm, &g_proc_id);
 #else
   g_proc_id = 0;
 #endif
diff --git a/src/lib/io/deri_write_stdout.c b/src/lib/io/deri_write_stdout.c
index 6b095d22e..db1c4f98e 100644
--- a/src/lib/io/deri_write_stdout.c
+++ b/src/lib/io/deri_write_stdout.c
@@ -66,7 +66,7 @@ void deri_write_stdout(su3adj** const df) {
             }
           }
 #ifdef TM_USE_MPI
-          MPI_Barrier(MPI_COMM_WORLD);
+          MPI_Barrier(app()->mpi.comm);
 #endif
         }
       }
diff --git a/src/lib/io/dml.c b/src/lib/io/dml.c
index 2650e108d..0465f1206 100644
--- a/src/lib/io/dml.c
+++ b/src/lib/io/dml.c
@@ -31,7 +31,7 @@ int DML_global_xor(uint32_t *x) {
   int status;
 
   status =
-      MPI_Allreduce((void *)&work, (void *)&dest, 1, MPI_UNSIGNED_LONG, MPI_BXOR, MPI_COMM_WORLD);
+      MPI_Allreduce((void *)&work, (void *)&dest, 1, MPI_UNSIGNED_LONG, MPI_BXOR, app()->mpi.comm);
 
   if (status == MPI_SUCCESS) {
     *x = (uint32_t)dest;
diff --git a/src/lib/io/eospinor_read.c b/src/lib/io/eospinor_read.c
index 2767d8b13..d23623595 100644
--- a/src/lib/io/eospinor_read.c
+++ b/src/lib/io/eospinor_read.c
@@ -99,7 +99,7 @@ int read_eospinor(spinor *const s, char *filename) {
                   "LIME read error occured with status = %d while reading file %s!\n Aborting...\n",
                   status, filename);
 #ifdef TM_USE_MPI
-              MPI_Abort(MPI_COMM_WORLD, 1);
+              MPI_Abort(app()->mpi.comm, 1);
               MPI_Finalize();
 #endif
               exit(500);
diff --git a/src/lib/io/eospinor_write.c b/src/lib/io/eospinor_write.c
index 619e11ac8..52c3e688d 100644
--- a/src/lib/io/eospinor_write.c
+++ b/src/lib/io/eospinor_write.c
@@ -67,7 +67,7 @@ int write_eospinor(spinor* const s, char* filename, const double evalue, const d
     if (limewriter == (LimeWriter*)NULL) {
       fprintf(stderr, "LIME error in file %s for writing!\n Aboring...\n", filename);
 #ifdef TM_USE_MPI
-      MPI_Abort(MPI_COMM_WORLD, 1);
+      MPI_Abort(app()->mpi.comm, 1);
       MPI_Finalize();
 #endif
       exit(500);
@@ -78,7 +78,7 @@ int write_eospinor(spinor* const s, char* filename, const double evalue, const d
     if (status < 0) {
       fprintf(stderr, "LIME write header (xlf-info) error %d\n", status);
 #ifdef TM_USE_MPI
-      MPI_Abort(MPI_COMM_WORLD, 1);
+      MPI_Abort(app()->mpi.comm, 1);
       MPI_Finalize();
 #endif
       exit(500);
@@ -94,7 +94,7 @@ int write_eospinor(spinor* const s, char* filename, const double evalue, const d
     if (status < 0) {
       fprintf(stderr, "LIME write header (eospinor-binary-data) error %d\n", status);
 #ifdef TM_USE_MPI
-      MPI_Abort(MPI_COMM_WORLD, 1);
+      MPI_Abort(app()->mpi.comm, 1);
       MPI_Finalize();
 #endif
       exit(500);
@@ -137,7 +137,7 @@ int write_eospinor(spinor* const s, char* filename, const double evalue, const d
               if (status < 0) {
                 fprintf(stderr, "LIME write error %d\n", status);
 #ifdef TM_USE_MPI
-                MPI_Abort(MPI_COMM_WORLD, 1);
+                MPI_Abort(app()->mpi.comm, 1);
                 MPI_Finalize();
 #endif
                 exit(500);
diff --git a/src/lib/io/gauge_read_binary.c b/src/lib/io/gauge_read_binary.c
index 473e4d9c7..28cebd7f7 100644
--- a/src/lib/io/gauge_read_binary.c
+++ b/src/lib/io/gauge_read_binary.c
@@ -219,7 +219,7 @@ int read_binary_gauge_data(LimeReader *limereader, DML_Checksum *checksum, param
                 "LIME read error occurred with status = %d while reading in gauge_read_binary.c!\n",
                 status);
 #ifdef TM_USE_MPI
-            MPI_Abort(MPI_COMM_WORLD, 1);
+            MPI_Abort(app()->mpi.comm, 1);
             MPI_Finalize();
 #endif
             return (-2);
diff --git a/src/lib/io/gauge_write.c b/src/lib/io/gauge_write.c
index 43ff7fd83..20605bc8a 100644
--- a/src/lib/io/gauge_write.c
+++ b/src/lib/io/gauge_write.c
@@ -50,7 +50,7 @@ int write_gauge_field(char *filename, const int prec, paramsXlfInfo const *xlfIn
     fflush(stdout);
   }
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
 #endif /* MPI */
 
   destruct_writer(writer);
diff --git a/src/lib/io/gauge_write_binary.c b/src/lib/io/gauge_write_binary.c
index ad3c7882e..64f504264 100644
--- a/src/lib/io/gauge_write_binary.c
+++ b/src/lib/io/gauge_write_binary.c
@@ -222,7 +222,7 @@ int write_binary_gauge_data(LimeWriter* limewriter, const int prec, DML_Checksum
               fprintf(stderr, "id = %d, bytes = %lu, size = %d\n", g_cart_id, bytes,
                       (int)(4 * sizeof(su3) / 8));
 #ifdef TM_USE_MPI
-              MPI_Abort(MPI_COMM_WORLD, 1);
+              MPI_Abort(app()->mpi.comm, 1);
               MPI_Finalize();
 #endif
               exit(500);
diff --git a/src/lib/io/io_cm.c b/src/lib/io/io_cm.c
index e08fd37c6..2f4bfa5a4 100644
--- a/src/lib/io/io_cm.c
+++ b/src/lib/io/io_cm.c
@@ -104,7 +104,7 @@ int read_spinorfield_cm_swap_single(spinor *const s, spinor *const r, char *file
   if (ifs == (FILE *)NULL) {
     fprintf(stderr, "Could not open file %s\n Aborting...\n", filename);
 #ifdef TM_USE_MPI
-    MPI_Abort(MPI_COMM_WORLD, 1);
+    MPI_Abort(app()->mpi.comm, 1);
     MPI_Finalize();
 #endif
     exit(500);
diff --git a/src/lib/io/spinor_read_binary.c b/src/lib/io/spinor_read_binary.c
index 81607a700..fbf97c32c 100644
--- a/src/lib/io/spinor_read_binary.c
+++ b/src/lib/io/spinor_read_binary.c
@@ -198,7 +198,7 @@ int read_binary_spinor_data(spinor *const s, spinor *const r, LimeReader *limere
                     "spinor_read_binary.c!\n",
                     status);
 #ifdef TM_USE_MPI
-            MPI_Abort(MPI_COMM_WORLD, 1);
+            MPI_Abort(app()->mpi.comm, 1);
             MPI_Finalize();
 #endif
             return (-2);
@@ -376,7 +376,7 @@ int read_binary_spinor_data_l(spinor *const s, LimeReader *limereader, DML_Check
                     "spinor_read_binary.c!\n",
                     status);
 #ifdef TM_USE_MPI
-            MPI_Abort(MPI_COMM_WORLD, 1);
+            MPI_Abort(app()->mpi.comm, 1);
             MPI_Finalize();
 #endif
             return (-2);
diff --git a/src/lib/io/spinor_write_binary.c b/src/lib/io/spinor_write_binary.c
index 560b5ce65..b8718cc9e 100644
--- a/src/lib/io/spinor_write_binary.c
+++ b/src/lib/io/spinor_write_binary.c
@@ -217,7 +217,7 @@ int write_binary_spinor_data(spinor *const s, spinor *const r, LimeWriter *limew
                       "write_binary_spinor_data (spinor_write_binary.c)!\n",
                       status);
 #ifdef TM_USE_MPI
-              MPI_Abort(MPI_COMM_WORLD, 1);
+              MPI_Abort(app()->mpi.comm, 1);
               MPI_Finalize();
 #endif
               exit(500);
@@ -458,7 +458,7 @@ int write_binary_spinor_data_l(spinor *const s, LimeWriter *limewriter, DML_Chec
                       "write_binary_spinor_data_l (spinor_write_binary.c)!\n",
                       status);
 #ifdef TM_USE_MPI
-              MPI_Abort(MPI_COMM_WORLD, 1);
+              MPI_Abort(app()->mpi.comm, 1);
               MPI_Finalize();
 #endif
               exit(500);
diff --git a/src/lib/io/spinor_write_stdout.c b/src/lib/io/spinor_write_stdout.c
index 9c13abf47..6c1fa1034 100644
--- a/src/lib/io/spinor_write_stdout.c
+++ b/src/lib/io/spinor_write_stdout.c
@@ -57,7 +57,7 @@ void spinor_write_stdout(spinor* const s) {
             fflush(stdout);
           }
 #ifdef TM_USE_MPI
-          MPI_Barrier(MPI_COMM_WORLD);
+          MPI_Barrier(app()->mpi.comm);
 #endif
         }
       }
diff --git a/src/lib/io/sw_write_stdout.c b/src/lib/io/sw_write_stdout.c
index 027edd172..b2cec0b29 100644
--- a/src/lib/io/sw_write_stdout.c
+++ b/src/lib/io/sw_write_stdout.c
@@ -66,7 +66,7 @@ void sw_write_stdout(su3** u) {
             }
           }
 #ifdef TM_USE_MPI
-          MPI_Barrier(MPI_COMM_WORLD);
+          MPI_Barrier(app()->mpi.comm);
 #endif
         }
       }
diff --git a/src/lib/io/utils_kill_with_error.c b/src/lib/io/utils_kill_with_error.c
index 322536bd7..282e5bc4b 100644
--- a/src/lib/io/utils_kill_with_error.c
+++ b/src/lib/io/utils_kill_with_error.c
@@ -14,7 +14,7 @@ void kill_with_error(LIME_FILE *fh, int const rank, char const *error) {
 #endif /* TM_USE_LEMON */
 
 #ifdef TM_USE_MPI
-  MPI_Abort(MPI_COMM_WORLD, 1);
+  MPI_Abort(app()->mpi.comm, 1);
   MPI_Finalize();
 #endif
   exit(500);
diff --git a/src/lib/linalg/assign_mul_add_r_and_square.c b/src/lib/linalg/assign_mul_add_r_and_square.c
index f9fb60f72..fc45c313f 100644
--- a/src/lib/linalg/assign_mul_add_r_and_square.c
+++ b/src/lib/linalg/assign_mul_add_r_and_square.c
@@ -99,7 +99,7 @@ double assign_mul_add_r_and_square(spinor *const R, const double c, const spinor
 
 #ifdef TM_USE_MPI
   if (parallel) {
-    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
     return (mres);
   }
 #endif
diff --git a/src/lib/linalg/diff_and_square_norm.c b/src/lib/linalg/diff_and_square_norm.c
index 01a134fc8..0b7d6612a 100644
--- a/src/lib/linalg/diff_and_square_norm.c
+++ b/src/lib/linalg/diff_and_square_norm.c
@@ -72,7 +72,7 @@ double diff_and_square_norm(spinor *const Q, spinor *const R, const int N) {
   }
   kc = ks + kc;
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&kc, &ks, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&kc, &ks, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   return ks;
 #else
   return kc;
diff --git a/src/lib/linalg/scalar_prod_body.inc b/src/lib/linalg/scalar_prod_body.inc
index 9a38cf5fb..3114d6119 100644
--- a/src/lib/linalg/scalar_prod_body.inc
+++ b/src/lib/linalg/scalar_prod_body.inc
@@ -71,7 +71,7 @@ _Complex double _PSWITCH(scalar_prod)(const _PTSWITCH(spinor) *const S,
 
 #ifdef TM_USE_MPI
   if (parallel == 1) {
-    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, app()->mpi.comm);
     return (mres);
   }
 #endif
@@ -113,7 +113,7 @@ _Complex double _PSWITCH(scalar_prod_ts)(const _PTSWITCH(spinor) *const S,
 
 #ifdef TM_USE_MPI
   if (parallel == 1) {
-    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, app()->mpi.comm);
     return (mres);
   }
 #endif
diff --git a/src/lib/linalg/scalar_prod_i.c b/src/lib/linalg/scalar_prod_i.c
index ca59c65d3..ca1d378b7 100644
--- a/src/lib/linalg/scalar_prod_i.c
+++ b/src/lib/linalg/scalar_prod_i.c
@@ -64,7 +64,7 @@ double scalar_prod_i(spinor *const S, spinor *const R, const int N, const int pa
 
 #if defined TM_USE_MPI
   if (parallel == 1) {
-    MPI_Allreduce(&kc, &ks, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&kc, &ks, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
     kc = ks;
   }
 #endif
diff --git a/src/lib/linalg/scalar_prod_r.c b/src/lib/linalg/scalar_prod_r.c
index c5288aa34..8e1c07af1 100644
--- a/src/lib/linalg/scalar_prod_r.c
+++ b/src/lib/linalg/scalar_prod_r.c
@@ -91,7 +91,7 @@ double scalar_prod_r(const spinor *const S, const spinor *const R, const int N,
 
 #if defined TM_USE_MPI
   if (parallel) {
-    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
     return mres;
   }
 #endif
diff --git a/src/lib/linalg/scalar_prod_r_32.c b/src/lib/linalg/scalar_prod_r_32.c
index 5bc512806..23f673ee2 100644
--- a/src/lib/linalg/scalar_prod_r_32.c
+++ b/src/lib/linalg/scalar_prod_r_32.c
@@ -64,7 +64,7 @@ float scalar_prod_r_32(const spinor32 *const S, const spinor32 *const R, const i
 
 #if defined TM_USE_MPI
   if (parallel) {
-    MPI_Allreduce(&res, &mres, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_FLOAT, MPI_SUM, app()->mpi.comm);
     return mres;
   }
 #endif
diff --git a/src/lib/linalg/square_and_minmax.c b/src/lib/linalg/square_and_minmax.c
index 3a97c1ede..6432d2602 100644
--- a/src/lib/linalg/square_and_minmax.c
+++ b/src/lib/linalg/square_and_minmax.c
@@ -75,12 +75,12 @@ void square_and_minmax(double *const sum, double *const min, double *const max,
 
 #if defined TM_USE_MPI
 
-  MPI_Allreduce(&kc, sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&kc, sum, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 
-  MPI_Allreduce(min, &kc, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+  MPI_Allreduce(min, &kc, 1, MPI_DOUBLE, MPI_MIN, app()->mpi.comm);
   *min = kc;
 
-  MPI_Allreduce(max, &kc, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+  MPI_Allreduce(max, &kc, 1, MPI_DOUBLE, MPI_MAX, app()->mpi.comm);
   *max = kc;
 
 #endif
@@ -129,12 +129,12 @@ void square_and_minmax_rel(double *const sum, double *const min, double *const m
 
 #if defined TM_USE_MPI
 
-  MPI_Allreduce(&kc, sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&kc, sum, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 
-  MPI_Allreduce(min, &kc, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+  MPI_Allreduce(min, &kc, 1, MPI_DOUBLE, MPI_MIN, app()->mpi.comm);
   *min = kc;
 
-  MPI_Allreduce(max, &kc, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+  MPI_Allreduce(max, &kc, 1, MPI_DOUBLE, MPI_MAX, app()->mpi.comm);
   *max = kc;
 
 #endif
@@ -256,18 +256,18 @@ void square_and_minmax_abs(double *const sum, double *const min, double *const m
 
 #if defined TM_USE_MPI
 
-  MPI_Allreduce(&kc, sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&kc, sum, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 
-  MPI_Allreduce(min, &kc, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+  MPI_Allreduce(min, &kc, 1, MPI_DOUBLE, MPI_MIN, app()->mpi.comm);
   *min = kc;
 
-  MPI_Allreduce(max, &kc, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+  MPI_Allreduce(max, &kc, 1, MPI_DOUBLE, MPI_MAX, app()->mpi.comm);
   *max = kc;
 
-  MPI_Allreduce(min_abs, &kc, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+  MPI_Allreduce(min_abs, &kc, 1, MPI_DOUBLE, MPI_MIN, app()->mpi.comm);
   *min_abs = kc;
 
-  MPI_Allreduce(max_abs, &kc, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+  MPI_Allreduce(max_abs, &kc, 1, MPI_DOUBLE, MPI_MAX, app()->mpi.comm);
   *max_abs = kc;
 
 #endif
@@ -428,18 +428,18 @@ void square_and_minmax_rel_abs(double *const sum, double *const min, double *con
 
 #if defined TM_USE_MPI
 
-  MPI_Allreduce(&kc, sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&kc, sum, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 
-  MPI_Allreduce(min, &kc, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+  MPI_Allreduce(min, &kc, 1, MPI_DOUBLE, MPI_MIN, app()->mpi.comm);
   *min = kc;
 
-  MPI_Allreduce(max, &kc, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+  MPI_Allreduce(max, &kc, 1, MPI_DOUBLE, MPI_MAX, app()->mpi.comm);
   *max = kc;
 
-  MPI_Allreduce(min_abs, &kc, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+  MPI_Allreduce(min_abs, &kc, 1, MPI_DOUBLE, MPI_MIN, app()->mpi.comm);
   *min_abs = kc;
 
-  MPI_Allreduce(max_abs, &kc, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+  MPI_Allreduce(max_abs, &kc, 1, MPI_DOUBLE, MPI_MAX, app()->mpi.comm);
   *max_abs = kc;
 
 #endif
diff --git a/src/lib/linalg/square_and_prod_r.c b/src/lib/linalg/square_and_prod_r.c
index 212b46593..b02babf42 100644
--- a/src/lib/linalg/square_and_prod_r.c
+++ b/src/lib/linalg/square_and_prod_r.c
@@ -80,7 +80,7 @@ void square_and_prod_r(double *const x1, double *const x2, spinor *const S, spin
 
 #if defined TM_USE_MPI
 
-  MPI_Allreduce(&xkc, x1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&xkc, x1, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 
 #endif
   kc = ks + kc;
@@ -88,7 +88,7 @@ void square_and_prod_r(double *const x1, double *const x2, spinor *const S, spin
 
 #if defined TM_USE_MPI
 
-  MPI_Allreduce(&kc, x2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&kc, x2, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 
 #endif
 }
diff --git a/src/lib/linalg/square_norm.c b/src/lib/linalg/square_norm.c
index fd4b77781..7c871f41b 100644
--- a/src/lib/linalg/square_norm.c
+++ b/src/lib/linalg/square_norm.c
@@ -91,7 +91,7 @@ double square_norm(const spinor *const P, const int N, const int parallel) {
 
 #ifdef TM_USE_MPI
   if (parallel) {
-    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
     return mres;
   }
 #endif
@@ -142,7 +142,7 @@ double square_norm_ts(const spinor *const P, const int N, const int parallel) {
 
 #ifdef TM_USE_MPI
   if (parallel) {
-    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
     return mres;
   }
 #endif
diff --git a/src/lib/linalg/square_norm_32.c b/src/lib/linalg/square_norm_32.c
index b331620cd..53207f454 100644
--- a/src/lib/linalg/square_norm_32.c
+++ b/src/lib/linalg/square_norm_32.c
@@ -65,7 +65,7 @@ float square_norm_32(const spinor32 *const P, const int N, const int parallel) {
 
 #ifdef TM_USE_MPI
   if (parallel) {
-    MPI_Allreduce(&res, &mres, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_FLOAT, MPI_SUM, app()->mpi.comm);
     return mres;
   }
 #endif
@@ -113,7 +113,7 @@ float square_norm_ts_32(const spinor32 *const P, const int N, const int parallel
 
 #ifdef TM_USE_MPI
   if (parallel) {
-    MPI_Allreduce(&res, &mres, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res, &mres, 1, MPI_FLOAT, MPI_SUM, app()->mpi.comm);
     return mres;
   }
 #endif
diff --git a/src/lib/meas/correlators.c b/src/lib/meas/correlators.c
index 36c8d2776..ed50eba3f 100644
--- a/src/lib/meas/correlators.c
+++ b/src/lib/meas/correlators.c
@@ -125,7 +125,7 @@ void correlators_measurement(const int traj, const int id, const int ieo) {
         t0 = (int)(measurement_list[id].max_source_slice * tmp);
       }
 #ifdef TM_USE_MPI
-      MPI_Bcast(&t0, 1, MPI_INT, 0, MPI_COMM_WORLD);
+      MPI_Bcast(&t0, 1, MPI_INT, 0, app()->mpi.comm);
 #endif
       if (g_debug_level > 1 && g_proc_id == 0) {
         printf("# timeslice set to %d (T=%d) for online measurement\n", t0, g_nproc_t * T);
@@ -252,7 +252,7 @@ void correlators_measurement(const int traj, const int id, const int ieo) {
       free(sCpp);
       free(sCpa);
       free(sCp4);
-      MPI_Barrier(MPI_COMM_WORLD);
+      MPI_Barrier(app()->mpi.comm);
 #else
       free(Cpp);
       free(Cpa);
diff --git a/src/lib/meas/measure_clover_field_strength_observables.c b/src/lib/meas/measure_clover_field_strength_observables.c
index 3d33ae999..b95a2390a 100644
--- a/src/lib/meas/measure_clover_field_strength_observables.c
+++ b/src/lib/meas/measure_clover_field_strength_observables.c
@@ -210,9 +210,9 @@ void measure_clover_field_strength_observables(const su3 **const gf,
 #endif
 
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&Eres, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&Eres, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   Eres = mres;
-  MPI_Allreduce(&Qres, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&Qres, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   Qres = mres;
 #endif
   fso->E = energy_density_normalization * Eres;
diff --git a/src/lib/meas/oriented_plaquettes.c b/src/lib/meas/oriented_plaquettes.c
index 51e73acf2..1fe9207fd 100644
--- a/src/lib/meas/oriented_plaquettes.c
+++ b/src/lib/meas/oriented_plaquettes.c
@@ -80,7 +80,7 @@ void measure_oriented_plaquettes(const su3 **const gf, double *plaq) {
   }
 
 #ifdef TM_USE_MPI
-  MPI_Allreduce(plaq, mplaq, 6, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(plaq, mplaq, 6, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   for (int j = 0; j < 6; j++) plaq[j] = mplaq[j];
 #endif
   tm_stopwatch_pop(&g_timers, 0, 2, "");
diff --git a/src/lib/meas/pion_norm.c b/src/lib/meas/pion_norm.c
index 1ee756145..38d13320b 100644
--- a/src/lib/meas/pion_norm.c
+++ b/src/lib/meas/pion_norm.c
@@ -68,7 +68,7 @@ void pion_norm_measurement(const int traj, const int id, const int ieo) {
   ranlxs(&tmp, 1);
   z0 = (int)(measurement_list[id].max_source_slice * tmp);
 #ifdef TM_USE_MPI
-  MPI_Bcast(&z0, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&z0, 1, MPI_INT, 0, app()->mpi.comm);
 #endif
 
   Cpp = (double *)calloc(g_nproc_z * LZ, sizeof(double));
diff --git a/src/lib/meas/polyakov_loop.c b/src/lib/meas/polyakov_loop.c
index 25deea402..14d7c9ba5 100644
--- a/src/lib/meas/polyakov_loop.c
+++ b/src/lib/meas/polyakov_loop.c
@@ -66,7 +66,7 @@ void polyakov_loop(_Complex double *pl_, const int mu) {
     fprintf(stderr, "Only direction %d and %d are allowed.\n", 2, 3);
     fprintf(stderr, "Actual value is %d! Aborting...\n", mu);
 #ifdef TM_USE_MPI
-    MPI_Abort(MPI_COMM_WORLD, 10);
+    MPI_Abort(app()->mpi.comm, 10);
     MPI_Finalize();
 #endif
     exit(0);
@@ -137,7 +137,7 @@ void polyakov_loop(_Complex double *pl_, const int mu) {
 
   /* Collect the results and return:*/
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&pl, &pls, 2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&pl, &pls, 2, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   pl = pls;
 #endif
 
diff --git a/src/lib/measure_gauge_action.c b/src/lib/measure_gauge_action.c
index ecbe7a888..654990166 100644
--- a/src/lib/measure_gauge_action.c
+++ b/src/lib/measure_gauge_action.c
@@ -96,7 +96,7 @@ double measure_plaquette(const su3 *const *const gf) {
   for (int i = 0; i < omp_num_threads; ++i) res += g_omp_acc_re[i];
 #endif
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   res = mres;
 #endif
   return res;
@@ -178,7 +178,7 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
   for (int i = 0; i < omp_num_threads; ++i) res += g_omp_acc_re[i];
 #endif
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   res = mres;
 #endif
   GaugeInfo.plaquetteEnergy = res;
diff --git a/src/lib/measure_rectangles.c b/src/lib/measure_rectangles.c
index e6f59b8b2..4f0ca07fd 100644
--- a/src/lib/measure_rectangles.c
+++ b/src/lib/measure_rectangles.c
@@ -129,7 +129,7 @@ double measure_rectangles(const su3 **const gf) {
 #else
 #endif
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   res = mres;
 #endif
 
diff --git a/src/lib/monomial/moment_energy.c b/src/lib/monomial/moment_energy.c
index 01776edfa..c9dfe16e6 100644
--- a/src/lib/monomial/moment_energy.c
+++ b/src/lib/monomial/moment_energy.c
@@ -69,7 +69,7 @@ double moment_energy(su3adj **const momenta) {
   kc = 0.5 * (ks + kc);
 #ifdef TM_USE_MPI
   ks = kc;
-  MPI_Allreduce(&ks, &kc, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&ks, &kc, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
 #endif
   etime = gettime();
   if (g_proc_id == 0) {
diff --git a/src/lib/monomial/monitor_forces.c b/src/lib/monomial/monitor_forces.c
index dc186ac05..ea64ddcfa 100644
--- a/src/lib/monomial/monitor_forces.c
+++ b/src/lib/monomial/monitor_forces.c
@@ -89,9 +89,9 @@ void monitor_forces(hamiltonian_field_t* const hf) {
 
       // output for force monitoring
 #ifdef TM_USE_MPI
-      MPI_Reduce(&sum, &sum2, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+      MPI_Reduce(&sum, &sum2, 1, MPI_DOUBLE, MPI_SUM, 0, app()->mpi.comm);
       sum = sum2;
-      MPI_Reduce(&max, &sum2, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
+      MPI_Reduce(&max, &sum2, 1, MPI_DOUBLE, MPI_MAX, 0, app()->mpi.comm);
       max = sum2;
 #endif
       if (g_proc_id == 0) {
diff --git a/src/lib/mpi_init.c b/src/lib/mpi_init.c
index f245f0556..25c79809a 100644
--- a/src/lib/mpi_init.c
+++ b/src/lib/mpi_init.c
@@ -159,7 +159,7 @@ void reduce_su3_ray(void *u_i /* in */, void *u_io /* in/out */, int *len /* in
 
   if (*dt != mpi_su3) {
     fprintf(stderr, "\nInvalid datatype for reduce_su3_ray(); abort.\n");
-    MPI_Abort(MPI_COMM_WORLD, 1);
+    MPI_Abort(app()->mpi.comm, 1);
   }
   for (n = 0; n < *len; n++) {
     _su3_times_su3(tmp, *(u + n), *(v + n)) _su3_assign(*(v + n), tmp)
@@ -189,7 +189,7 @@ void tmlqcd_mpi_init(int argc, char *argv[]) {
 
 #ifdef TM_USE_MPI
 #ifdef TM_USE_SHMEM
-  /* we need that the PE number in MPI_COMM_WORL  */
+  /* we need that the PE number in MPI_COMM_WORLD  */
   /* exactly correspond to the one in g_cart_grid */
   reorder = 0;
 #endif
@@ -249,8 +249,8 @@ void tmlqcd_mpi_init(int argc, char *argv[]) {
   dims[2] = N_PROC_Y;
   dims[3] = N_PROC_Z;
 
-  MPI_Comm_size(MPI_COMM_WORLD, &g_nproc);
-  MPI_Comm_rank(MPI_COMM_WORLD, &g_proc_id);
+  MPI_Comm_size(app()->mpi.comm, &g_nproc);
+  MPI_Comm_rank(app()->mpi.comm, &g_proc_id);
   MPI_Get_processor_name(processor_name, &namelen);
   MPI_Dims_create(g_nproc, nalldims, dims);
   if (g_proc_id == 0) {
@@ -273,7 +273,7 @@ void tmlqcd_mpi_init(int argc, char *argv[]) {
       fprintf(stderr, "Please check your number of processors and the Nr?Procs input variables\n");
       fprintf(stderr, "Aborting...!\n");
     }
-    MPI_Abort(MPI_COMM_WORLD, 1);
+    MPI_Abort(app()->mpi.comm, 1);
     MPI_Finalize();
     exit(-1);
   }
@@ -332,7 +332,7 @@ void tmlqcd_mpi_init(int argc, char *argv[]) {
   halffield_buffer_z2 = (halfspinor *)malloc(T * LX * LY / 2 * sizeof(halfspinor));
 #endif
 
-  MPI_Cart_create(MPI_COMM_WORLD, nalldims, dims, periods, reorder, &g_cart_grid);
+  MPI_Cart_create(app()->mpi.comm, nalldims, dims, periods, reorder, &g_cart_grid);
   MPI_Comm_rank(g_cart_grid, &g_cart_id);
   MPI_Cart_coords(g_cart_grid, g_cart_id, nalldims, g_proc_coords);
   if (g_debug_level > 1) {
diff --git a/src/lib/operator/clover_det.c b/src/lib/operator/clover_det.c
index 65545d547..3bc057997 100644
--- a/src/lib/operator/clover_det.c
+++ b/src/lib/operator/clover_det.c
@@ -171,7 +171,7 @@ double sw_trace(const int ieo, const double mu, const int single_flavor) {
 #endif
 
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   tm_stopwatch_pop(&g_timers, 0, 1, "");
   return (fac * mres);
 #else
@@ -261,7 +261,7 @@ double sw_trace_nd(const int ieo, const double mu, const double eps) {
 #endif
 
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&res, &mres, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   tm_stopwatch_pop(&g_timers, 0, 1, "");
   return (mres);
 #else
diff --git a/src/lib/prepare_source.c b/src/lib/prepare_source.c
index 3b45c35a1..8f15c9c7b 100644
--- a/src/lib/prepare_source.c
+++ b/src/lib/prepare_source.c
@@ -91,7 +91,7 @@ void prepare_source(const int nstore, const int isample, const int ix, const int
               }
             }
 #ifdef TM_USE_MPI
-            MPI_Bcast(&t, 1, MPI_INT, 0, MPI_COMM_WORLD);
+            MPI_Bcast(&t, 1, MPI_INT, 0, app()->mpi.comm);
 #endif
             SourceInfo.t = t;
           }
@@ -180,7 +180,7 @@ void prepare_source(const int nstore, const int isample, const int ix, const int
           t = (int)(u * g_nproc_t * T);
         }
 #ifdef TM_USE_MPI
-        MPI_Bcast(&t, 1, MPI_INT, 0, MPI_COMM_WORLD);
+        MPI_Bcast(&t, 1, MPI_INT, 0, app()->mpi.comm);
 #endif
         SourceInfo.t = t;
       }
@@ -206,7 +206,7 @@ void prepare_source(const int nstore, const int isample, const int ix, const int
           }
         }
 #ifdef TM_USE_MPI
-        MPI_Bcast(&t, 1, MPI_INT, 0, MPI_COMM_WORLD);
+        MPI_Bcast(&t, 1, MPI_INT, 0, app()->mpi.comm);
 #endif
         SourceInfo.t = t;
       }
@@ -315,7 +315,7 @@ void prepare_source(const int nstore, const int isample, const int ix, const int
         if (read_spinor(g_spinor_field[2], g_spinor_field[3], source_filename, 0) != 0) {
           fprintf(stderr, "Error reading source! Aborting...\n");
 #ifdef TM_USE_MPI
-          MPI_Abort(MPI_COMM_WORLD, 1);
+          MPI_Abort(app()->mpi.comm, 1);
           MPI_Finalize();
 #endif
           exit(-1);
diff --git a/src/lib/qphix/qphix_interface.cpp b/src/lib/qphix/qphix_interface.cpp
index 2c61427dd..4469d3109 100644
--- a/src/lib/qphix/qphix_interface.cpp
+++ b/src/lib/qphix/qphix_interface.cpp
@@ -224,9 +224,9 @@ void _initQphix(int argc, char **argv, tm_QPhiXParams_t params, int c12, QphixPr
                  g_proc_coords[2], g_proc_coords[3], g_proc_coords[0]);
           free(qmp_coords);
           fflush(stdout);
-          MPI_Barrier(MPI_COMM_WORLD);
+          MPI_Barrier(app()->mpi.comm);
         } else {
-          MPI_Barrier(MPI_COMM_WORLD);
+          MPI_Barrier(app()->mpi.comm);
         }
       }
     }
diff --git a/src/lib/quda_interface.c b/src/lib/quda_interface.c
index fea10947a..e86a1aba3 100644
--- a/src/lib/quda_interface.c
+++ b/src/lib/quda_interface.c
@@ -2929,20 +2929,20 @@ void quda_mg_tune_params(void *spinorOut, void *spinorIn, const int max_iter) {
     copy_quda_mg_tunable_params(&tunable_params[0], &cur_params);
     print_tunable_params_pair(&cur_params, &tunable_params[0], mg_n_level);
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     tm_stopwatch_push(&g_timers, "updateMultigridQuda", "");
     updateMultigridQuda(quda_mg_preconditioner, &quda_mg_param);
     tm_stopwatch_pop(&g_timers, 0, 1, "TM_QUDA");
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
   }
 
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
   tm_stopwatch_push(&g_timers, "invertQuda", "");
   invertQuda(spinorOut, spinorIn, &inv_param);
   tunable_params[0].tts = inv_param.secs;
   tunable_params[0].iter = inv_param.iter;
   tm_stopwatch_pop(&g_timers, 0, 1, "TM_QUDA");
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
 
   for (i = 1; i < quda_mg_tuning_plan.mg_tuning_iterations; i++) {
     // the best params from all previous iterations
@@ -2987,16 +2987,16 @@ void quda_mg_tune_params(void *spinorOut, void *spinorIn, const int max_iter) {
 
     print_tunable_params_pair(&cur_params, &tunable_params[i], mg_n_level);
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     tm_stopwatch_push(&g_timers, "updateMultigridQuda", "");
     updateMultigridQuda(quda_mg_preconditioner, &quda_mg_param);
     tm_stopwatch_pop(&g_timers, 0, 1, "TM_QUDA");
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
     tm_stopwatch_push(&g_timers, "invertQuda", "");
     invertQuda(spinorOut, spinorIn, &inv_param);
     tm_stopwatch_pop(&g_timers, 0, 1, "TM_QUDA");
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
     tunable_params[i].tts = inv_param.secs;
     tunable_params[i].iter = inv_param.iter;
diff --git a/src/lib/reweighting_factor.c b/src/lib/reweighting_factor.c
index 417a23f9d..553116c15 100644
--- a/src/lib/reweighting_factor.c
+++ b/src/lib/reweighting_factor.c
@@ -163,6 +163,6 @@ void reweighting_factor(const int N, const int nstore) {
   free(data);
   free(trlog);
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
 #endif  // TM_USE_MPI
 }
diff --git a/src/lib/sighandler.c b/src/lib/sighandler.c
index 7bcb36d15..1447a678e 100644
--- a/src/lib/sighandler.c
+++ b/src/lib/sighandler.c
@@ -49,7 +49,7 @@ void catch_ill_inst(int s) {
   fprintf(stderr, "Aborting...\n");
   fflush(stdout);
 #ifdef TM_USE_MPI
-  MPI_Abort(MPI_COMM_WORLD, 1);
+  MPI_Abort(app()->mpi.comm, 1);
   MPI_Finalize();
 #endif
   exit(0);
diff --git a/src/lib/solver/dfl_projector.c b/src/lib/solver/dfl_projector.c
index 5f0ce2026..e165dfc15 100644
--- a/src/lib/solver/dfl_projector.c
+++ b/src/lib/solver/dfl_projector.c
@@ -428,7 +428,7 @@ void little_project(_Complex double *const out, _Complex double *const in, const
   }
 
 #ifdef TM_USE_MPI
-  MPI_Allreduce(phi, psi, N, MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(phi, psi, N, MPI_DOUBLE_COMPLEX, MPI_SUM, app()->mpi.comm);
 #else
   memcpy(psi, phi, N * sizeof(_Complex double));
 #endif
@@ -487,7 +487,7 @@ void little_project2(_Complex double *const out, _Complex double *const in, cons
     phi[i] = lscalar_prod(little_dfl_fields[i], in, nb_blocks * N, 0);
   }
 #ifdef TM_USE_MPI
-  MPI_Allreduce(phi, psi, g_N_s, MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(phi, psi, g_N_s, MPI_DOUBLE_COMPLEX, MPI_SUM, app()->mpi.comm);
 #else
   memcpy(psi, phi, g_N_s * sizeof(_Complex double));
 #endif
@@ -1023,7 +1023,7 @@ void check_little_D_inversion(const int repro) {
   little_D(result, invvec); /* This should be a proper inverse now */
 
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
 #endif
 
   ldiff(invvec, result, inprod, nb_blocks * g_N_s);
diff --git a/src/lib/solver/gcr4complex_body.inc b/src/lib/solver/gcr4complex_body.inc
index 29351d87f..556838c83 100644
--- a/src/lib/solver/gcr4complex_body.inc
+++ b/src/lib/solver/gcr4complex_body.inc
@@ -218,7 +218,7 @@ _F_TYPE _PSWITCH(lsquare_norm)(_C_TYPE *const Q, const int N, const int parallel
 #ifdef TM_USE_MPI
   if (parallel) {
     double nrm2 = nrm;
-    MPI_Allreduce(&nrm2, &nrm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&nrm2, &nrm, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   }
 #endif
 
@@ -253,7 +253,7 @@ _C_TYPE _PSWITCH(lscalar_prod)(_C_TYPE *const R, _C_TYPE *const S, const int N,
 #ifdef TM_USE_MPI
   if (parallel) {
     _Complex double res2 = res;
-    MPI_Allreduce(&res2, &res, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res2, &res, 1, MPI_DOUBLE_COMPLEX, MPI_SUM, app()->mpi.comm);
   }
 #endif
 
@@ -289,7 +289,7 @@ _F_TYPE _PSWITCH(lscalar_prod_r)(_C_TYPE *const R, _C_TYPE *const S, const int N
 #ifdef TM_USE_MPI
   if (parallel) {
     double res2 = res;
-    MPI_Allreduce(&res2, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    MPI_Allreduce(&res2, &res, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   }
 #endif
 
diff --git a/src/lib/solver/little_project_eo_body.inc b/src/lib/solver/little_project_eo_body.inc
index eaca17537..7d43100eb 100644
--- a/src/lib/solver/little_project_eo_body.inc
+++ b/src/lib/solver/little_project_eo_body.inc
@@ -37,7 +37,7 @@ void _PSWITCH(little_project_eo)(_Complex _F_TYPE *const out, _Complex _F_TYPE *
   }
 
 #ifdef TM_USE_MPI
-  MPI_Allreduce(phi, psi, N, _MPI_C_TYPE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(phi, psi, N, _MPI_C_TYPE, MPI_SUM, app()->mpi.comm);
 #else
   memcpy(psi, phi, N * sizeof(_Complex _F_TYPE));
 #endif
diff --git a/src/lib/start.c b/src/lib/start.c
index 7316ec9bd..5504a83a2 100644
--- a/src/lib/start.c
+++ b/src/lib/start.c
@@ -218,7 +218,7 @@ void random_spinor_field_lexic(spinor *const k, const int repro, const enum RN_T
     } else if (g_proc_id == 0) {
       rlxd_get(rlxd_state);
     }
-    MPI_Bcast(rlxd_state, 105, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Bcast(rlxd_state, 105, MPI_INT, 0, app()->mpi.comm);
     if (g_proc_id != 0) {
       rlxd_reset(rlxd_state);
     }
@@ -289,7 +289,7 @@ void random_spinor_field_eo(spinor *const k, const int repro, const enum RN_TYPE
     } else if (g_proc_id == 0) {
       rlxd_get(rlxd_state);
     }
-    MPI_Bcast(rlxd_state, 105, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Bcast(rlxd_state, 105, MPI_INT, 0, app()->mpi.comm);
     if (g_proc_id != 0) {
       rlxd_reset(rlxd_state);
     }
@@ -432,7 +432,7 @@ void random_gauge_field(const int repro, su3 **const gf) {
     } else if (g_proc_id == 0) {
       rlxd_get(rlxd_state);
     }
-    MPI_Bcast(rlxd_state, 105, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Bcast(rlxd_state, 105, MPI_INT, 0, app()->mpi.comm);
     rlxd_reset(rlxd_state);
 #endif
     for (t0 = 0; t0 < g_nproc_t * T; t0++) {
@@ -506,7 +506,7 @@ double random_su3adj_field(const int repro, su3adj **const momenta) {
     } else if (g_proc_id == 0) {
       rlxd_get(rlxd_state);
     }
-    MPI_Bcast(rlxd_state, 105, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Bcast(rlxd_state, 105, MPI_INT, 0, app()->mpi.comm);
     rlxd_reset(rlxd_state);
 #endif
     for (int t0 = 0; t0 < g_nproc_t * T; t0++) {
@@ -589,7 +589,7 @@ double random_su3adj_field(const int repro, su3adj **const momenta) {
     kc = 0.5 * (ks + kc);
   }
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&kc, &ks, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&kc, &ks, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   return ks;
 #endif
   return kc;
@@ -838,7 +838,7 @@ void start_ranlux(int level, int seed) {
 #ifdef TM_USE_MPI
   unsigned int *seeds = calloc(g_nproc, sizeof(unsigned int));
   if (seeds == NULL) fatal_error("Memory allocation for seeds buffer failed!", "start_ranlux");
-  MPI_Gather(&loc_seed, 1, MPI_UNSIGNED, seeds, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD);
+  MPI_Gather(&loc_seed, 1, MPI_UNSIGNED, seeds, 1, MPI_UNSIGNED, 0, app()->mpi.comm);
   if (g_proc_id == 0) {
     for (int i = 0; i < g_nproc; ++i) {
       for (int j = i + 1; j < g_nproc; ++j) {
diff --git a/src/lib/test/check_xchange.c b/src/lib/test/check_xchange.c
index a20f86df4..7b4c2067b 100644
--- a/src/lib/test/check_xchange.c
+++ b/src/lib/test/check_xchange.c
@@ -85,9 +85,9 @@ int check_xchange() {
     }
 #endif
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     xchange_field(g_spinor_field[0], 0);
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
     x = (double*)&g_spinor_field[0][VOLUME / 2];
     for (i = 0; i < LX * LY * LZ / 2 * 24; i++, x++) {
@@ -95,7 +95,7 @@ int check_xchange() {
         printf("The exchange up of fields in time direction\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_up);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -107,7 +107,7 @@ int check_xchange() {
         printf("The exchange down of fields in time direction\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_dn);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -120,7 +120,7 @@ int check_xchange() {
         printf("The exchange up of fields in x direction\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_x_up);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -132,7 +132,7 @@ int check_xchange() {
         printf("The exchange down of fields in x direction\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_x_dn);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -146,7 +146,7 @@ int check_xchange() {
         printf("The exchange up of fields in y direction\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_y_up);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -159,7 +159,7 @@ int check_xchange() {
         printf("The exchange down of fields in y direction\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_y_dn);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -178,9 +178,9 @@ int check_xchange() {
       }
     }
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     xchange_field(g_spinor_field[0], 1);
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
     x = (double*)&g_spinor_field[0][VOLUME / 2 + 2 * LX * LY * LZ / 2 + 2 * T * LY * LZ / 2 +
                                     2 * T * LX * LZ / 2];
@@ -189,7 +189,7 @@ int check_xchange() {
         printf("The exchange up of fields in z (1) direction up\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_z_up);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -202,7 +202,7 @@ int check_xchange() {
         printf("The exchange down of fields in z (1) direction down\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_z_dn);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -219,9 +219,9 @@ int check_xchange() {
       }
     }
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     xchange_field(g_spinor_field[0], 1);
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
     x = (double*)&g_spinor_field[0][VOLUME / 2 + 2 * LX * LY * LZ / 2 + 2 * T * LY * LZ / 2 +
                                     2 * T * LX * LZ / 2];
@@ -230,7 +230,7 @@ int check_xchange() {
         printf("The exchange up of fields in z (0) direction up\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_z_up);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -243,7 +243,7 @@ int check_xchange() {
         printf("The exchange down of fields in z (0) direction down\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_z_dn);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -312,9 +312,9 @@ int check_xchange() {
     }
 #endif
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     xchange_gauge(g_gauge_field);
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
     x = (double*)&g_gauge_field[T * LX * LY * LZ][0];
     for (i = 0; i < LX * LY * LZ * 72; i++, x++) {
@@ -322,7 +322,7 @@ int check_xchange() {
         printf("The exchange up of gaugefields in time direction\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_up);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -334,7 +334,7 @@ int check_xchange() {
         printf("The exchange down of gaugefields in time direction\n");
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_dn);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -348,7 +348,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_x_up);
         printf("%d %d %d\n", g_cart_id, i, (int)(*x));
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -361,7 +361,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_x_dn);
         printf("%d %d %d\n", g_cart_id, i, (int)(*x));
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -376,7 +376,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_y_up);
         printf("%d %d %d\n", g_cart_id, i, (int)(*x));
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -389,7 +389,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, g_nb_y_dn);
         printf("%d %d %d\n", g_cart_id, i, (int)(*x));
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -404,7 +404,7 @@ int check_xchange() {
         printf("between %d and %d is not correct, down is %d\n", g_cart_id, g_nb_z_up, g_nb_z_dn);
         printf("%d %d %d\n", g_cart_id, i, (int)(*x));
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -418,7 +418,7 @@ int check_xchange() {
         printf("between %d and %d is not correct, up is %d\n", g_cart_id, g_nb_z_dn, g_nb_z_up);
         printf("%d %d %d\n", g_cart_id, i, (int)(*x));
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -499,9 +499,9 @@ int check_xchange() {
       }
     }
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     xchange_gauge(g_gauge_field);
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
     /* The edges */
 #if (defined TM_PARALLELXT || defined TM_PARALLELXYT || defined TM_PARALLELXYZT)
@@ -531,7 +531,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -544,7 +544,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -557,7 +557,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -570,7 +570,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -600,7 +600,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -613,7 +613,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -626,7 +626,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -639,7 +639,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -667,7 +667,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -680,7 +680,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -693,7 +693,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -706,7 +706,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -735,7 +735,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -748,7 +748,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -761,7 +761,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -774,7 +774,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -802,7 +802,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -816,7 +816,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -830,7 +830,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -844,7 +844,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -873,7 +873,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -887,7 +887,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mp);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -901,7 +901,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, pm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -915,7 +915,7 @@ int check_xchange() {
         printf("between %d and %d is not correct\n", g_cart_id, mm);
         printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
         printf("Program aborted\n");
-        MPI_Abort(MPI_COMM_WORLD, 5);
+        MPI_Abort(app()->mpi.comm, 5);
         MPI_Finalize();
         exit(0);
       }
@@ -973,9 +973,9 @@ int check_xchange() {
         }
       }
 
-      MPI_Barrier(MPI_COMM_WORLD);
+      MPI_Barrier(app()->mpi.comm);
       xchange_gauge(g_gauge_field);
-      MPI_Barrier(MPI_COMM_WORLD);
+      MPI_Barrier(app()->mpi.comm);
 
       x = (double*)&g_gauge_field[VOLUMEPLUSRAND][0];
       for (i = 0; i < LX * LY * LZ * 72; i++, x++) {
@@ -983,7 +983,7 @@ int check_xchange() {
           printf("The exchange up of gaugefields in time direction\n");
           printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_up);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -995,7 +995,7 @@ int check_xchange() {
           printf("The exchange up of gaugefields in time direction\n");
           printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_up);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1009,7 +1009,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, g_nb_x_up);
           printf("%d %d %d\n", g_cart_id, i, (int)(*x));
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1022,7 +1022,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, g_nb_x_dn);
           printf("%d %d %d\n", g_cart_id, i, (int)(*x));
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1037,7 +1037,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, g_nb_y_up);
           printf("%d %d %d\n", g_cart_id, i, (int)(*x));
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1051,7 +1051,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, g_nb_y_dn);
           printf("%d %d %d\n", g_cart_id, i, (int)(*x));
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1067,7 +1067,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, g_nb_z_up);
           printf("%d %d %d\n", g_cart_id, i, (int)(*x));
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1081,7 +1081,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, g_nb_z_dn);
           printf("%d %d %d\n", g_cart_id, i, (int)(*x));
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1108,9 +1108,9 @@ int check_xchange() {
         }
       }
 
-      MPI_Barrier(MPI_COMM_WORLD);
+      MPI_Barrier(app()->mpi.comm);
       xchange_gauge(g_gauge_field);
-      MPI_Barrier(MPI_COMM_WORLD);
+      MPI_Barrier(app()->mpi.comm);
 
       /* Now there should be in the t and t2 Rand certain values set */
 
@@ -1128,7 +1128,7 @@ int check_xchange() {
               printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_up);
               printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), g_nb_t_up);
               printf("Program aborted\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -1141,7 +1141,7 @@ int check_xchange() {
               printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_up);
               printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), g_nb_t_up);
               printf("Program aborted\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -1154,7 +1154,7 @@ int check_xchange() {
               printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_dn);
               printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), g_nb_t_dn);
               printf("Program aborted\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -1167,7 +1167,7 @@ int check_xchange() {
               printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_dn);
               printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), g_nb_t_dn);
               printf("Program aborted\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -1181,7 +1181,7 @@ int check_xchange() {
               printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_up);
               printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), g_nb_t_up);
               printf("Program aborted\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -1194,7 +1194,7 @@ int check_xchange() {
               printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_up);
               printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), g_nb_t_up);
               printf("Program aborted\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -1207,7 +1207,7 @@ int check_xchange() {
               printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_dn);
               printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), g_nb_t_dn);
               printf("Program aborted\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -1220,7 +1220,7 @@ int check_xchange() {
               printf("between %d and %d is not correct\n", g_cart_id, g_nb_t_dn);
               printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), g_nb_t_dn);
               printf("Program aborted\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -1328,9 +1328,9 @@ int check_xchange() {
         }
       }
 #endif
-      MPI_Barrier(MPI_COMM_WORLD);
+      MPI_Barrier(app()->mpi.comm);
       xchange_gauge(g_gauge_field);
-      MPI_Barrier(MPI_COMM_WORLD);
+      MPI_Barrier(app()->mpi.comm);
 
 #if (defined TM_PARALLELXT || defined TM_PARALLELXYT || defined TM_PARALLELXYZT)
       di[0] = (g_proc_coords[0] - 1) % g_nproc_t;
@@ -1355,7 +1355,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1368,7 +1368,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1381,7 +1381,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1394,7 +1394,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1407,7 +1407,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1420,7 +1420,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1433,7 +1433,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1446,7 +1446,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1477,7 +1477,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1490,7 +1490,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1503,7 +1503,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1516,7 +1516,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1529,7 +1529,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1542,7 +1542,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1555,7 +1555,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1568,7 +1568,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1596,7 +1596,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1609,7 +1609,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1622,7 +1622,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1635,7 +1635,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1648,7 +1648,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1661,7 +1661,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1674,7 +1674,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1687,7 +1687,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1717,7 +1717,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1730,7 +1730,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1744,7 +1744,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1758,7 +1758,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1772,7 +1772,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1786,7 +1786,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1800,7 +1800,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1814,7 +1814,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1844,7 +1844,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1858,7 +1858,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1872,7 +1872,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1886,7 +1886,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1900,7 +1900,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1914,7 +1914,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1928,7 +1928,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1942,7 +1942,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1972,7 +1972,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -1986,7 +1986,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -2000,7 +2000,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -2014,7 +2014,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -2028,7 +2028,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -2042,7 +2042,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, pm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), pm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -2056,7 +2056,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mp);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mp);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -2070,7 +2070,7 @@ int check_xchange() {
           printf("between %d and %d is not correct\n", g_cart_id, mm);
           printf("%d %d (%d != %d)\n", g_cart_id, i, (int)(*x), mm);
           printf("Program aborted\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -2190,9 +2190,9 @@ int check_xchange() {
     }
 #endif
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     xchange_deri(df0);
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
 #if defined TM_PARALLELT
     for (x1 = 0; x1 < LX; x1++) {
@@ -2205,7 +2205,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_t_up) {
                 printf("Exchange of derivatives is working not correctly (1u)!\n");
                 printf("Aborting program!");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2218,7 +2218,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_t_dn) {
                 printf("Exchange of derivatives is working not correctly (1d)!\n");
                 printf("Aborting program!");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2239,7 +2239,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_t_up) {
                 printf("Exchange of derivatives is working not correctly (2u)!\n");
                 printf("Aborting program!");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2252,7 +2252,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_t_dn) {
                 printf("Exchange of derivatives is working not correctly (2d)!\n");
                 printf("Aborting program!");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2261,7 +2261,7 @@ int check_xchange() {
         }
       }
     }
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     for (x0 = 1; x0 < T - 1; x0++) {
       for (x2 = 0; x2 < LY; x2++) {
         for (x3 = 0; x3 < LZ; x3++) {
@@ -2272,7 +2272,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_x_up) {
                 printf("Exchange of derivatives is working not correctly (3u)!\n");
                 printf("Aborting program!");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2285,7 +2285,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_x_dn) {
                 printf("Exchange of derivatives is working not correctly (3d)!\n");
                 printf("Aborting program!");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2303,7 +2303,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_up + g_nb_t_up) {
               printf("Exchange of derivatives is working not correctly (4uu)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2316,7 +2316,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_up + g_nb_t_dn) {
               printf("Exchange of derivatives is working not correctly (4ud)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2329,7 +2329,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_dn + g_nb_t_up) {
               printf("Exchange of derivatives is working not correctly (4du)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2342,7 +2342,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_dn + g_nb_t_dn) {
               printf("Exchange of derivatives is working not correctly (4dd)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2362,7 +2362,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_t_up) {
                 printf("Exchange of derivatives is working not correctly (5u)!\n");
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2375,7 +2375,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_t_dn) {
                 printf("Exchange of derivatives is working not correctly (5d)!\n");
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2394,7 +2394,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_x_up) {
                 printf("Exchange of derivatives is working not correctly (6u)!\n");
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2407,7 +2407,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_x_dn) {
                 printf("Exchange of derivatives is working not correctly (6d)!\n");
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2426,7 +2426,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_y_up) {
                 printf("Exchange of derivatives is working not correctly (7u)!\n");
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2439,7 +2439,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_y_dn) {
                 printf("Exchange of derivatives is working not correctly (7d)!\n");
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2457,7 +2457,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_up + g_nb_t_up) {
               printf("Exchange of derivatives is working not correctly (8uu)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2470,7 +2470,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_up + g_nb_t_dn) {
               printf("Exchange of derivatives is working not correctly (8ud)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2483,7 +2483,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_dn + g_nb_t_up) {
               printf("Exchange of derivatives is working not correctly (8du)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2496,7 +2496,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_dn + g_nb_t_dn) {
               printf("Exchange of derivatives is working not correctly (8dd)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2513,7 +2513,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_t_up + g_nb_y_up) {
               printf("Exchange of derivatives is working not correctly (9uu)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2526,7 +2526,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_t_up + g_nb_y_dn) {
               printf("Exchange of derivatives is working not correctly (9ud)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2541,7 +2541,7 @@ int check_xchange() {
               printf("%d %d %d %d %d %d %d\n", (int)x[j], g_nb_t_dn, g_nb_t_up, g_nb_y_dn,
                      g_nb_y_up, x1, x3);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2554,7 +2554,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_t_dn + g_nb_y_dn) {
               printf("Exchange of derivatives is working not correctly (9dd)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2571,7 +2571,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_up + g_nb_y_up) {
               printf("Exchange of derivatives is working not correctly (10uu)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2584,7 +2584,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_up + g_nb_y_dn) {
               printf("Exchange of derivatives is working not correctly (10ud)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2597,7 +2597,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_dn + g_nb_y_up) {
               printf("Exchange of derivatives is working not correctly (10du)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2610,7 +2610,7 @@ int check_xchange() {
             if ((int)x[j] != g_nb_x_dn + g_nb_y_dn) {
               printf("Exchange of derivatives is working not correctly (10dd)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2626,7 +2626,7 @@ int check_xchange() {
           if ((int)x[j] != g_nb_x_up + g_nb_y_up + g_nb_t_up) {
             printf("Exchange of derivatives is working not correctly (11uuu)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2639,7 +2639,7 @@ int check_xchange() {
           if ((int)x[j] != g_nb_x_dn + g_nb_y_up + g_nb_t_up) {
             printf("Exchange of derivatives is working not correctly (11duu)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2652,7 +2652,7 @@ int check_xchange() {
           if ((int)x[j] != g_nb_x_dn + g_nb_y_up + g_nb_t_dn) {
             printf("Exchange of derivatives is working not correctly (11dud)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2665,7 +2665,7 @@ int check_xchange() {
           if ((int)x[j] != g_nb_x_dn + g_nb_y_dn + g_nb_t_up) {
             printf("Exchange of derivatives is working not correctly (11ddu)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2678,7 +2678,7 @@ int check_xchange() {
           if ((int)x[j] != g_nb_x_up + g_nb_y_dn + g_nb_t_dn) {
             printf("Exchange of derivatives is working not correctly (11udd)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2691,7 +2691,7 @@ int check_xchange() {
           if ((int)x[j] != g_nb_x_up + g_nb_y_up + g_nb_t_dn) {
             printf("Exchange of derivatives is working not correctly (11uud)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2704,7 +2704,7 @@ int check_xchange() {
           if ((int)x[j] != g_nb_x_up + g_nb_y_dn + g_nb_t_up) {
             printf("Exchange of derivatives is working not correctly (11udu)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2717,7 +2717,7 @@ int check_xchange() {
           if ((int)x[j] != g_nb_x_dn + g_nb_y_dn + g_nb_t_dn) {
             printf("Exchange of derivatives is working not correctly (11ddd)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2735,7 +2735,7 @@ int check_xchange() {
                 if ((int)x[j] != 0) {
                   printf("Exchange of derivatives is working not correctly (bulk XYT)!\n");
                   printf("Aborting program!\n");
-                  MPI_Abort(MPI_COMM_WORLD, 5);
+                  MPI_Abort(app()->mpi.comm, 5);
                   MPI_Finalize();
                   exit(0);
                 }
@@ -2761,7 +2761,7 @@ int check_xchange() {
                 printf("%d %d %d %d %d\n", x1, x2, x3, ix, g_proc_id);
                 printf("%f %d %d\n", df0[ix][mu].d8, g_nb_t_up, g_nb_t_dn);
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2780,7 +2780,7 @@ int check_xchange() {
               if ((int)x[j] != g_nb_x_up) {
                 printf("Exchange of derivatives is working not correctly (13)!\n");
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2800,7 +2800,7 @@ int check_xchange() {
                 printf("Exchange of derivatives is working not correctly (14)!\n");
                 printf("%d %d %d %d %d\n", x0, x1, x3, ix, g_proc_id);
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2820,7 +2820,7 @@ int check_xchange() {
                 printf("Exchange of derivatives is working not correctly (15)!\n");
                 printf("%d %d %d %d %d\n", x0, x1, x3, ix, g_proc_id);
                 printf("Aborting program!\n");
-                MPI_Abort(MPI_COMM_WORLD, 5);
+                MPI_Abort(app()->mpi.comm, 5);
                 MPI_Finalize();
                 exit(0);
               }
@@ -2838,7 +2838,7 @@ int check_xchange() {
             if ((int)x[j] != (g_nb_x_up + g_nb_t_up)) {
               printf("Exchange of derivatives is working not correctly (16)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2855,7 +2855,7 @@ int check_xchange() {
             if ((int)x[j] != (g_nb_y_up + g_nb_t_up)) {
               printf("Exchange of derivatives is working not correctly (17)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2872,7 +2872,7 @@ int check_xchange() {
             if ((int)x[j] != (g_nb_y_up + g_nb_x_up)) {
               printf("Exchange of derivatives is working not correctly (18)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2890,7 +2890,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (19)!\n");
               printf("%f %d %d %d\n", df0[ix][mu].d1, g_nb_x_up + g_nb_z_up, g_nb_x_up, g_nb_z_up);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2907,7 +2907,7 @@ int check_xchange() {
             if ((int)x[j] != (g_nb_y_up + g_nb_z_up)) {
               printf("Exchange of derivatives is working not correctly (20)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2924,7 +2924,7 @@ int check_xchange() {
             if ((int)x[j] != (g_nb_t_up + g_nb_z_up)) {
               printf("Exchange of derivatives is working not correctly (21)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -2940,7 +2940,7 @@ int check_xchange() {
           if ((int)x[j] != (g_nb_t_up + g_nb_x_up + g_nb_y_up)) {
             printf("Exchange of derivatives is working not correctly (22)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2955,7 +2955,7 @@ int check_xchange() {
           if ((int)x[j] != (g_nb_t_up + g_nb_x_up + g_nb_z_up)) {
             printf("Exchange of derivatives is working not correctly (23)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2970,7 +2970,7 @@ int check_xchange() {
           if ((int)x[j] != (g_nb_t_up + g_nb_z_up + g_nb_y_up)) {
             printf("Exchange of derivatives is working not correctly (24)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2985,7 +2985,7 @@ int check_xchange() {
           if ((int)x[j] != (g_nb_z_up + g_nb_x_up + g_nb_y_up)) {
             printf("Exchange of derivatives is working not correctly (25)!\n");
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -2999,7 +2999,7 @@ int check_xchange() {
         if ((int)x[j] != (g_nb_z_up + g_nb_x_up + g_nb_y_up + g_nb_t_up)) {
           printf("Exchange of derivatives is working not correctly (26)!\n");
           printf("Aborting program!\n");
-          MPI_Abort(MPI_COMM_WORLD, 5);
+          MPI_Abort(app()->mpi.comm, 5);
           MPI_Finalize();
           exit(0);
         }
@@ -3135,9 +3135,9 @@ int check_xchange() {
 
 #endif
 
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
     xchange_deri(df0);
-    MPI_Barrier(MPI_COMM_WORLD);
+    MPI_Barrier(app()->mpi.comm);
 
 #if (defined TM_PARALLELXT || defined TM_PARALLELXYT || defined TM_PARALLELXYZT)
 
@@ -3171,7 +3171,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e5mm)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_t_up, g_nb_x_up, pp, mm);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3185,7 +3185,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e5mp)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_t_up, g_nb_x_up, pm, mp);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3199,7 +3199,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e5pm)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_t_up, g_nb_x_up, pm, mp);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3213,7 +3213,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e5pp)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_t_up, g_nb_x_up, pp, mm);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3251,7 +3251,7 @@ int check_xchange() {
             if ((int)x[j] != mm) {
               printf("Exchange of derivatives is working not correctly (e6mm)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3265,7 +3265,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e6pm)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_x_up, g_nb_y_up, pm, mp);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3279,7 +3279,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e6mp)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_x_up, g_nb_y_up, pm, mp);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3293,7 +3293,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e6pp)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_x_up, g_nb_y_up, pp, mm);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3326,7 +3326,7 @@ int check_xchange() {
             if ((int)x[j] != mm) {
               printf("Exchange of derivatives is working not correctly (e7mm)!\n");
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3340,7 +3340,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e7pm)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_t_up, g_nb_y_up, pm, pm);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3354,7 +3354,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e7mp)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_t_up, g_nb_y_up, pm, mp);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3368,7 +3368,7 @@ int check_xchange() {
               printf("Exchange of derivatives is working not correctly (e7pp)!\n");
               printf("%f %d %d %d %d\n", x[j], g_nb_t_up, g_nb_y_up, pp, mm);
               printf("Aborting program!\n");
-              MPI_Abort(MPI_COMM_WORLD, 5);
+              MPI_Abort(app()->mpi.comm, 5);
               MPI_Finalize();
               exit(0);
             }
@@ -3402,7 +3402,7 @@ int check_xchange() {
             printf("Exchange of derivatives is working not correctly (e8mmm)!\n");
             printf("%d %d %d %d %d\n", (int)x[j], mm, mp, pm, pp);
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -3435,7 +3435,7 @@ int check_xchange() {
             printf("Exchange of derivatives is working not correctly (e8pmm)!\n");
             printf("%d %d %d %d %d\n", (int)x[j], mm, mp, pm, pp);
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -3468,7 +3468,7 @@ int check_xchange() {
             printf("Exchange of derivatives is working not correctly (e8pmp)!\n");
             printf("%d %d %d %d %d\n", (int)x[j], mm, mp, pm, pp);
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -3501,7 +3501,7 @@ int check_xchange() {
             printf("Exchange of derivatives is working not correctly (e8ppp)!\n");
             printf("%d %d %d %d %d\n", (int)x[j], mm, mp, pm, pp);
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -3534,7 +3534,7 @@ int check_xchange() {
             printf("Exchange of derivatives is working not correctly (e8mpm)!\n");
             printf("%d %d %d %d %d\n", (int)x[j], mm, mp, pm, pp);
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -3567,7 +3567,7 @@ int check_xchange() {
             printf("Exchange of derivatives is working not correctly (e8mmp)!\n");
             printf("%d %d %d %d %d\n", (int)x[j], mm, mp, pm, pp);
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -3600,7 +3600,7 @@ int check_xchange() {
             printf("Exchange of derivatives is working not correctly (e8mpp)!\n");
             printf("%d %d %d %d %d\n", (int)x[j], mm, mp, pm, pp);
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -3633,7 +3633,7 @@ int check_xchange() {
             printf("Exchange of derivatives is working not correctly (e8ppm)!\n");
             printf("%d %d %d %d %d\n", (int)x[j], mm, mp, pm, pp);
             printf("Aborting program!\n");
-            MPI_Abort(MPI_COMM_WORLD, 5);
+            MPI_Abort(app()->mpi.comm, 5);
             MPI_Finalize();
             exit(0);
           }
@@ -3651,7 +3651,7 @@ int check_xchange() {
                 if ((int)x[j] != 0) {
                   printf("Exchange of derivatives is working not correctly (ebulk XYT)!\n");
                   printf("Aborting program!\n");
-                  MPI_Abort(MPI_COMM_WORLD, 5);
+                  MPI_Abort(app()->mpi.comm, 5);
                   MPI_Finalize();
                   exit(0);
                 }
diff --git a/src/lib/test/measure_rectangles.debug.c b/src/lib/test/measure_rectangles.debug.c
index 422f681b2..9c4f51215 100644
--- a/src/lib/test/measure_rectangles.debug.c
+++ b/src/lib/test/measure_rectangles.debug.c
@@ -131,7 +131,7 @@ double measure_rectangles() {
   fclose(debugfile);
   ga = (kc + ks) / 3.0;
 #ifdef TM_USE_MPI
-  MPI_Allreduce(&ga, &gas, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&ga, &gas, 1, MPI_DOUBLE, MPI_SUM, app()->mpi.comm);
   return gas;
 #else
   return ga;
diff --git a/src/lib/update_tm.c b/src/lib/update_tm.c
index 3f1cdc5d5..df0e0fd54 100644
--- a/src/lib/update_tm.c
+++ b/src/lib/update_tm.c
@@ -174,7 +174,7 @@ int update_tm(double *plaquette_energy, double *rectangle_energy, char *filename
      the other sites */
   ranlxd(yy, 1);
 #ifdef TM_USE_MPI
-  MPI_Bcast(&yy[0], 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&yy[0], 1, MPI_DOUBLE, 0, app()->mpi.comm);
 #endif
 
   /* when acctest is 0 (i.e. do not perform acceptance test), the trajectory is accepted whatever
@@ -277,7 +277,7 @@ int update_tm(double *plaquette_energy, double *rectangle_energy, char *filename
 
 #ifdef TM_USE_MPI
     tmp = ret_gauge_diff;
-    MPI_Reduce(&tmp, &ret_gauge_diff, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+    MPI_Reduce(&tmp, &ret_gauge_diff, 1, MPI_DOUBLE, MPI_SUM, 0, app()->mpi.comm);
 #endif
     /* compute the total H */
     tmp = enep;
diff --git a/src/lib/wrapper/lib_wrapper.c b/src/lib/wrapper/lib_wrapper.c
index 19d36ddc6..71d78bed4 100644
--- a/src/lib/wrapper/lib_wrapper.c
+++ b/src/lib/wrapper/lib_wrapper.c
@@ -98,7 +98,7 @@ int tmLQCD_invert_init(int argc, char* argv[], const int _verbose, const int ext
   g_use_clover_flag = 0;
 
 #ifdef TM_USE_MPI
-  MPI_Comm_rank(MPI_COMM_WORLD, &g_proc_id);
+  MPI_Comm_rank(app()->mpi.comm, &g_proc_id);
 #else
   g_proc_id = 0;
 #endif
@@ -341,7 +341,7 @@ int tmLQCD_finalise() {
   free_moment_field();
   free_chi_spinor_field();
 #ifdef TM_USE_MPI
-  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Barrier(app()->mpi.comm);
 #endif
   return (0);
 }
diff --git a/tests/test_buffers.c b/tests/test_buffers.c
index d398c9768..7fedf2fe7 100644
--- a/tests/test_buffers.c
+++ b/tests/test_buffers.c
@@ -13,7 +13,7 @@ TEST_SUITES{TEST_SUITE_ADD(BUFFERS_GAUGE), TEST_SUITES_CLOSURE};
 int main(int argc, char *argv[]) {
 #ifdef TM_USE_MPI
   MPI_Init(&argc, &argv);
-  MPI_Comm_rank(MPI_COMM_WORLD, &g_proc_id);
+  MPI_Comm_rank(app()->mpi.comm, &g_proc_id);
 #else
   g_proc_id = 0;
 #endif

From cf852e2f82fe5480fc02289c862a23032711ea9a Mon Sep 17 00:00:00 2001
From: JingJing Li <lijingjing416@gmail.com>
Date: Mon, 20 Apr 2026 17:09:29 +0200
Subject: [PATCH 2/8] update cmake

---
 src/lib/CMakeLists.txt       |   7 ++
 src/lib/include/ptbc.h       |   6 +
 src/lib/init/init_parallel.c |   2 +
 src/lib/ptbc/defect.c        |  62 ++++++++++
 src/lib/read_input.l         | 225 ++++++++++++++++++++++++++++++++++-
 5 files changed, 301 insertions(+), 1 deletion(-)
 create mode 100644 src/lib/include/ptbc.h
 create mode 100644 src/lib/ptbc/defect.c

diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt
index 9c5322997..2793fa91c 100644
--- a/src/lib/CMakeLists.txt
+++ b/src/lib/CMakeLists.txt
@@ -78,6 +78,12 @@ list(
   ALLOC_SRC_C
   alloc/alloc.c)
 
+list(
+  APPEND
+  PTBC_SRC_C
+  ptbc/defect.c
+)
+
 list(
   APPEND
   SOLVER_SRC_C
@@ -380,6 +386,7 @@ list(
   ${INIT_SRC_C}
   ${APP_CONTEXT_SRC_C}
   ${ALLOC_SRC_C}
+  ${PTBC_SRC_C}
   ${SOLVER_SRC_C}
   ${TEST_SRC_C}
   ${MEAS_SRC_C}
diff --git a/src/lib/include/ptbc.h b/src/lib/include/ptbc.h
new file mode 100644
index 000000000..42b434cf6
--- /dev/null
+++ b/src/lib/include/ptbc.h
@@ -0,0 +1,6 @@
+
+#include <mpi.h>
+
+
+bool is_defect(PTBCDefect *def, int const ix, int const mu);
+double get_ptbc_coeff(int const ix, int const mu);
\ No newline at end of file
diff --git a/src/lib/init/init_parallel.c b/src/lib/init/init_parallel.c
index 9dfdbb0c5..7d830e8bb 100644
--- a/src/lib/init/init_parallel.c
+++ b/src/lib/init/init_parallel.c
@@ -85,6 +85,8 @@ void init_parallel_and_read_input(int argc, char *argv[], const char input_filen
     exit(-1);
   }
 
+  app()->ptbc.initialize();
+
 #ifdef TM_USE_OMP
   init_openmp();
 #endif
diff --git a/src/lib/ptbc/defect.c b/src/lib/ptbc/defect.c
new file mode 100644
index 000000000..176dd7241
--- /dev/null
+++ b/src/lib/ptbc/defect.c
@@ -0,0 +1,62 @@
+#include <stdbool.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <mpi.h>
+#include <ptbc.h>
+#include "global.h"
+
+
+// check if a link lie in defect region
+// 1. last slice of along 2. 3d cube with 0 origin
+bool is_defect(PTBCDefect *def, int const ix, int const mu) {
+  if (ix >= VOLUME) return false;
+  
+  int const global_dim[4] = {T*g_nproc_t, LX*g_nproc_t, LY*g_nproc_x, LZ*g_nproc_z};
+
+  // check start point is at second last slice and link stays on the slice
+  if (g_coord[ix][def->along] == global_dim[def->along] - 2 && mu != def->along){
+    // xf = end point of link 
+    int xf[4] = {g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3]};
+    xf[mu]  = xf[mu] + 1;
+    int dim3[3];
+    int count = 0;
+    for (int d=0; d<4; d++) {
+      if (d != def->along) {
+        dim3[count] = d;
+        count++;
+      }
+    }
+    
+    for (int d=0; d<3; d++) {
+      if (g_coord[ix][dim3[d]] >= def->Ld[d] || xf[dim3[d]] >= def->Ld[d]) {
+        return false;
+      }
+    }
+    if (g_proc_id==0) printf("defect coordinate %d %d %d %d mu = %d\n", g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3], mu);
+    return true;
+  }
+  else {
+    return false;
+  }
+}
+
+// get multiplying factor of parallel tempering locally (assumne no overlapping defects!)
+double get_ptbc_coeff(int const ix, int const mu) {
+  int const inst = app()->ptbc.instance_id;
+  const PTBCInstance *instance = &(app()->ptbc.instances[inst]);
+  
+  // loop over defects
+  for (int i=0; i<instance->n_coeffs; i++) {
+    PTBCDefect *def = instance->defects[i];
+    // apply coeff if within defect
+    if (is_defect(def, ix, mu)) {
+      return instance->coefficients[i];
+    }
+  }
+
+  return 1.;
+}
\ No newline at end of file
diff --git a/src/lib/read_input.l b/src/lib/read_input.l
index 5eb542f87..8b2f068eb 100644
--- a/src/lib/read_input.l
+++ b/src/lib/read_input.l
@@ -52,6 +52,7 @@ EQL {SPC}*={SPC}*
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
+#include <stdarg.h>
 #define INIT_GLOBALS
 #include "global.h"
 #undef INIT_GLOBALS
@@ -61,12 +62,14 @@ EQL {SPC}*={SPC}*
 #include "solver/solver_types.h"
 #include "meas/measurements.h"
 #include "integrator.h"
+#include "include/app.h"
 #include "operator.h"
 #include "phmc.h"
 #include <io/params.h>
 #include "qphix_types.h"
 #include "quda_types.h"
 #include "misc_types.h"
+#include "alloc.h"
 
 #include <ctype.h>
 
@@ -87,6 +90,28 @@ static inline void rmQuotes(char *str){
   *strsave='\0';
 }
 
+
+static inline int strlist_count_tokens(char const * const input)
+{
+  int n = 0;
+  for (int i = 0; i < strlen(input); i++)
+    if (input[i] == ',') n++;
+  return n+1;
+}
+
+static inline void fail(const bool test, const char* fmt, ...)
+{
+  if (test) {
+    va_list args;
+    char msg[1024];
+    va_start(args, fmt);
+    vsnprintf(msg, 1024, fmt, args);
+    va_end(args);
+    yy_fatal_error(msg);
+  }
+}
+
+
 /* tokenize the comma-delimited list 'input' of the form
     'list = token1, token2, ...'
    and return the first token, which is the name of the list
@@ -473,6 +498,60 @@ static inline double fltlist_next_token(int * const list_end){
     free(input_copy);
   }
   
+  static inline void parse_int_par_array(char const * const input, int * par_array, const int max_size) {
+    char paramname[100];
+    char error_message[ERR_MSG_LEN];
+    int list_end = 0;
+    int element = 0;
+
+    char * input_copy = (char*)NULL;
+    strlist_tokenize(input, &input_copy, paramname, 100);
+    int parval = (int) fltlist_next_token(&list_end);
+    while( list_end != 1 ){
+      if( element >= max_size ){
+        snprintf(error_message, ERR_MSG_LEN, "Exceeded maximum number of elements (%d) parsing %s!\n", max_size, paramname);
+        yy_fatal_error(error_message);
+      }
+
+      par_array[element] = parval;
+      if(myverbose){
+        printf("  %s, element %d set to %d line %d\n", paramname,
+                element, par_array[element], line_of_file);
+      }
+
+      element++;
+      parval = (int) fltlist_next_token(&list_end);
+    }
+    free(input_copy);
+  }
+
+  static inline void parse_dbl_par_array(char const * const input, double * par_array, const int max_size) {
+    char paramname[100];
+    char error_message[ERR_MSG_LEN];
+    int list_end = 0;
+    int element = 0;
+
+    char * input_copy = (char*)NULL;
+    strlist_tokenize(input, &input_copy, paramname, 100);
+    double parval = fltlist_next_token(&list_end);
+    while( list_end != 1 ){
+      if( element >= max_size ){
+        snprintf(error_message, ERR_MSG_LEN, "Exceeded maximum number of elements (%d) parsing %s!\n", max_size, paramname);
+        yy_fatal_error(error_message);
+      }
+
+      par_array[element] = parval;
+      if(myverbose){
+        printf("  %s, element %d set to %e line %d\n", paramname,
+                element, par_array[element], line_of_file);
+      }
+
+      element++;
+      parval = fltlist_next_token(&list_end);
+    }
+    free(input_copy);
+  }
+
 %}
 
 %option never-interactive
@@ -556,6 +635,11 @@ static inline double fltlist_next_token(int * const list_end){
 %x INITINTEGRATOR
 %x INTEGRATOR
 
+%x INITPTBC
+%x PTBC
+%x PTBCDEFECT
+%x PTBCINSTANCE
+
 %x DEFLATION
 %x INITDEFLATION
 
@@ -712,6 +796,9 @@ static inline double fltlist_next_token(int * const list_end){
 ^BeginInt                          BEGIN(INITINTEGRATOR);
 ^BeginOperator{SPC}+               BEGIN(INITOPERATOR);
 
+
+^BeginPTBC                         BEGIN(INITPTBC);
+
 ^BeginExternalInverter{SPC}+       BEGIN(INITEXTERNALINVERTER);
 
 ^BeginTuneMGParams{SPC}+           BEGIN(TUNEMGPARAMS);
@@ -3071,6 +3158,142 @@ static inline double fltlist_next_token(int * const list_end){
   }
 }
 
+<INITPTBC>{SPC}* {
+  appm()->ptbc.active = false;
+  appm()->ptbc.n_instances = 0;
+  appm()->ptbc.n_defects = 0;
+  appm()->ptbc.instance_id = 0;
+  if(myverbose) printf("Initialising PTBC line %d\n", line_of_file);
+  BEGIN(PTBC);
+}
+<PTBC>{
+  {SPC}*BeginPTBCDefect{SPC}+{DIGIT}+ {
+    sscanf(yytext, " %[a-zA-Z] %d", name, &b);
+    fail(b<0 || b>MAX_N_DEFECTS, "PTBC Defect id = %d out of bounds!", b);
+    appm()->ptbc.defects[b].active = true;
+    appm()->ptbc.n_defects++;
+    fail(appm()->ptbc.n_defects>=MAX_N_DEFECTS, "To many PTBC defects! Limit is %d", MAX_N_DEFECTS);
+    if(myverbose) printf("  Initialising PTBC defect (index = %d) in line %d\n", b, line_of_file);
+    BEGIN(PTBCDEFECT);
+  }
+  <PTBCDEFECT>{
+    {SPC}*Extents{EQL}{STRLIST} {
+
+      PTBCDefect* cd = appm()->ptbc.defects + b;
+
+      int n_parameters = strlist_count_tokens(yytext);
+      fail(n_parameters != 3, "Coefficients must provide exactly 3 numbers, not %d!\n", n_parameters);
+      sscanf(yytext, " %[a-zA-Z] = %d, %d, %d", name, cd->Ld, cd->Ld+1, cd->Ld+2);
+      if(myverbose) printf("    Ld = [%d, %d, %d] line %d\n", cd->Ld[0], cd->Ld[1], cd->Ld[2], line_of_file);
+      BEGIN(PTBCDEFECT);
+    }
+    {SPC}*Along{EQL}T {
+      appm()->ptbc.defects[b].along = DIRECTION_T;
+      if(myverbose) printf("    PTBC defect along T direction line %d\n", line_of_file);
+      BEGIN(PTBCDEFECT);
+    }
+    {SPC}*Along{EQL}X {
+      appm()->ptbc.defects[b].along = DIRECTION_X;
+      if(myverbose) printf("    PTBC defect along X direction line %d\n", line_of_file);
+      BEGIN(PTBCDEFECT);
+    }
+    {SPC}*Along{EQL}Y {
+      appm()->ptbc.defects[b].along = DIRECTION_Y;
+      if(myverbose) printf("    PTBC defect along Y direction line %d\n", line_of_file);
+      BEGIN(PTBCDEFECT);
+    }
+    {SPC}*Along{EQL}Z {
+      appm()->ptbc.defects[b].along = DIRECTION_Z;
+      if(myverbose) printf("    PTBC defect along Z direction line %d\n", line_of_file);
+      BEGIN(PTBCDEFECT);
+    }
+    {SPC}*EndPTBCDefect{SPC}* {
+      if(myverbose) printf("  PTBC defect parsed line %d\n", line_of_file);
+      BEGIN(PTBC);
+    }
+  }
+  {SPC}*BeginPTBCInstance{SPC}+{DIGIT}+ {
+    sscanf(yytext, " %[a-zA-Z] %d", name, &b);
+    fail(b<0 || b>MAX_N_INSTANCES, "PTBC Instance id = %d out of bounds!", b);
+    appm()->ptbc.instances[b].active = true;
+    appm()->ptbc.instances[b].n_coeffs = -1;
+    appm()->ptbc.n_instances++;
+    fail(appm()->ptbc.n_instances>=MAX_N_INSTANCES, "To many PTBC instances! Limit is %d", MAX_N_INSTANCES);
+    if(myverbose) printf("  Initialising PTBC instance (index = %d) in line %d\n", b, line_of_file);
+    BEGIN(PTBCINSTANCE);
+  }
+  <PTBCINSTANCE>{
+    {SPC}*Coefficients{EQL}{STRLIST} {
+
+      PTBCInstance* ci = appm()->ptbc.instances + b;
+
+      int n_parameters = strlist_count_tokens(yytext);
+      if (ci->coefficients == NULL)
+        ci->coefficients = safe_malloc(n_parameters*sizeof(double*));
+
+      parse_dbl_par_array(yytext, ci->coefficients, n_parameters);
+
+      for (int i = 0; i < n_parameters; ++i) {
+        fail(ci->coefficients[i]<0 || ci->coefficients[i] > 1,
+          "Coefficient %d with value %e is out of bounds", i, ci->coefficients[i]);
+      }
+
+      if (ci->n_coeffs == -1) ci->n_coeffs = n_parameters;
+      fail(ci->n_coeffs != n_parameters, "Number of PTBC coefficients must be equal to number of defects!\n");
+
+      for (int i = 0; i < n_parameters; ++i) {
+        if(myverbose) printf("    Coefficients[%d] = %f line %d\n", i, ci->coefficients[i], line_of_file);
+      }
+
+      BEGIN(PTBCINSTANCE);
+    }
+    {SPC}*Defects{EQL}{STRLIST} {
+
+      PTBCInstance* ci = appm()->ptbc.instances + b;
+
+      int n_parameters = strlist_count_tokens(yytext);
+      if (ci->defects == NULL) {
+        ci->defects = safe_malloc(n_parameters*sizeof(PTBCDefect*));
+      }
+
+      int* ids = safe_malloc(n_parameters*sizeof(int*));
+      parse_int_par_array(yytext, ids, n_parameters);
+
+      for (int i = 0; i < n_parameters; ++i) {
+        if(myverbose) printf("    Defects = %d line %d\n", ids[i], line_of_file);
+        fail(ids[i]<0 || ids[i]>MAX_N_DEFECTS, "PTBC Defect id = %d out of bounds!", ids[i]);
+        ci->defects[i] = &(appm()->ptbc.defects[ids[i]]);
+      }
+
+      if (ci->n_coeffs == -1) ci->n_coeffs = n_parameters;
+
+      fail(ci->n_coeffs != n_parameters,
+        "Number of PTBC defects must be equal to number of coefficients!\n");
+
+      free(ids);
+      BEGIN(PTBCINSTANCE);
+    }
+    {SPC}*EndPTBCInstance{SPC}* {
+      if(myverbose) printf("  PTBC instance parsed line %d\n", line_of_file);
+      BEGIN(PTBC);
+    }
+  }
+  EndPTBC{SPC}* {
+    fail(app()->ptbc.n_instances == 1, "Number of PTBC chains must be larger than 1!\n");
+
+    for (int i = 0; i < app()->ptbc.n_instances; ++i) {
+      fail(app()->ptbc.instances[i].active == false, "PTBC instance %d is not active!\n", i);
+
+      for (int j = 0; j < app()->ptbc.instances[i].n_coeffs; ++j) {
+        fail(app()->ptbc.instances[i].defects[j]->active == false, "PTBC instance %d does not refer to a valid defect!\n", i);
+      }
+    }
+
+    if(myverbose) printf("PTBC parsed line %d\n\n", line_of_file);
+    BEGIN(0);
+  }
+}
+
 <SOURCETYPE>{
   Point {
     SourceInfo.type = SRC_TYPE_POINT;
@@ -3729,7 +3952,7 @@ static inline double fltlist_next_token(int * const list_end){
   BEGIN(comment_caller);
 }
 
-<INITMONOMIAL,DETMONOMIAL,CLDETMONOMIAL,CLDETRATMONOMIAL,CLDETRATRWMONOMIAL,NDPOLYMONOMIAL,NDRATMONOMIAL,NDRATCORMONOMIAL,NDCLRATMONOMIAL,NDCLRATCORMONOMIAL,CLPOLYMONOMIAL,GAUGEMONOMIAL,INTEGRATOR,INITINTEGRATOR,INITMEASUREMENT,PIONNORMMEAS,ONLINEMEAS,ORIENTEDPLAQUETTESMEAS,GRADIENTFLOWMEAS,INITOPERATOR,TMOP,DBTMOP,OVERLAPOP,WILSONOP,CLOVEROP,DBCLOVEROP,POLYMONOMIAL,PLOOP,RATMONOMIAL,RATCORMONOMIAL,CLRATMONOMIAL,CLRATCORMONOMIAL,INITDEFLATION,DEFLATION,INITMULTIGRID,MULTIGRID,INITEXTERNALINVERTER,QUDAINVERTER,QPHIXINVERTER,NDDETRATMONOMIAL,NDCLDETRATMONOMIAL,TUNEQUDAMGPARAMS>{SPC}*\n   {
+<INITMONOMIAL,DETMONOMIAL,CLDETMONOMIAL,CLDETRATMONOMIAL,CLDETRATRWMONOMIAL,NDPOLYMONOMIAL,NDRATMONOMIAL,NDRATCORMONOMIAL,NDCLRATMONOMIAL,NDCLRATCORMONOMIAL,CLPOLYMONOMIAL,GAUGEMONOMIAL,INTEGRATOR,INITINTEGRATOR,PTBC,INITPTBC,PTBCDEFECT,PTBCINSTANCE,INITMEASUREMENT,PIONNORMMEAS,ONLINEMEAS,ORIENTEDPLAQUETTESMEAS,GRADIENTFLOWMEAS,INITOPERATOR,TMOP,DBTMOP,OVERLAPOP,WILSONOP,CLOVEROP,DBCLOVEROP,POLYMONOMIAL,PLOOP,RATMONOMIAL,RATCORMONOMIAL,CLRATMONOMIAL,CLRATCORMONOMIAL,INITDEFLATION,DEFLATION,INITMULTIGRID,MULTIGRID,INITEXTERNALINVERTER,QUDAINVERTER,QPHIXINVERTER,NDDETRATMONOMIAL,NDCLDETRATMONOMIAL,TUNEQUDAMGPARAMS>{SPC}*\n   {
   line_of_file++;
 }
 <*>{SPC}*\n                       {

From 655052f29a12b3485c72f943c63a24a557dbb40b Mon Sep 17 00:00:00 2001
From: JingJing Li <lijingjing416@gmail.com>
Date: Mon, 20 Apr 2026 18:12:00 +0200
Subject: [PATCH 3/8] defect added

---
 src/lib/get_rectangle_staples.c | 31 +++++++++++++++++++++++++------
 src/lib/get_staples.c           | 25 +++++++++++++++++++------
 src/lib/measure_gauge_action.c  | 18 ++++++++++++++----
 src/lib/measure_rectangles.c    |  6 +++++-
 src/lib/su3.h                   | 33 +++++++++++++++++++++++++++++++++
 5 files changed, 96 insertions(+), 17 deletions(-)

diff --git a/src/lib/get_rectangle_staples.c b/src/lib/get_rectangle_staples.c
index c8f69596b..78ab70b4b 100644
--- a/src/lib/get_rectangle_staples.c
+++ b/src/lib/get_rectangle_staples.c
@@ -25,6 +25,7 @@
 #include "get_rectangle_staples.h"
 #include "global.h"
 #include "su3.h"
+#include "ptbc.h"
 
 void get_rectangle_staples(su3 *const v, const int x, const int mu) {
   get_rectangle_staples_general(v, x, mu, (const su3 *const *const)g_gauge_field);
@@ -51,14 +52,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
       _su3_times_su3(tmp1, *a, *b);
       z = g_iup[y][nu];
       c = &gf[z][mu];
+      double const ptbc_fac0 = get_ptbc_coeff(x, nu) * get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, mu);
       _su3_times_su3(tmp2, tmp1, *c);
 
       y = g_iup[x][mu];
       d = &gf[y][nu];
       z = g_iup[y][nu];
       e = &gf[z][nu];
+      double const ptbc_fac1 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, nu);
       _su3_times_su3(tmp1, *d, *e);
-      _su3_times_su3d_acc((*v), tmp2, tmp1);
+      //_su3_times_su3d_acc((*v), tmp2, tmp1);
+      _real_times_su3_times_su3d_acc((*v), tmp2, tmp1, ptbc_fac0 * ptbc_fac1);
 
       /* 1 contr. starting idn[idn[x][nu]][nu]
        * e^+ d^+ a b c
@@ -73,14 +77,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
       a = &gf[z][mu];
       _su3d_times_su3(tmp1, *d, *a);
       e = &gf[y][nu];
+      double const ptbc_fac2 = get_ptbc_coeff(z, nu) * get_ptbc_coeff(z, mu) * get_ptbc_coeff(y, nu);
       _su3d_times_su3(tmp2, *e, tmp1);
 
       y = g_iup[z][mu];
       b = &gf[y][nu];
       z = g_iup[y][nu];
       c = &gf[z][nu];
+      double const ptbc_fac3 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, nu);
       _su3_times_su3(tmp1, *b, *c);
-      _su3_times_su3_acc((*v), tmp2, tmp1);
+      //_su3_times_su3_acc((*v), tmp2, tmp1);
+      _real_times_su3_times_su3_acc((*v), tmp2, tmp1, ptbc_fac2 * ptbc_fac3);
 
       /* second contr. starting from x
        * a b c e^+ d^+
@@ -96,14 +103,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
       _su3_times_su3(tmp1, *a, *b);
       z = g_iup[y][mu];
       c = &gf[z][mu];
+      double const ptbc_fac4 = get_ptbc_coeff(x, nu) * get_ptbc_coeff(y, mu) * get_ptbc_coeff(z, mu);
       _su3_times_su3(tmp2, tmp1, *c);
 
       y = g_iup[x][mu];
       d = &gf[y][mu];
       z = g_iup[y][mu];
       e = &gf[z][nu];
+      double const ptbc_fac5 = get_ptbc_coeff(y, mu) * get_ptbc_coeff(z, nu);
       _su3_times_su3(tmp1, *d, *e);
-      _su3_times_su3d_acc((*v), tmp2, tmp1);
+      //_su3_times_su3d_acc((*v), tmp2, tmp1);
+      _real_times_su3_times_su3d_acc((*v), tmp2, tmp1, ptbc_fac4 * ptbc_fac5);
 
       /* 1 contr. starting idn[x][nu]
        * d^+ a b c e^+
@@ -119,14 +129,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
       _su3d_times_su3(tmp1, *d, *a);
       z = g_iup[y][mu];
       b = &gf[z][mu];
+      double const ptbc_fac6 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(y, mu) * get_ptbc_coeff(z, mu);
       _su3_times_su3(tmp2, tmp1, *b);
 
       y = g_iup[z][mu];
       c = &gf[y][nu];
       z = g_iup[x][mu];
       e = &gf[z][mu];
+      double const ptbc_fac7 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, mu);
       _su3_times_su3d(tmp1, *c, *e);
-      _su3_times_su3_acc((*v), tmp2, tmp1);
+      //_su3_times_su3_acc((*v), tmp2, tmp1);
+      _real_times_su3_times_su3_acc((*v), tmp2, tmp1, ptbc_fac6 * ptbc_fac7);
 
       /* 1 contr. starting idn[idn[x][mu]][nu]
        *  e^+ d^+ a b c
@@ -142,14 +155,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
       a = &gf[z][mu];
       _su3d_times_su3(tmp1, *d, *a);
       e = &gf[y][mu];
+      double const ptbc_fac8 = get_ptbc_coeff(z, nu) * get_ptbc_coeff(z, mu) * get_ptbc_coeff(y, mu);
       _su3d_times_su3(tmp2, *e, tmp1);
 
       y = g_idn[x][nu];
       b = &gf[y][mu];
       z = g_iup[y][mu];
       c = &gf[z][nu];
+      double const ptbc_fac9 = get_ptbc_coeff(y, mu) * get_ptbc_coeff(z, nu);
       _su3_times_su3(tmp1, *b, *c);
-      _su3_times_su3_acc((*v), tmp2, tmp1);
+      //_su3_times_su3_acc((*v), tmp2, tmp1);
+      _real_times_su3_times_su3_acc((*v), tmp2, tmp1, ptbc_fac8 * ptbc_fac9);
 
       /* 1 contr. starting idn[x][mu]
        * d^+ a b c e^+
@@ -165,14 +181,17 @@ void get_rectangle_staples_general(su3 *const v, const int x, const int mu,
       a = &gf[y][nu];
       _su3d_times_su3(tmp1, *d, *a);
       b = &gf[z][mu];
+      double const ptbc_fac10 = get_ptbc_coeff(y, mu) * get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, mu);
       _su3_times_su3(tmp2, tmp1, *b);
 
       y = g_iup[x][mu];
       e = &gf[y][nu];
       z = g_iup[x][nu];
       c = &gf[z][mu];
+      double const ptbc_fac11 = get_ptbc_coeff(y, nu) * get_ptbc_coeff(z, mu);
       _su3_times_su3d(tmp1, *c, *e);
-      _su3_times_su3_acc((*v), tmp2, tmp1);
+      //_su3_times_su3_acc((*v), tmp2, tmp1);
+      _real_times_su3_times_su3_acc((*v), tmp2, tmp1, ptbc_fac10 * ptbc_fac11);
     }
   }
 }
diff --git a/src/lib/get_staples.c b/src/lib/get_staples.c
index ae7f19d09..0de554f56 100644
--- a/src/lib/get_staples.c
+++ b/src/lib/get_staples.c
@@ -29,6 +29,7 @@
 #include "start.h"
 #include "su3.h"
 #include "su3adj.h"
+#include "ptbc.h"
 
 void get_staples(su3* const staple, const int x, const int mu, const su3** in_gauge_field) {
   int iy;
@@ -41,20 +42,24 @@ void get_staples(su3* const staple, const int x, const int mu, const su3** in_ga
       w1 = &in_gauge_field[x][k];
       w2 = &in_gauge_field[g_iup[x][k]][mu];
       w3 = &in_gauge_field[g_iup[x][mu]][k];
+      double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
 
       /* st = w2 * w3^d */
       _su3_times_su3d(st, *w2, *w3);
       /* v = v + w1 * st */
-      _su3_times_su3_acc(*staple, *w1, st);
+      //_su3_times_su3_acc(*staple, *w1, st);
+      _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
 
       iy = g_idn[x][k];
       w1 = &in_gauge_field[iy][k];
       w2 = &in_gauge_field[iy][mu];
       w3 = &in_gauge_field[g_iup[iy][mu]][k];
+      double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
       /* st = w2 * w3 */
       _su3_times_su3(st, *w2, *w3);
       /* v = v + w1^d * st */
-      _su3d_times_su3_acc(*staple, *w1, st);
+      //_su3d_times_su3_acc(*staple, *w1, st);
+      _real_times_su3d_times_su3_acc(*staple, *w1, st, ptbc_fac1);
     }
   }
 }
@@ -71,20 +76,24 @@ void get_spacelike_staples(su3* const staple, const int x, const int mu,
       w1 = &in_gauge_field[x][k];
       w2 = &in_gauge_field[g_iup[x][k]][mu];
       w3 = &in_gauge_field[g_iup[x][mu]][k];
+       double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
 
       /* st = w2 * w3^d */
       _su3_times_su3d(st, *w2, *w3);
       /* v = v + w1 * st */
-      _su3_times_su3_acc(*staple, *w1, st);
+      //_su3_times_su3_acc(*staple, *w1, st);
+      _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
 
       iy = g_idn[x][k];
       w1 = &in_gauge_field[iy][k];
       w2 = &in_gauge_field[iy][mu];
       w3 = &in_gauge_field[g_iup[iy][mu]][k];
+      double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
       /* st = w2 * w3 */
       _su3_times_su3(st, *w2, *w3);
       /* v = v + w1^d * st */
-      _su3d_times_su3_acc(*staple, *w1, st);
+      //_su3d_times_su3_acc(*staple, *w1, st);
+      _real_times_su3d_times_su3_acc(*staple, *w1, st, ptbc_fac1);
     }
   }
 }
@@ -101,19 +110,23 @@ void get_timelike_staples(su3* const staple, const int x, const int mu,
     w1 = &in_gauge_field[x][k];
     w2 = &in_gauge_field[g_iup[x][k]][mu];
     w3 = &in_gauge_field[g_iup[x][mu]][k];
+    double const ptbc_fac0 = get_ptbc_coeff(x, k) * get_ptbc_coeff(g_iup[x][k], mu) * get_ptbc_coeff(g_iup[x][mu], k);
 
     /* st = w2 * w3^d */
     _su3_times_su3d(st, *w2, *w3);
     /* v = v + w1 * st */
-    _su3_times_su3_acc(*staple, *w1, st);
+    //_su3_times_su3_acc(*staple, *w1, st);
+    _real_times_su3_times_su3_acc(*staple, *w1, st, ptbc_fac0);
 
     iy = g_idn[x][k];
     w1 = &in_gauge_field[iy][k];
     w2 = &in_gauge_field[iy][mu];
     w3 = &in_gauge_field[g_iup[iy][mu]][k];
+    double const ptbc_fac1 = get_ptbc_coeff(iy, k) * get_ptbc_coeff(iy, mu) * get_ptbc_coeff(g_iup[iy][mu], k);
     /* st = w2 * w3 */
     _su3_times_su3(st, *w2, *w3);
     /* v = v + w1^d * st */
-    _su3d_times_su3_acc(*staple, *w1, st);
+    //_su3d_times_su3_acc(*staple, *w1, st);
+    _real_times_su3d_times_su3_acc(*staple, *w1, st, ptbc_fac1);
   }
 }
diff --git a/src/lib/measure_gauge_action.c b/src/lib/measure_gauge_action.c
index 654990166..f503b45f5 100644
--- a/src/lib/measure_gauge_action.c
+++ b/src/lib/measure_gauge_action.c
@@ -39,6 +39,7 @@
 #include "measure_gauge_action.h"
 #include "su3.h"
 #include "su3adj.h"
+#include "ptbc.h"
 
 double measure_plaquette(const su3 *const *const gf) {
   static double res;
@@ -73,8 +74,11 @@ double measure_plaquette(const su3 *const *const gf) {
           v = &gf[ix][mu2];
           w = &gf[ix2][mu1];
           _su3_times_su3(pr2, *v, *w);
+          // compute local parallel tempering factor: fac_1 * fac_2 * fac_3 *  fac_4
+          double const ptbc_fac = get_ptbc_coeff(ix, mu1) * get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, mu1);
+          //double const ptbc_fac = 1.;
           _trace_su3_times_su3d(ac, pr1, pr2);
-          tr = ac + kc;
+          tr = ac*ptbc_fac + kc;
           ts = tr + ks;
           tt = ts - ks;
           ks = ts;
@@ -135,9 +139,12 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
         v = &gf[ix][mu2];
         w = &gf[ix2][0];
         _su3_times_su3(pr2, *v, *w);
+        // parallel tempering factor = fac_1 * fac_2 * fac_3 *  fac_4 computed in 2 parts
+        double const ptbc_fac = get_ptbc_coeff(ix, 0) * get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, 0);
         _trace_su3_times_su3d(ac, pr1, pr2);
         ac *= (1 + lambda);
-        tr = ac + kc;
+        tr = ac*ptbc_fac + kc;
+        //tr = ac + kc;
         ts = tr + ks;
         tt = ts - ks;
         ks = ts;
@@ -146,6 +153,7 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
       // magnetic part
       for (int mu1 = 1; mu1 < 3; mu1++) {
         ix1 = g_iup[ix][mu1];
+        double const ptbc_fac1 = get_ptbc_coeff(ix, mu1);
         for (int mu2 = mu1 + 1; mu2 < 4; mu2++) {
           ix2 = g_iup[ix][mu2];
           v = &gf[ix][mu1];
@@ -154,9 +162,11 @@ double measure_gauge_action(const su3 *const *const gf, const double lambda) {
           v = &gf[ix][mu2];
           w = &gf[ix2][mu1];
           _su3_times_su3(pr2, *v, *w);
+          double const ptbc_fac2 = get_ptbc_coeff(ix1, mu2) * get_ptbc_coeff(ix, mu2) * get_ptbc_coeff(ix2, mu1);
           _trace_su3_times_su3d(ac, pr1, pr2);
-          ac *= (1 - lambda);
-          tr = ac + kc;
+          //ac *= (1 - lambda);
+          tr = ac*(1 - lambda)*(ptbc_fac2*ptbc_fac1) + kc;
+          //tr = ac + kc;
           ts = tr + ks;
           tt = ts - ks;
           ks = ts;
diff --git a/src/lib/measure_rectangles.c b/src/lib/measure_rectangles.c
index 4f0ca07fd..ba8636e7d 100644
--- a/src/lib/measure_rectangles.c
+++ b/src/lib/measure_rectangles.c
@@ -45,6 +45,7 @@
 #include "measure_rectangles.h"
 #include "su3.h"
 #include "su3adj.h"
+#include "ptbc.h"
 
 double measure_rectangles(const su3 **const gf) {
   static double res;
@@ -86,6 +87,7 @@ double measure_rectangles(const su3 **const gf) {
             _su3_times_su3(tmp, *v, *w);
             v = &gf[k][nu];
             _su3_times_su3(pr1, tmp, *v);
+            double const ptbc_fac0 = get_ptbc_coeff(i, mu) * get_ptbc_coeff(j, nu) * get_ptbc_coeff(k, nu);
             /*
               ->
               ^
@@ -100,12 +102,14 @@ double measure_rectangles(const su3 **const gf) {
             _su3_times_su3(tmp, *v, *w);
             v = &gf[k][mu];
             _su3_times_su3(pr2, tmp, *v);
+            double const ptbc_fac1 = get_ptbc_coeff(i, nu) * get_ptbc_coeff(j, nu) * get_ptbc_coeff(k, mu);
 
             /* Trace it */
             _trace_su3_times_su3d(ac, pr1, pr2);
             /* 	  printf("i mu nu: %d %d %d, ac = %e\n", i, mu, nu, ac); */
             /* Kahan summation */
-            tr = ac + kc;
+            tr = ac*(ptbc_fac0 * ptbc_fac1) + kc;
+            //tr = ac + kc;
             ts = tr + ks;
             tt = ts - ks;
             ks = ts;
diff --git a/src/lib/su3.h b/src/lib/su3.h
index 4dd851811..2e53b2210 100644
--- a/src/lib/su3.h
+++ b/src/lib/su3.h
@@ -490,6 +490,17 @@ typedef struct {
   (u).c21 += (v).c20 * (w).c01 + (v).c21 * (w).c11 + (v).c22 * (w).c21; \
   (u).c22 += (v).c20 * (w).c02 + (v).c21 * (w).c12 + (v).c22 * (w).c22;
 
+#define _real_times_su3_times_su3_acc(u, v, w, r)                       \
+  (u).c00 += (r) * ((v).c00 * (w).c00 + (v).c01 * (w).c10 + (v).c02 * (w).c20); \
+  (u).c01 += (r) * ((v).c00 * (w).c01 + (v).c01 * (w).c11 + (v).c02 * (w).c21); \
+  (u).c02 += (r) * ((v).c00 * (w).c02 + (v).c01 * (w).c12 + (v).c02 * (w).c22); \
+  (u).c10 += (r) * ((v).c10 * (w).c00 + (v).c11 * (w).c10 + (v).c12 * (w).c20); \
+  (u).c11 += (r) * ((v).c10 * (w).c01 + (v).c11 * (w).c11 + (v).c12 * (w).c21); \
+  (u).c12 += (r) * ((v).c10 * (w).c02 + (v).c11 * (w).c12 + (v).c12 * (w).c22); \
+  (u).c20 += (r) * ((v).c20 * (w).c00 + (v).c21 * (w).c10 + (v).c22 * (w).c20); \
+  (u).c21 += (r) * ((v).c20 * (w).c01 + (v).c21 * (w).c11 + (v).c22 * (w).c21); \
+  (u).c22 += (r) * ((v).c20 * (w).c02 + (v).c21 * (w).c12 + (v).c22 * (w).c22);
+
 #define _su3_times_su3d(u, v, w)                                                         \
   (u).c00 = (v).c00 * conj((w).c00) + (v).c01 * conj((w).c01) + (v).c02 * conj((w).c02); \
   (u).c01 = (v).c00 * conj((w).c10) + (v).c01 * conj((w).c11) + (v).c02 * conj((w).c12); \
@@ -512,6 +523,17 @@ typedef struct {
   (u).c21 += (v).c20 * conj((w).c10) + (v).c21 * conj((w).c11) + (v).c22 * conj((w).c12); \
   (u).c22 += (v).c20 * conj((w).c20) + (v).c21 * conj((w).c21) + (v).c22 * conj((w).c22);
 
+#define _real_times_su3_times_su3d_acc(u, v, w, r)                                                      \
+  (u).c00 += (r) * ((v).c00 * conj((w).c00) + (v).c01 * conj((w).c01) + (v).c02 * conj((w).c02)); \
+  (u).c01 += (r) * ((v).c00 * conj((w).c10) + (v).c01 * conj((w).c11) + (v).c02 * conj((w).c12)); \
+  (u).c02 += (r) * ((v).c00 * conj((w).c20) + (v).c01 * conj((w).c21) + (v).c02 * conj((w).c22)); \
+  (u).c10 += (r) * ((v).c10 * conj((w).c00) + (v).c11 * conj((w).c01) + (v).c12 * conj((w).c02)); \
+  (u).c11 += (r) * ((v).c10 * conj((w).c10) + (v).c11 * conj((w).c11) + (v).c12 * conj((w).c12)); \
+  (u).c12 += (r) * ((v).c10 * conj((w).c20) + (v).c11 * conj((w).c21) + (v).c12 * conj((w).c22)); \
+  (u).c20 += (r) * ((v).c20 * conj((w).c00) + (v).c21 * conj((w).c01) + (v).c22 * conj((w).c02)); \
+  (u).c21 += (r) * ((v).c20 * conj((w).c10) + (v).c21 * conj((w).c11) + (v).c22 * conj((w).c12)); \
+  (u).c22 += (r) * ((v).c20 * conj((w).c20) + (v).c21 * conj((w).c21) + (v).c22 * conj((w).c22));
+
 #define _su3d_times_su3(u, v, w)                                                         \
   (u).c00 = conj((v).c00) * (w).c00 + conj((v).c10) * (w).c10 + conj((v).c20) * (w).c20; \
   (u).c01 = conj((v).c00) * (w).c01 + conj((v).c10) * (w).c11 + conj((v).c20) * (w).c21; \
@@ -534,6 +556,17 @@ typedef struct {
   (u).c21 += conj((v).c02) * (w).c01 + conj((v).c12) * (w).c11 + conj((v).c22) * (w).c21; \
   (u).c22 += conj((v).c02) * (w).c02 + conj((v).c12) * (w).c12 + conj((v).c22) * (w).c22;
 
+#define _real_times_su3d_times_su3_acc(u, v, w, r)                                                      \
+  (u).c00 += (r) * (conj((v).c00) * (w).c00 + conj((v).c10) * (w).c10 + conj((v).c20) * (w).c20); \
+  (u).c01 += (r) * (conj((v).c00) * (w).c01 + conj((v).c10) * (w).c11 + conj((v).c20) * (w).c21); \
+  (u).c02 += (r) * (conj((v).c00) * (w).c02 + conj((v).c10) * (w).c12 + conj((v).c20) * (w).c22); \
+  (u).c10 += (r) * (conj((v).c01) * (w).c00 + conj((v).c11) * (w).c10 + conj((v).c21) * (w).c20); \
+  (u).c11 += (r) * (conj((v).c01) * (w).c01 + conj((v).c11) * (w).c11 + conj((v).c21) * (w).c21); \
+  (u).c12 += (r) * (conj((v).c01) * (w).c02 + conj((v).c11) * (w).c12 + conj((v).c21) * (w).c22); \
+  (u).c20 += (r) * (conj((v).c02) * (w).c00 + conj((v).c12) * (w).c10 + conj((v).c22) * (w).c20); \
+  (u).c21 += (r) * (conj((v).c02) * (w).c01 + conj((v).c12) * (w).c11 + conj((v).c22) * (w).c21); \
+  (u).c22 += (r) * (conj((v).c02) * (w).c02 + conj((v).c12) * (w).c12 + conj((v).c22) * (w).c22);
+
 #define _su3_minus_const_times_im_trace_su3(w, c, v)                     \
   (w).c00 -= I * c * (cimag((v).c00) + cimag((v).c11) + cimag((v).c22)); \
   (w).c11 -= I * c * (cimag((v).c00) + cimag((v).c11) + cimag((v).c22)); \

From 4f33ca572460c8a4261817b215b1726fac04f41d Mon Sep 17 00:00:00 2001
From: JingJing Li <lijingjing416@gmail.com>
Date: Fri, 8 May 2026 10:09:37 +0200
Subject: [PATCH 4/8] rebase again +  defect check

---
 src/lib/ptbc/defect.c | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/lib/ptbc/defect.c b/src/lib/ptbc/defect.c
index 176dd7241..121f1f57d 100644
--- a/src/lib/ptbc/defect.c
+++ b/src/lib/ptbc/defect.c
@@ -36,7 +36,7 @@ bool is_defect(PTBCDefect *def, int const ix, int const mu) {
         return false;
       }
     }
-    if (g_proc_id==0) printf("defect coordinate %d %d %d %d mu = %d\n", g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3], mu);
+    //if (g_proc_id==0) printf("defect coordinate %d %d %d %d mu = %d\n", g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3], mu);
     return true;
   }
   else {
@@ -48,6 +48,9 @@ bool is_defect(PTBCDefect *def, int const ix, int const mu) {
 double get_ptbc_coeff(int const ix, int const mu) {
   int const inst = app()->ptbc.instance_id;
   const PTBCInstance *instance = &(app()->ptbc.instances[inst]);
+
+  // if instance is not active, return 1
+  if (!instance->active) return 1.;
   
   // loop over defects
   for (int i=0; i<instance->n_coeffs; i++) {

From adb5bc2db141007e2ddbb5926313a6cb1ce51cd0 Mon Sep 17 00:00:00 2001
From: JingJing Li <lijingjing416@gmail.com>
Date: Thu, 21 May 2026 17:18:11 +0200
Subject: [PATCH 5/8] changes to defect region. include surface links

---
 src/lib/include/app.h |  6 ++++--
 src/lib/ptbc/defect.c | 49 +++++++++++++++++++++++++++++++++++++++----
 src/lib/read_input.l  | 32 ++++++++++------------------
 3 files changed, 60 insertions(+), 27 deletions(-)

diff --git a/src/lib/include/app.h b/src/lib/include/app.h
index 169f48f22..97ccac83a 100644
--- a/src/lib/include/app.h
+++ b/src/lib/include/app.h
@@ -67,8 +67,10 @@ typedef struct {
 
 typedef struct {
     bool active;        // Whether the defect is active or not
-    int Ld[3];          // Extents of the defect
-    direction_t along;  // Along which dimension
+    //int Ld[3];          // Extents of the defect
+    //direction_t along;  // Along which dimension
+    int pos[4];          // Position of the defect
+    int Ld[4];           // Extent of the defect
 } PTBCDefect;
 
 
diff --git a/src/lib/ptbc/defect.c b/src/lib/ptbc/defect.c
index 121f1f57d..d08d7abeb 100644
--- a/src/lib/ptbc/defect.c
+++ b/src/lib/ptbc/defect.c
@@ -10,10 +10,52 @@
 #include "global.h"
 
 
-// check if a link lie in defect region
-// 1. last slice of along 2. 3d cube with 0 origin
+// check if a link lie in defect region, cuts at pos + 1/2 and pos + Ld + 1/2
+// if either or end of the link is in the defect region, return true 
+// include both case when link is in the defect and case link is crossing defect boundary
 bool is_defect(PTBCDefect *def, int const ix, int const mu) {
   if (ix >= VOLUME) return false;
+  int const global_dim[4] = {T*g_nproc_t, LX*g_nproc_t, LY*g_nproc_x, LZ*g_nproc_z};
+  int* coords = g_coord[ix];
+
+  // find three non-mu directions
+  int dim3[3];
+  int count = 0;
+  for (int d=0; d<4; d++) {
+    if (d != mu) {
+      dim3[count] = d;
+      count++;
+    }
+  }
+  // check if mu direction is in the defect or crossing the cut
+  if (coords[mu] >= def->pos[mu] && coords[mu] + 1 <= def->pos[mu] + def->Ld[mu] + 1) {
+    // find three non-mu directions
+    int dim3[3];
+    int count = 0;
+    for (int d=0; d<4; d++) {
+      if (d != mu) {
+        dim3[count] = d;
+        count++;
+      }
+    }
+    // check if the other three directions are within the defect region
+    for (int d=0; d<3; d++) {
+      if (coords[dim3[d]] > def->pos[dim3[d]] && coords[dim3[d]] <= def->pos[dim3[d]] + def->Ld[dim3[d]]){
+        continue;
+      }
+      else{
+        return false;
+      }
+    }
+    return true;
+  }
+  else
+    return false;  
+}
+
+
+/* bool is_defect(PTBCDefect *def, int const ix, int const mu) {
+  if (ix >= VOLUME) return false;
   
   int const global_dim[4] = {T*g_nproc_t, LX*g_nproc_t, LY*g_nproc_x, LZ*g_nproc_z};
 
@@ -36,13 +78,12 @@ bool is_defect(PTBCDefect *def, int const ix, int const mu) {
         return false;
       }
     }
-    //if (g_proc_id==0) printf("defect coordinate %d %d %d %d mu = %d\n", g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3], mu);
     return true;
   }
   else {
     return false;
   }
-}
+} */
 
 // get multiplying factor of parallel tempering locally (assumne no overlapping defects!)
 double get_ptbc_coeff(int const ix, int const mu) {
diff --git a/src/lib/read_input.l b/src/lib/read_input.l
index 8b2f068eb..f2ea925ef 100644
--- a/src/lib/read_input.l
+++ b/src/lib/read_input.l
@@ -3182,29 +3182,19 @@ static inline double fltlist_next_token(int * const list_end){
       PTBCDefect* cd = appm()->ptbc.defects + b;
 
       int n_parameters = strlist_count_tokens(yytext);
-      fail(n_parameters != 3, "Coefficients must provide exactly 3 numbers, not %d!\n", n_parameters);
-      sscanf(yytext, " %[a-zA-Z] = %d, %d, %d", name, cd->Ld, cd->Ld+1, cd->Ld+2);
-      if(myverbose) printf("    Ld = [%d, %d, %d] line %d\n", cd->Ld[0], cd->Ld[1], cd->Ld[2], line_of_file);
+      fail(n_parameters != 4, "Coefficients must provide exactly 4 numbers, not %d!\n", n_parameters);
+      sscanf(yytext, " %[a-zA-Z] = %d, %d, %d, %d", name, cd->Ld, cd->Ld+1, cd->Ld+2, cd->Ld+3);
+      if(myverbose) printf("    Ld = [%d, %d, %d, %d] line %d\n", cd->Ld[0], cd->Ld[1], cd->Ld[2], cd->Ld[3], line_of_file);
       BEGIN(PTBCDEFECT);
     }
-    {SPC}*Along{EQL}T {
-      appm()->ptbc.defects[b].along = DIRECTION_T;
-      if(myverbose) printf("    PTBC defect along T direction line %d\n", line_of_file);
-      BEGIN(PTBCDEFECT);
-    }
-    {SPC}*Along{EQL}X {
-      appm()->ptbc.defects[b].along = DIRECTION_X;
-      if(myverbose) printf("    PTBC defect along X direction line %d\n", line_of_file);
-      BEGIN(PTBCDEFECT);
-    }
-    {SPC}*Along{EQL}Y {
-      appm()->ptbc.defects[b].along = DIRECTION_Y;
-      if(myverbose) printf("    PTBC defect along Y direction line %d\n", line_of_file);
-      BEGIN(PTBCDEFECT);
-    }
-    {SPC}*Along{EQL}Z {
-      appm()->ptbc.defects[b].along = DIRECTION_Z;
-      if(myverbose) printf("    PTBC defect along Z direction line %d\n", line_of_file);
+    {SPC}*Position{EQL}{STRLIST} {
+      
+      PTBCDefect* cd = appm()->ptbc.defects + b;
+
+      int n_parameters = strlist_count_tokens(yytext);
+      fail(n_parameters != 4, "Position must provide exactly 4 numbers, not %d!\n", n_parameters);
+      sscanf(yytext, " %[a-zA-Z] = %d, %d, %d, %d", name, cd->pos, cd->pos+1, cd->pos+2, cd->pos+3);
+      if(myverbose) printf("    pos = [%d, %d, %d, %d] line %d\n", cd->pos[0], cd->pos[1], cd->pos[2], cd->pos[3], line_of_file);
       BEGIN(PTBCDEFECT);
     }
     {SPC}*EndPTBCDefect{SPC}* {

From 6fe7bbd974dd8aaa41013f6d48c0726deda307ae Mon Sep 17 00:00:00 2001
From: JingJing Li <lijingjing416@gmail.com>
Date: Wed, 3 Jun 2026 19:29:48 +0200
Subject: [PATCH 6/8] ptbc topology initialiser & utility functions

---
 ...-ptbc-hmc-rgmixedcg-tmcloverdetratio.input |  36 +-
 src/bin/hmc_tm.c                              |   8 +
 src/lib/app_context/app_context.c             |  10 +
 src/lib/include/ptbc.h                        |  43 +-
 src/lib/ptbc/defect.c                         | 399 +++++++++++++++---
 src/lib/start.c                               |   5 +
 6 files changed, 445 insertions(+), 56 deletions(-)

diff --git a/doc/sample-input/sample-ptbc-hmc-rgmixedcg-tmcloverdetratio.input b/doc/sample-input/sample-ptbc-hmc-rgmixedcg-tmcloverdetratio.input
index 89ff0511f..55f963189 100644
--- a/doc/sample-input/sample-ptbc-hmc-rgmixedcg-tmcloverdetratio.input
+++ b/doc/sample-input/sample-ptbc-hmc-rgmixedcg-tmcloverdetratio.input
@@ -16,16 +16,35 @@ ReversibilityCheckIntervall = 4
 InitialStoreCounter = 0
 DebugLevel = 1
 
+
+# complex PTBC topo example, 2 periodic, 2 tentacles
+# 2 periodic instances connected by dummy defect 2
+
 BeginPTBC
 
   BeginPTBCDefect 0
-    Extents = 2, 2, 2
-    Along = Y
+    Extents = 2, 2, 2, 2
+    Position = 1, 1, 1, 1
+  EndPTBCDefect
+
+  BeginPTBCDefect 1
+    Extents = 0, 2, 2, 2
+    Position = 4, 4, 4, 4
+  EndPTBCDefect
+
+  BeginPTBCDefect 2
+    Extents = 0, 0, 0, 1
+    Position = 0, 0, 3, 0
+  EndPTBCDefect
+
+  BeginPTBCDefect 3
+    Extents = 0, 1, 1, 1
+    Position = 6, 6, 6, 6
   EndPTBCDefect
 
   BeginPTBCInstance 0
-    Coefficients = 1.0
-    Defects = 0
+    Coefficients = 1.0, 1.0
+    Defects = 0, 2
   EndPTBCInstance
 
   BeginPTBCInstance 1
@@ -34,19 +53,20 @@ BeginPTBC
   EndPTBCInstance
 
   BeginPTBCInstance 2
-    Coefficients = 0.75
-    Defects = 0
+    Coefficients = 1.0, 1.0, 1.0
+    Defects = 1, 2, 3
   EndPTBCInstance
 
   BeginPTBCInstance 3
-    Coefficients = 0.5
-    Defects = 0
+    Coefficients = 0.1, 0.0
+    Defects = 1, 3
   EndPTBCInstance
 
 EndPTBC
 
 
 
+
 # since this is a test file, we employ the reproducible random numbers mode
 ReproduceRandomNumbers = yes
 Seed = 127782
diff --git a/src/bin/hmc_tm.c b/src/bin/hmc_tm.c
index 6e3654c2c..2d2b8165d 100644
--- a/src/bin/hmc_tm.c
+++ b/src/bin/hmc_tm.c
@@ -73,6 +73,7 @@
 #ifdef TM_USE_QUDA
 #include "quda_interface.h"
 #endif
+#include "ptbc.h"
 
 extern int nstore;
 
@@ -142,6 +143,13 @@ int main(int argc, char *argv[]) {
   NO_OF_SPINORFIELDS_32 = 6;
 
   tmlqcd_mpi_init(argc, argv);
+
+  // initialise ptbc
+  if (app()->ptbc.active) {
+    init_ptbc_tree();
+    if (g_proc_id == 0) print_ptbc_topo();
+  }
+
   tm_stopwatch_push(&g_timers, "HMC", "");
 
   if (nstore == -1) {
diff --git a/src/lib/app_context/app_context.c b/src/lib/app_context/app_context.c
index 0d2e514fa..78ad3495c 100644
--- a/src/lib/app_context/app_context.c
+++ b/src/lib/app_context/app_context.c
@@ -215,6 +215,16 @@ static void initialize(void)
         app_instance.ptbc.instance_id, app_instance.ptbc.n_instances,
         instance_rank, instance_size);
 
+    
+    // verify that local ranks are organised regularly
+    int expected_world_rank = (app_instance.ptbc.instance_id * instance_size) + instance_rank;
+
+    // If this assertion fails, the user launched the MPI program with a cyclic layout
+    if (app_instance.mpi.world_rank != expected_world_rank) {
+        fprintf(stderr, "Error: MPI rank topology is not block-contiguous! "
+                        "Please check your mpirun/srun pinning settings.\n");
+        MPI_Abort(app_instance.mpi.world_comm, 1);
+    }
 
     //err(true, "bailing out");
     /*if (app_instance.ptbc.instance_id != 0) {
diff --git a/src/lib/include/ptbc.h b/src/lib/include/ptbc.h
index 42b434cf6..9652693cb 100644
--- a/src/lib/include/ptbc.h
+++ b/src/lib/include/ptbc.h
@@ -1,6 +1,47 @@
 
 #include <mpi.h>
+#include <app.h>
 
 
+/* for computing action with PTBC */
 bool is_defect(PTBCDefect *def, int const ix, int const mu);
-double get_ptbc_coeff(int const ix, int const mu);
\ No newline at end of file
+double get_ptbc_coeff(int const ix, int const mu);
+
+
+/* for swapping rng */
+typedef struct{
+    int state_s[105];   // state of rlxs
+    int state_d[105];   // state of rlxd
+} SwapRNG;
+
+void swap_rng(int const dest_inst);
+
+/* Graph of instances */
+typedef struct {
+  int parent; // parent node
+  int children[MAX_N_DEFECTS]; // children nodes
+  int n_children; // number of children nodes  
+} Node;
+
+typedef struct {
+    int root; // root node of the graph 
+} Tree;
+
+// setters
+void set_tree_root(int const root);
+void set_node_parent(int const node_id, int const parent_id);
+void add_node_child(int const node_id, int const child_id);
+void set_edge(int const parent_id, int const child_id);
+void init_node(int const node_id);
+
+
+// getters
+int const get_tree_root();
+int const get_node_parent(int const node_id);
+int const* get_node_children(int const node_id);
+int const get_node_n_children(int const node_id);
+
+
+// initialiser
+void init_ptbc_tree();
+void print_ptbc_topo();
\ No newline at end of file
diff --git a/src/lib/ptbc/defect.c b/src/lib/ptbc/defect.c
index d08d7abeb..6a8d8c8e4 100644
--- a/src/lib/ptbc/defect.c
+++ b/src/lib/ptbc/defect.c
@@ -8,27 +8,55 @@
 #include <mpi.h>
 #include <ptbc.h>
 #include "global.h"
+#include "ranlxs.h"
+#include "ranlxd.h"
 
+MPI_Datatype swap_info_type= MPI_DATATYPE_NULL; // MPI datatype for swapping info
+static Tree tree;
+static Node nodes[MAX_N_INSTANCES]; // each node corresponds to an instance
 
-// check if a link lie in defect region, cuts at pos + 1/2 and pos + Ld + 1/2
-// if either or end of the link is in the defect region, return true 
-// include both case when link is in the defect and case link is crossing defect boundary
+#define err(test, ...) err_impl(test, __func__, __FILE__, __LINE__, __VA_ARGS__)
+static void err_impl(const bool test, const char* func, const char* file, const int line, const char* format, ...)
+{
+    if (test) {
+        va_list args;
+        char message[1024];
+        va_start(args, format);
+        vsnprintf(message, 1024, format, args);
+        va_end(args);
+        char location[1024];
+        snprintf(location, 1024, "%s:%d %s", file, line, func);
+        fatal_error(message, location);
+    }
+}
+
+// find distance end-start corrected for peridic bc 
+int const static dist(int start, int end, int const g_length) {
+  int const d = end - start;
+  int const dist = (d>=0) ? d : d + g_length; // assuming periodic bc
+
+  return dist;
+}
+
+/**
+ * @brief      check if a link lie in defect region, cuts at pos + 1/2 and pos + Ld + 1/2
+ *             if either or end of the link is in the defect region, return true 
+ *             include both case when link is in the defect and case link is crossing defect boundary
+ *
+ * @param      def  defect region
+ * @param      ix   link start point
+ * @param      mu   link direction
+ */
 bool is_defect(PTBCDefect *def, int const ix, int const mu) {
   if (ix >= VOLUME) return false;
   int const global_dim[4] = {T*g_nproc_t, LX*g_nproc_t, LY*g_nproc_x, LZ*g_nproc_z};
   int* coords = g_coord[ix];
 
-  // find three non-mu directions
-  int dim3[3];
-  int count = 0;
-  for (int d=0; d<4; d++) {
-    if (d != mu) {
-      dim3[count] = d;
-      count++;
-    }
-  }
   // check if mu direction is in the defect or crossing the cut
-  if (coords[mu] >= def->pos[mu] && coords[mu] + 1 <= def->pos[mu] + def->Ld[mu] + 1) {
+  // i.e. start >= pos[mu] and end <= pos[mu] + Ld[mu] + 1
+  // dist >= 0 or dist + 1 <= Ld + 1
+  int const dist_mu = dist(def->pos[mu], coords[mu], global_dim[mu]);
+  if (dist_mu >= 0 && dist_mu <= def->Ld[mu]) {
     // find three non-mu directions
     int dim3[3];
     int count = 0;
@@ -39,8 +67,11 @@ bool is_defect(PTBCDefect *def, int const ix, int const mu) {
       }
     }
     // check if the other three directions are within the defect region
+    // i.e. start at > pos end at <= pos + Ld
+    // dist > 0 or dist + 1 <= Ld
     for (int d=0; d<3; d++) {
-      if (coords[dim3[d]] > def->pos[dim3[d]] && coords[dim3[d]] <= def->pos[dim3[d]] + def->Ld[dim3[d]]){
+      int const dist_d = dist(def->pos[dim3[d]], coords[dim3[d]], global_dim[dim3[d]]);
+      if (dist_d > 0 && dist_d < def->Ld[dim3[d]]){
         continue;
       }
       else{
@@ -53,39 +84,12 @@ bool is_defect(PTBCDefect *def, int const ix, int const mu) {
     return false;  
 }
 
-
-/* bool is_defect(PTBCDefect *def, int const ix, int const mu) {
-  if (ix >= VOLUME) return false;
-  
-  int const global_dim[4] = {T*g_nproc_t, LX*g_nproc_t, LY*g_nproc_x, LZ*g_nproc_z};
-
-  // check start point is at second last slice and link stays on the slice
-  if (g_coord[ix][def->along] == global_dim[def->along] - 2 && mu != def->along){
-    // xf = end point of link 
-    int xf[4] = {g_coord[ix][0], g_coord[ix][1], g_coord[ix][2], g_coord[ix][3]};
-    xf[mu]  = xf[mu] + 1;
-    int dim3[3];
-    int count = 0;
-    for (int d=0; d<4; d++) {
-      if (d != def->along) {
-        dim3[count] = d;
-        count++;
-      }
-    }
-    
-    for (int d=0; d<3; d++) {
-      if (g_coord[ix][dim3[d]] >= def->Ld[d] || xf[dim3[d]] >= def->Ld[d]) {
-        return false;
-      }
-    }
-    return true;
-  }
-  else {
-    return false;
-  }
-} */
-
-// get multiplying factor of parallel tempering locally (assumne no overlapping defects!)
+/**
+ * @brief      get multiplying factor of parallel tempering locally (assumne no overlapping defects!)
+ *
+ * @param      ix   starting point of link
+ * @param      mu   direction of link
+ */
 double get_ptbc_coeff(int const ix, int const mu) {
   int const inst = app()->ptbc.instance_id;
   const PTBCInstance *instance = &(app()->ptbc.instances[inst]);
@@ -103,4 +107,305 @@ double get_ptbc_coeff(int const ix, int const mu) {
   }
 
   return 1.;
+}
+
+
+/**
+ * @brief      Swap RNG state of with another instance
+ *
+ * @param      dest_inst  Destination instance ID
+ */
+void swap_rng(int const dest_inst) {
+  SwapRNG info;
+  
+  if (dest_inst == app()->ptbc.instance_id) return;
+
+  // find dest rank and src rank
+  int n_ranks, local_rank;
+  MPI_Comm_rank(app()->mpi.comm, &local_rank);
+  MPI_Comm_size(app()->mpi.world_comm, &n_ranks);
+  int const dest_rank = n_ranks / app()->ptbc.n_instances * dest_inst + local_rank;  
+  int const src_rank = app()->mpi.world_rank;
+
+  // get states
+  rlxs_get(info.state_s);
+  rlxd_get(info.state_d);
+
+  // define MPI datatype for swapping info if not defined
+  if (swap_info_type == MPI_DATATYPE_NULL) {
+    PTBCInstance const *instance = &(app()->ptbc.instances[app()->ptbc.instance_id]);
+    int const blocklengths[3] = {105, 105, instance->n_coeffs};
+    MPI_Aint const displacements[3] = {0, 105 * sizeof(int), 2 * 105 * sizeof(int)};
+    MPI_Datatype const types[3] = {MPI_INT, MPI_INT, MPI_DOUBLE};
+    MPI_Type_create_struct(3, blocklengths, displacements, types, &swap_info_type);
+    MPI_Type_commit(&swap_info_type);
+  }
+
+  // swap rng and coeff info
+  MPI_Status status;
+  MPI_Sendrecv_replace(&info, 1, swap_info_type, dest_rank, 123, src_rank, 123, app()->mpi.world_comm, &status);
+  if (status.MPI_ERROR != MPI_SUCCESS) {
+    fprintf(stderr, "Error in MPI_Sendrecv_replace in swap_rng_coeff\n");
+    MPI_Abort(app()->mpi.world_comm, status.MPI_ERROR);
+  }
+
+  // set state
+  rlxs_reset(info.state_s);
+  rlxd_reset(info.state_d);
+
+  return;
+}
+
+typedef struct {
+  unsigned int n_coeffs;
+  double coeffs[MAX_N_INSTANCES];
+  unsigned int inst_id;
+} InstanceInfo;
+
+int compare_coeff(void const *ia, void const *ib) {
+  unsigned int count1 = 0;
+  unsigned int count2 = 0;
+
+  const InstanceInfo * a = (const InstanceInfo *) ia;
+  const InstanceInfo * b = (const InstanceInfo *) ib;
+  
+  for (unsigned int i=0; i<a->n_coeffs; i++) {
+    if (a->coeffs[i] < b->coeffs[i]) count1++;
+    else if (a->coeffs[i] > b->coeffs[i]) count2++;
+  }
+  if (count1 && count2==0) return 1;  // ib has bigger coeffs
+  else if (count2 && count1==0) return -1;  // ia has bigger coeffs
+  else if (count1==0 && count2==0) return 0; // two instances with exactly the same coeffs
+  
+  err((count1 && count2), "Error in compare_coeff: cannot sort instance coefficients in a strictly descending order!");
+  return 0;
+}
+
+static int if_periodic(PTBCInstance const *instance) {
+  for (int i=0; i<instance->n_coeffs; i++) {
+    if (instance->coefficients[i] != 1) {
+      return false;
+    }
+  }
+  return true;
+}
+
+
+void set_tree_root(int const root) {tree.root = root;};
+void set_node_parent(int const node_id, int const parent_id) {nodes[node_id].parent = parent_id;};
+void add_node_child(int const node_id, int const child_id) {
+    nodes[node_id].children[nodes[node_id].n_children] = child_id;
+    nodes[node_id].n_children++;
+};
+void set_edge(int const parent_id, int const child_id) {
+    set_node_parent(child_id, parent_id);
+    add_node_child(parent_id, child_id);
+}
+void init_node(int const node_id) {nodes[node_id].parent=-1; nodes[node_id].n_children=0;};
+
+int const get_tree_root() {return tree.root;};
+int const get_node_parent(int const node_id) {return nodes[node_id].parent;};
+int const* get_node_children(int const node_id) {return nodes[node_id].children;};
+int const get_node_n_children(int const node_id) {return nodes[node_id].n_children;};
+
+
+/**
+ * @brief      Initialise PTBC instance connection graph.
+ */
+void init_ptbc_tree() {
+  PTBCContext const *ptbc_ctx = &(app()->ptbc);
+  InstanceInfo info[MAX_N_DEFECTS][MAX_N_INSTANCES]; // instance info n_tentacles x instance per tentacle
+  int n_tentacles = 0;
+  int tentacle_type[MAX_N_DEFECTS][MAX_N_DEFECTS]; // store defect id for each tentacle
+  int tentacle_n_defect[MAX_N_DEFECTS]; // number of defects in each tentacle
+  int tentacle_length[MAX_N_DEFECTS] = {0}; // count how many instances in each tentacle
+  int n_periodic = 0; // number of periodic instances
+  int periodic_id[MAX_N_INSTANCES]; // store periodic instance id, may be multiple periodic instances
+
+  // initialise
+  for (int i=0; i<ptbc_ctx->n_instances; i++) {
+    init_node(i);
+  }
+
+  // log tentacles excluding periodic instances
+  for (int id=0; id<ptbc_ctx->n_instances; id++) {
+    PTBCInstance const * ins = &(ptbc_ctx->instances[id]);
+
+    if (!if_periodic(ins)) {
+      // check if defect already recorded, if not add to list, first exclude periodic instances
+      int flag_new_tentacle = 1;
+      for (int t=0; t<n_tentacles && flag_new_tentacle; t++) {
+        int defect_match = 0;
+        for (int nc=0; nc<ins->n_coeffs; nc++) {
+          if (ins->defects[nc] - ptbc_ctx->defects == tentacle_type[t][nc]) {
+            defect_match++;
+          }
+        }
+
+        // if match, record instance info in the tentacle
+        if (defect_match == ins->n_coeffs) {
+          flag_new_tentacle = 0;
+          info[t][tentacle_length[t]].n_coeffs = ins->n_coeffs;
+          info[t][tentacle_length[t]].inst_id = id;
+          for (int nc=0; nc<ins->n_coeffs; nc++) {
+            info[t][tentacle_length[t]].coeffs[nc] = ins->coefficients[nc];
+          }
+          tentacle_length[t]++;
+        }
+      }
+
+      // add new tentacle if no match
+      if (flag_new_tentacle) {
+        for (int nc=0; nc<ins->n_coeffs; nc++) {
+          tentacle_type[n_tentacles][nc] = ins->defects[nc] - ptbc_ctx->defects;
+        }
+        tentacle_n_defect[n_tentacles] = ins->n_coeffs;
+
+        // record first instance info in the new tentacle
+        info[n_tentacles][0].n_coeffs = ins->n_coeffs;
+        info[n_tentacles][0].inst_id = id;
+        for (int nc=0; nc<ins->n_coeffs; nc++) {
+          info[n_tentacles][0].coeffs[nc] = ins->coefficients[nc];
+        }
+        tentacle_length[n_tentacles]++;
+        n_tentacles++;
+      }
+    }
+    else {
+      periodic_id[n_periodic] = id;
+      n_periodic++;
+    }
+  }
+
+  // order instances in each tentacle by coeff 1->0
+  for (int t=0; t<n_tentacles; t++) {
+    qsort(info[t], tentacle_length[t], sizeof(InstanceInfo), compare_coeff);
+  }
+
+  // log tree structure
+  for (int t=0; t<n_tentacles; t++) {
+    // connect instances in the same tentacle, ordered by coeff 1->0
+    for (int n=0; n<tentacle_length[t]-1; n++) {
+      set_edge(info[t][n].inst_id, info[t][n+1].inst_id);
+    }
+  }
+
+  // set tree root
+  if (n_periodic) {
+    set_tree_root(periodic_id[0]);
+  }
+  else{
+    err(n_periodic, "Error in init_ptbc_tree: no periodic instances found!");
+  }
+
+  // Sort periodic instances
+  int connected_tentacles=0;
+  for (int p=0; p<n_periodic; p++) {
+    int periodic_edge_type[MAX_N_DEFECTS]; // store defect id for each periodic edge
+    PTBCInstance const *period_instance = &(ptbc_ctx->instances[periodic_id[p]]);
+    int n_periodic_edge=0;
+    int matches[MAX_N_DEFECTS] = {0}; // count how many matches per tentacle
+    
+    // attach tentacles to periodic instance(s)
+    for (int nc=0; nc<period_instance->n_coeffs; nc++) {
+      int flag_match = 0;
+      for (int t=0; t<n_tentacles; t++) {
+        for (int nc_t=0; nc_t<tentacle_n_defect[t]; nc_t++) {
+          if (period_instance->defects[nc] - ptbc_ctx->defects == tentacle_type[t][nc_t]) {
+            matches[t]++;
+            flag_match = 1;
+          }
+        }
+      }
+      if (!flag_match) {
+        // must be a periodic-periodic edge
+        periodic_edge_type[n_periodic_edge] = period_instance->defects[nc] - ptbc_ctx->defects;
+        n_periodic_edge++;
+      }
+    }
+
+    // connect periodic instance to tentacles
+    for (int t=0; t<n_tentacles; t++) {
+      if (matches[t] == tentacle_n_defect[t]) {
+        // this tentacle is connected to the periodic root
+        set_edge(periodic_id[p], info[t][0].inst_id);
+        connected_tentacles++;
+      }
+      else if (matches[t] > 0) {
+        err(matches[t], "Error in init_ptbc_tree: incomplete overlap of periodic instance with tentacles!");
+      }
+    }
+
+    // connect between periodic instances
+    // skip the last periodic as it has been checked already
+    if (p==n_periodic-1) continue;
+    int const cnt = n_periodic_edge;
+    // loop over other periodic instances to find match
+    for (int p2=p+1; p2<n_periodic; p2++) {
+      for (int type=0; type<cnt; type++) {
+        PTBCInstance const *instance_p2 = &(ptbc_ctx->instances[periodic_id[p2]]);
+        for (int nc=0; nc<instance_p2->n_coeffs; nc++) {
+          if (instance_p2->defects[nc] - ptbc_ctx->defects == periodic_edge_type[type]) {
+            // connect two periodic instances
+            set_edge(periodic_id[p], periodic_id[p2]);
+            n_periodic_edge--;
+          }
+        }
+      }
+    }
+    err(n_periodic_edge!=0, "Error in init_ptbc_tree: invalid connections between periodic instances!");
+  }
+
+  // err if tentacles connected to multiple periodic or loose tentacles exist
+  err(connected_tentacles>n_tentacles, "Error in init_ptbc_tree: a tentacle is connected to multiple periodic instances!");
+  err(connected_tentacles<n_tentacles, "Error in init_ptbc_tree: a tentacle is not connected to any periodic instance!");
+
+
+  // check for loose connections 
+  for (int id=0; id<ptbc_ctx->n_instances; id++) {
+    // a lost node that is not connected to anything
+    if (get_node_n_children(id) == 0 && get_node_parent(id) == -1) {
+      err(1, "Error in init_ptbc_tree: disconnected instance(s) excist(s)!");
+    }
+
+    // a non-root periodic node that is not connected to another periodic node
+    if (id != get_tree_root() && if_periodic(&(ptbc_ctx->instances[id])) && get_node_parent(id) == -1){
+      err(1, "Error in init_ptbc_tree: disconnected PTBC chains!");
+    }
+  }
+
+  if(g_proc_id == 0) printf("PTBC graph initialised: %d periodic instance(s), %d tentacle(s). \n", n_periodic, n_tentacles);
+  return;
+}
+
+
+void print_ptbc_topo() {
+  // Iterate through the nodes
+  const char* typeLabels[] = {
+    "open",         // index 0
+    "periodic",     // index 1
+    "intermediate"  // index 2
+  };
+
+  for (int id=0; id<app()->ptbc.n_instances; id++) {
+    int type;
+    if (if_periodic(&(app()->ptbc.instances[id]))) {
+      type = 1; // periodic
+    }
+    else if (get_node_n_children(id) == 0) {
+      type = 0; // open
+    }
+    else {
+      type = 2; // intermediate
+    }
+    
+    printf("Instance ID: %d\tType: %s\t", id, typeLabels[type]);
+    printf("Upstream ID: ");
+    if (get_node_parent(id) >= 0) printf("%d", get_node_parent(id)); // if parent is -1 the node is the root
+    printf("\tDownstream ID: ");
+    for (int c=0; c<get_node_n_children(id); c++) {
+      printf("%d, ", get_node_children(id)[c]);
+    }
+    printf("\n\n");
+  }
 }
\ No newline at end of file
diff --git a/src/lib/start.c b/src/lib/start.c
index 5504a83a2..1605d5ba0 100644
--- a/src/lib/start.c
+++ b/src/lib/start.c
@@ -831,6 +831,11 @@ void start_ranlux(int level, int seed) {
                       g_proc_coords[3];
 
   max_seed = 2147483647 / g_nproc;
+  
+  // offset for ptbc. If PTBC not active, instance_id=0, n_instances=1 ptbc_offset=0
+  unsigned int ptbc_offset = app()->ptbc.instance_id * max_seed / app()->ptbc.n_instances;
+  seed += ptbc_offset;
+
   loc_seed = (seed + step * max_seed) % 2147483647;
 
   if (loc_seed == 0) loc_seed++;

From 73ae240479c954c91e949df81f16a3bd54ca9257 Mon Sep 17 00:00:00 2001
From: JingJing Li <lijingjing416@gmail.com>
Date: Thu, 4 Jun 2026 11:47:17 +0200
Subject: [PATCH 7/8] update PTBC swap MPI utilities and add instance base rank
 tracking

---
 src/lib/include/ptbc.h | 10 +++++-
 src/lib/ptbc/defect.c  | 77 +++++++++++++++++++++++++++++++++++-------
 2 files changed, 73 insertions(+), 14 deletions(-)

diff --git a/src/lib/include/ptbc.h b/src/lib/include/ptbc.h
index 9652693cb..bdc8053d9 100644
--- a/src/lib/include/ptbc.h
+++ b/src/lib/include/ptbc.h
@@ -1,3 +1,5 @@
+#ifndef PTBC_H
+#define PTBC_H
 
 #include <mpi.h>
 #include <app.h>
@@ -8,6 +10,10 @@ bool is_defect(PTBCDefect *def, int const ix, int const mu);
 double get_ptbc_coeff(int const ix, int const mu);
 
 
+void mpi_gather_base_rank();
+void mpi_bcast_base_rank();
+void mpi_base_rank_update_fini();
+
 /* for swapping rng */
 typedef struct{
     int state_s[105];   // state of rlxs
@@ -44,4 +50,6 @@ int const get_node_n_children(int const node_id);
 
 // initialiser
 void init_ptbc_tree();
-void print_ptbc_topo();
\ No newline at end of file
+void print_ptbc_topo();
+
+#endif
\ No newline at end of file
diff --git a/src/lib/ptbc/defect.c b/src/lib/ptbc/defect.c
index 6a8d8c8e4..9310e072b 100644
--- a/src/lib/ptbc/defect.c
+++ b/src/lib/ptbc/defect.c
@@ -11,10 +11,6 @@
 #include "ranlxs.h"
 #include "ranlxd.h"
 
-MPI_Datatype swap_info_type= MPI_DATATYPE_NULL; // MPI datatype for swapping info
-static Tree tree;
-static Node nodes[MAX_N_INSTANCES]; // each node corresponds to an instance
-
 #define err(test, ...) err_impl(test, __func__, __FILE__, __LINE__, __VA_ARGS__)
 static void err_impl(const bool test, const char* func, const char* file, const int line, const char* format, ...)
 {
@@ -30,6 +26,15 @@ static void err_impl(const bool test, const char* func, const char* file, const
     }
 }
 
+static int base_rank[MAX_N_INSTANCES];  // MPI rank offset of each instance
+static MPI_Comm leader_comm;
+static bool leader_comm_initialised = false;
+static MPI_Request stats[2];
+
+MPI_Datatype swap_info_type= MPI_DATATYPE_NULL; // MPI datatype for swapping info
+static Tree tree;
+static Node nodes[MAX_N_INSTANCES]; // each node corresponds to an instance
+
 // find distance end-start corrected for peridic bc 
 int const static dist(int start, int end, int const g_length) {
   int const d = end - start;
@@ -121,10 +126,9 @@ void swap_rng(int const dest_inst) {
   if (dest_inst == app()->ptbc.instance_id) return;
 
   // find dest rank and src rank
-  int n_ranks, local_rank;
+  int local_rank;
   MPI_Comm_rank(app()->mpi.comm, &local_rank);
-  MPI_Comm_size(app()->mpi.world_comm, &n_ranks);
-  int const dest_rank = n_ranks / app()->ptbc.n_instances * dest_inst + local_rank;  
+  int const dest_rank = base_rank[dest_inst] + local_rank;  
   int const src_rank = app()->mpi.world_rank;
 
   // get states
@@ -133,11 +137,10 @@ void swap_rng(int const dest_inst) {
 
   // define MPI datatype for swapping info if not defined
   if (swap_info_type == MPI_DATATYPE_NULL) {
-    PTBCInstance const *instance = &(app()->ptbc.instances[app()->ptbc.instance_id]);
-    int const blocklengths[3] = {105, 105, instance->n_coeffs};
-    MPI_Aint const displacements[3] = {0, 105 * sizeof(int), 2 * 105 * sizeof(int)};
-    MPI_Datatype const types[3] = {MPI_INT, MPI_INT, MPI_DOUBLE};
-    MPI_Type_create_struct(3, blocklengths, displacements, types, &swap_info_type);
+    int const blocklengths[2] = {105, 105};
+    MPI_Aint const displacements[2] = {0, 105 * sizeof(int)};
+    MPI_Datatype const types[2] = {MPI_INT, MPI_INT};
+    MPI_Type_create_struct(2, blocklengths, displacements, types, &swap_info_type);
     MPI_Type_commit(&swap_info_type);
   }
 
@@ -156,6 +159,51 @@ void swap_rng(int const dest_inst) {
   return;
 }
 
+
+/**
+ * @brief Gather base rank offset from leader rank of all instances. Must use with the other two!
+ * 
+ */
+void mpi_gather_base_rank() {
+  int my_rank;
+  MPI_Comm_rank(app()->mpi.comm, &my_rank);
+
+  if (my_rank != 0) return; // only local rank 0 needs to send information, avoid redundancy
+
+  // set communicator if not yet
+  if (!leader_comm_initialised) {
+    int leader_color = (my_rank == 0) ? 0 : MPI_UNDEFINED;
+    
+    // leader rank == instance_id
+    MPI_Comm_split(app()->mpi.world_comm, leader_color, app()->ptbc.instance_id, &leader_comm);
+    leader_comm_initialised = true;
+  }
+
+  // collect base rank offsets from leader ranks
+  MPI_Iallgather(&(app()->mpi.world_rank), 1, MPI_INT, base_rank, 1, MPI_INT, leader_comm, stats);
+}
+
+/**
+ * @brief Broadcast base rank offset to all local ranks from leader rank. Must use with the other two!
+ * 
+ */
+void mpi_bcast_base_rank() {
+  // check if gather complete
+  MPI_Wait(stats, MPI_STATUS_IGNORE);
+
+  // broadcast from leader rank (rank 0) to all other local ranks
+  MPI_Ibcast(base_rank, app()->ptbc.n_instances, MPI_INT, 0, app()->mpi.comm, stats+1);
+}
+
+/**
+ * @brief Finalise base rank update.
+ * 
+ */
+void mpi_base_rank_update_fini(){
+  // finalise
+  MPI_Wait(stats+1, MPI_STATUS_IGNORE);
+}
+
 typedef struct {
   unsigned int n_coeffs;
   double coeffs[MAX_N_INSTANCES];
@@ -378,7 +426,10 @@ void init_ptbc_tree() {
   return;
 }
 
-
+/**
+ * @brief Print information of instances, including instance type upstream IDs and downstream IDs 
+ * 
+ */
 void print_ptbc_topo() {
   // Iterate through the nodes
   const char* typeLabels[] = {

From a498ca4e2094c9a501a78d62d8af09f756ec37f8 Mon Sep 17 00:00:00 2001
From: JingJing Li <lijingjing416@gmail.com>
Date: Fri, 12 Jun 2026 18:28:37 +0200
Subject: [PATCH 8/8] Add defect overlap check

---
 src/lib/include/ptbc.h |  20 +++++++
 src/lib/ptbc/defect.c  | 118 ++++++++++++++++++++++++++++++-----------
 2 files changed, 107 insertions(+), 31 deletions(-)

diff --git a/src/lib/include/ptbc.h b/src/lib/include/ptbc.h
index bdc8053d9..0af457753 100644
--- a/src/lib/include/ptbc.h
+++ b/src/lib/include/ptbc.h
@@ -1,3 +1,23 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2026 JingJing Li
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD 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 tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ *******************************************************************************/
+
 #ifndef PTBC_H
 #define PTBC_H
 
diff --git a/src/lib/ptbc/defect.c b/src/lib/ptbc/defect.c
index 9310e072b..7cd12b381 100644
--- a/src/lib/ptbc/defect.c
+++ b/src/lib/ptbc/defect.c
@@ -1,3 +1,23 @@
+/***********************************************************************
+ *
+ * Copyright (C) 2026 JingJing Li
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD 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 3 of the License, or
+ * (at your option) any later version.
+ *
+ * tmLQCD 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 tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ *******************************************************************************/
+
 #include <stdbool.h>
 #include <sys/types.h>
 #include <sys/stat.h>
@@ -43,6 +63,34 @@ int const static dist(int start, int end, int const g_length) {
   return dist;
 }
 
+
+/** 
+ * @brief           check a point's relative position to a 4D hypercube
+ * 
+ * @param coords    coordinate of the point, 4D array
+ * @param pos       position of hypercube, 4D array
+ * @param Ld        extent of hypercube, 4D array  
+ * 
+ * @returns         1 for in hypercubic, 2 for potentially cutting on the left, 0 for not in hypercube
+ */
+int static in_hypercube(int const * coords, int const *pos, int const *Ld) {
+  if (dist(pos[0], coords[0], T*g_nproc_t)<=Ld[0] && dist(pos[1], coords[1], LX*g_nproc_x)<=Ld[1]
+      && dist(pos[2], coords[2], LY*g_nproc_y)<=Ld[2] && dist(pos[3], coords[3], LZ*g_nproc_z)<=Ld[3]
+      && dist(pos[0]+Ld[0], coords[0], T*g_nproc_t)<=Ld[0] && dist(pos[1]+Ld[1], coords[1], LX*g_nproc_x)<=Ld[1]
+      && dist(pos[2]+Ld[2], coords[2], LY*g_nproc_y)<=Ld[2] && dist(pos[3]+Ld[3], coords[3], LZ*g_nproc_z)<=Ld[3]){
+    if(dist(pos[0], coords[0], T*g_nproc_t)>0 && dist(pos[1], coords[1], LX*g_nproc_x)>0
+      && dist(pos[2], coords[2], LY*g_nproc_y)>0 && dist(pos[3], coords[3], LZ*g_nproc_z)>0) {
+        return 1; // in hypercube
+    }
+    else if (!(dist(pos[0], coords[0], T*g_nproc_t)==0 && dist(pos[1], coords[1], LX*g_nproc_x)==0
+      && dist(pos[2], coords[2], LY*g_nproc_y)==0 && dist(pos[3], coords[3], LZ*g_nproc_z)==0)){
+      return 2; // slightly to the left of hypercube, vector starting from these points in +ve direction can cut hypercube
+    }
+  }
+  return 0; // not in hypercube
+
+}
+
 /**
  * @brief      check if a link lie in defect region, cuts at pos + 1/2 and pos + Ld + 1/2
  *             if either or end of the link is in the defect region, return true 
@@ -54,39 +102,20 @@ int const static dist(int start, int end, int const g_length) {
  */
 bool is_defect(PTBCDefect *def, int const ix, int const mu) {
   if (ix >= VOLUME) return false;
-  int const global_dim[4] = {T*g_nproc_t, LX*g_nproc_t, LY*g_nproc_x, LZ*g_nproc_z};
   int* coords = g_coord[ix];
-
-  // check if mu direction is in the defect or crossing the cut
-  // i.e. start >= pos[mu] and end <= pos[mu] + Ld[mu] + 1
-  // dist >= 0 or dist + 1 <= Ld + 1
-  int const dist_mu = dist(def->pos[mu], coords[mu], global_dim[mu]);
-  if (dist_mu >= 0 && dist_mu <= def->Ld[mu]) {
-    // find three non-mu directions
-    int dim3[3];
-    int count = 0;
-    for (int d=0; d<4; d++) {
-      if (d != mu) {
-        dim3[count] = d;
-        count++;
-      }
-    }
-    // check if the other three directions are within the defect region
-    // i.e. start at > pos end at <= pos + Ld
-    // dist > 0 or dist + 1 <= Ld
-    for (int d=0; d<3; d++) {
-      int const dist_d = dist(def->pos[dim3[d]], coords[dim3[d]], global_dim[dim3[d]]);
-      if (dist_d > 0 && dist_d < def->Ld[dim3[d]]){
-        continue;
-      }
-      else{
-        return false;
-      }
-    }
-    return true;
+  
+  int const point_stat = in_hypercube(coords, def->pos, def->Ld);
+
+  if (point_stat == 0) return false;  // start point outside
+  else if (point_stat == 1) return true;  // start point inside
+  else if (point_stat == 2) {
+    int end_coord[] = {coords[0], coords[1], coords[2], coords[3]};
+    end_coord[mu] += 1;
+    if (in_hypercube(end_coord, def->pos, def->Ld) == 1) return true;  // start point outside but ends inside
+    else return false;  // start/end point outside
   }
-  else
-    return false;  
+  
+  return false;
 }
 
 /**
@@ -270,6 +299,33 @@ void init_ptbc_tree() {
   int n_periodic = 0; // number of periodic instances
   int periodic_id[MAX_N_INSTANCES]; // store periodic instance id, may be multiple periodic instances
 
+  // check if defects are valid / no overlap
+  for (int i=0; i<ptbc_ctx->n_defects-1; i++) {
+    // define diagonal of the hypercubic
+    PTBCDefect const* def_ref = &(ptbc_ctx->defects[i]);
+    
+    // loop over compare defects
+    for (int j=i+1; i<ptbc_ctx->n_defects; j++) {
+      PTBCDefect const* def = &(ptbc_ctx->defects[j]);
+      int const start[] = {def->Ld[0], def->Ld[1], def->Ld[2], def->Ld[3]};
+      int const end[] = {def->Ld[0] + def->pos[0], def->Ld[1] + def->pos[1], 
+                        def->Ld[2] + def->pos[2], def->Ld[3] + def->pos[3]};
+
+      bool def_overlap=false;
+      // start point must be outside
+      if (in_hypercube(start, def_ref->Ld, def_ref->Ld) ){
+        def_overlap = true;
+        break;
+      }
+      else if(in_hypercube(end, def_ref->Ld, def_ref->Ld)){
+        def_overlap=true;
+        break;
+      }
+      err((def_overlap), "Some defects overlap!");
+    }
+  }
+
+
   // initialise
   for (int i=0; i<ptbc_ctx->n_instances; i++) {
     init_node(i);