diff --git a/bfps/DNS.py b/bfps/DNS.py
index 886f973cfc7de463d05761cac0ea20bfee39355b..44c1cca1366dca33c846f7bdafef9d1aef2444b0 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -611,35 +611,23 @@ class DNS(_code):
         parser_NSVEparticles_no_output = subparsers.add_parser(
                 'NSVEparticles_no_output',
                 help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers, checkpoints are NOT SAVED')
-        self.simulation_parser_arguments(parser_NSVEparticles_no_output)
-        self.job_parser_arguments(parser_NSVEparticles_no_output)
-        self.particle_parser_arguments(parser_NSVEparticles_no_output)
-        self.parameters_to_parser_arguments(parser_NSVEparticles_no_output)
-        self.parameters_to_parser_arguments(
-                parser_NSVEparticles_no_output,
-                self.NSVEp_extra_parameters)
 
         parser_NSVEp2 = subparsers.add_parser(
                 'NSVEparticles',
                 help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
-        self.simulation_parser_arguments(parser_NSVEp2)
-        self.job_parser_arguments(parser_NSVEp2)
-        self.particle_parser_arguments(parser_NSVEp2)
-        self.parameters_to_parser_arguments(parser_NSVEp2)
-        self.parameters_to_parser_arguments(
-                parser_NSVEp2,
-                self.NSVEp_extra_parameters)
 
         parser_NSVEp2p = subparsers.add_parser(
                 'NSVEparticlesP2P',
                 help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
-        self.simulation_parser_arguments(parser_NSVEp2p)
-        self.job_parser_arguments(parser_NSVEp2p)
-        self.particle_parser_arguments(parser_NSVEp2p)
-        self.parameters_to_parser_arguments(parser_NSVEp2p)
-        self.parameters_to_parser_arguments(
-                parser_NSVEp2p,
-                self.NSVEp_extra_parameters)
+
+        for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p']:
+            eval('self.simulation_parser_arguments({0})'.format('parser_' + parser))
+            eval('self.job_parser_arguments({0})'.format('parser_' + parser))
+            eval('self.particle_parser_arguments({0})'.format('parser_' + parser))
+            eval('self.parameters_to_parser_arguments({0})'.format('parser_' + parser))
+            eval('self.parameters_to_parser_arguments('
+                    'parser_{0},'
+                    'self.NSVEp_extra_parameters)'.format(parser))
         return None
     def prepare_launch(
             self,
diff --git a/bfps/cpp/full_code/NSVEparticlesP2P.cpp b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
index c858b38c7bb5c3148a443342c518917d9331b25c..ff1a594ffcd450a975ab3cc0627a3d1182d11995 100644
--- a/bfps/cpp/full_code/NSVEparticlesP2P.cpp
+++ b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
@@ -14,8 +14,8 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
     p2p_computer<double, long long int> current_p2p_computer;
     // TODO: particle interactions are switched off manually for testing purposes.
     // this needs to be fixed once particle interactions can be properly resolved.
-    //current_p2p_computer.setEnable(enable_p2p);
-    current_p2p_computer.setEnable(false);
+    current_p2p_computer.setEnable(enable_p2p);
+    //current_p2p_computer.setEnable(false);
 
     particles_inner_computer<double, long long int> current_particles_inner_computer(inner_v0);
     current_particles_inner_computer.setEnable(enable_inner);
diff --git a/bfps/cpp/particles/particles_inner_computer.hpp b/bfps/cpp/particles/particles_inner_computer.hpp
index 442e5220279830448af6a103e29a5b74a4aa8778..5c855cf5d4f9c9ce69200212e4929420a0c37228 100644
--- a/bfps/cpp/particles/particles_inner_computer.hpp
+++ b/bfps/cpp/particles/particles_inner_computer.hpp
@@ -41,18 +41,11 @@ public:
         for(partsize_t idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
             const partsize_t idx0 = idx_part*size_particle_positions + 3;
             const partsize_t idx1 = idx_part*size_particle_rhs + 3;
-            //real_number alpha[3] = {0, 0, 0};
             // check that orientation is unit vector:
             real_number orientation_size = sqrt(
                     pos_part[idx0+IDX_X]*pos_part[idx0+IDX_X] +
                     pos_part[idx0+IDX_Y]*pos_part[idx0+IDX_Y] +
                     pos_part[idx0+IDX_Z]*pos_part[idx0+IDX_Z]);
-            //DEBUG_MSG("particle ID %d\n"
-            //          "pos_part[%d] = %g, pos_part[%d] = %g, pos_part[%d] = %g\n",
-            //          idx_part,
-            //          IDX_X, pos_part[idx0 + IDX_X],
-            //          IDX_Y, pos_part[idx0 + IDX_Y],
-            //          IDX_Z, pos_part[idx0 + IDX_Z]);
             assert(orientation_size > 0.99);
             assert(orientation_size < 1.01);
             // I call "rotation" to be the right hand side of the orientation part of the ODE
@@ -61,14 +54,25 @@ public:
                     pos_part[idx0+IDX_X]*rhs_part[idx1+IDX_X] +
                     pos_part[idx0+IDX_Y]*rhs_part[idx1+IDX_Y] +
                     pos_part[idx0+IDX_Z]*rhs_part[idx1+IDX_Z]);
-            //// alpha is the vector that makes rotation perpendicular to orientation.
-            //// note that the following three lines assume the current orientation is a unit vector.
-            //alpha[IDX_X] = -;
-            //alpha[IDX_Y] = -;
-            //alpha[IDX_Z] = -;
-            //    DEBUG_MSG("projection = %g\n"
+
+            // now remove parallel bit.
+            rhs_part[idx1+IDX_X] -= pos_part[idx0+IDX_X]*projection;
+            rhs_part[idx1+IDX_Y] -= pos_part[idx0+IDX_Y]*projection;
+            rhs_part[idx1+IDX_Z] -= pos_part[idx0+IDX_Z]*projection;
+
+            // DEBUG
+            // sanity check, for debugging purposes
+            // compute dot product between orientation and orientation change
+            //real_number dotproduct = (
+            //        rhs_part[idx1 + IDX_X]*pos_part[idx0 + IDX_X] +
+            //        rhs_part[idx1 + IDX_Y]*pos_part[idx0 + IDX_Y] +
+            //        rhs_part[idx1 + IDX_Z]*pos_part[idx0 + IDX_Z]);
+            //if (dotproduct > 0.1)
+            //{
+            //    DEBUG_MSG("dotproduct = %g, projection = %g\n"
             //              "pos_part[%d] = %g, pos_part[%d] = %g, pos_part[%d] = %g\n"
             //              "rhs_part[%d] = %g, rhs_part[%d] = %g, rhs_part[%d] = %g\n",
+            //            dotproduct,
             //            projection,
             //            IDX_X, pos_part[idx0 + IDX_X],
             //            IDX_Y, pos_part[idx0 + IDX_Y],
@@ -76,32 +80,8 @@ public:
             //            IDX_X, rhs_part[idx1 + IDX_X],
             //            IDX_Y, rhs_part[idx1 + IDX_Y],
             //            IDX_Z, rhs_part[idx1 + IDX_Z]);
-
-            // now remove parallel bit.
-            rhs_part[idx1+IDX_X] -= pos_part[idx0+IDX_X]*projection;
-            rhs_part[idx1+IDX_Y] -= pos_part[idx0+IDX_Y]*projection;
-            rhs_part[idx1+IDX_Z] -= pos_part[idx0+IDX_Z]*projection;
-
-            // compute dot product between orientation and orientation change
-            real_number dotproduct = (
-                    rhs_part[idx1 + IDX_X]*pos_part[idx0 + IDX_X] +
-                    rhs_part[idx1 + IDX_Y]*pos_part[idx0 + IDX_Y] +
-                    rhs_part[idx1 + IDX_Z]*pos_part[idx0 + IDX_Z]);
-            if (dotproduct > 0.1)
-            {
-                DEBUG_MSG("dotproduct = %g, projection = %g\n"
-                          "pos_part[%d] = %g, pos_part[%d] = %g, pos_part[%d] = %g\n"
-                          "rhs_part[%d] = %g, rhs_part[%d] = %g, rhs_part[%d] = %g\n",
-                        dotproduct,
-                        projection,
-                        IDX_X, pos_part[idx0 + IDX_X],
-                        IDX_Y, pos_part[idx0 + IDX_Y],
-                        IDX_Z, pos_part[idx0 + IDX_Z],
-                        IDX_X, rhs_part[idx1 + IDX_X],
-                        IDX_Y, rhs_part[idx1 + IDX_Y],
-                        IDX_Z, rhs_part[idx1 + IDX_Z]);
-                assert(false);
-            }
+            //    assert(false);
+            //}
             //assert(dotproduct <= 0.1);
         }
     }