From b8eec341c1baa1e3c32ca6002d86b0efab3f33ce Mon Sep 17 00:00:00 2001
From: Olivier Coulaud <Olivier.Coulaud@inria.fr>
Date: Wed, 25 Apr 2018 15:17:26 +0200
Subject: [PATCH] Add periodic boundary conditions in Exemples ; Check periodic
 boundary conditions.

---
 Examples/MPIChebyshevInterpolationFMM.cpp     |  33 +
 Examples/MPIInterpolationFMM.hpp              | 319 ++++++
 Examples/MPIUniformInterpolationFMM.cpp       |  34 +
 .../GroupTree/testFMMInterpolationStarPU.hpp  |   2 +-
 .../Kernels/testSphericalDlpolyAlgorithm.cpp  | 978 +++++++++---------
 Tests/Utils/testFmmAlgorithmPeriodic.cpp      |  26 +-
 6 files changed, 893 insertions(+), 499 deletions(-)
 create mode 100644 Examples/MPIChebyshevInterpolationFMM.cpp
 create mode 100644 Examples/MPIInterpolationFMM.hpp
 create mode 100644 Examples/MPIUniformInterpolationFMM.cpp

diff --git a/Examples/MPIChebyshevInterpolationFMM.cpp b/Examples/MPIChebyshevInterpolationFMM.cpp
new file mode 100644
index 000000000..0b1054cb7
--- /dev/null
+++ b/Examples/MPIChebyshevInterpolationFMM.cpp
@@ -0,0 +1,33 @@
+// ==== CMAKE =====
+// @FUSE_BLAS
+// ================
+// Keep in private GIT
+//
+//
+/** \brief Chebyshev FMM example
+ *
+ * \file
+ * \authors O. Coulaud
+ *
+ * This program runs the FMM Algorithm with the interpolation kernel based on
+ * Chebyshev (grid points) interpolation (1/r kernel). It then compares the
+ * results with a direct computation.
+ */
+#include <string> 
+#include "Kernels/Chebyshev/FChebCell.hpp"
+#include "Kernels/Chebyshev/FChebSymKernel.hpp"
+//
+template<typename FReal, int ORDER> 
+using FInterpolationCell =  FChebCell<FReal, ORDER>;
+
+template<typename FReal, typename GroupCellClass,
+	 typename GroupContainerClass,
+	 typename MatrixKernelClass, int ORDER>  
+using FInterpolationKernel = FChebSymKernel<FReal,
+					    GroupCellClass,
+					    GroupContainerClass,
+					    MatrixKernelClass,
+					    ORDER> ;
+const std::string interpolationType("Chebyshev interpolation");
+
+#include "MPIInterpolationFMM.hpp" 
diff --git a/Examples/MPIInterpolationFMM.hpp b/Examples/MPIInterpolationFMM.hpp
new file mode 100644
index 000000000..fdd5a4cde
--- /dev/null
+++ b/Examples/MPIInterpolationFMM.hpp
@@ -0,0 +1,319 @@
+// See LICENCE file at project root
+
+// ==== CMAKE =====
+// @FUSE_MPI
+// @FUSE_BLAS
+// ================
+
+#include <iostream>
+#include <stdexcept>
+#include <cstdio>
+#include <cstdlib>
+
+
+#include "ScalFmmConfig.h"
+#include "Containers/FOctree.hpp"
+#include "Utils/FMpi.hpp"
+#include "Core/FFmmAlgorithmThreadProc.hpp"
+#include "Core/FFmmAlgorithmThreadProcPeriodic.hpp"
+
+#include "Files/FFmaGenericLoader.hpp"
+#include "Files/FMpiFmaGenericLoader.hpp"
+#include "Files/FMpiTreeBuilder.hpp"
+
+#include "Utils/FLeafBalance.hpp"
+
+#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
+#include "Kernels/Chebyshev/FChebSymKernel.hpp"
+#include "Kernels/Chebyshev/FChebCell.hpp"
+
+#include "Components/FSimpleLeaf.hpp"
+#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
+
+#include "Utils/FParameters.hpp"
+#include "Utils/FParameterNames.hpp"
+
+//
+// Order of the Interpolation approximation
+static constexpr unsigned ORDER = 6 ;
+using FReal                 = double;
+//   1/r kernel
+//
+using MatrixKernelClass     = FInterpMatrixKernelR<FReal> ;
+//
+// Typedefs
+using ContainerClass = FP2PParticleContainerIndexed<FReal>;
+using LeafClass      = FSimpleLeaf<FReal, ContainerClass>;
+
+using CellClass      = FInterpolationCell<FReal, ORDER>;
+
+
+using OctreeClass    = FOctree<FReal,CellClass,ContainerClass,LeafClass>;
+
+using MatrixKernelClass = FInterpMatrixKernelR<FReal>;
+const MatrixKernelClass MatrixKernel;
+
+using KernelClass    = FInterpolationKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> ;
+
+using FmmClassProc     = FFmmAlgorithmThreadProc<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
+using FmmClassProcPER  = FFmmAlgorithmThreadProcPeriodic<FReal,OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
+
+/// \file
+//!
+//! \brief This program runs the MPI FMM with Chebyshev/Lagrange interpolation of 1/r kernel
+//!  \authors B. Bramas, O. Coulaud
+//!
+//!  This code is a short example to use the FMM Algorithm Proc with Chebyshev or equispaced grid points Interpolation for the 1/r kernel
+
+
+// Simply create particles and try the kernels
+int main(int argc, char* argv[])
+{
+  ///////// PARAMETERS HANDLING //////////////////////////////////////
+  //  const FParameterNames LocalOptionPeriodicDeep { {"-periodic","-per"}, "Perdiodic boundary condition (-per 5) the box grow by a facor (3L)^5"};
+  FHelpDescribeAndExit(argc, argv,
+                       "Driver for Chebyshev Interpolation kernel using MPI  (1/r kernel). "
+                       "Usully run using : mpirun -np nb_proc_needed ./ChebyshevInterpolationAlgorithm [params].",
+                       FParameterDefinitions::OctreeHeight,
+                       FParameterDefinitions::OctreeSubHeight,
+                       FParameterDefinitions::InputFile,
+                       FParameterDefinitions::OutputFile,
+                       FParameterDefinitions::NbThreads,
+                       FParameterDefinitions::PeriodicityNbLevels);
+
+  const std::string defaultFile("../Data/test20k.fma");
+  const std::string  filename      = FParameters::getStr(argc,argv,    FParameterDefinitions::InputFile.options, defaultFile.c_str());
+  const unsigned int TreeHeight    = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
+  const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
+  const unsigned int NbThreads     = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
+  bool periodicCondition = false ;
+  if(FParameters::existParameter(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options)){
+      periodicCondition = true;
+    }
+  const unsigned int aboveTree = FParameters::getValue(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options, 5);
+
+  omp_set_num_threads(NbThreads);
+  std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
+
+  //
+  std::cout << "Parameters"<< std::endl
+            << "      Octree Depth      " << TreeHeight    << std::endl
+            << "      SubOctree depth   " << SubTreeHeight << std::endl;
+  if(periodicCondition){
+      std::cout << "      AboveTree    "<< aboveTree <<std::endl;
+
+    }
+  std::cout    << "      Input file  name: " << filename      << std::endl
+               << "      Thread count :    " << NbThreads     << std::endl
+               << std::endl;
+
+
+  ///////// VAR INIT /////////////////////////////////////////////////
+
+  // Initialize values for MPI
+  FMpi app(argc,argv);
+  //
+  // Initialize timer
+  FTic time;
+
+  // Creation of the particle loader
+  FMpiFmaGenericLoader<FReal> loader(filename,app.global());
+  if(!loader.isOpen()) {
+      throw std::runtime_error("Particle file couldn't be opened!") ;
+    }
+
+  // Initialize empty oct-tree
+  OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
+
+  FSize localParticlesNumber = 0 ;
+
+  { // -----------------------------------------------------
+    if(app.global().processId() == 0){
+        std::cout << "Loading & Inserting " << loader.getNumberOfParticles()
+                  << " particles ..." << std::endl;
+        std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
+      }
+    time.tic();
+
+    /* Mock particle structure to balance the tree over the processes. */
+    struct TestParticle{
+      FSize index;             // Index of the particle in the original file.
+      FPoint<FReal> position;  // Spatial position of the particle.
+      FReal physicalValue;     // Physical value of the particle.
+      /* Returns the particle position. */
+      const FPoint<FReal>& getPosition(){
+        return position;
+      }
+    };
+
+    // Temporary array of particles read by this process.
+    TestParticle* particles = new TestParticle[loader.getMyNumberOfParticles()];
+    memset(particles, 0, (sizeof(TestParticle) * loader.getMyNumberOfParticles()));
+
+    // Index (in file) of the first particle that will be read by this process.
+    FSize idxStart = loader.getStart();
+    std::cout << "Proc:" << app.global().processId() << " start-index: " << idxStart << std::endl;
+
+    // Read particles from parts.
+    for(FSize idxPart = 0 ; idxPart < loader.getMyNumberOfParticles() ; ++idxPart){
+        // Store the index (in the original file) the particle.
+        particles[idxPart].index = idxPart + idxStart;
+        // Read particle from file
+        loader.fillParticle(&particles[idxPart].position,
+                            &particles[idxPart].physicalValue);
+      }
+
+    // Final vector of particles
+    FVector<TestParticle> finalParticles;
+    FLeafBalance balancer;
+    // Redistribute particules between processes
+    FMpiTreeBuilder< FReal, TestParticle >::
+        DistributeArrayToContainer(app.global(),
+                                   particles,
+                                   loader.getMyNumberOfParticles(),
+                                   tree.getBoxCenter(),
+                                   tree.getBoxWidth(),
+                                   tree.getHeight(),
+                                   &finalParticles,
+                                   &balancer);
+
+    // Free temporary array memory.
+    delete[] particles;
+
+    // Insert final particles into tree.
+
+    for(FSize idx = 0 ; idx < finalParticles.getSize(); ++idx){
+        tree.insert(finalParticles[idx].position,
+                    finalParticles[idx].index,
+                    finalParticles[idx].physicalValue);
+      }
+
+    time.tac();
+
+    localParticlesNumber = finalParticles.getSize() ;
+
+    double timeUsed = time.elapsed();
+    double minTime,maxTime;
+    std::cout << "Proc:" << app.global().processId()
+              << " "     << finalParticles.getSize()
+              << " particles have been inserted in the tree. (@Reading and Inserting Particles = "
+              << time.elapsed() << " s)."
+              << std::endl;
+
+    MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
+    MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
+    if(app.global().processId() == 0){
+        std::cout << "readinsert-time-min:" << minTime
+                  << " readinsert-time-max:" << maxTime
+                  << std::endl;
+      }
+  } // -----------------------------------------------------
+  std::vector<MortonIndex> mortonLeafDistribution(2*app.global().processCount());
+  { // -----------------------------------------------------
+    std::cout << "\n"<<interpolationType<<" FMM Proc (ORDER="<< ORDER << ") ... " << std::endl;
+
+    time.tic();
+
+    // Kernels to use (pointer because of the limited size of the stack)
+
+    FAbstractAlgorithm * algorithm  = nullptr;
+    FAlgorithmTimers   * timer      = nullptr;
+    // non periodic FMM algorithm
+    std::unique_ptr<KernelClass> kernelsNoPer(new KernelClass(TreeHeight, loader.getBoxWidth(),
+                                                              loader.getCenterOfBox(),&MatrixKernel));
+    FmmClassProc    algoNoPer(app.global(),&tree, kernelsNoPer.get());
+    //
+    // periodic FMM algorithm
+    FmmClassProcPER algoPer(app.global(),&tree, aboveTree);
+    KernelClass kernelsPer(algoPer.extendedTreeHeight(), algoPer.extendedBoxWidth(),
+                           algoPer.extendedBoxCenter(),&MatrixKernel);
+    algoPer.setKernel(&kernelsPer);
+    ///////////////////////////////////////////////////////////////////////////////////////////////////
+    if(! periodicCondition) {// Non periodic case
+        algorithm  = &algoNoPer ;
+        timer      = &algoNoPer ;
+      }
+    else {  // Periodic case
+        algorithm  = &algoPer ;
+        timer      = &algoPer ;
+      }
+    //
+    // FMM exectution  FFmmFarField
+    algorithm->execute();
+    //
+
+    time.tac();
+    double timeUsed = time.elapsed();
+    double minTime,maxTime;
+    std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s)." << std::endl;
+    MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
+    MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
+    if(app.global().processId() == 0){
+        std::cout << "exec-time-min:" << minTime
+                  << " exec-time-max:" << maxTime
+                  << std::endl;
+  }
+  // -----------------------------------------------------
+  //
+  // Some output
+  //
+  //
+  { // -----------------------------------------------------
+    FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= (loader.getNumberOfParticles()-1); ;
+    FReal energy =0.0 ;
+    //
+    //   Loop over all leaves
+    //
+    std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
+    std::cout << std::scientific;
+    std::cout.precision(10) ;
+
+    tree.forEachLeaf([&](LeafClass* leaf){
+      const FReal*const posX = leaf->getTargets()->getPositions()[0];
+      const FReal*const posY = leaf->getTargets()->getPositions()[1];
+      const FReal*const posZ = leaf->getTargets()->getPositions()[2];
+
+      const FReal*const potentials = leaf->getTargets()->getPotentials();
+      const FReal*const forcesX = leaf->getTargets()->getForcesX();
+      const FReal*const forcesY = leaf->getTargets()->getForcesY();
+      const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
+      const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
+      const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
+
+      const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
+
+      for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
+          const FSize indexPartOrig = indexes[idxPart];
+          if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3)  ) {
+              std::cout << "Proc "<< app.global().processId() << " Index "<< indexPartOrig <<"  potential  " << potentials[idxPart]
+                           << " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
+                              << "   Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
+            }
+          energy += potentials[idxPart]*physicalValues[idxPart] ;
+        }
+    });
+    FReal gloEnergy = app.global().reduceSum(energy);
+    if(0 == app.global().processId()){
+        std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
+        std::cout <<std::endl <<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
+      }
+  }
+  // -----------------------------------------------------
+  if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
+      algorithm->getAndBuildMortonIndexAtLeaf(mortonLeafDistribution);
+      //
+      if(app.global().processId() == 0)
+        {
+          std::cout << " Morton distribution "<< std::endl ;
+          for(auto v : mortonLeafDistribution)
+            std::cout << v << " ";
+
+          std::cout <<  std::endl;
+        }
+      std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.fma"));
+      FMpiFmaGenericWriter<FReal> paraWriter(name,app);
+      paraWriter.writeDistributionOfParticlesFromOctree(tree,loader.getNumberOfParticles(),localParticlesNumber,mortonLeafDistribution);
+
+    }
+  return 0;
+}
diff --git a/Examples/MPIUniformInterpolationFMM.cpp b/Examples/MPIUniformInterpolationFMM.cpp
new file mode 100644
index 000000000..c464f2891
--- /dev/null
+++ b/Examples/MPIUniformInterpolationFMM.cpp
@@ -0,0 +1,34 @@
+// ==== CMAKE =====
+// @FUSE_BLAS
+// @FUSE_FFT
+// ================
+// Keep in private GIT
+//
+//
+/** \brief Uniform FMM example
+ *
+ * \file
+ * \authors B. Bramas, O. Coulaud
+ *
+ * This program runs the FMM Algorithm with the interpolation kernel based on
+ * uniform (grid points) interpolation (1/r kernel). It then compares the
+ * results with a direct computation.
+ */
+#include <string> 
+#include "Kernels/Uniform/FUnifCell.hpp"
+#include "Kernels/Uniform/FUnifKernel.hpp"
+//
+template<typename FReal, int ORDER> 
+using FInterpolationCell =  FUnifCell<FReal, ORDER>;
+
+template<typename FReal, typename GroupCellClass,
+	 typename GroupContainerClass,
+	 typename MatrixKernelClass, int ORDER>  
+using FInterpolationKernel = FUnifKernel<FReal,
+					    GroupCellClass,
+					    GroupContainerClass,
+					    MatrixKernelClass,
+					    ORDER> ;
+const std::string interpolationType("Uniform interpolation");
+
+#include "MPIInterpolationFMM.hpp" 
diff --git a/Tests/GroupTree/testFMMInterpolationStarPU.hpp b/Tests/GroupTree/testFMMInterpolationStarPU.hpp
index 443aa54df..7c44bd0b0 100644
--- a/Tests/GroupTree/testFMMInterpolationStarPU.hpp
+++ b/Tests/GroupTree/testFMMInterpolationStarPU.hpp
@@ -280,7 +280,7 @@ int main(int argc, char *argv[]) {
       std::string outputFile(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options,   "outputMPIFMA"));
 
       FMpiFmaGenericWriter<FReal> paraWriter(outputFile,FMpiComm);
-      paraWriter.writeDistributionOfParticlesFromOctree(*computeOnGroupTree,loader.getNumberOfParticles(),mortonCellDistribution);
+      paraWriter.writeDistributionOfParticlesFromGroupedOctree(*computeOnGroupTree,loader.getNumberOfParticles(),mortonCellDistribution);
     }
   return 0;
 }
diff --git a/Tests/Kernels/testSphericalDlpolyAlgorithm.cpp b/Tests/Kernels/testSphericalDlpolyAlgorithm.cpp
index 3e6169de1..51fccdcd1 100644
--- a/Tests/Kernels/testSphericalDlpolyAlgorithm.cpp
+++ b/Tests/Kernels/testSphericalDlpolyAlgorithm.cpp
@@ -51,512 +51,520 @@
 
 template <class FReal>
 struct EwalParticle {
-    FPoint<FReal> position;
-	FReal forces[3];
-	FReal physicalValue;
-	FReal potential;
-	int index;
+  FPoint<FReal> position;
+  FReal forces[3];
+  FReal physicalValue;
+  FReal potential;
+  int index;
 };
 
 // Simply create particles and try the kernels
 int main(int argc, char ** argv){
-    FHelpDescribeAndExit(argc, argv, "Please read the code to know more, sorry");
-
-    typedef double FReal;
-	typedef FP2PParticleContainerIndexed<FReal>                    ContainerClass;
-    typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
+  FHelpDescribeAndExit(argc, argv, "Please read the code to know more, sorry",
+                       FParameterDefinitions::OctreeHeight,
+                       FParameterDefinitions::OctreeSubHeight,
+                       FParameterDefinitions::InputFile,
+                       FParameterDefinitions::OutputFile,
+                       FParameterDefinitions::NbThreads,
+                       FParameterDefinitions::EnabledVerbose,
+                       FParameterDefinitions::PeriodicityNbLevels);
+
+  typedef double FReal;
+  typedef FP2PParticleContainerIndexed<FReal>                    ContainerClass;
+  typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
 
 #ifdef  SCALFMM_USE_BLAS
-	// begin Chebyshev kernel
-	// accuracy
-	const unsigned int ORDER = 13;
-	// typedefs
-    typedef FInterpMatrixKernelR<FReal>                                              MatrixKernelClass;
-    typedef FChebCell<FReal,ORDER>                                                  CellClass;
-    typedef FOctree<FReal, CellClass,ContainerClass,LeafClass>                       OctreeClass;
-    typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER>  KernelClass;
+  // begin Chebyshev kernel
+  // accuracy
+  const unsigned int ORDER = 6;
+  // typedefs
+  typedef FInterpMatrixKernelR<FReal>                                              MatrixKernelClass;
+  typedef FChebCell<FReal,ORDER>                                                  CellClass;
+  typedef FOctree<FReal, CellClass,ContainerClass,LeafClass>                       OctreeClass;
+  typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER>  KernelClass;
 
 #else
-    typedef FSphericalCell<FReal>                                    CellClass;
-    typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
-    typedef FSphericalKernel<FReal, CellClass, ContainerClass >     KernelClass;
-    const int DevP          = FParameters::getValue(argc,argv,"-P", 9);
+  typedef FSphericalCell<FReal>                                    CellClass;
+  typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
+  typedef FSphericalKernel<FReal, CellClass, ContainerClass >     KernelClass;
+  const int DevP          = FParameters::getValue(argc,argv,"-P", 9);
 #endif
 
-    using FmmClass     = FFmmAlgorithmPeriodic<FReal,OctreeClass,  CellClass, ContainerClass, 
-					       KernelClass, LeafClass > ;
-    using FmmClassNoPer= FFmmAlgorithm<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass >       ;
-	///////////////////////What we do/////////////////////////////
-	if( FParameters::existParameter(argc, argv, "-help")){
-		std::cout << ">> This executable has to be used to compute direct interaction either for periodic or non periodic system.\n";
-		std::cout << ">> options are -h H -sh SH  [-per perdeep,  -noper] -fin filenameIN (-bin)  -[no]scale \n";
-		std::cout << ">> Recommended files : ../Data/EwalTest_Periodic.run ../Data/EwalTest_NoPeriodic.run\n";
-		std::cout << " Options " << std::endl;
-		std::cout << "     -per perDeep    " << std::endl;
-		std::cout << "     -noper no periodic boundary conditions   " << std::endl;
-		std::cout << "     -verbose : print index x y z fx fy fy Q and V" << std::endl;
-		std::cout << "     -noscale no scaling and we don't remove the dipole term " << std::endl;
-		exit(-1);
-
-	}	//////////////////////////////////////////////////////////////
-
-    const int NbLevels         = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options,   4);
-    const int SizeSubLevels    = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options,  2);
-	const int PeriodicDeep     = FParameters::getValue(argc,argv,"-per", 3);
-    const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/EwalTest_Periodic.run");
-	//  file for -saveError option
-	std::ofstream errorfile("outputEwaldError.txt",std::ios::out);
-
-	FTic counter;
-	//    c     conversion factor for coulombic terms in internal units
-	//    c     i.e. (unit(charge)**2/(4 pi eps0 unit(length))/unit(energy)
-	const FReal r4pie0 = FReal(138935.4835);
-	FReal scaleEnergy, scaleForce  ;
-	// -----------------------------------------------------
-	//  LOADER
-	//  -----------------------------------------------------
-	std::cout << "Opening : " << filename << "\n";
-    FDlpolyLoader<FReal>  *loader = nullptr ;
-	if(FParameters::existParameter(argc, argv, "-bin")){
-        loader  = new FDlpolyBinLoader<FReal>(filename);
-	}
-	else {
-        loader  = new FDlpolyAsciiLoader<FReal>(filename);
-	}
-
-	if(! loader->isOpen()){
-		std::cout << "Loader Error, " << filename << " is missing\n";
-		return 1;
-	}
-	// ---------------------------------------------------
-	//        DL_POLY CONSTANT
-	//  ---------------------------------------------------
-	bool scale = true ;
-	if(FParameters::existParameter(argc, argv, "-noscale")){
-		scale = false ;
-		scaleEnergy =  1.0;   // kcal mol^{-1}
-		scaleForce  = 1.0 ;           // 10 J mol^{-1} A^{-1}
-	}
-	else {
-		scaleEnergy =  r4pie0 / 418.4 ;   // kcal mol^{-1}
-		scaleForce  = -r4pie0 ;           // 10 J mol^{-1} A^{-1}
-	}
-
-
-	//
+  using FmmClass     = FFmmAlgorithmPeriodic<FReal,OctreeClass,  CellClass, ContainerClass,
+  KernelClass, LeafClass > ;
+  using FmmClassNoPer= FFmmAlgorithm<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass >       ;
+  ///////////////////////What we do/////////////////////////////
+  if( FParameters::existParameter(argc, argv, "-help")){
+      std::cout << ">> This executable has to be used to compute direct interaction either for periodic or non periodic system.\n";
+      std::cout << ">> options are -h H -sh SH  [-per perdeep,  -noper] -fin filenameIN (-bin)  -[no]scale \n";
+      std::cout << ">> Recommended files : ../Data/EwalTest_Periodic.run ../Data/EwalTest_NoPeriodic.run\n";
+      std::cout << " Options " << std::endl;
+      std::cout << "     -per perDeep    " << std::endl;
+      std::cout << "     -noper no periodic boundary conditions   " << std::endl;
+      std::cout << "     -verbose : print index x y z fx fy fy Q and V" << std::endl;
+      std::cout << "     -noscale no scaling and we don't remove the dipole term " << std::endl;
+      exit(-1);
+
+    }	//////////////////////////////////////////////////////////////
+
+  const int NbLevels         = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeHeight.options,   4);
+  const int SizeSubLevels    = FParameters::getValue(argc,argv,FParameterDefinitions::OctreeSubHeight.options,  2);
+  const int PeriodicDeep     = FParameters::getValue(argc,argv,FParameterDefinitions::PeriodicityNbLevels.options, 3);
+  const std::string filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/EwalTest_Periodic.run");
+  //  file for -saveError option
+  std::ofstream errorfile("outputEwaldError.txt",std::ios::out);
+
+  FTic counter;
+  //    c     conversion factor for coulombic terms in internal units
+  //    c     i.e. (unit(charge)**2/(4 pi eps0 unit(length))/unit(energy)
+  const FReal r4pie0 = FReal(138935.4835);
+  FReal scaleEnergy, scaleForce  ;
+  // -----------------------------------------------------
+  //  LOADER
+  //  -----------------------------------------------------
+  std::cout << "Opening : " << filename << "\n";
+  FDlpolyLoader<FReal>  *loader = nullptr ;
+  if( filename.find(".bin") != std::string::npos ){
+      loader  = new FDlpolyBinLoader<FReal>(filename.c_str());
+    }
+  else {
+      loader  = new FDlpolyAsciiLoader<FReal>(filename.c_str());
+    }
+
+  if(! loader->isOpen()){
+      std::cout << "Loader Error, " << filename << " is missing\n";
+      return 1;
+    }
+  // ---------------------------------------------------
+  //        DL_POLY CONSTANT
+  //  ---------------------------------------------------
+  bool scale = true ;
+  if(FParameters::existParameter(argc, argv, "-noscale")){
+      scale = false ;
+      scaleEnergy =  1.0;   // kcal mol^{-1}
+      scaleForce  = 1.0 ;           // 10 J mol^{-1} A^{-1}
+    }
+  else {
+      scaleEnergy =  r4pie0 / 418.4 ;   // kcal mol^{-1}
+      scaleForce  = -r4pie0 ;           // 10 J mol^{-1} A^{-1}
+    }
+
+
+  //
 #ifndef  SCALFMM_USE_BLAS
-	CellClass::Init(DevP);
-	std::cout << " $$$$$$$$$$$$$$$  SPHERICAL VERSION $$$$$$$$$$$$"<<std::endl;
-	std::cout << " $$$$$$$$$$$$$$$  Order "<< DevP <<" $$$$$$$$$$$$"<<std::endl;
+  CellClass::Init(DevP);
+  std::cout << " $$$$$$$$$$$$$$$  SPHERICAL VERSION $$$$$$$$$$$$"<<std::endl;
+  std::cout << " $$$$$$$$$$$$$$$  Order "<< DevP <<" $$$$$$$$$$$$"<<std::endl;
 #else
-	std::cout << " $$$$$$$$$$$$$$$  CHEBYCHEV VERSION $$$$$$$$$$$$" <<std::endl;
-	std::cout << " $$$$$$$$$$$$$$$  Order "<<ORDER <<" $$$$$$$$$$$$"<<std::endl;
+  std::cout << " $$$$$$$$$$$$$$$  CHEBYCHEV VERSION $$$$$$$$$$$$" <<std::endl;
+  std::cout << " $$$$$$$$$$$$$$$  Order "<<ORDER <<" $$$$$$$$$$$$"<<std::endl;
 #endif
-	OctreeClass tree(NbLevels, SizeSubLevels, loader->getBoxWidth(), loader->getCenterOfBox());
-	// ---------------------------------------------------------------------------------
-	// Insert particles in the Octree
-	// ---------------------------------------------------------------------------------   std::cout << "Creating & Inserting " << loader->getNumberOfParticles() << " particles ..." << std::endl;
-	std::cout << "\tWidth : " << loader->getBoxWidth() << " \t center x : " << loader->getCenterOfBox().getX()
-	    				<< " y : " << loader->getCenterOfBox().getY() << " z : " << loader->getCenterOfBox().getZ() << std::endl;
-	std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
-
-	counter.tic();
-    FPoint<FReal> electricMoment(0.0,0.0,0.0) ;
-    EwalParticle<FReal> * const particles = new EwalParticle<FReal>[loader->getNumberOfParticles()];
-    memset(particles, 0, sizeof(EwalParticle<FReal>) * loader->getNumberOfParticles());
-	double totalCharge = 0.0;
-	for(FSize idxPart = 0 ; idxPart < loader->getNumberOfParticles() ; ++idxPart){
-		//
-		loader->fillParticle(&particles[idxPart].position, particles[idxPart].forces,
-				&particles[idxPart].physicalValue,&particles[idxPart].index);
-		//
-		totalCharge += particles[idxPart].physicalValue ;
-		electricMoment.incX(particles[idxPart].physicalValue*particles[idxPart].position.getX() );
-		electricMoment.incY(particles[idxPart].physicalValue*particles[idxPart].position.getY() );
-		electricMoment.incZ(particles[idxPart].physicalValue*particles[idxPart].position.getZ() );
-		// reset forces and insert in the tree
-		tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
-	}
-
-	counter.tac();
-	double dipoleNorm = electricMoment.norm2() ;
-	double volume     =  loader->getBoxWidth()*loader->getBoxWidth()*loader->getBoxWidth() ;
-	double coeffCorrectionDLPOLY = 2.0*FMath::FPi<FReal>()/volume/3.0 ;
-
-	std::cout << std::endl;
-	std::cout << "Total Charge         = "<< totalCharge <<std::endl;
-	std::cout << "Electric Moment      = "<< electricMoment <<std::endl;
-	std::cout << "Electric Moment norm = "<< dipoleNorm  <<std::endl;
-	std::cout << std::endl;
-	std::cout << std::endl;
-
-	std::cout << "Done  " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl;
-	//
-	// ---------------------------------------------------------------------------------
-	// write Octree information in octreeData.txt file
-	// ---------------------------------------------------------------------------------
-	std::ofstream octreeData("octreeData.txt",std::ios::out);
-	typename  OctreeClass::Iterator  octreeIterator(&tree);
-	octreeIterator.gotoBottomLeft();
-	int inTreeHeight = NbLevels ;
-	double  inBoxWidth = loader->getBoxWidth() ;
-    FPoint<FReal> inBoxCenter(loader->getCenterOfBox()) ;
-	double  widthAtLeafLevel(inBoxWidth/FReal(1 << (inTreeHeight-1))) , widthAtLeafLevelDiv2 = widthAtLeafLevel/2;
-    FPoint<FReal>  boxCorner(inBoxCenter.getX()-(inBoxWidth/2),inBoxCenter.getY()-(inBoxWidth/2),
-			inBoxCenter.getZ()-(inBoxWidth/2));
-	octreeData << "Box Width  " << inBoxWidth << std::endl ;
-	octreeData << "Leaf width " << widthAtLeafLevel << std::endl ;
-	octreeData << "Box Corner "<< boxCorner << std::endl<< std::endl ;
-
-	do{
-		auto * const FRestrict cell = octreeIterator.getCurrentCell();
-		FTreeCoordinate coordinate  = cell->getCoordinate() ;
-        FPoint<FReal> leafCenter(FReal(coordinate.getX()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
-				FReal(coordinate.getY()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
-				FReal(coordinate.getZ()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX());
-		octreeData << "Leaf " << cell->getMortonIndex() << std::endl
-			   << "    Center   "<<  coordinate <<    std::endl
-			   << "    Center   "<<  leafCenter
-			   << std::endl ;
-	} while(octreeIterator.moveRight());
-	//
-	FIOVtk<FReal> vtkfile ;
-	vtkfile.writeOctree("octreeFile.vtk","Octree ", tree) ;
-	//
-	// ---------------------------------------------------------------------------------
-
-	std::cout << "Create kernel & run simu ..." << std::endl;
-	counter.tic();
-
-	const FInterpMatrixKernelR<FReal> MatrixKernel;
-
-	FTreeCoordinate min{0,0,0} , max{0,0,0};
-
-	if( FParameters::existParameter(argc, argv, "-noper") ){
+  OctreeClass tree(NbLevels, SizeSubLevels, loader->getBoxWidth(), loader->getCenterOfBox());
+  // ---------------------------------------------------------------------------------
+  // Insert particles in the Octree
+  // ---------------------------------------------------------------------------------   std::cout << "Creating & Inserting " << loader->getNumberOfParticles() << " particles ..." << std::endl;
+  std::cout << "\tWidth : " << loader->getBoxWidth() << " \t center x : " << loader->getCenterOfBox().getX()
+            << " y : " << loader->getCenterOfBox().getY() << " z : " << loader->getCenterOfBox().getZ() << std::endl;
+  std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
+
+  counter.tic();
+  FPoint<FReal> electricMoment(0.0,0.0,0.0) ;
+  EwalParticle<FReal> * const particles = new EwalParticle<FReal>[loader->getNumberOfParticles()];
+  memset(particles, 0, sizeof(EwalParticle<FReal>) * loader->getNumberOfParticles());
+  double totalCharge = 0.0;
+  for(FSize idxPart = 0 ; idxPart < loader->getNumberOfParticles() ; ++idxPart){
+      //
+      loader->fillParticle(&particles[idxPart].position, particles[idxPart].forces,
+                           &particles[idxPart].physicalValue,&particles[idxPart].index);
+      //
+      totalCharge += particles[idxPart].physicalValue ;
+      electricMoment.incX(particles[idxPart].physicalValue*particles[idxPart].position.getX() );
+      electricMoment.incY(particles[idxPart].physicalValue*particles[idxPart].position.getY() );
+      electricMoment.incZ(particles[idxPart].physicalValue*particles[idxPart].position.getZ() );
+      // reset forces and insert in the tree
+      tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
+    }
+
+  counter.tac();
+  double dipoleNorm = electricMoment.norm2() ;
+  double volume     =  loader->getBoxWidth()*loader->getBoxWidth()*loader->getBoxWidth() ;
+  double coeffCorrectionDLPOLY = 2.0*FMath::FPi<FReal>()/volume/3.0 ;
+
+  std::cout << std::endl;
+  std::cout << "Total Charge         = "<< totalCharge <<std::endl;
+  std::cout << "Electric Moment      = "<< electricMoment <<std::endl;
+  std::cout << "Electric Moment norm = "<< dipoleNorm  <<std::endl;
+  std::cout << std::endl;
+  std::cout << std::endl;
+
+  std::cout << "Done  " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl;
+  //
+  // ---------------------------------------------------------------------------------
+  // write Octree information in octreeData.txt file
+  // ---------------------------------------------------------------------------------
+  std::ofstream octreeData("octreeData.txt",std::ios::out);
+  typename  OctreeClass::Iterator  octreeIterator(&tree);
+  octreeIterator.gotoBottomLeft();
+  int inTreeHeight = NbLevels ;
+  double  inBoxWidth = loader->getBoxWidth() ;
+  FPoint<FReal> inBoxCenter(loader->getCenterOfBox()) ;
+  double  widthAtLeafLevel(inBoxWidth/FReal(1 << (inTreeHeight-1))) , widthAtLeafLevelDiv2 = widthAtLeafLevel/2;
+  FPoint<FReal>  boxCorner(inBoxCenter.getX()-(inBoxWidth/2),inBoxCenter.getY()-(inBoxWidth/2),
+                           inBoxCenter.getZ()-(inBoxWidth/2));
+  octreeData << "Box Width  " << inBoxWidth << std::endl ;
+  octreeData << "Leaf width " << widthAtLeafLevel << std::endl ;
+  octreeData << "Box Corner "<< boxCorner << std::endl<< std::endl ;
+
+  do{
+      auto * const FRestrict cell = octreeIterator.getCurrentCell();
+      FTreeCoordinate coordinate  = cell->getCoordinate() ;
+      FPoint<FReal> leafCenter(FReal(coordinate.getX()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
+                               FReal(coordinate.getY()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
+                               FReal(coordinate.getZ()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX());
+      octreeData << "Leaf " << cell->getMortonIndex() << std::endl
+                 << "    Center   "<<  coordinate <<    std::endl
+                 << "    Center   "<<  leafCenter
+                 << std::endl ;
+    } while(octreeIterator.moveRight());
+  //
+  FIOVtk<FReal> vtkfile ;
+  vtkfile.writeOctree("octreeFile.vtk","Octree ", tree) ;
+  //
+  // ---------------------------------------------------------------------------------
+
+  std::cout << "Create kernel & run simu ..." << std::endl;
+  counter.tic();
+
+  const FInterpMatrixKernelR<FReal> MatrixKernel;
+
+  FTreeCoordinate min{0,0,0} , max{0,0,0};
+
+  if( ! FParameters::existParameter(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options) ){
 #ifndef  SCALFMM_USE_BLAS
-	  KernelClass kernels( DevP, NbLevels, loader->getBoxWidth(), loader->getCenterOfBox());
+      KernelClass kernels( DevP, NbLevels, loader->getBoxWidth(), loader->getCenterOfBox());
 #else
-	  KernelClass kernels( NbLevels, loader->getBoxWidth(), loader->getCenterOfBox(),&MatrixKernel);
+      KernelClass kernels( NbLevels, loader->getBoxWidth(), loader->getCenterOfBox(),&MatrixKernel);
 #endif
-	  FmmClassNoPer algo(&tree,&kernels);
-	  algo.execute();
-	}
-	else{
-	  FmmClass algo(&tree,PeriodicDeep);
-	  std::cout << "The simulation box is repeated " << algo.theoricalRepetition() << " in X/Y/Z" << std::endl;
-	  std::cout << "Simulated box: " << algo.extendedBoxWidth()<<std::endl;
+      FmmClassNoPer algo(&tree,&kernels);
+      algo.execute();
+    }
+  else{
+      FmmClass algo(&tree,PeriodicDeep);
+      std::cout << "The simulation box is repeated " << algo.theoricalRepetition() << " in X/Y/Z" << std::endl;
+      std::cout << "Simulated box: " << algo.extendedBoxWidth()<<std::endl;
 
 #ifndef  SCALFMM_USE_BLAS
-	  KernelClass kernels( DevP, algo.extendedTreeHeight(), algo.extendedBoxWidth(),algo.extendedBoxCenter());
+      KernelClass kernels( DevP, algo.extendedTreeHeight(), algo.extendedBoxWidth(),algo.extendedBoxCenter());
 #else
-	  KernelClass kernels(algo.extendedTreeHeight(), algo.extendedBoxWidth(),algo.extendedBoxCenter(),&MatrixKernel);
+      KernelClass kernels(algo.extendedTreeHeight(), algo.extendedBoxWidth(),algo.extendedBoxCenter(),&MatrixKernel);
 #endif
-	  algo.setKernel(&kernels);
-	  algo.execute();
-	  algo.repetitionsIntervals(&min, &max);
-	}
-
-	counter.tac();
-
-	std::cout << "Done  " << "(@FMM Algorithm = " << counter.elapsed() << " s)." << std::endl;
-	std::cout << "\n"<< "END FMM "
-		  << "-------------------------------------------------------------------------" << std::endl << std::endl ;
-
-	// ----------------------------------------------------------------------------------------------------------
-	//                                  DIRECT COMPUTATION
-	// ----------------------------------------------------------------------------------------------------------
-	EwalParticle<FReal>* particlesDirect = nullptr;
-	FReal directEnergy = 0.0;
+      algo.setKernel(&kernels);
+      algo.execute();
+      algo.repetitionsIntervals(&min, &max);
+    }
+
+  counter.tac();
+
+  std::cout << "Done  " << "(@FMM Algorithm = " << counter.elapsed() << " s)." << std::endl;
+  std::cout << "\n"<< "END FMM "
+            << "-------------------------------------------------------------------------" << std::endl << std::endl ;
+
+  // ----------------------------------------------------------------------------------------------------------
+  //                                  DIRECT COMPUTATION
+  // ----------------------------------------------------------------------------------------------------------
+  EwalParticle<FReal>* particlesDirect = nullptr;
+  FReal directEnergy = 0.0;
   
-	//
-	// direct computation
-	//
-	bool direct = false ;
-	if(FParameters::existParameter(argc, argv, "-direct")){
-	  direct = true ;
-	  printf("Compute direct:\n");
-	  printf("Box [%d;%d][%d;%d][%d;%d]\n", min.getX(), max.getX(), min.getY(),
-		 max.getY(), min.getZ(), max.getZ());
-    
-	  particlesDirect = new EwalParticle<FReal>[loader->getNumberOfParticles()];
-    
-	  FReal denergy = 0.0;
-	  FMath::FAccurater<FReal> dfx, dfy, dfz ;
-	  //
-	  //
-	  // particles       : Ewald results
-	  // particlesDirect : Direct results
-	  //
-	  counter.tic();
-	  for(FSize idxTarget = 0 ; idxTarget < loader->getNumberOfParticles() ; ++idxTarget){
-	    particlesDirect[idxTarget] = particles[idxTarget];
-	    EwalParticle<FReal> & part        = particlesDirect[idxTarget];
-	    part.forces[0] = part.forces[1] = part.forces[2] = 0.0;
-	    part.potential = 0.0;
-	    //
-	    // compute with all other except itself
-	    //
-	    // Compute force and potential between  particles[idxTarget] and particles inside the box
-	    //
-	    for(FSize idxOther = 0; idxOther < loader->getNumberOfParticles() ; ++idxOther){
-	      if( idxOther != idxTarget ){
-		FP2P::NonMutualParticles(
-					 part.position.getX(), part.position.getY(),
-					 part.position.getZ(),part.physicalValue,
-					 &part.forces[0],&part.forces[1],
-					 &part.forces[2],&part.potential,
-					 particles[idxOther].position.getX(), particles[idxOther].position.getY(),
-					 particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
-					 &MatrixKernel);
-	      }
-	    }
-	    //
-	    // compute with image boxes
-	    //
-	    for(int idxX = min.getX() ; idxX <= max.getX() ; ++idxX){
-	      for(int idxY = min.getY() ; idxY <= max.getY() ; ++idxY){
-		for(int idxZ = min.getZ() ; idxZ <= max.getZ() ; ++idxZ){
-	    
-		  {
-		    //int nL = PeriodicDeep ;
-		    //	for(int idxX = -nL ; idxX <= nL -1; ++idxX){
-		    //	  for(int idxY = -nL ; idxY <=nL-1 ; ++idxY){
-		    //	    for(int idxZ = -nL ; idxZ <= nL-1; ++idxZ){
-		    if(idxX == 0 && idxY == 0 && idxZ == 0) continue;
-	      
-		    const FPoint<FReal> offset(loader->getBoxWidth() * FReal(idxX),
-					       loader->getBoxWidth() * FReal(idxY),
-					       loader->getBoxWidth() * FReal(idxZ));
-		    //							std::cout <<" ( "<< idxX<<" , "<<idxY << " , "<< idxZ << " ) "<< offset <<std::endl;
-		    for(FSize idxSource = 0 ; idxSource < loader->getNumberOfParticles() ; ++idxSource){
-		      EwalParticle<FReal> source = particles[idxSource];
-		      source.position += offset;
-		      //								std::cout << "Part "<<idxSource<< " " <<source.position.getX()<< " " << source.position.getY()<< " " <<source.position.getZ()<< " " <<source.physicalValue <<std::endl ;
-		      FP2P::NonMutualParticles(
-					       part.position.getX(), part.position.getY(),part.position.getZ(),part.physicalValue,
-					       &part.forces[0],&part.forces[1],&part.forces[2],&part.potential,
-					       source.position.getX(), source.position.getY(),source.position.getZ(),source.physicalValue,
-					       &MatrixKernel);
-		    }
-		    //						std::cout <<std::endl<<std::endl<<std::endl;
-		  }
-		}
-	      } // END
-	    }
-	    if(scale){
-	      //
-	      // remove dipole correction for DL_POLY
-	      //
-	      part.forces[0] -= 2.0*part.physicalValue*coeffCorrectionDLPOLY*electricMoment.getX() ;
-	      part.forces[1] -= 2.0*part.physicalValue*coeffCorrectionDLPOLY*electricMoment.getY() ;
-	      part.forces[2] -= 2.0*part.physicalValue*coeffCorrectionDLPOLY*electricMoment.getZ() ;
-	      //
-	      //  Scaling
-	      //
-	      part.forces[0] *= scaleForce ;
-	      part.forces[1] *= scaleForce ;
-	      part.forces[2] *= scaleForce ;
-	    }
-	    if(FParameters::existParameter(argc, argv, FParameterDefinitions::EnabledVerbose.options)){
-	      std::cout << ">> index " << particles[idxTarget].index << std::endl;
-	      std::cout << "Good    x " << particles[idxTarget].position.getX() << " y " << particles[idxTarget].position.getY() << " z " << particles[idxTarget].position.getZ() << std::endl;
-	      std::cout << "Good    fx " <<particles[idxTarget].forces[0] << " fy " << particles[idxTarget].forces[1] << " fz " << particles[idxTarget].forces[2] << std::endl;
-	      std::cout << "DIRECT  fx " << part.forces[0]<< " fy " <<part.forces[1] << " fz " << part.forces[2] << std::endl;
-	      std::cout << "ratio  fx " << particles[idxTarget].forces[0]/part.forces[0] << " fy " <<particles[idxTarget].forces[1]/part.forces[1] << " fz " << particles[idxTarget].forces[2]/part.forces[2] << std::endl;
-	      std::cout << "DIRECT  physical value " << part.physicalValue << " potential " << part.potential << std::endl;
-	
-	      std::cout << "\n";
-	    }
-	    //
-	    dfx.add(particles[idxTarget].forces[0],part.forces[0]);
-	    dfy.add(particles[idxTarget].forces[1],part.forces[1]);
-	    dfz.add(particles[idxTarget].forces[2],part.forces[2]);
-      
-	    denergy += part.potential * part.physicalValue;
-      
-	    //	particlesDirect[idxTarget] = part;
-	  }
-	  counter.tac();
-	  denergy *= 0.5*scaleEnergy ;
-	  if(scale){
-	    denergy -= coeffCorrectionDLPOLY*dipoleNorm*scaleEnergy ;
-	  }
-	  directEnergy	= denergy ;
-	  printf("Difference between direct and Ewald (DL_POLY)\n");
-	  printf("DirectEwald Fx diff is = \n");
-	  printf("DirectEwald   L2Norm  %e\n",dfx.getRelativeL2Norm() );
-	  printf("DirectEwald   InfNorm %e\n",dfx.getRelativeInfNorm());
-	  printf("DirectEwaldFy diff is = \n");
-	  printf("DirectEwald   L2Norm  %e\n",dfx.getRelativeL2Norm());
-	  printf("DirectEwald   InfNorm %e\n",dfy.getRelativeInfNorm());
-	  printf("DirectEwaldFz diff is = \n");
-	  printf("DirectEwald   L2Norm  %e\n",dfx.getRelativeL2Norm());
-	  printf("DirectEwald   InfNorm %e\n",dfz.getRelativeInfNorm());
-	  double L2error = (dfx.getRelativeL2Norm()*dfx.getRelativeL2Norm() + dfy.getRelativeL2Norm()*dfy.getRelativeL2Norm()  + dfz.getRelativeL2Norm() *dfz.getRelativeL2Norm()  );
-	  printf("DirectEwald L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
-	  printf("DirectEwald Energy DIRECT=   %.12e\n",denergy);
-	  printf("DirectEwald Energy EWALD =   %.12e\n",loader->getEnergy());
-	  printf("DirectEwald Energy EWALD - Energy DIRECT              = %e\n", loader->getEnergy()-denergy );
-	  printf("DirectEwald |Energy EWALD - Energy DIRECT|/directEnergy= %e\n",FMath::Abs((loader->getEnergy()-denergy)/loader->getEnergy()));
-
-	  //
-	  if(FParameters::existParameter(argc, argv, "-saveError")){
-	    errorfile << std::endl << "      END DIRECT " << std::endl ;
-	  }
-	  std::cout << "Done  " << "(@DIRECT Algorithm = " << counter.elapsed() << " s)." << std::endl;
-
-	  std::cout << "\n"<< "END DIRECT "
-		    << "-------------------------------------------------------------------------" << std::endl << std::endl ;
-	}
-	// ----------------------------------------------------------------------------------------------------------
-	//                                  DIRECT -- FMM  Comparisons
-	//                                  Ewald  -- FMM  Comparisons
-	// ----------------------------------------------------------------------------------------------------------
-	{
-	  // particles       : Ewald results
-	  // particlesDirect : Direct results
-	  //
-
-	  FReal energy = 0.0;
-	  FMath::FAccurater<FReal> fx, fy, fz, fmmdfx, fmmdfy, fmmdfz,fmmpot;
-
-	  tree.forEachLeaf([&](LeafClass* leaf){
-	      const FReal*const potentials = leaf->getTargets()->getPotentials();
-	      FReal*const forcesX = leaf->getTargets()->getForcesX();
-	      FReal*const forcesY = leaf->getTargets()->getForcesY();
-	      FReal*const forcesZ = leaf->getTargets()->getForcesZ();
-	      const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
-	      const FReal*const positionsX = leaf->getTargets()->getPositions()[0];
-	      const FReal*const positionsY = leaf->getTargets()->getPositions()[1];
-	      const FReal*const positionsZ = leaf->getTargets()->getPositions()[2];
-	      const FSize nbParticlesInLeaf  = leaf->getTargets()->getNbParticles();
-	      const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
-	      //
-	      for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
-		const FSize indexPartOrig = indexes[idxPart];
-		if(scale){
-		  //
-		  // remove dipole correction for DL_POLY
-		  //
-		  forcesX[idxPart] -= 2.0*physicalValues[idxPart]*coeffCorrectionDLPOLY*electricMoment.getX() ;
-		  forcesY[idxPart] -= 2.0*physicalValues[idxPart]*coeffCorrectionDLPOLY*electricMoment.getY() ;
-		  forcesZ[idxPart] -= 2.0*physicalValues[idxPart]*coeffCorrectionDLPOLY*electricMoment.getZ() ;
-		  //
-		  forcesX[idxPart] *= scaleForce;
-		  forcesY[idxPart] *= scaleForce;
-		  forcesZ[idxPart] *= scaleForce;
-		}
-		if(direct) {
-		  //					std::cout << "Good    x " << particlesDirect[idxPart].position.getX() << " y " << particlesDirect[idxPart].position.getY() << " z " << particlesDirect[idxPart].position.getZ() << std::endl;
-		  //					std::cout << "Direct fx " <<particlesDirect[idxPart].forces[0]<< " fy " << particlesDirect[idxPart].forces[1] << " fz " << particlesDirect[idxPart].forces[2] << std::endl;
-		  //					std::cout << "FMM  fx " << forcesX[idxPart] << " fy " <<forcesY[idxPart] << " fz " << forcesZ[idxPart] << std::endl;
-
-		  fmmdfx.add(particlesDirect[indexPartOrig].forces[0],forcesX[idxPart]);
-		  fmmdfy.add(particlesDirect[indexPartOrig].forces[1],forcesY[idxPart]);
-		  fmmdfz.add(particlesDirect[indexPartOrig].forces[2],forcesZ[idxPart]);
-		  fmmpot.add(particlesDirect[indexPartOrig].potential,potentials[idxPart]);
-		} else {
-		  fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]); // Ewald - FMM
-		  fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
-		  fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
-		}
-		energy +=  potentials[idxPart]*physicalValues[idxPart];
-		//
-                if(FParameters::existParameter(argc, argv, FParameterDefinitions::EnabledVerbose.options)){
-		  std::cout << ">> index " << particles[indexPartOrig].index << std::endl;
-		  std::cout << "Good x " << particles[indexPartOrig].position.getX() << " y " << particles[indexPartOrig].position.getY() << " z " << particles[indexPartOrig].position.getZ() << std::endl;
-		  std::cout << std::fixed  << std::setprecision(5) ;
-
-		  std::cout << "FMM  x " << positionsX[idxPart] << " y " << positionsY[idxPart] << " z " << positionsZ[idxPart] << std::endl;
-		  std::cout << "Good fx " <<particles[indexPartOrig].forces[0] << " fy " << particles[indexPartOrig].forces[1] << " fz " << particles[indexPartOrig].forces[2] << std::endl;
-		  std::cout << "FMM  fx " << forcesX[idxPart] << " fy " <<forcesY[idxPart] << " fz " << forcesZ[idxPart] << std::endl;
-		  std::cout << std::scientific  << std::setprecision(5) ;
-		  std::cout << "Diff  fx " << particles[indexPartOrig].forces[0]-forcesX[idxPart] << " fy " <<particles[indexPartOrig].forces[1]-forcesY[idxPart] << " fz " << particles[indexPartOrig].forces[2]-forcesZ[idxPart] << std::endl;
-		  //					std::cout << "GOOD physical value " << particles[indexPartOrig].physicalValue << " potential " << particles[indexPartOrig].potential << std::endl;
-		  //					std::cout << "FMM  physical value " << physicalValues[idxPart] << " potential " << potentials[idxPart] <<  " energy cumul " << energy<<std::endl;
-		  std::cout << std::fixed  << std::setprecision(5) ;
-		  std::cout << "\n";
-		}
-		if(FParameters::existParameter(argc, argv, "-saveError")){
-		  double ratio,tmp ;
-		  ratio = std::fabs(1-particles[indexPartOrig].forces[0]/forcesX[idxPart]) ;
-		  tmp    = std::fabs(1-particles[indexPartOrig].forces[1]/forcesY[idxPart]) ;
-		  ratio = std::max(tmp, ratio) ;
-		  ratio = std::max(std::fabs(1-particles[indexPartOrig].forces[2]/forcesZ[idxPart]) , ratio) ;
-		  if(ratio >=0.25) {
-		    errorfile << ">> index " << particles[indexPartOrig].index << std::endl;
-		    errorfile << "Good x " << particles[indexPartOrig].position.getX() << " y " << particles[indexPartOrig].position.getY() << " z " << particles[indexPartOrig].position.getZ() << std::endl;
-		    errorfile << "FMM  x " << positionsX[idxPart] << " y " << positionsY[idxPart] << " z " << positionsZ[idxPart] << std::endl;
-		    errorfile << "Good fx " <<particles[indexPartOrig].forces[0] << " fy " << particles[indexPartOrig].forces[1] << " fz " << particles[indexPartOrig].forces[2] << std::endl;
-		    errorfile << "FMM  fx " << forcesX[idxPart] << " fy " <<forcesY[idxPart] << " fz " << forcesZ[idxPart] << std::endl;
-		    errorfile << "ratio  fx " << particles[indexPartOrig].forces[0]/forcesX[idxPart] << " fy " <<particles[indexPartOrig].forces[1]/forcesY[idxPart] << " fz " << particles[indexPartOrig].forces[2]/forcesZ[idxPart] << std::endl;
-		    errorfile << "GOOD physical value " << particles[indexPartOrig].physicalValue << " potential " << particles[indexPartOrig].potential << std::endl;
-		    errorfile << "FMM  physical value " << physicalValues[idxPart] << " potential " << potentials[idxPart] << std::endl;
-		    errorfile << "\n";
-		  }
-		}
-	      }
-	    });
-	  energy *= 0.5*scaleEnergy ;
-	  if(scale){
-	    energy -= coeffCorrectionDLPOLY*dipoleNorm*scaleEnergy ;
-	  }
-	  //    printf("(Energy EWALD - Energy FMM)/dipoleNorm = %e\n",(loader->getEnergy()-energy)*volume/(dipoleNorm));
-	  //    printf("(dipoleNorm /Volume = %e\n",correctEnergy);
-	  //
-	  if(direct) {
-	    printf("Difference between FMM and DiRECT:\n");
-	    printf("DirectFmm Fx diff is = \n");
-	    printf("DirectFmm   L2Norm  %e\n",fmmdfx.getRelativeL2Norm());
-	    printf("DirectFmm   InfNorm %e\n",fmmdfx.getRelativeInfNorm());
-	    printf("DirectFmm Fy diff is = \n");
-	    printf("DirectFmm   L2Norm  %e\n",fmmdfy.getRelativeL2Norm());
-	    printf("DirectFmm   InfNorm %e\n",fmmdfy.getRelativeInfNorm());
-	    printf("DirectFmm Fz diff is = \n");
-	    printf("DirectFmm   L2Norm  %e\n",fmmdfz.getRelativeL2Norm());
-	    printf("DirectFmm   InfNorm %e\n",fmmdfz.getRelativeInfNorm());
-	    double L2error = (fmmdfx.getRelativeL2Norm()*fmmdfx.getRelativeL2Norm() + fmmdfy.getRelativeL2Norm()*fmmdfy.getRelativeL2Norm()  + fmmdfz.getRelativeL2Norm() *fmmdfz.getRelativeL2Norm()  );
-	    printf("DirectFmm L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
-	    printf("DirectFmm L2 potential Error= %e\n",fmmpot.getRelativeL2Norm()) ;
-
-	    //
-	    printf("DirectFmm Energy FMM=    %.12e\n",energy);
-	    printf("DirectFmm Energy DIRECT= %.12e\n",directEnergy);
-	    printf("DirectFmm |Energy DIRECT - Energy DIRECT|/directEnergy= %e\n",FMath::Abs((directEnergy-energy)/directEnergy));
-	  }
-	  else {
-	    printf("FmmEwald Difference between FMM and Ewald DL_POLY:\n");
-	    printf("FmmEwald Fx diff is = \n");
-	    printf("FmmEwald   L2Norm  %e\n",fx.getRelativeL2Norm());
-	    printf("FmmEwald   InfNorm %e\n",fx.getRelativeInfNorm());
-	    printf("FmmEwald Fy diff is = \n");
-	    printf("FmmEwald   L2Norm  %e\n",fy.getRelativeL2Norm());
-	    printf("FmmEwald   InfNorm %e\n",fy.getInfNorm());
-	    printf("FmmEwald Fz diff is = \n");
-	    printf("FmmEwald   L2Norm  %e\n",fz.getRelativeL2Norm());
-	    printf("FmmEwald   InfNorm %e\n",fz.getRelativeInfNorm());
-	    //
-	    double L2error = (fx.getL2Norm()*fx.getL2Norm() + fy.getL2Norm()*fy.getL2Norm()  + fz.getL2Norm() *fz.getL2Norm()  );
-	    printf("FmmEwald RMS Force Error= %e\n",FMath::Sqrt(L2error/static_cast<double>(loader->getNumberOfParticles()))) ;
-	    //
-	    printf("FmmEwald Energy FMM=   %.12e\n",energy);
-	    printf("FmmEwald Energy EWALD= %.12e\n",loader->getEnergy());
-	    printf("FmmEwald Energy EWALD - Energy FMM = %e\n",loader->getEnergy()-energy);
-	    printf("FmmEwald |Energy EWALD -Energy FMM|/Energy EWALD= %e\n",FMath::Abs((loader->getEnergy()-energy)/loader->getEnergy()));
-
-	  }
-	  printf("vtkParticles\n");
-
-	  vtkfile.writeParticules(tree, true,1,1,"vtkParticles","FMM results");
-
-	}
-	//
-	// ----------------------------------------------------------------------------
-	//
-	delete[] particles;
-	delete[] particlesDirect;
-
-	return 0;
+  //
+  // direct computation
+  //
+  bool direct = false ;
+  if(FParameters::existParameter(argc, argv, "-direct")){
+      direct = true ;
+      printf("Compute direct:\n");
+      printf("Box [%d;%d][%d;%d][%d;%d]\n", min.getX(), max.getX(), min.getY(),
+             max.getY(), min.getZ(), max.getZ());
+
+      particlesDirect = new EwalParticle<FReal>[loader->getNumberOfParticles()];
+
+      FReal denergy = 0.0;
+      FMath::FAccurater<FReal> dfx, dfy, dfz ;
+      //
+      //
+      // particles       : Ewald results
+      // particlesDirect : Direct results
+      //
+      counter.tic();
+      for(FSize idxTarget = 0 ; idxTarget < loader->getNumberOfParticles() ; ++idxTarget){
+          particlesDirect[idxTarget] = particles[idxTarget];
+          EwalParticle<FReal> & part        = particlesDirect[idxTarget];
+          part.forces[0] = part.forces[1] = part.forces[2] = 0.0;
+          part.potential = 0.0;
+          //
+          // compute with all other except itself
+          //
+          // Compute force and potential between  particles[idxTarget] and particles inside the box
+          //
+          for(FSize idxOther = 0; idxOther < loader->getNumberOfParticles() ; ++idxOther){
+              if( idxOther != idxTarget ){
+                  FP2P::NonMutualParticles(
+                        part.position.getX(), part.position.getY(),
+                        part.position.getZ(),part.physicalValue,
+                        &part.forces[0],&part.forces[1],
+                      &part.forces[2],&part.potential,
+                      particles[idxOther].position.getX(), particles[idxOther].position.getY(),
+                      particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
+                      &MatrixKernel);
+                }
+            }
+          //
+          // compute with image boxes
+          //
+          for(int idxX = min.getX() ; idxX <= max.getX() ; ++idxX){
+              for(int idxY = min.getY() ; idxY <= max.getY() ; ++idxY){
+                  for(int idxZ = min.getZ() ; idxZ <= max.getZ() ; ++idxZ){
+
+                      {
+                        //int nL = PeriodicDeep ;
+                        //	for(int idxX = -nL ; idxX <= nL -1; ++idxX){
+                        //	  for(int idxY = -nL ; idxY <=nL-1 ; ++idxY){
+                        //	    for(int idxZ = -nL ; idxZ <= nL-1; ++idxZ){
+                        if(idxX == 0 && idxY == 0 && idxZ == 0) continue;
+
+                        const FPoint<FReal> offset(loader->getBoxWidth() * FReal(idxX),
+                                                   loader->getBoxWidth() * FReal(idxY),
+                                                   loader->getBoxWidth() * FReal(idxZ));
+                        //							std::cout <<" ( "<< idxX<<" , "<<idxY << " , "<< idxZ << " ) "<< offset <<std::endl;
+                        for(FSize idxSource = 0 ; idxSource < loader->getNumberOfParticles() ; ++idxSource){
+                            EwalParticle<FReal> source = particles[idxSource];
+                            source.position += offset;
+                            //								std::cout << "Part "<<idxSource<< " " <<source.position.getX()<< " " << source.position.getY()<< " " <<source.position.getZ()<< " " <<source.physicalValue <<std::endl ;
+                            FP2P::NonMutualParticles(
+                                  part.position.getX(), part.position.getY(),part.position.getZ(),part.physicalValue,
+                                  &part.forces[0],&part.forces[1],&part.forces[2],&part.potential,
+                                source.position.getX(), source.position.getY(),source.position.getZ(),source.physicalValue,
+                                &MatrixKernel);
+                          }
+                        //						std::cout <<std::endl<<std::endl<<std::endl;
+                      }
+                    }
+                } // END
+            }
+          if(scale){
+              //
+              // remove dipole correction for DL_POLY
+              //
+              part.forces[0] -= 2.0*part.physicalValue*coeffCorrectionDLPOLY*electricMoment.getX() ;
+              part.forces[1] -= 2.0*part.physicalValue*coeffCorrectionDLPOLY*electricMoment.getY() ;
+              part.forces[2] -= 2.0*part.physicalValue*coeffCorrectionDLPOLY*electricMoment.getZ() ;
+              //
+              //  Scaling
+              //
+              part.forces[0] *= scaleForce ;
+              part.forces[1] *= scaleForce ;
+              part.forces[2] *= scaleForce ;
+            }
+          if(FParameters::existParameter(argc, argv, FParameterDefinitions::EnabledVerbose.options)){
+              std::cout << ">> index " << particles[idxTarget].index << std::endl;
+              std::cout << "Good    x " << particles[idxTarget].position.getX() << " y " << particles[idxTarget].position.getY() << " z " << particles[idxTarget].position.getZ() << std::endl;
+              std::cout << "Good    fx " <<particles[idxTarget].forces[0] << " fy " << particles[idxTarget].forces[1] << " fz " << particles[idxTarget].forces[2] << std::endl;
+              std::cout << "DIRECT  fx " << part.forces[0]<< " fy " <<part.forces[1] << " fz " << part.forces[2] << std::endl;
+              std::cout << "ratio  fx " << particles[idxTarget].forces[0]/part.forces[0] << " fy " <<particles[idxTarget].forces[1]/part.forces[1] << " fz " << particles[idxTarget].forces[2]/part.forces[2] << std::endl;
+              std::cout << "DIRECT  physical value " << part.physicalValue << " potential " << part.potential << std::endl;
+
+              std::cout << "\n";
+            }
+          //
+          dfx.add(particles[idxTarget].forces[0],part.forces[0]);
+          dfy.add(particles[idxTarget].forces[1],part.forces[1]);
+          dfz.add(particles[idxTarget].forces[2],part.forces[2]);
+
+          denergy += part.potential * part.physicalValue;
+
+          //	particlesDirect[idxTarget] = part;
+        }
+      counter.tac();
+      denergy *= 0.5*scaleEnergy ;
+      if(scale){
+          denergy -= coeffCorrectionDLPOLY*dipoleNorm*scaleEnergy ;
+        }
+      directEnergy	= denergy ;
+      printf("Difference between direct and Ewald (DL_POLY)\n");
+      printf("DirectEwald Fx diff is = \n");
+      printf("DirectEwald   L2Norm  %e\n",dfx.getRelativeL2Norm() );
+      printf("DirectEwald   InfNorm %e\n",dfx.getRelativeInfNorm());
+      printf("DirectEwaldFy diff is = \n");
+      printf("DirectEwald   L2Norm  %e\n",dfx.getRelativeL2Norm());
+      printf("DirectEwald   InfNorm %e\n",dfy.getRelativeInfNorm());
+      printf("DirectEwaldFz diff is = \n");
+      printf("DirectEwald   L2Norm  %e\n",dfx.getRelativeL2Norm());
+      printf("DirectEwald   InfNorm %e\n",dfz.getRelativeInfNorm());
+      double L2error = (dfx.getRelativeL2Norm()*dfx.getRelativeL2Norm() + dfy.getRelativeL2Norm()*dfy.getRelativeL2Norm()  + dfz.getRelativeL2Norm() *dfz.getRelativeL2Norm()  );
+      printf("DirectEwald L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
+      printf("DirectEwald Energy DIRECT=   %.12e\n",denergy);
+      printf("DirectEwald Energy EWALD =   %.12e\n",loader->getEnergy());
+      printf("DirectEwald Energy EWALD - Energy DIRECT              = %e\n", loader->getEnergy()-denergy );
+      printf("DirectEwald |Energy EWALD - Energy DIRECT|/directEnergy= %e\n",FMath::Abs((loader->getEnergy()-denergy)/loader->getEnergy()));
+
+      //
+      if(FParameters::existParameter(argc, argv, "-saveError")){
+          errorfile << std::endl << "      END DIRECT " << std::endl ;
+        }
+      std::cout << "Done  " << "(@DIRECT Algorithm = " << counter.elapsed() << " s)." << std::endl;
+
+      std::cout << "\n"<< "END DIRECT "
+                << "-------------------------------------------------------------------------" << std::endl << std::endl ;
+    }
+  // ----------------------------------------------------------------------------------------------------------
+  //                                  DIRECT -- FMM  Comparisons
+  //                                  Ewald  -- FMM  Comparisons
+  // ----------------------------------------------------------------------------------------------------------
+  {
+    // particles       : Ewald results
+    // particlesDirect : Direct results
+    //
+
+    FReal energy = 0.0;
+    FMath::FAccurater<FReal> fx, fy, fz, fmmdfx, fmmdfy, fmmdfz,fmmpot;
+
+    tree.forEachLeaf([&](LeafClass* leaf){
+      const FReal*const potentials = leaf->getTargets()->getPotentials();
+      FReal*const forcesX = leaf->getTargets()->getForcesX();
+      FReal*const forcesY = leaf->getTargets()->getForcesY();
+      FReal*const forcesZ = leaf->getTargets()->getForcesZ();
+      const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
+      const FReal*const positionsX = leaf->getTargets()->getPositions()[0];
+      const FReal*const positionsY = leaf->getTargets()->getPositions()[1];
+      const FReal*const positionsZ = leaf->getTargets()->getPositions()[2];
+      const FSize nbParticlesInLeaf  = leaf->getTargets()->getNbParticles();
+      const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
+      //
+      for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
+          const FSize indexPartOrig = indexes[idxPart];
+          if(scale){
+              //
+              // remove dipole correction for DL_POLY
+              //
+              forcesX[idxPart] -= 2.0*physicalValues[idxPart]*coeffCorrectionDLPOLY*electricMoment.getX() ;
+              forcesY[idxPart] -= 2.0*physicalValues[idxPart]*coeffCorrectionDLPOLY*electricMoment.getY() ;
+              forcesZ[idxPart] -= 2.0*physicalValues[idxPart]*coeffCorrectionDLPOLY*electricMoment.getZ() ;
+              //
+              forcesX[idxPart] *= scaleForce;
+              forcesY[idxPart] *= scaleForce;
+              forcesZ[idxPart] *= scaleForce;
+            }
+          if(direct) {
+              //					std::cout << "Good    x " << particlesDirect[idxPart].position.getX() << " y " << particlesDirect[idxPart].position.getY() << " z " << particlesDirect[idxPart].position.getZ() << std::endl;
+              //					std::cout << "Direct fx " <<particlesDirect[idxPart].forces[0]<< " fy " << particlesDirect[idxPart].forces[1] << " fz " << particlesDirect[idxPart].forces[2] << std::endl;
+              //					std::cout << "FMM  fx " << forcesX[idxPart] << " fy " <<forcesY[idxPart] << " fz " << forcesZ[idxPart] << std::endl;
+
+              fmmdfx.add(particlesDirect[indexPartOrig].forces[0],forcesX[idxPart]);
+              fmmdfy.add(particlesDirect[indexPartOrig].forces[1],forcesY[idxPart]);
+              fmmdfz.add(particlesDirect[indexPartOrig].forces[2],forcesZ[idxPart]);
+              fmmpot.add(particlesDirect[indexPartOrig].potential,potentials[idxPart]);
+            } else {
+              fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]); // Ewald - FMM
+              fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
+              fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
+            }
+          energy +=  potentials[idxPart]*physicalValues[idxPart];
+          //
+          if(FParameters::existParameter(argc, argv, FParameterDefinitions::EnabledVerbose.options)){
+              std::cout << ">> index " << particles[indexPartOrig].index << std::endl;
+              std::cout << "Good x " << particles[indexPartOrig].position.getX() << " y " << particles[indexPartOrig].position.getY() << " z " << particles[indexPartOrig].position.getZ() << std::endl;
+              std::cout << std::fixed  << std::setprecision(5) ;
+
+              std::cout << "FMM  x " << positionsX[idxPart] << " y " << positionsY[idxPart] << " z " << positionsZ[idxPart] << std::endl;
+              std::cout << "Good fx " <<particles[indexPartOrig].forces[0] << " fy " << particles[indexPartOrig].forces[1] << " fz " << particles[indexPartOrig].forces[2] << std::endl;
+              std::cout << "FMM  fx " << forcesX[idxPart] << " fy " <<forcesY[idxPart] << " fz " << forcesZ[idxPart] << std::endl;
+              std::cout << std::scientific  << std::setprecision(5) ;
+              std::cout << "Diff  fx " << particles[indexPartOrig].forces[0]-forcesX[idxPart] << " fy " <<particles[indexPartOrig].forces[1]-forcesY[idxPart] << " fz " << particles[indexPartOrig].forces[2]-forcesZ[idxPart] << std::endl;
+              //					std::cout << "GOOD physical value " << particles[indexPartOrig].physicalValue << " potential " << particles[indexPartOrig].potential << std::endl;
+              //					std::cout << "FMM  physical value " << physicalValues[idxPart] << " potential " << potentials[idxPart] <<  " energy cumul " << energy<<std::endl;
+              std::cout << std::fixed  << std::setprecision(5) ;
+              std::cout << "\n";
+            }
+          if(FParameters::existParameter(argc, argv, "-saveError")){
+              double ratio,tmp ;
+              ratio = std::fabs(1-particles[indexPartOrig].forces[0]/forcesX[idxPart]) ;
+              tmp    = std::fabs(1-particles[indexPartOrig].forces[1]/forcesY[idxPart]) ;
+              ratio = std::max(tmp, ratio) ;
+              ratio = std::max(std::fabs(1-particles[indexPartOrig].forces[2]/forcesZ[idxPart]) , ratio) ;
+              if(ratio >=0.25) {
+                  errorfile << ">> index " << particles[indexPartOrig].index << std::endl;
+                  errorfile << "Good x " << particles[indexPartOrig].position.getX() << " y " << particles[indexPartOrig].position.getY() << " z " << particles[indexPartOrig].position.getZ() << std::endl;
+                  errorfile << "FMM  x " << positionsX[idxPart] << " y " << positionsY[idxPart] << " z " << positionsZ[idxPart] << std::endl;
+                  errorfile << "Good fx " <<particles[indexPartOrig].forces[0] << " fy " << particles[indexPartOrig].forces[1] << " fz " << particles[indexPartOrig].forces[2] << std::endl;
+                  errorfile << "FMM  fx " << forcesX[idxPart] << " fy " <<forcesY[idxPart] << " fz " << forcesZ[idxPart] << std::endl;
+                  errorfile << "ratio  fx " << particles[indexPartOrig].forces[0]/forcesX[idxPart] << " fy " <<particles[indexPartOrig].forces[1]/forcesY[idxPart] << " fz " << particles[indexPartOrig].forces[2]/forcesZ[idxPart] << std::endl;
+                  errorfile << "GOOD physical value " << particles[indexPartOrig].physicalValue << " potential " << particles[indexPartOrig].potential << std::endl;
+                  errorfile << "FMM  physical value " << physicalValues[idxPart] << " potential " << potentials[idxPart] << std::endl;
+                  errorfile << "\n";
+                }
+            }
+        }
+    });
+    energy *= 0.5*scaleEnergy ;
+    std::cout << " Energy before correction " << energy << " Correction: " <<coeffCorrectionDLPOLY*dipoleNorm*scaleEnergy <<std::endl;
+    if(scale){
+        energy -= coeffCorrectionDLPOLY*dipoleNorm*scaleEnergy ;
+      }
+    //    printf("(Energy EWALD - Energy FMM)/dipoleNorm = %e\n",(loader->getEnergy()-energy)*volume/(dipoleNorm));
+    //    printf("(dipoleNorm /Volume = %e\n",correctEnergy);
+    //
+    if(direct) {
+        printf("Difference between FMM and DiRECT:\n");
+        printf("DirectFmm Fx diff is = \n");
+        printf("DirectFmm   L2Norm  %e\n",fmmdfx.getRelativeL2Norm());
+        printf("DirectFmm   InfNorm %e\n",fmmdfx.getRelativeInfNorm());
+        printf("DirectFmm Fy diff is = \n");
+        printf("DirectFmm   L2Norm  %e\n",fmmdfy.getRelativeL2Norm());
+        printf("DirectFmm   InfNorm %e\n",fmmdfy.getRelativeInfNorm());
+        printf("DirectFmm Fz diff is = \n");
+        printf("DirectFmm   L2Norm  %e\n",fmmdfz.getRelativeL2Norm());
+        printf("DirectFmm   InfNorm %e\n",fmmdfz.getRelativeInfNorm());
+        double L2error = (fmmdfx.getRelativeL2Norm()*fmmdfx.getRelativeL2Norm() + fmmdfy.getRelativeL2Norm()*fmmdfy.getRelativeL2Norm()  + fmmdfz.getRelativeL2Norm() *fmmdfz.getRelativeL2Norm()  );
+        printf("DirectFmm L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
+        printf("DirectFmm L2 potential Error= %e\n",fmmpot.getRelativeL2Norm()) ;
+
+        //
+        printf("DirectFmm Energy FMM=    %.12e\n",energy);
+        printf("DirectFmm Energy DIRECT= %.12e\n",directEnergy);
+        printf("DirectFmm |Energy DIRECT - Energy DIRECT|/directEnergy= %e\n",FMath::Abs((directEnergy-energy)/directEnergy));
+      }
+    else {
+        printf("FmmEwald Difference between FMM and Ewald DL_POLY:\n");
+        printf("FmmEwald Fx diff is = \n");
+        printf("FmmEwald   L2Norm  %e\n",fx.getRelativeL2Norm());
+        printf("FmmEwald   InfNorm %e\n",fx.getRelativeInfNorm());
+        printf("FmmEwald Fy diff is = \n");
+        printf("FmmEwald   L2Norm  %e\n",fy.getRelativeL2Norm());
+        printf("FmmEwald   InfNorm %e\n",fy.getInfNorm());
+        printf("FmmEwald Fz diff is = \n");
+        printf("FmmEwald   L2Norm  %e\n",fz.getRelativeL2Norm());
+        printf("FmmEwald   InfNorm %e\n",fz.getRelativeInfNorm());
+        //
+        double L2error = (fx.getL2Norm()*fx.getL2Norm() + fy.getL2Norm()*fy.getL2Norm()  + fz.getL2Norm() *fz.getL2Norm()  );
+        printf("FmmEwald RMS Force Error= %e\n",FMath::Sqrt(L2error/static_cast<double>(loader->getNumberOfParticles()))) ;
+        //
+        printf("FmmEwald Energy FMM=   %.12e\n",energy);
+        printf("FmmEwald Energy EWALD= %.12e\n",loader->getEnergy());
+        printf("FmmEwald Energy EWALD - Energy FMM = %e\n",loader->getEnergy()-energy);
+        printf("FmmEwald |Energy EWALD -Energy FMM|/Energy EWALD= %e\n",FMath::Abs((loader->getEnergy()-energy)/loader->getEnergy()));
+
+      }
+    printf("vtkParticles\n");
+
+    vtkfile.writeParticules(tree, true,1,1,"vtkParticles","FMM results");
+
+  }
+  //
+  // ----------------------------------------------------------------------------
+  //
+  delete[] particles;
+  delete[] particlesDirect;
+
+  return 0;
 }
 
 
diff --git a/Tests/Utils/testFmmAlgorithmPeriodic.cpp b/Tests/Utils/testFmmAlgorithmPeriodic.cpp
index 34e91209d..710e9a21d 100644
--- a/Tests/Utils/testFmmAlgorithmPeriodic.cpp
+++ b/Tests/Utils/testFmmAlgorithmPeriodic.cpp
@@ -6,28 +6,28 @@
 #include <cstdio>
 #include <cstdlib>
 
-#include "../../Src/Utils/FParameters.hpp"
-#include "../../Src/Utils/FTic.hpp"
+#include "Utils/FParameters.hpp"
+#include "Utils/FTic.hpp"
 
-#include "../../Src/Files/FRandomLoader.hpp"
+#include "Files/FRandomLoader.hpp"
 
-#include "../../Src/Files/FPerLeafLoader.hpp"
+#include "Files/FPerLeafLoader.hpp"
 
-#include "../../Src/Containers/FOctree.hpp"
-#include "../../Src/Containers/FVector.hpp"
+#include "Containers/FOctree.hpp"
+#include "Containers/FVector.hpp"
 
-#include "../../Src/Components/FSimpleLeaf.hpp"
+#include "Components/FSimpleLeaf.hpp"
 
-#include "../../Src/Components/FTestParticleContainer.hpp"
+#include "Components/FTestParticleContainer.hpp"
 
-#include "../../Src/Utils/FPoint.hpp"
+#include "Utils/FPoint.hpp"
 
-#include "../../Src/Components/FTestCell.hpp"
-#include "../../Src/Components/FTestKernels.hpp"
+#include "Components/FTestCell.hpp"
+#include "Components/FTestKernels.hpp"
 
-#include "../../Src/Core/FFmmAlgorithmPeriodic.hpp"
+#include "Core/FFmmAlgorithmPeriodic.hpp"
 
-#include "../../Src/Utils/FParameterNames.hpp"
+#include "Utils/FParameterNames.hpp"
 
 /** This program show an example of use of
   * the fmm basic algo
-- 
GitLab