diff --git a/ippl/test/particle/p3m3dHeating.cpp b/ippl/test/particle/p3m3dHeating.cpp
index 217251023159cc1d3c540b08ca552429d8832e5c..ca6209f4c63dccb9a3a7779a11a74a313a2920ed 100644
--- a/ippl/test/particle/p3m3dHeating.cpp
+++ b/ippl/test/particle/p3m3dHeating.cpp
@@ -36,10 +36,10 @@
 #include "Particle/PairBuilder/HashPairBuilderPeriodic.h"
 #include "Particle/PairBuilder/HashPairBuilderPeriodicParallel.h"
 #include "Particle/PairBuilder/PairConditions.h"
+#include "Utility/PAssert.h"
 #include "math.h"
 
 #include <random>
-#include <cassert>
 
 #include "VTKFieldWriterParallel.hpp"
 #include "ChargedParticleFactory.hpp"
@@ -477,9 +477,9 @@ public:
 #if defined (NDEBUG)
         (void)h5err;
 #endif
-        assert (h5err != H5_ERR);
+        PAssert (h5err != H5_ERR);
         H5f_m = H5OpenFile (fn.c_str(), H5_O_RDONLY, props);
-        assert (H5f_m != (h5_file_t)H5_ERR);
+        PAssert (H5f_m != (h5_file_t)H5_ERR);
         H5CloseProp (props);
     }
 
diff --git a/ippl/test/particle/p3m3dMicrobunching.cpp b/ippl/test/particle/p3m3dMicrobunching.cpp
index 5251efda61f1c63739cdd381425409c35185bf5e..13cafb213c1fce9088474d4341fa1184ff2ef29e 100644
--- a/ippl/test/particle/p3m3dMicrobunching.cpp
+++ b/ippl/test/particle/p3m3dMicrobunching.cpp
@@ -22,7 +22,6 @@
 // along with OPAL. If not, see <https://www.gnu.org/licenses/>.
 //
 #include "Ippl.h"
-#include <cassert>
 #include <string>
 #include <vector>
 #include <iostream>
@@ -35,6 +34,7 @@
 #include "Particle/PairBuilder/HashPairBuilderPeriodic.h"
 #include "Particle/PairBuilder/HashPairBuilderPeriodicParallel.h"
 #include "Particle/PairBuilder/PairConditions.h"
+#include "Utility/PAssert.h"
 #include "math.h"
 //#include "FixedAlgebra/FMatrix.h"
 
@@ -61,8 +61,8 @@ typedef Field<std::complex<double>, Dim, Mesh_t, Center_t>            CxField_t;
 typedef FFT<CCTransform, Dim, double>                                 FFT_t;
 
 typedef IntCIC                                                        IntrplCIC_t;
-typedef IntNGP                                                        IntrplNGP_t;
-typedef IntTSC                                                        IntrplTSC_t;
+//typedef IntNGP                                                        IntrplNGP_t;
+//typedef IntTSC                                                        IntrplTSC_t;
 
 typedef UniformCartesian<2, double>                                   Mesh2d_t;
 typedef CenteredFieldLayout<2, Mesh2d_t, Center_t>                    FieldLayout2d_t;
@@ -464,7 +464,7 @@ class ChargedParticles : public IpplParticleBase<PL> {
                     h5_prop_t props = H5CreateFileProp ();
                     MPI_Comm comm = Ippl::getComm();
                     h5_err_t h5err = H5SetPropFileMPIOCollective (props, &comm);
-                    assert (h5err != H5_ERR);
+                    PAssert (h5err != H5_ERR);
                     H5f_m = H5OpenFile(fn.c_str(), H5_O_WRONLY, props);
                 }
 
diff --git a/src/Amr/AmrBoxLib.cpp b/src/Amr/AmrBoxLib.cpp
index f5a985c7395f0fd3484a47164edb49108fa545f5..06d1e292e956e33e2742a2a10bb64d8462a8d414 100644
--- a/src/Amr/AmrBoxLib.cpp
+++ b/src/Amr/AmrBoxLib.cpp
@@ -27,9 +27,10 @@
 #include "Algorithms/AmrPartBunch.h"
 #include "Structure/FieldSolver.h"
 #include "Solvers/PoissonSolver.h"
+#include "Utility/PAssert.h"
 
 #include "Amr/AmrYtWriter.h"
-    
+
 #include <AMReX_MultiFabUtil.H>
 
 #include <AMReX_ParmParse.H> // used in initialize function
@@ -721,7 +722,7 @@ void AmrBoxLib::doRegrid_m(int lbase, double time) {
     
     MakeNewGrids(lbase, time, new_finest, new_grids);
 
-    BL_ASSERT(new_finest <= finest_level+1);
+    PAssert(new_finest <= finest_level+1);
     
     for (int lev = lbase+1; lev <= new_finest; ++lev)
     {
diff --git a/src/Amr/BoxLibLayout.hpp b/src/Amr/BoxLibLayout.hpp
index c898753d4f07a595418b9644f66103a487dcaeca..82e4f28fc0d09a9e010f02bed640eb3c9745b121 100644
--- a/src/Amr/BoxLibLayout.hpp
+++ b/src/Amr/BoxLibLayout.hpp
@@ -39,9 +39,9 @@
 #include "BoxLibLayout.h"
 
 #include "Message/Formatter.h"
+#include "Utility/PAssert.h"
 #include "Utilities/OpalException.h"
 
-#include <cassert>
 #include <cmath>
 
 template <class T, unsigned Dim>
@@ -473,7 +473,7 @@ void BoxLibLayout<T, Dim>::buildLevelMask(int lev, const int ncells) {
 
 template <class T, unsigned Dim>
 void BoxLibLayout<T, Dim>::clearLevelMask(int lev) {
-    assert(lev < (int)masks_m.size());
+    PAssert(lev < (int)masks_m.size());
     masks_m[lev].reset(nullptr);
 }
 
@@ -548,9 +548,9 @@ bool BoxLibLayout<T, Dim>::Where(AmrParticleBase< BoxLibLayout<T,Dim> >& p,
     if (lev_max == -1)
         lev_max = finestLevel();
 
-    BL_ASSERT(lev_max <= finestLevel());
+    PAssert(lev_max <= finestLevel());
 
-    BL_ASSERT(nGrow == 0 || (nGrow >= 0 && lev_min == lev_max));
+    PAssert(nGrow == 0 || (nGrow >= 0 && lev_min == lev_max));
 
     std::vector< std::pair<int, AmrBox_t> > isects;
 
@@ -558,7 +558,7 @@ bool BoxLibLayout<T, Dim>::Where(AmrParticleBase< BoxLibLayout<T,Dim> >& p,
     {
         const AmrIntVect_t& iv = Index(p, ip, lev);
         const AmrGrid_t& ba = ParticleBoxArray(lev);
-        BL_ASSERT(ba.ixType().cellCentered());
+        PAssert(ba.ixType().cellCentered());
 
         if (lev == (int)p.Level[ip]) {
             // The fact that we are here means this particle does not belong to any finer grids.
@@ -605,7 +605,7 @@ bool BoxLibLayout<T, Dim>::EnforcePeriodicWhere (AmrParticleBase< BoxLibLayout<T
     if (lev_max == -1)
         lev_max = finestLevel();
 
-    BL_ASSERT(lev_max <= finestLevel());
+    PAssert(lev_max <= finestLevel());
     //
     // Create a copy "dummy" particle to check for periodic outs.
     //
@@ -678,7 +678,7 @@ bool BoxLibLayout<T, Dim>::PeriodicShift (SingleParticlePos_t R) const
                 //
                 R[i] += .125*geom.CellSize(i);
 
-            BL_ASSERT(R[i] >= geom.ProbLo(i));
+            PAssert(R[i] >= geom.ProbLo(i));
 
             shifted = true;
         }
@@ -700,7 +700,7 @@ bool BoxLibLayout<T, Dim>::PeriodicShift (SingleParticlePos_t R) const
                 //
                 R[i] -= .125*geom.CellSize(i);
 
-            BL_ASSERT(R[i] <= geom.ProbHi(i));
+            PAssert(R[i] <= geom.ProbHi(i));
 
             shifted = true;
         }
@@ -810,8 +810,8 @@ void BoxLibLayout<T, Dim>::initBaseBox_m(int nGridPoints,
     AmrDomain_t real_box;
     for (int d = 0; d < AMREX_SPACEDIM; ++d) {
 
-        assert(lowerBound[d] < 0);
-        assert(upperBound[d] > 0);
+        PAssert(lowerBound[d] < 0);
+        PAssert(upperBound[d] > 0);
 
         real_box.setLo(d, lowerBound[d] * (1.0 + dh));
         real_box.setHi(d, upperBound[d] * (1.0 + dh));
diff --git a/src/Classic/Algorithms/PartBins.cpp b/src/Classic/Algorithms/PartBins.cpp
index fd6652b615d389ecf1beec6a8941aa980079f469..1cc6efe7a0980f4f1fa0075854b1137789d19f5c 100644
--- a/src/Classic/Algorithms/PartBins.cpp
+++ b/src/Classic/Algorithms/PartBins.cpp
@@ -1,12 +1,30 @@
+//
+// Class PartBins
+//   Defines a structure to hold energy bins and their associated data
+//
+// Copyright (c) 2007-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
+// All rights reserved
+//
+// This file is part of OPAL.
+//
+// OPAL 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
+
 #include "Algorithms/PartBins.h"
 #include "Algorithms/PBunchDefs.h"
 #include "Physics/Physics.h"
+#include "Utility/Inform.h"
 #include <cfloat>
 #include <limits>
 #include <vector>
 #include <cmath>
 
-
 PartBins::PartBins(int bins, int sbins) :
     gamma_m(1.0),
     bins_m(bins),
@@ -14,17 +32,16 @@ PartBins::PartBins(int bins, int sbins) :
     xmin_m(0.0),
     xmax_m(0.0),
     hBin_m(0.0),
-    nemittedBins_m(0),
-    nemittedSBins_m(0) {
-    
+    nemittedBins_m(0) {
+
     // number of particles in the bins on the local node
     nBin_m        = std::unique_ptr<size_t[]>(new size_t[bins_m]);
     xbinmin_m     = std::unique_ptr<double[]>(new double[bins_m]);
     xbinmax_m     = std::unique_ptr<double[]>(new double[bins_m]);
-    
+
     // flag whether the bin contain particles or not
     binsEmitted_m = std::unique_ptr<bool[]>(new bool[bins_m]);
-    
+
     nDelBin_m     = std::unique_ptr<size_t[]>(new size_t[bins_m]);
 
     for(int i = 0; i < bins_m; i++) {
@@ -57,46 +74,6 @@ size_t PartBins::getTotalNumPerBin(int b) {
     return s;
 }
 
-void PartBins::updateStatus(const int bunchCount, const size_t partInBin) {
-    // array index of binsEmitted_m[] starts from 0
-    // nemittedBins_m and bins_m index starts from 1
-    binsEmitted_m[bunchCount - 1] = true;
-    size_t NpartInBin = partInBin;
-    reduce(NpartInBin, NpartInBin, OpAddAssign());
-    nBin_m[bunchCount - 1] = NpartInBin;
-    nemittedBins_m++;
-}
-
-void PartBins::updateDeletedPartsInBin(size_t countLost[]) {
-    Inform msg("updateDeletedPartsInBin ");
-
-    const int lastEmittedBin = getLastemittedBin();
-    reduce(&(countLost[0]), &(countLost[0]) + lastEmittedBin, &(countLost[0]), OpAddAssign());
-    for(int ii = 0; ii < lastEmittedBin; ii++) {
-        if(countLost[ii] > 0) {
-            nDelBin_m[ii] = countLost[ii];
-            msg << "In Bin: " << ii << ", " << nDelBin_m[ii] << " particle(s) lost" << endl;
-        }
-    }
-}
-
-
-void PartBins::updatePartInBin(size_t countLost[]) {
-
-    Inform msg0("updatePartInBin ");
-
-    for(int ii = 0; ii < nemittedBins_m; ii++) {
-        msg0 << "In Bin: " << ii << ", " << nBin_m[ii] << " particles " << endl;
-    }
-
-    reduce(&(countLost[0]), &(countLost[0]) + nemittedBins_m, &(countLost[0]), OpAddAssign());
-    for(int ii = 0; ii < nemittedBins_m; ii++) {
-        if(countLost[ii] > 0) {
-            nBin_m[ii] -= countLost[ii];
-            msg0 << "In Bin: " << ii << ", " << countLost[ii] << " particle(s) lost" << endl;
-        }
-    }
-}
 
 void PartBins::updatePartInBin_cyc(size_t countLost[]) {
 
@@ -106,16 +83,6 @@ void PartBins::updatePartInBin_cyc(size_t countLost[]) {
   }
 }
 
-
-void PartBins::resetPartInBin(size_t newPartNum[]) {
-    reduce(&(newPartNum[0]), &(newPartNum[0]) + nemittedBins_m, &(newPartNum[0]), OpAddAssign());
-    for(int ii = 0; ii < nemittedBins_m; ii++) {
-        nBin_m[ii] = newPartNum[ii];
-        INFOMSG("After reset Bin: " << ii << ", particle(s): " << newPartNum[ii] << endl);
-    }
-}
-
-
 void PartBins::resetPartInBin_cyc(size_t newPartNum[], int maxbinIndex) {
     reduce(maxbinIndex, maxbinIndex, OpMaxAssign());
     nemittedBins_m =  maxbinIndex + 1;
@@ -127,7 +94,6 @@ void PartBins::resetPartInBin_cyc(size_t newPartNum[], int maxbinIndex) {
 }
 
 
-
 PartBins::~PartBins() {
     tmppart_m.clear();
     isEmitted_m.clear();
@@ -171,10 +137,6 @@ void PartBins::sortArray() {
 }
 
 
-void PartBins::sortArrayT() {
-    setActualemittedBin(0);
-}
-
 void PartBins::calcHBins() {
 
     for(unsigned int n = 0; n < tmppart_m.size(); n++)
@@ -182,27 +144,6 @@ void PartBins::calcHBins() {
     calcExtrema();
 }
 
-size_t PartBins::getSum() {
-    size_t s = 0;
-    for(int n = 0; n < bins_m; n++)
-        s += nBin_m[n];
-    return s;
-}
-
-void PartBins::calcGlobalExtrema() {
-    xmin_m = std::numeric_limits<double>::max();
-    xmax_m = -xmin_m;
-    for(unsigned int n = 0; n < tmppart_m.size(); n++) {
-        if(tmppart_m[n][2] <= xmin_m)
-            xmin_m = tmppart_m[n][2];
-        if(tmppart_m[n][2] >= xmax_m)
-            xmax_m = tmppart_m[n][2];
-    }
-    double xdiff = 0.01 * (xmax_m - xmin_m);
-    xmin_m -= xdiff;
-    xmax_m += xdiff;
-}
-
 void PartBins::calcExtrema() {
     for(unsigned int n = 0; n < tmppart_m.size(); n++) {
         if(xbinmin_m[(int)tmppart_m[n][6]] >= tmppart_m[n][2])
diff --git a/src/Classic/Algorithms/PartBins.h b/src/Classic/Algorithms/PartBins.h
index 67637a9b516f3c6516594130065671b51506696e..7d77d894193ccb4027c74f25c362952db9e8e3ba 100644
--- a/src/Classic/Algorithms/PartBins.h
+++ b/src/Classic/Algorithms/PartBins.h
@@ -1,33 +1,33 @@
-/** \file
- *  \brief     Defines a structure to hold energy bins and their
- *             associated data
- *
- *
- *
- *  \author    Andreas Adelmann
- *  \date      xx. November 2007
- *
- *  \warning   None.
- *  \attention
- *  \bug       sure no bug in this code :=)
- *  \todo
- */
+//
+// Class PartBins
+//   Defines a structure to hold energy bins and their associated data
+//
+// Copyright (c) 2007-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
+// All rights reserved
+//
+// This file is part of OPAL.
+//
+// OPAL 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
+//
 
 #ifndef OPAL_Bins_HH
 #define OPAL_Bins_HH
 
-#ifndef PartBinTest
-#include "Algorithms/PBunchDefs.h"
-#else
-#include "ranlib.h"
-#define Inform ostream
-#endif
-
-#include <gsl/gsl_rng.h>
 #include <gsl/gsl_histogram.h>
-#include <gsl/gsl_cdf.h>
-#include <gsl/gsl_randist.h>
 
+#include <functional>
+#include <memory>
+#include <vector>
+
+#include "Algorithms/PBunchDefs.h"
+
+class Inform;
 
 class PartBins {
 
@@ -48,14 +48,6 @@ public:
     virtual ~PartBins();
 
 
-    /** \brief How many deleted particles are on one bin */
-    inline int getDelBinCont(int bin) {
-        int a = nDelBin_m[bin];
-        reduce(a, a, OpAddAssign());
-        return a;
-    }
-
-
     /** \brief Add a particle to the temporary container */
     void fill(std::vector<double> &p) {
         tmppart_m.push_back(p);
@@ -73,11 +65,6 @@ public:
     /** set particles number in given bin */
     void setPartNum(int bin, long long num) {nBin_m[bin] = num;}
 
-    /** assume we emmit in monotinic increasing order */
-    void setBinEmitted(int bin) {binsEmitted_m[bin] = true;}
-
-    bool getBinEmitted(int bin) {return binsEmitted_m[bin];}
-
     /** \brief Is true if we still have particles to emit */
     bool doEmission() {return getNp() > 0;}
 
@@ -85,19 +72,6 @@ public:
         return isEmitted_m[n]; //(isEmitted_m[n][0]==1) && (isEmitted_m[n][1] == bin);
     }
 
-    void setEmitted(int n, int /*bin*/) {
-        isEmitted_m[n] = true;
-    }
-
-    void updatePartPos(int n, int /*bin*/, double z) {
-        tmppart_m[n][2] = z;
-    }
-
-    void updateExtramePos(int bin, double dz) {
-        xbinmax_m[bin] += dz;
-        xbinmin_m[bin] += dz;
-    }
-
     /** assigns the proper position of particle n if it belongs to bin 'bin' */
     bool getPart(size_t n, int bin, std::vector<double> &p);
 
@@ -106,29 +80,19 @@ public:
     */
     void sortArray();
 
+ private:
     /** assigns the particles to the bins */
     void calcHBins();
-    size_t getSum();
-    void calcGlobalExtrema();
     void calcExtrema();
-    void getExtrema(double &min, double &max) {
-        min = xmin_m;
-        max = xmax_m;
-    }
+    /** assume we emit in monotonic increasing order */
+    void setBinEmitted(int bin) {binsEmitted_m[bin] = true;}
 
-    /** update global bin parameters after inject a new bunch */
-    void updateStatus(int bunchCount, size_t nPartInBin);
-    /** update particles number in bin after reset Bin ID of PartBunch  */
-    void resetPartInBin(size_t newPartNum[]);
+ public:
     /** update local particles number in bin after reset Bin ID of PartBunch  */
     void resetPartInBin_cyc(size_t newPartNum[], int binID);
-    /** update particles number in bin after particle deletion */
-    void updatePartInBin(size_t countLost[]);
     /** update local particles number in bin after particle deletion */
     void updatePartInBin_cyc(size_t countLost[]);
 
-    void updateDeletedPartsInBin(size_t countLost[]) ;
-
     void setGamma(double gamma) { gamma_m = gamma;}
     double getGamma() {return gamma_m;}
 
@@ -159,13 +123,8 @@ protected:
     std::vector< std::vector<double> > tmppart_m;
     std::vector< bool > isEmitted_m;
     /** holds information whether all particles of a bin are emitted */
-    //  std::vector< bool > binsEmitted_m;
     std::unique_ptr<bool[]> binsEmitted_m;
 
-    /**
-        Here comes the new stuff, t-binning
-    */
-
 public:
 
     Inform &print(Inform &os);
@@ -175,35 +134,9 @@ public:
     /** get the number of used bin */
     virtual int getNBins() {return gsl_histogram_bins(h_m.get()) / sBins_m; }
 
-    /** Get the total number of sampled bins */
-    virtual int getNSBins() {return gsl_histogram_bins(h_m.get()); }
-
-    int getBinToEmit() {
-        int save;
-        if(nemittedBins_m < getNBins()) {
-            save =  nemittedBins_m;
-            nemittedBins_m++;
-            return save;
-        } else
-            return -1;
-    }
-
-    int getSBinToEmit() {
-        int save;
-        if(nemittedSBins_m < getNSBins()) {
-            save = nemittedSBins_m;
-            nemittedSBins_m++;
-            return save;
-        } else
-            return -1;
-    }
-
     /** the last emitted bin is always smaller or equal getNbins */
     int getLastemittedBin() {return nemittedBins_m; }
 
-    /** set the actual emitted bib */
-    void setActualemittedBin(int bin) {nemittedBins_m = bin; }
-
     /** \brief If the bunch object rebins we need to call resetBins() */
     void resetBins() {
         h_m.reset(nullptr);
@@ -213,31 +146,6 @@ public:
         return h_m != nullptr;
     }
 
-    /** sort the vector of particles according to the bin number */
-    void sortArrayT();
-
-    inline void setHistogram(gsl_histogram *h) { h_m.reset(h);}
-
-    /** \brief How many particles are on one bin */
-    virtual inline size_t getGlobalBinCount(int bin) {
-        size_t a = gsl_histogram_get(h_m.get(), bin);
-        reduce(a, a, OpAddAssign());
-        return a;
-    }
-
-    /** \brief How many particles are on one energy bin */
-    virtual inline size_t getLocalBinCount(int bin) {
-        size_t ret = 0;
-        for(int i = sBins_m * bin; i < sBins_m * (bin + 1); i++) {
-            ret += gsl_histogram_get(h_m.get(), i);
-        }
-        return ret;
-    }
-
-    /** \brief How many particles are in one sampling bin */
-    inline size_t getLocalSBinCount(int bin) { return gsl_histogram_get(h_m.get(), bin);}
-
-
     /** \brief How many particles are in all the bins */
     size_t getTotalNum();
 
@@ -250,9 +158,6 @@ protected:
     /** number of emitted bins */
     int nemittedBins_m;
 
-    /** Number of total sampled bins emitted */
-    int nemittedSBins_m;
-
     /** number of particles in the bins, the sum of all the nodes */
     std::unique_ptr<size_t[]> nBin_m;
 
diff --git a/src/Classic/Algorithms/PartBinsCyc.cpp b/src/Classic/Algorithms/PartBinsCyc.cpp
index 9b047c2c4c68d45f3d881006f728b24d33f5b5ad..dfc15aab271e1fda7e0d061283fa72c5ddbc057c 100644
--- a/src/Classic/Algorithms/PartBinsCyc.cpp
+++ b/src/Classic/Algorithms/PartBinsCyc.cpp
@@ -26,6 +26,7 @@
 //
 #include "Algorithms/PartBinsCyc.h"
 #include "Physics/Physics.h"
+#include "Utility/Inform.h"
 #include <cfloat>
 #include <vector>
 extern Inform *gmsg;
@@ -36,7 +37,7 @@ PartBinsCyc::PartBinsCyc(int specifiedNumBins, int bins, size_t  partInBin[])
 
     bins_m = specifiedNumBins;        // max bin number
     nemittedBins_m = bins;            // the bin number with particles
-    
+
     for(int i = 0; i < nemittedBins_m; i++) {
         nBin_m[i] = partInBin[i];
 
@@ -51,8 +52,8 @@ PartBinsCyc::PartBinsCyc(int specifiedNumBins, int bins)
 
     bins_m = specifiedNumBins;        // max bin number
     nemittedBins_m = bins;            // the bin number with particles
-    
+
     for(int i = 0; i < nemittedBins_m; i++) {
       binsEmitted_m[i] = true;
     }
-}
+}
\ No newline at end of file
diff --git a/src/Classic/Algorithms/PartBinsCyc.h b/src/Classic/Algorithms/PartBinsCyc.h
index 2dd9233139998633658579028efeb77c0b39bdb1..021dcef123a0c6967d1e92675009b19f2a7489ec 100644
--- a/src/Classic/Algorithms/PartBinsCyc.h
+++ b/src/Classic/Algorithms/PartBinsCyc.h
@@ -4,7 +4,7 @@
 //   associated data for multi-bunch tracking in cyclotrons
 //
 // Copyright (c) 2010, Jianjun Yang, Paul Scherrer Institut, Villigen PSI, Switzerland
-// Copyright (c) 2017-2019, Paul Scherrer Institut, Villigen PSI, Switzerland
+// Copyright (c) 2017-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
 // All rights reserved
 //
 // Implemented as part of the PhD thesis
@@ -27,47 +27,22 @@
 #ifndef OPAL_BinsCyc_HH
 #define OPAL_BinsCyc_HH
 
-#ifndef PartBinTest
-#else
-#include "ranlib.h"
-#define Inform ostream
-#endif
-
 #include "Algorithms/PartBins.h"
-#include <gsl/gsl_rng.h>
-#include <gsl/gsl_histogram.h>
-#include <gsl/gsl_cdf.h>
-#include <gsl/gsl_randist.h>
-
-#include <cassert>
 
 class PartBinsCyc: public PartBins {
 
 public:
 
 
-    /** constructer function for cyclotron*/
+    /** constructor function for cyclotron*/
     PartBinsCyc(int bunches, int bins, size_t  partInBin[]);
     PartBinsCyc(int specifiedNumBins, int bins);
 
     /** get the number of used bin */
-    int getNBins() {return bins_m; }
-
-    /** \brief How many particles are on one bin */
-    inline size_t getGlobalBinCount(int bin) {
-      size_t a = nBin_m[bin];
-      reduce(a, a, OpAddAssign());
-      return a;
-    }
-
-    /** \brief How many particles are on one energy bin */
-    inline size_t getLocalBinCount(int bin) {
-        assert(bin < bins_m);
-        return nBin_m[bin];
-    }
+    virtual int getNBins() {return bins_m; }
 
-    bool weHaveBins() {
-      return  ( nemittedBins_m > 0 );
+    virtual bool weHaveBins() {
+        return ( nemittedBins_m > 0 );
     }
 };
 
diff --git a/src/Classic/Fields/Fieldmap.cpp b/src/Classic/Fields/Fieldmap.cpp
index abaac55f96c8350746584159b130b7c5e8ed8af2..b01f97905c9abe24eadff81d4bcd35bd9de9d93b 100644
--- a/src/Classic/Fields/Fieldmap.cpp
+++ b/src/Classic/Fields/Fieldmap.cpp
@@ -1,4 +1,7 @@
 #include "Fields/Fieldmap.h"
+
+#include "Utility/PAssert.h"
+
 #include "Fields/FM3DDynamic.h"
 #include "Fields/FM3DH5BlockBase.h"
 #include "Fields/FM3DH5Block.h"
@@ -38,7 +41,6 @@
 #include <iostream>
 #include <fstream>
 #include <ios>
-#include <cassert>
 
 namespace fs = boost::filesystem;
 
@@ -366,22 +368,22 @@ MapType Fieldmap::readHeader(std::string Filename) {
         h5_prop_t props = H5CreateFileProp ();
         MPI_Comm comm = Ippl::getComm();
         h5err = H5SetPropFileMPIOCollective (props, &comm);
-        assert (h5err != H5_ERR);
+        PAssert (h5err != H5_ERR);
         h5_file_t file = H5OpenFile (Filename.c_str(), H5_O_RDONLY, props);
-        assert (file != (h5_file_t)H5_ERR);
+        PAssert (file != (h5_file_t)H5_ERR);
         H5CloseProp (props);
 
         h5err = H5SetStep(file, 0);
-        assert (h5err != H5_ERR);
+        PAssert (h5err != H5_ERR);
 
         h5_int64_t num_fields = H5BlockGetNumFields(file);
-        assert (num_fields != H5_ERR);
+        PAssert (num_fields != H5_ERR);
         MapType maptype = UNKNOWN;
 
         for (h5_ssize_t i = 0; i < num_fields; ++ i) {
             h5err = H5BlockGetFieldInfo(
                 file, (h5_size_t)i, name, len_name, NULL, NULL, NULL, NULL);
-            assert (h5err != H5_ERR);
+            PAssert (h5err != H5_ERR);
             // using field name "Bfield" and "Hfield" to distinguish the type
             if (std::strcmp(name, "Bfield") == 0) {
                 maptype = T3DMagnetoStaticH5Block;
@@ -392,7 +394,7 @@ MapType Fieldmap::readHeader(std::string Filename) {
             }
         }
         h5err = H5CloseFile(file);
-        assert (h5err != H5_ERR);
+        PAssert (h5err != H5_ERR);
         if (maptype != UNKNOWN)
             return maptype;
     }
@@ -747,7 +749,7 @@ void Fieldmap::write3DField(unsigned int nx,
     });
     std::ofstream of;
     of.open (fname);
-    assert (of.is_open ());
+    PAssert (of.is_open ());
     of.precision (6);
 
     const double hx = (xrange.second - xrange.first) / (nx - 1);
diff --git a/src/Classic/Solvers/CollimatorPhysics.cpp b/src/Classic/Solvers/CollimatorPhysics.cpp
index b2ea0faa5aeb4fba8f6d1906134e6984849177ca..11e529447c39268fb911327c03d783e0a73576d2 100644
--- a/src/Classic/Solvers/CollimatorPhysics.cpp
+++ b/src/Classic/Solvers/CollimatorPhysics.cpp
@@ -37,6 +37,8 @@
 
 #include "Ippl.h"
 
+#include <gsl/gsl_randist.h>
+
 #include <cmath>
 #include <iostream>
 #include <fstream>
diff --git a/src/Classic/Solvers/GreenWakeFunction.cpp b/src/Classic/Solvers/GreenWakeFunction.cpp
index e429f4a5e402481e5f4ec028eb91b2d53aa91b66..2f08a8d096eacc769c28610ca01d0a700e91ad5d 100644
--- a/src/Classic/Solvers/GreenWakeFunction.cpp
+++ b/src/Classic/Solvers/GreenWakeFunction.cpp
@@ -136,9 +136,9 @@ void GreenWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
             spacing = rmax(0) * rmax(0) + rmax(1) * rmax(1);
             break;
         default:
-            throw GeneralClassicException("GreenWakeFunction", "invalid direction specified");
+            throw GeneralClassicException("GreenWakeFunction::apply", "invalid direction specified");
     }
-    assert(NBin_m > 0);
+    PAssert(NBin_m > 0);
     spacing /= (NBin_m - 1); //FIXME: why -1? CKR: because grid spacings = grid points - 1
 
     // Calculate the Wakefield if needed
@@ -187,7 +187,7 @@ void GreenWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
                 int idx = (int)(floor((bunch->R[i](2) - mindist) / spacing));
                 //IFF: should be ok
                 if(idx == NBin_m) idx--;
-                assert(idx >= 0 && idx < NBin_m);
+                PAssert(idx >= 0 && idx < NBin_m);
                 double dE = OutEnergy[idx];
                 bunch->Ef[i](2) += dE;
 
@@ -201,12 +201,12 @@ void GreenWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
                 int idx = (int)(floor((bunch->R[i](2) - mindist) / spacing));
                 //IFF: should be ok
                 if(idx == NBin_m) idx--;
-                assert(idx >= 0 && idx < NBin_m);
+                PAssert(idx >= 0 && idx < NBin_m);
                 double dE = OutEnergy[idx];
 
                 // ACHTUNG spacing auch in transversal richtung
                 double dist = sqrt(bunch->R[i](0) * bunch->R[i](0) + bunch->R[i](1) * bunch->R[i](1));
-                assert(dist > 0);
+                PAssert(dist > 0);
                 bunch->Ef[i](0) += dE * bunch->R[i](0) / dist;
                 bunch->Ef[i](1) += dE * bunch->R[i](1) / dist;
 
@@ -214,7 +214,7 @@ void GreenWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
             break;
 
         default:
-            throw GeneralClassicException("GreenWakeFunction", "invalid direction specified");
+            throw GeneralClassicException("GreenWakeFunction::apply", "invalid direction specified");
     }
 
 #ifdef ENABLE_WAKE_DUMP
diff --git a/src/Classic/Solvers/GreenWakeFunction.hh b/src/Classic/Solvers/GreenWakeFunction.hh
index 52f6e0748bfdcff22e16c72e4fa673b892a9ce9c..a5da6e742ff89829e059eeeec1a3e0ce1a3d2337 100644
--- a/src/Classic/Solvers/GreenWakeFunction.hh
+++ b/src/Classic/Solvers/GreenWakeFunction.hh
@@ -21,9 +21,9 @@
 #include "Solvers/WakeFunction.hh"
 #include "Physics/Physics.h"
 #include "Utility/IpplInfo.h"
+#include "Utility/PAssert.h"
 
 #include <vector>
-#include <cassert>
 #include <map>
 #include <string>
 #include <complex>
@@ -138,8 +138,8 @@ private:
      *
      */
     template<class F> double simpson(F &f, double a, double b, unsigned int N) {
-        assert(b > a);
-        assert(N > 0);
+        PAssert(b > a);
+        PAssert(N > 0);
 
         double result = 0;
         double h = (b - a) / N;
diff --git a/src/Classic/Structure/LossDataSink.cpp b/src/Classic/Structure/LossDataSink.cpp
index 866a39b3f62bb7a7a416c1bf6c8a21f3ada5f123..fe465ead4b393cd8849f9bf2c0aa41c4afcfb874 100644
--- a/src/Classic/Structure/LossDataSink.cpp
+++ b/src/Classic/Structure/LossDataSink.cpp
@@ -9,7 +9,9 @@
 #include <boost/filesystem.hpp>
 #include "OPALconfig.h"
 
-#include <cassert>
+#include "Message/GlobalComm.h"
+#include "Utility/IpplInfo.h"
+
 #include <cmath>
 
 #define ADD_ATTACHMENT( fname ) {             \
@@ -479,7 +481,7 @@ void LossDataSink::saveASCII() {
         }
         bool res = Ippl::Comm->send(smsg, 0, tag);
         if(! res)
-            ERRORMSG("LossDataSink Ippl::Comm->send(smsg, 0, tag) failed " << endl;);
+            ERRORMSG("LossDataSink Ippl::Comm->send(smsg, 0, tag) failed " << endl);
     }
 }
 
diff --git a/src/Expressions/Expressions.cpp b/src/Expressions/Expressions.cpp
index 47a8044129bb79312dbfa510a4bd618bdac26de9..8f0170bb4fd8b6ef7e4c44ea20b1d0e5bb51ce8d 100644
--- a/src/Expressions/Expressions.cpp
+++ b/src/Expressions/Expressions.cpp
@@ -53,7 +53,6 @@
 #include "Utilities/ParseError.h"
 #include "ValueDefinitions/BoolConstant.h"
 #include <algorithm>
-#include <cassert>
 #include <cerrno>
 #include <cmath>
 #include <list>
diff --git a/src/OpalParser/OpalParser.cpp b/src/OpalParser/OpalParser.cpp
index 4e5e24f9ff31e7a6ec3ab409ba62400143b96925..8cd8ecbc022814922c8a689a9b78d639cafda26e 100644
--- a/src/OpalParser/OpalParser.cpp
+++ b/src/OpalParser/OpalParser.cpp
@@ -34,7 +34,6 @@
 #include "Utilities/OpalException.h"
 #include "Utilities/ParseError.h"
 #include "Utilities/Options.h"
-#include <cassert>
 #include <cmath>
 #include <ctime>
 #include <exception>
@@ -178,7 +177,7 @@ void OpalParser::parseAction(Statement &stat) const {
     } else {
         std::string hint = getHint(cmdName, "command");
         if (hint != "") {
-            throw ParseError("OpalParser::parse()",
+            throw ParseError("OpalParser::parseAction()",
                              "Syntax error, " + hint);
         }
 
diff --git a/src/Solvers/ArbitraryDomain.cpp b/src/Solvers/ArbitraryDomain.cpp
index d878d959405d8652f9b928a6452767065e9d0725..26fc09696474c51b84a4fe2bc8049afa51c62a78 100644
--- a/src/Solvers/ArbitraryDomain.cpp
+++ b/src/Solvers/ArbitraryDomain.cpp
@@ -34,7 +34,6 @@
 #include <cmath>
 #include <iostream>
 #include <tuple>
-#include <cassert>
 #include "Utilities/OpalException.h"
 
 ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
diff --git a/src/Solvers/BoxCornerDomain.cpp b/src/Solvers/BoxCornerDomain.cpp
index 0980a52a91eb32dbf0c73056baa9ce04f3cc969c..8f0bd5d32453f85d0398c52e28f2e05b742f0be8 100644
--- a/src/Solvers/BoxCornerDomain.cpp
+++ b/src/Solvers/BoxCornerDomain.cpp
@@ -30,7 +30,6 @@
 #include <string>
 #include <cmath>
 #include <iostream>
-#include <cassert>
 
 //FIXME: ORDER HOW TO TRAVERSE NODES IS FIXED, THIS SHOULD BE MORE GENERIC! (PLACES MARKED)
 
diff --git a/src/Solvers/BoxLibSolvers/FMGPoissonSolver.cpp b/src/Solvers/BoxLibSolvers/FMGPoissonSolver.cpp
index 8e7cda8ac756990ab09e0e34a1a577fdd6eb7d0d..8db36e7236309b65aa0bbe32da3a39c23e0ba498 100644
--- a/src/Solvers/BoxLibSolvers/FMGPoissonSolver.cpp
+++ b/src/Solvers/BoxLibSolvers/FMGPoissonSolver.cpp
@@ -22,13 +22,13 @@
 //
 #include "FMGPoissonSolver.h"
 
+#include "Utility/PAssert.h"
+
 #include "Utilities/OpalException.h"
 
 #include <AMReX_ParmParse.H>
 #include <AMReX_Interpolater.H>
 
-#include <cassert>
-
 FMGPoissonSolver::FMGPoissonSolver(AmrBoxLib* itsAmrObject_p)
     : AmrPoissonSolver<AmrBoxLib>(itsAmrObject_p),
       reltol_m(1.0e-15),
@@ -70,14 +70,14 @@ void FMGPoissonSolver::solve(AmrScalarFieldContainer_t& rho,
 
     amrex::Vector< AmrField_t > tmp(rho.size());
 
-    assert(baseLevel <= finestLevel);
-    assert(finestLevel < (int)rho.size());
+    PAssert(baseLevel <= finestLevel);
+    PAssert(finestLevel < (int)rho.size());
 
     for (int lev = baseLevel; lev <= finestLevel ; ++lev) {
         const AmrProcMap_t& dmap = rho[lev]->DistributionMap();
         grad_phi_edge[lev].resize(AMREX_SPACEDIM);
         tmp[lev].define(rho[lev]->boxArray(), dmap, AMREX_SPACEDIM, 1);
-        
+
         for (int n = 0; n < AMREX_SPACEDIM; ++n) {
             AmrGrid_t ba = rho[lev]->boxArray();
             grad_phi_edge[lev][n].reset(new AmrField_t(ba.surroundingNodes(n), dmap, 1, 1));
diff --git a/src/Solvers/EllipticDomain.cpp b/src/Solvers/EllipticDomain.cpp
index 0b7be840f88885c8b7000bbfd15eeda3858c5c7d..ac9d2a03798c4461ec3b0401e96a5f9ee9c7ecf6 100644
--- a/src/Solvers/EllipticDomain.cpp
+++ b/src/Solvers/EllipticDomain.cpp
@@ -30,7 +30,6 @@
 #include <map>
 #include <cmath>
 #include <iostream>
-#include <cassert>
 
 //FIXME: ORDER HOW TO TRAVERSE NODES IS FIXED, THIS SHOULD BE MORE GENERIC! (PLACES MARKED)
 
diff --git a/src/Solvers/IrregularDomain.cpp b/src/Solvers/IrregularDomain.cpp
index 0597bcc8406e38c509b23f5861168afa54b55fcc..f9e1cc6a2e26704b3ba718099b22c8ee99cf177e 100644
--- a/src/Solvers/IrregularDomain.cpp
+++ b/src/Solvers/IrregularDomain.cpp
@@ -25,9 +25,9 @@
 //
 
 #include "Solvers/IrregularDomain.h"
-#include "Utilities/OpalException.h"
 
-#include <cassert>
+#include "Utility/PAssert.h"
+#include "Utilities/OpalException.h"
 
 IrregularDomain::IrregularDomain(const IntVector_t& nr, const Vector_t& hr,
                                  const std::string& interpl)
@@ -114,7 +114,7 @@ void IrregularDomain::getBoundaryStencil(int x, int y, int z, StencilValue_t& va
     }
 
     // stencil center value has to be positive!
-    assert(value.center > 0);
+    PAssert(value.center > 0);
 }
 
 
diff --git a/src/Structure/BoundaryGeometry.cpp b/src/Structure/BoundaryGeometry.cpp
index 8c59c6e2f421c9324cd40fade1e929e72de321ef..4fd93edb990830bfcc52128d661bdb17dbe35e15 100644
--- a/src/Structure/BoundaryGeometry.cpp
+++ b/src/Structure/BoundaryGeometry.cpp
@@ -13,15 +13,16 @@
 // (at your option) any later version.
 //
 // You should have received a copy of the GNU General Public License
-// along with OPAL.  If not, see <https://www.gnu.org/licenses/>.
+// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
 //
 
 //#define   ENABLE_DEBUG
 
 #include "Structure/BoundaryGeometry.h"
 
-#include <fstream>
 #include <cmath>
+#include <fstream>
+#include <string>
 
 #include "H5hut.h"
 
@@ -685,7 +686,7 @@ static inline Vector_t normalVector (
     ) {
     const Vector_t N = cross (B - A, C - A);
     const double magnitude = std::sqrt (SQR (N (0)) + SQR (N (1)) + SQR (N (2)));
-    assert (gsl_fcmp (magnitude, 0.0, EPS) > 0); // in case we have degenerted triangles
+    PAssert (gsl_fcmp (magnitude, 0.0, EPS) > 0); // in case we have degenerated triangles
     return N / magnitude;
 }
 
@@ -746,17 +747,17 @@ BoundaryGeometry::BoundaryGeometry() :
 
     itsAttr[L1] = Attributes::makeReal
         ("L1",
-         "In case of BOXCORNER Specifies first part with hight == B [m]",
+         "In case of BOXCORNER Specifies first part with height == B [m]",
          0.5);
 
     itsAttr[L2] = Attributes::makeReal
         ("L2",
-         "In case of BOXCORNER Specifies first second with hight == B-C [m]",
+         "In case of BOXCORNER Specifies first second with height == B-C [m]",
          0.2);
 
     itsAttr[C] = Attributes::makeReal
         ("C",
-         "In case of BOXCORNER Specifies hight of corner C [m]",
+         "In case of BOXCORNER Specifies height of corner C [m]",
          0.01);
 
     itsAttr[XYZSCALE] = Attributes::makeReal
@@ -1589,7 +1590,7 @@ Change orientation if diff is:
             for (unsigned int triangle_id = 0; triangle_id < bg->Triangles_m.size(); triangle_id++) {
                 for (unsigned int j = 1; j <= 3; j++) {
                     auto pt_id = bg->PointID (triangle_id, j);
-                    assert (pt_id < bg->Points_m.size ());
+                    PAssert (pt_id < bg->Points_m.size ());
                     adjacencies_to_pt [pt_id].insert (triangle_id);
                 }
             }
@@ -1708,7 +1709,7 @@ Change orientation if diff is:
                     }
                 }
             }
-            assert (n == 2);
+            PAssert (n == 2);
         edge_found:
             int diff = id[1] - id[0];
             if ((((ic[1] - ic[0]) == 1) && ((diff == 1) || (diff == -2))) ||
@@ -1798,7 +1799,7 @@ Change orientation if diff is:
     (void)rc;
 #endif
     rc = H5SetErrorHandler (H5AbortErrorhandler);
-    assert (rc != H5_ERR);
+    PAssert (rc != H5_ERR);
     H5SetVerbosityLevel (1);
 
     h5_prop_t props = H5CreateFileProp ();
@@ -2052,7 +2053,7 @@ BoundaryGeometry::intersectTinyLineSegmentBoundary (
             }
             break;
         case -1:                    // triangle is degenerated
-            assert (tmp_intersect_result != -1);
+            PAssert (tmp_intersect_result != -1);
             exit (42);              // terminate even if NDEBUG is set
         }
     }                   // end for all triangles
@@ -2192,7 +2193,7 @@ void
 BoundaryGeometry::writeGeomToVtk (std::string fn) {
     std::ofstream of;
     of.open (fn.c_str ());
-    assert (of.is_open ());
+    PAssert (of.is_open ());
     of.precision (6);
     of << "# vtk DataFile Version 2.0" << std::endl;
     of << "generated using DataSink::writeGeoToVtk" << std::endl;
diff --git a/src/Structure/BoundaryGeometry.h b/src/Structure/BoundaryGeometry.h
index 549b2e45ea0754c3632409bcff13b73dbc0878bb..2ec308c6307b3e1ffb257f9a9116c6035fb54e0e 100644
--- a/src/Structure/BoundaryGeometry.h
+++ b/src/Structure/BoundaryGeometry.h
@@ -36,20 +36,19 @@
 class OpalBeamline;
 class ElementBase;
 
-#include <cassert>
 #include <unordered_map>
 #include <unordered_set>
 #include <array>
+#include <vector>
 
 #include "AbstractObjects/Definition.h"
 #include "Attributes/Attributes.h"
 #include "Utilities/Util.h"
 #include "Utility/IpplTimings.h"
+#include "Utility/PAssert.h"
 
 #include <gsl/gsl_rng.h>
 
-extern Inform* gmsg;
-
 class BoundaryGeometry : public Definition {
 
 public:
@@ -286,7 +285,7 @@ private:
     }
 
     inline const Vector_t& getPoint (const int triangle_id, const int vertex_id) {
-        assert (1 <= vertex_id && vertex_id <=3);
+        PAssert (1 <= vertex_id && vertex_id <=3);
         return Points_m[Triangles_m[triangle_id][vertex_id]];
     }
 
@@ -312,12 +311,12 @@ private:
         FGEOM,    // file holding the geometry
         LENGTH,   // length of elliptic tube or boxcorner
         S,        // start of the geometry
-        L1,       // in case of BOXCORNER first part of geometry with hight B
-        L2,       // in case of BOXCORNER second part of geometry with hight B-C
+        L1,       // in case of BOXCORNER first part of geometry with height B
+        L2,       // in case of BOXCORNER second part of geometry with height B-C
         A,        // major semi-axis of elliptic tube
         B,        // minor semi-axis of ellitpic tube
-        C,        // in case of BOXCORNER hight of corner
-        TOPO,     // RECTANGULAR, BOXCORNER, ELLIPTIC if FGEOM is selected topo is over-written
+        C,        // in case of BOXCORNER height of corner
+        TOPO,     // RECTANGULAR, BOXCORNER, ELLIPTIC if FGEOM is selected TOPO is over-written
         ZSHIFT,   // Shift in z direction
         XYZSCALE, // Multiplicative scaling factor for coordinates
         XSCALE,   // Multiplicative scaling factor for x-coordinates
@@ -339,4 +338,4 @@ inline Inform &operator<< (Inform& os, const BoundaryGeometry& b) {
 // c-basic-offset: 4
 // indent-tabs-mode: nil
 // require-final-newline: nil
-// End:
+// End:
\ No newline at end of file
diff --git a/src/Structure/H5PartWrapper.cpp b/src/Structure/H5PartWrapper.cpp
index 109ef19225d6963ed66a089614a8b338f0c06435..a8184b2f0053b8531869a60d579cdac752e5d21a 100644
--- a/src/Structure/H5PartWrapper.cpp
+++ b/src/Structure/H5PartWrapper.cpp
@@ -6,8 +6,13 @@
 
 #include "OPALconfig.h"
 #include "AbstractObjects/OpalData.h"
-#include "Utilities/Options.h"
 #include "Physics/Physics.h"
+#include "Utilities/Options.h"
+#include "Utilities/OpalException.h"
+
+#include "Message/Communicate.h"
+#include "Message/Message.h"
+#include "Utility/PAssert.h"
 
 #include <boost/filesystem.hpp>
 
@@ -40,8 +45,6 @@ H5PartWrapper::H5PartWrapper(const std::string &fileName, int restartStep, std::
     numSteps_m(0),
     startedFromExistingFile_m(true)
 {
-    namespace fs = boost::filesystem;
-
     if (sourceFile == "") sourceFile = fileName_m;
 
     copyFile(sourceFile, restartStep, flags);
@@ -72,9 +75,9 @@ void H5PartWrapper::open(h5_int32_t flags) {
 #if defined (NDEBUG)
     (void)h5err;
 #endif
-    assert (h5err != H5_ERR);
+    PAssert (h5err != H5_ERR);
     file_m = H5OpenFile (fileName_m.c_str(), flags, props);
-    assert (file_m != (h5_file_t)H5_ERR);
+    PAssert (file_m != (h5_file_t)H5_ERR);
     H5CloseProp (props);
 }
 
@@ -135,9 +138,9 @@ void H5PartWrapper::copyFile(const std::string &sourceFile, int lastStep, h5_int
 #if defined (NDEBUG)
         (void)h5err;
 #endif
-        assert (h5err != H5_ERR);
+        PAssert (h5err != H5_ERR);
         h5_file_t source = H5OpenFile (sourceFile.c_str(), H5_O_RDONLY, props);
-        assert (source != (h5_file_t)H5_ERR);
+        PAssert (source != (h5_file_t)H5_ERR);
         H5CloseProp (props);
         h5_ssize_t numStepsInSource = H5GetNumSteps(source);
 
@@ -169,9 +172,9 @@ void H5PartWrapper::copyFile(const std::string &sourceFile, int lastStep, h5_int
         props = H5CreateFileProp ();
         comm = Ippl::getComm();
         h5err = H5SetPropFileMPIOCollective (props, &comm);
-        assert (h5err != H5_ERR);
+        PAssert (h5err != H5_ERR);
         source = H5OpenFile (sourceFileName.c_str(), H5_O_RDONLY, props);
-        assert (source != (h5_file_t)H5_ERR);
+        PAssert (source != (h5_file_t)H5_ERR);
         H5CloseProp (props);
         copyHeader(source);
 
@@ -204,9 +207,9 @@ void H5PartWrapper::copyFile(const std::string &sourceFile, int lastStep, h5_int
 #if defined (NDEBUG)
         (void)h5err;
 #endif
-        assert (h5err != H5_ERR);
+        PAssert (h5err != H5_ERR);
         h5_file_t source = H5OpenFile (sourceFile.c_str(), H5_O_RDONLY, props);
-        assert (source != (h5_file_t)H5_ERR);
+        PAssert (source != (h5_file_t)H5_ERR);
         H5CloseProp (props);
         h5_ssize_t numStepsInSource = H5GetNumSteps(source);