diff --git a/bfps/cpp/distributed_particles.cpp b/bfps/cpp/distributed_particles.cpp
index f486568160bdbd7ba6b8b9d7d0c66bc32ae678c9..3773b5f4e0c84811f456a7cef0cac7b750847381 100644
--- a/bfps/cpp/distributed_particles.cpp
+++ b/bfps/cpp/distributed_particles.cpp
@@ -174,7 +174,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
     buffer_size = (buffer_size > npr[0])? buffer_size : npr[0];
     buffer_size = (buffer_size > npr[1])? buffer_size : npr[1];
     //DEBUG_MSG("buffer size is %d\n", buffer_size);
-    double *buffer = new double[buffer_size*this->ncomponents*(1+vals.size())];
+    double *buffer = new double[buffer_size*state_dimension(particle_type)*(1+vals.size())];
     for (rsrc = 0; rsrc<this->nprocs; rsrc++)
         for (int i=0; i<2; i++)
         {
@@ -192,21 +192,21 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
                 for (int p: ps[i])
                 {
                     std::copy(x[p].data,
-                              x[p].data + this->ncomponents,
-                              buffer + pcounter*(1+vals.size())*this->ncomponents);
+                              x[p].data + state_dimension(particle_type),
+                              buffer + pcounter*(1+vals.size())*state_dimension(particle_type));
                     x.erase(p);
                     for (int tindex=0; tindex<vals.size(); tindex++)
                     {
                         std::copy(vals[tindex][p].data,
-                                  vals[tindex][p].data + this->ncomponents,
-                                  buffer + (pcounter*(1+vals.size()) + tindex+1)*this->ncomponents);
+                                  vals[tindex][p].data + state_dimension(particle_type),
+                                  buffer + (pcounter*(1+vals.size()) + tindex+1)*state_dimension(particle_type));
                         vals[tindex].erase(p);
                     }
                     pcounter++;
                 }
                 MPI_Send(
                         buffer,
-                        nps[i]*(1+vals.size())*this->ncomponents,
+                        nps[i]*(1+vals.size())*state_dimension(particle_type),
                         MPI_DOUBLE,
                         rdst,
                         2*(rsrc*this->nprocs + rdst)+1,
@@ -224,7 +224,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
                         MPI_STATUS_IGNORE);
                 MPI_Recv(
                         buffer,
-                        npr[1-i]*(1+vals.size())*this->ncomponents,
+                        npr[1-i]*(1+vals.size())*state_dimension(particle_type),
                         MPI_DOUBLE,
                         rsrc,
                         2*(rsrc*this->nprocs + rdst)+1,
@@ -233,10 +233,10 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::redistrib
                 int pcounter = 0;
                 for (int p: pr[1-i])
                 {
-                    x[p] = buffer + (pcounter*(1+vals.size()))*this->ncomponents;
+                    x[p] = buffer + (pcounter*(1+vals.size()))*state_dimension(particle_type);
                     for (int tindex=0; tindex<vals.size(); tindex++)
                     {
-                        vals[tindex][p] = buffer + (pcounter*(1+vals.size()) + tindex+1)*this->ncomponents;
+                        vals[tindex][p] = buffer + (pcounter*(1+vals.size()) + tindex+1)*state_dimension(particle_type);
                     }
                     pcounter++;
                 }
@@ -267,7 +267,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBash
 {
     this->get_rhs(this->state, this->rhs[0]);
     for (auto &pp: this->state)
-        for (int i=0; i<this->ncomponents; i++)
+        for (int i=0; i<state_dimension(particle_type); i++)
             switch(nsteps)
             {
                 case 1:
@@ -322,7 +322,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::step()
 template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
 {
-    double *temp = new double[this->chunk_size*this->ncomponents];
+    double *temp = new double[this->chunk_size*state_dimension(particle_type)];
     for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
     {
         //read state
@@ -330,14 +330,14 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
             this->read_state_chunk(cindex, temp);
         MPI_Bcast(
                 temp,
-                this->chunk_size*this->ncomponents,
+                this->chunk_size*state_dimension(particle_type),
                 MPI_DOUBLE,
                 0,
                 this->comm);
         for (int p=0; p<this->chunk_size; p++)
         {
-            if (this->vel->get_rank(temp[this->ncomponents*p+2]) == this->myrank)
-                this->state[p+cindex*this->chunk_size] = temp + this->ncomponents*p;
+            if (this->vel->get_rank(temp[state_dimension(particle_type)*p+2]) == this->myrank)
+                this->state[p+cindex*this->chunk_size] = temp + state_dimension(particle_type)*p;
         }
         //read rhs
         if (this->iteration > 0)
@@ -347,7 +347,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
                     this->read_rhs_chunk(cindex, i, temp);
                 MPI_Bcast(
                         temp,
-                        this->chunk_size*this->ncomponents,
+                        this->chunk_size*state_dimension(particle_type),
                         MPI_DOUBLE,
                         0,
                         this->comm);
@@ -355,7 +355,7 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
                 {
                     auto pp = this->state.find(p+cindex*this->chunk_size);
                     if (pp != this->state.end())
-                        this->rhs[i][p+cindex*this->chunk_size] = temp + this->ncomponents*p;
+                        this->rhs[i][p+cindex*this->chunk_size] = temp + state_dimension(particle_type)*p;
                 }
             }
     }
@@ -399,24 +399,24 @@ template <particle_types particle_type, class rnumber, int interp_neighbours>
 void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
         const bool write_rhs)
 {
-    double *temp0 = new double[this->chunk_size*this->ncomponents];
-    double *temp1 = new double[this->chunk_size*this->ncomponents];
+    double *temp0 = new double[this->chunk_size*state_dimension(particle_type)];
+    double *temp1 = new double[this->chunk_size*state_dimension(particle_type)];
     for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
     {
         //write state
-        std::fill_n(temp0, this->ncomponents*this->chunk_size, 0);
+        std::fill_n(temp0, state_dimension(particle_type)*this->chunk_size, 0);
         for (int p=0; p<this->chunk_size; p++)
         {
             auto pp = this->state.find(p + cindex*this->chunk_size);
             if (pp != this->state.end())
                 std::copy(pp->second.data,
-                          pp->second.data + this->ncomponents,
-                          temp0 + pp->first*this->ncomponents);
+                          pp->second.data + state_dimension(particle_type),
+                          temp0 + pp->first*state_dimension(particle_type));
         }
         MPI_Allreduce(
                 temp0,
                 temp1,
-                this->ncomponents*this->chunk_size,
+                state_dimension(particle_type)*this->chunk_size,
                 MPI_DOUBLE,
                 MPI_SUM,
                 this->comm);
@@ -426,19 +426,19 @@ void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
         if (write_rhs)
             for (int i=0; i<this->integration_steps; i++)
             {
-                std::fill_n(temp0, this->ncomponents*this->chunk_size, 0);
+                std::fill_n(temp0, state_dimension(particle_type)*this->chunk_size, 0);
                 for (int p=0; p<this->chunk_size; p++)
                 {
                     auto pp = this->rhs[i].find(p + cindex*this->chunk_size);
                     if (pp != this->rhs[i].end())
                         std::copy(pp->second.data,
-                                  pp->second.data + this->ncomponents,
-                                  temp0 + pp->first*this->ncomponents);
+                                  pp->second.data + state_dimension(particle_type),
+                                  temp0 + pp->first*state_dimension(particle_type));
                 }
                 MPI_Allreduce(
                         temp0,
                         temp1,
-                        this->ncomponents*this->chunk_size,
+                        state_dimension(particle_type)*this->chunk_size,
                         MPI_DOUBLE,
                         MPI_SUM,
                         this->comm);
diff --git a/bfps/cpp/particles.cpp b/bfps/cpp/particles.cpp
index 9f55f65db12d05f12df9d7daa452c1d79f9b6fcf..459ea7a37847a6ce4f217bc937130f9304466ed1 100644
--- a/bfps/cpp/particles.cpp
+++ b/bfps/cpp/particles.cpp
@@ -55,7 +55,7 @@ particles<particle_type, rnumber, interp_neighbours>::particles(
            (INTEGRATION_STEPS >= 1));
     this->vel = FIELD;
     this->integration_steps = INTEGRATION_STEPS;
-    this->array_size = this->nparticles * this->ncomponents;
+    this->array_size = this->nparticles * state_dimension(particle_type);
     this->state = new double[this->array_size];
     std::fill_n(this->state, this->array_size, 0.0);
     for (int i=0; i < this->integration_steps; i++)
@@ -81,7 +81,7 @@ void particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, do
     switch(particle_type)
     {
         case VELOCITY_TRACER:
-            this->vel->sample(this->nparticles, this->ncomponents, x, y);
+            this->vel->sample(this->nparticles, state_dimension(particle_type), x, y);
             break;
     }
 }
@@ -107,26 +107,26 @@ void particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
     {
         case 1:
             for (int p=0; p<this->nparticles; p++)
-                for (int i=0; i<this->ncomponents; i++)
+                for (int i=0; i<state_dimension(particle_type); i++)
                 {
-                    ii = p*this->ncomponents+i;
+                    ii = p*state_dimension(particle_type)+i;
                     this->state[ii] += this->dt*this->rhs[0][ii];
                 }
             break;
         case 2:
             for (int p=0; p<this->nparticles; p++)
-                for (int i=0; i<this->ncomponents; i++)
+                for (int i=0; i<state_dimension(particle_type); i++)
                 {
-                    ii = p*this->ncomponents+i;
+                    ii = p*state_dimension(particle_type)+i;
                     this->state[ii] += this->dt*(3*this->rhs[0][ii]
                                                -   this->rhs[1][ii])/2;
                 }
             break;
         case 3:
             for (int p=0; p<this->nparticles; p++)
-                for (int i=0; i<this->ncomponents; i++)
+                for (int i=0; i<state_dimension(particle_type); i++)
                 {
-                    ii = p*this->ncomponents+i;
+                    ii = p*state_dimension(particle_type)+i;
                     this->state[ii] += this->dt*(23*this->rhs[0][ii]
                                                - 16*this->rhs[1][ii]
                                                +  5*this->rhs[2][ii])/12;
@@ -134,9 +134,9 @@ void particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
             break;
         case 4:
             for (int p=0; p<this->nparticles; p++)
-                for (int i=0; i<this->ncomponents; i++)
+                for (int i=0; i<state_dimension(particle_type); i++)
                 {
-                    ii = p*this->ncomponents+i;
+                    ii = p*state_dimension(particle_type)+i;
                     this->state[ii] += this->dt*(55*this->rhs[0][ii]
                                                - 59*this->rhs[1][ii]
                                                + 37*this->rhs[2][ii]
@@ -145,9 +145,9 @@ void particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
             break;
         case 5:
             for (int p=0; p<this->nparticles; p++)
-                for (int i=0; i<this->ncomponents; i++)
+                for (int i=0; i<state_dimension(particle_type); i++)
                 {
-                    ii = p*this->ncomponents+i;
+                    ii = p*state_dimension(particle_type)+i;
                     this->state[ii] += this->dt*(1901*this->rhs[0][ii]
                                                - 2774*this->rhs[1][ii]
                                                + 2616*this->rhs[2][ii]
@@ -157,9 +157,9 @@ void particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
             break;
         case 6:
             for (int p=0; p<this->nparticles; p++)
-                for (int i=0; i<this->ncomponents; i++)
+                for (int i=0; i<state_dimension(particle_type); i++)
                 {
-                    ii = p*this->ncomponents+i;
+                    ii = p*state_dimension(particle_type)+i;
                     this->state[ii] += this->dt*(4277*this->rhs[0][ii]
                                                - 7923*this->rhs[1][ii]
                                                + 9982*this->rhs[2][ii]
@@ -189,10 +189,10 @@ void particles<particle_type, rnumber, interp_neighbours>::read()
     if (this->myrank == 0)
         for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
         {
-            this->read_state_chunk(cindex, this->state+cindex*this->chunk_size*this->ncomponents);
+            this->read_state_chunk(cindex, this->state+cindex*this->chunk_size*state_dimension(particle_type));
             if (this->iteration > 0)
                 for (int i=0; i<this->integration_steps; i++)
-                    this->read_rhs_chunk(cindex, i, this->rhs[i]+cindex*this->chunk_size*this->ncomponents);
+                    this->read_rhs_chunk(cindex, i, this->rhs[i]+cindex*this->chunk_size*state_dimension(particle_type));
         }
     MPI_Bcast(
             this->state,
@@ -217,10 +217,10 @@ void particles<particle_type, rnumber, interp_neighbours>::write(
     if (this->myrank == 0)
         for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
         {
-            this->write_state_chunk(cindex, this->state+cindex*this->chunk_size*this->ncomponents);
+            this->write_state_chunk(cindex, this->state+cindex*this->chunk_size*state_dimension(particle_type));
             if (write_rhs)
                 for (int i=0; i<this->integration_steps; i++)
-                    this->write_rhs_chunk(cindex, i, this->rhs[i]+cindex*this->chunk_size*this->ncomponents);
+                    this->write_rhs_chunk(cindex, i, this->rhs[i]+cindex*this->chunk_size*state_dimension(particle_type));
         }
 }
 
@@ -230,7 +230,7 @@ void particles<particle_type, rnumber, interp_neighbours>::sample(
         const char *dset_name)
 {
     double *y = new double[this->nparticles*3];
-    field->sample(this->nparticles, this->ncomponents, this->state, y);
+    field->sample(this->nparticles, state_dimension(particle_type), this->state, y);
     if (this->myrank == 0)
         for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
             this->write_point3D_chunk(dset_name, cindex, y+cindex*this->chunk_size*3);
diff --git a/bfps/cpp/particles.hpp b/bfps/cpp/particles.hpp
index 06caefdd5e28881bf9a8e2dd80730645d5294583..03daf3e3fc866ac485b3649a28dfb13cf1b50ff1 100644
--- a/bfps/cpp/particles.hpp
+++ b/bfps/cpp/particles.hpp
@@ -75,7 +75,7 @@ class particles: public particles_io_base<particle_type>
                 interpolator_base<rnumber, interp_neighbours> *field,
                 double *y)
         {
-            field->sample(this->nparticles, this->ncomponents, this->state, y);
+            field->sample(this->nparticles, state_dimension(particle_type), this->state, y);
         }
 
         void get_rhs(
diff --git a/bfps/cpp/particles_base.cpp b/bfps/cpp/particles_base.cpp
index 5fd9e79c1dd535f552b88555c5836bd3b6f19b87..59d6d078dd63a319f34d3ae03acccc1f259714be 100644
--- a/bfps/cpp/particles_base.cpp
+++ b/bfps/cpp/particles_base.cpp
@@ -119,7 +119,6 @@ particles_io_base<particle_type>::particles_io_base(
         const hid_t data_file_id,
         MPI_Comm COMM)
 {
-    this->ncomponents = state_dimension(particle_type);
     this->name = std::string(NAME);
     this->traj_skip = TRAJ_SKIP;
     this->comm = COMM;
@@ -134,7 +133,7 @@ particles_io_base<particle_type>::particles_io_base(
         dspace = H5Dget_space(dset);
         this->hdf5_state_dims.resize(H5Sget_simple_extent_ndims(dspace));
         H5Sget_simple_extent_dims(dspace, &this->hdf5_state_dims.front(), NULL);
-        assert(this->hdf5_state_dims[this->hdf5_state_dims.size()-1] == this->ncomponents);
+        assert(this->hdf5_state_dims[this->hdf5_state_dims.size()-1] == state_dimension(particle_type));
         this->nparticles = 1;
         for (int i=1; i<this->hdf5_state_dims.size()-1; i++)
             this->nparticles *= this->hdf5_state_dims[i];
diff --git a/bfps/cpp/particles_base.hpp b/bfps/cpp/particles_base.hpp
index a207e6b41c0670a7bf1db2946610ef9945bf5ea3..114661fcc4c8f8e7d1a38bbb6d646e5745b9415f 100644
--- a/bfps/cpp/particles_base.hpp
+++ b/bfps/cpp/particles_base.hpp
@@ -78,7 +78,6 @@ class particles_io_base
         MPI_Comm comm;
 
         int nparticles;
-        int ncomponents;
 
         std::string name;
         int chunk_size;
diff --git a/bfps/cpp/rFFTW_distributed_particles.cpp b/bfps/cpp/rFFTW_distributed_particles.cpp
index 32aefaf3b7f40e1931d2ff5e66e7bb7d61655289..412671b89ea000efadd713e9cf67221009b0d2cb 100644
--- a/bfps/cpp/rFFTW_distributed_particles.cpp
+++ b/bfps/cpp/rFFTW_distributed_particles.cpp
@@ -336,7 +336,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::red
     buffer_size = (buffer_size > npr[0])? buffer_size : npr[0];
     buffer_size = (buffer_size > npr[1])? buffer_size : npr[1];
     //DEBUG_MSG("buffer size is %d\n", buffer_size);
-    double *buffer = new double[buffer_size*this->ncomponents*(1+vals.size())];
+    double *buffer = new double[buffer_size*state_dimension(particle_type)*(1+vals.size())];
     for (rsrc = 0; rsrc<this->nprocs; rsrc++)
         for (int i=0; i<2; i++)
         {
@@ -354,19 +354,19 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::red
                 for (int p: ps[i])
                 {
                     std::copy(x[p].data,
-                              x[p].data + this->ncomponents,
-                              buffer + pcounter*(1+vals.size())*this->ncomponents);
+                              x[p].data + state_dimension(particle_type),
+                              buffer + pcounter*(1+vals.size())*state_dimension(particle_type));
                     for (int tindex=0; tindex<vals.size(); tindex++)
                     {
                         std::copy(vals[tindex][p].data,
-                                  vals[tindex][p].data + this->ncomponents,
-                                  buffer + (pcounter*(1+vals.size()) + tindex+1)*this->ncomponents);
+                                  vals[tindex][p].data + state_dimension(particle_type),
+                                  buffer + (pcounter*(1+vals.size()) + tindex+1)*state_dimension(particle_type));
                     }
                     pcounter++;
                 }
                 MPI_Send(
                         buffer,
-                        nps[i]*(1+vals.size())*this->ncomponents,
+                        nps[i]*(1+vals.size())*state_dimension(particle_type),
                         MPI_DOUBLE,
                         rdst,
                         2*(rsrc*this->nprocs + rdst)+1,
@@ -384,7 +384,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::red
                         MPI_STATUS_IGNORE);
                 MPI_Recv(
                         buffer,
-                        npr[1-i]*(1+vals.size())*this->ncomponents,
+                        npr[1-i]*(1+vals.size())*state_dimension(particle_type),
                         MPI_DOUBLE,
                         rsrc,
                         2*(rsrc*this->nprocs + rdst)+1,
@@ -393,11 +393,11 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::red
                 int pcounter = 0;
                 for (int p: pr[1-i])
                 {
-                    x[p] = buffer + (pcounter*(1+vals.size()))*this->ncomponents;
+                    x[p] = buffer + (pcounter*(1+vals.size()))*state_dimension(particle_type);
                     newdp[1-i].insert(p);
                     for (int tindex=0; tindex<vals.size(); tindex++)
                     {
-                        vals[tindex][p] = buffer + (pcounter*(1+vals.size()) + tindex+1)*this->ncomponents;
+                        vals[tindex][p] = buffer + (pcounter*(1+vals.size()) + tindex+1)*state_dimension(particle_type);
                     }
                     pcounter++;
                 }
@@ -431,7 +431,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::Ada
 {
     this->get_rhs(this->state, this->domain_particles, this->rhs[0]);
     for (auto &pp: this->state)
-        for (int i=0; i<this->ncomponents; i++)
+        for (int i=0; i<state_dimension(particle_type); i++)
             switch(nsteps)
             {
                 case 1:
@@ -526,7 +526,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::sor
 template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::read()
 {
-    double *temp = new double[this->chunk_size*this->ncomponents];
+    double *temp = new double[this->chunk_size*state_dimension(particle_type)];
     int tmpint1, tmpint2;
     for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
     {
@@ -535,15 +535,15 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::rea
             this->read_state_chunk(cindex, temp);
         MPI_Bcast(
                 temp,
-                this->chunk_size*this->ncomponents,
+                this->chunk_size*state_dimension(particle_type),
                 MPI_DOUBLE,
                 0,
                 this->comm);
         for (int p=0; p<this->chunk_size; p++)
         {
-            if (this->vel->get_rank_info(temp[this->ncomponents*p+2], tmpint1, tmpint2))
+            if (this->vel->get_rank_info(temp[state_dimension(particle_type)*p+2], tmpint1, tmpint2))
             {
-                this->state[p+cindex*this->chunk_size] = temp + this->ncomponents*p;
+                this->state[p+cindex*this->chunk_size] = temp + state_dimension(particle_type)*p;
             }
         }
         //read rhs
@@ -554,7 +554,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::rea
                     this->read_rhs_chunk(cindex, i, temp);
                 MPI_Bcast(
                         temp,
-                        this->chunk_size*this->ncomponents,
+                        this->chunk_size*state_dimension(particle_type),
                         MPI_DOUBLE,
                         0,
                         this->comm);
@@ -562,7 +562,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::rea
                 {
                     auto pp = this->state.find(p+cindex*this->chunk_size);
                     if (pp != this->state.end())
-                        this->rhs[i][p+cindex*this->chunk_size] = temp + this->ncomponents*p;
+                        this->rhs[i][p+cindex*this->chunk_size] = temp + state_dimension(particle_type)*p;
                 }
             }
     }
@@ -615,14 +615,14 @@ template <particle_types particle_type, class rnumber, int interp_neighbours>
 void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::write(
         const bool write_rhs)
 {
-    double *temp0 = new double[this->chunk_size*this->ncomponents];
-    double *temp1 = new double[this->chunk_size*this->ncomponents];
+    double *temp0 = new double[this->chunk_size*state_dimension(particle_type)];
+    double *temp1 = new double[this->chunk_size*state_dimension(particle_type)];
     int zmin_rank, zmax_rank;
     int pindex = 0;
     for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
     {
         //write state
-        std::fill_n(temp0, this->ncomponents*this->chunk_size, 0);
+        std::fill_n(temp0, state_dimension(particle_type)*this->chunk_size, 0);
         pindex = cindex*this->chunk_size;
         for (int p=0; p<this->chunk_size; p++, pindex++)
         {
@@ -630,14 +630,14 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::wri
                 this->domain_particles[ 0].find(pindex) != this->domain_particles[ 0].end())
             {
                 std::copy(this->state[pindex].data,
-                          this->state[pindex].data + this->ncomponents,
-                          temp0 + p*this->ncomponents);
+                          this->state[pindex].data + state_dimension(particle_type),
+                          temp0 + p*state_dimension(particle_type));
             }
         }
         MPI_Allreduce(
                 temp0,
                 temp1,
-                this->ncomponents*this->chunk_size,
+                state_dimension(particle_type)*this->chunk_size,
                 MPI_DOUBLE,
                 MPI_SUM,
                 this->comm);
@@ -647,7 +647,7 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::wri
         if (write_rhs)
             for (int i=0; i<this->integration_steps; i++)
             {
-                std::fill_n(temp0, this->ncomponents*this->chunk_size, 0);
+                std::fill_n(temp0, state_dimension(particle_type)*this->chunk_size, 0);
                 pindex = cindex*this->chunk_size;
                 for (int p=0; p<this->chunk_size; p++, pindex++)
                 {
@@ -655,14 +655,14 @@ void rFFTW_distributed_particles<particle_type, rnumber, interp_neighbours>::wri
                         this->domain_particles[ 0].find(pindex) != this->domain_particles[ 0].end())
                     {
                         std::copy(this->rhs[i][pindex].data,
-                                  this->rhs[i][pindex].data + this->ncomponents,
-                                  temp0 + p*this->ncomponents);
+                                  this->rhs[i][pindex].data + state_dimension(particle_type),
+                                  temp0 + p*state_dimension(particle_type));
                     }
                 }
                 MPI_Allreduce(
                         temp0,
                         temp1,
-                        this->ncomponents*this->chunk_size,
+                        state_dimension(particle_type)*this->chunk_size,
                         MPI_DOUBLE,
                         MPI_SUM,
                         this->comm);