diff --git a/data/see/constant_yield.dat b/data/see/constant_yield.dat
new file mode 100644
index 0000000..060c8ad
--- /dev/null
+++ b/data/see/constant_yield.dat
@@ -0,0 +1,4 @@
+#Relative energy (eV) yield ()
+0.000 1.000E-01
+1.000 1.000E-01
+
diff --git a/doc/user-manual/.gitignore b/doc/user-manual/.gitignore
index b2ba4d5..61ad7e0 100644
--- a/doc/user-manual/.gitignore
+++ b/doc/user-manual/.gitignore
@@ -19,3 +19,4 @@ fpakc_UserManual-blx.bib
 *.gls
 *.ist
 fpakc_UserManual.run.xml
+*.bib.sav
diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90
index cd553a8..657ca28 100644
--- a/src/modules/common/moduleRandom.f90
+++ b/src/modules/common/moduleRandom.f90
@@ -66,6 +66,24 @@ MODULE moduleRandom
 
     END FUNCTION randomMaxwellian
 
+    FUNCTION randomHalfMaxwellian() RESULT(rnd)
+      USE moduleConstParam, ONLY: PI
+      IMPLICIT NONE
+
+      REAL(8):: rnd
+      REAL(8):: x, y
+
+      rnd = 0.D0
+      x = 0.D0
+      DO WHILE (x == 0.D0)
+        CALL RANDOM_NUMBER(x)
+      END DO
+      y = random(-PI, PI)
+
+      rnd = DSQRT(-2.D0*DLOG(x))*DCOS(y)
+
+    END FUNCTION randomHalfMaxwellian
+
     !Returns a random number weighted with the cumWeight array
     FUNCTION randomWeighted(cumWeight, sumWeight) RESULT(rnd)
       IMPLICIT NONE
diff --git a/src/modules/common/moduleTable.f90 b/src/modules/common/moduleTable.f90
index e4da335..ac5cbc7 100644
--- a/src/modules/common/moduleTable.f90
+++ b/src/modules/common/moduleTable.f90
@@ -30,6 +30,8 @@ MODULE moduleTable
         READ(id, '(A)', iostat = stat) dummy
         !If EOF or error, exit file
         IF (stat /= 0) EXIT
+        !If empty line, exit file
+        IF (dummy == "") EXIT
         !Skip comment
         IF (INDEX(dummy,'#') /= 0) CYCLE
         !Add row
@@ -55,6 +57,7 @@ MODULE moduleTable
         !TODO: Make this a function
         IF (INDEX(dummy,'#') /= 0) CYCLE
         IF (stat /= 0) EXIT
+        IF (dummy == "") EXIT
         !Add data
         !TODO: substitute with extracting information from dummy
         BACKSPACE(id)
diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90
index 9891806..9f16195 100644
--- a/src/modules/init/moduleInput.f90
+++ b/src/modules/init/moduleInput.f90
@@ -339,10 +339,9 @@ MODULE moduleInput
       !Mean velocity and temperature at particle position
       REAL(8):: velocityXi(1:3), temperatureXi
       INTEGER:: nNewPart = 0.D0
-      CLASS(meshCell), POINTER:: cell
       TYPE(particle), POINTER:: partNew
+      CLASS(meshCell), POINTER:: cell
       REAL(8):: vTh
-      TYPE(lNode), POINTER:: partCurr, partNext
 
       CALL config%info('solver.initial', found, n_children = nInitial)
 
@@ -371,9 +370,12 @@ MODULE moduleInput
             !Calculate number of particles
             nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
 
-            !Allocate new particles
+            !Allocate array of particles for this species
+            ALLOCATE(partOld(sp)%p(1:nNewPart))
+
+            !Create new particles
             DO p = 1, nNewPart
-              ALLOCATE(partNew)
+              partNew = partOld(sp)%p(p)
               partNew%species => species(sp)%obj
               partNew%r       =  mesh%cells(e)%obj%randPos()
               partNew%Xi      =  mesh%cells(e)%obj%phy2log(partNew%r)
@@ -408,9 +410,6 @@ MODULE moduleInput
 
               partNew%weight = species(sp)%obj%weight
 
-              !Assign particle to temporal list of particles
-              CALL partInitial%add(partNew)
-
               !Assign particle to list in volume
               IF (listInCells) THEN
                 cell => meshforMCC%cells(partNew%cellColl)%obj
@@ -430,30 +429,10 @@ MODULE moduleInput
 
           END DO
 
+          nPartOld(sp) = SIZE(partOld(sp)%p)
+
         END DO
 
-        !Convert temporal list of particles into initial partOld array
-        !Deallocate the list of initial particles
-        nNewPart = partInitial%amount
-        IF (nNewPart > 0) THEN
-          ALLOCATE(partOld(1:nNewPart))
-          partCurr => partInitial%head
-          DO p = 1, nNewPart
-            partNext => partCurr%next
-            partOld(p) = partCurr%part
-            DEALLOCATE(partCurr)
-            partCurr => partNext
-
-          END DO
-
-          IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head)
-          IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail)
-          partInitial%amount = 0
-
-        END IF
-
-        nPartOld = SIZE(partOld)
-
       END IF
 
     END SUBROUTINE readInitial
@@ -600,7 +579,10 @@ MODULE moduleInput
 
       END DO
 
+      !Allocate the wrapper array that contains particles
+      ALLOCATE(partOld(1:nSpecies))
       !Set number of particles to 0 for init state
+      ALLOCATE(nPartOld(1:nSpecies))
       nPartOld = 0
 
       !Initialize the lock for the non-analogue (NA) list of particles
@@ -806,7 +788,7 @@ MODULE moduleInput
       REAL(8):: effTime
       REAL(8):: eThreshold !Energy threshold
       INTEGER:: speciesID
-      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection
+      CHARACTER(:), ALLOCATABLE:: speciesName, crossSection, yield
       LOGICAL:: found
       INTEGER:: nTypes
 
@@ -872,6 +854,15 @@ MODULE moduleInput
 
             CALL initWallTemperature(boundary(i)%bTypes(s)%obj, Tw, cw)
 
+          CASE('secondaryEmission')
+            CALL config%get(object // '.yield', yield, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'yield' for secondary emission", 'readBoundary')
+            CALL config%get(object // '.electron', speciesName, found)
+            IF (.NOT. found) CALL criticalError("missing parameter 'electron' for secondary emission", 'readBoundary')
+            speciesID = speciesName2Index(speciesName)
+
+            CALL initSEE(boundary(i)%bTypes(s)%obj, yield, speciesID)
+
           CASE('axis')
             ALLOCATE(boundaryAxis:: boundary(i)%bTypes(s)%obj)
 
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index 7e30326..203047f 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -816,9 +816,8 @@ MODULE moduleMesh
 
       IF (MOD(t, everyColl) == 0) THEN
         !Collisions need to be performed in this iteration
-        !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
+        !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j, cell)
         DO e=1, self%numCells
-
           cell => self%cells(e)%obj
 
           !TODO: Simplify this, to many sublevels
@@ -885,7 +884,7 @@ MODULE moduleMesh
                   !Obtain the cross sections for the different processes
                   !TODO: From here it might be a procedure in interactionMatrix
                   vRel = NORM2(part_i%v-part_j%v)
-                  rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
+                  rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass)
                   eRel = rMass*vRel**2
                   CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
 
@@ -1080,7 +1079,7 @@ MODULE moduleMesh
 
             !Compute changes in velocity for each particle
             deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
-            mass_ij(p)       = pair%sp_i%m*partTemp%part%weight
+            mass_ij(p)       = pair%sp_i%mass*partTemp%part%weight
             p_ij(p,1:3)      = mass_ij(p)*partTemp%part%v
 
             !Move to the next particle in the list
@@ -1163,7 +1162,7 @@ MODULE moduleMesh
 
               !Compute changes in velocity for each particle
               deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
-              mass_ji(p)       = pair%sp_j%m*partTemp%part%weight
+              mass_ji(p)       = pair%sp_j%mass*partTemp%part%weight
               p_ji(p,1:3)      = mass_ji(p)*partTemp%part%v
 
               !Move to the next particle in the list
diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90
index 4f72e10..e7e4d65 100644
--- a/src/modules/mesh/moduleMeshBoundary.f90
+++ b/src/modules/mesh/moduleMeshBoundary.f90
@@ -125,7 +125,7 @@ MODULE moduleMeshBoundary
 
       SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
       TYPE IS(boundaryIonization)
-        mRel = reducedMass(bound%m0, part%species%m)
+        mRel = reducedMass(bound%m0, part%species%mass)
         vRel = SUM(DABS(part%v-bound%v0))
         eRel = mRel*vRel**2*5.D-1
 
@@ -193,7 +193,7 @@ MODULE moduleMeshBoundary
 
       END SELECT
 
-      !Removes ionizing electron regardless the number of pair created
+      !Removes ionizing electron regardless the number of pairs created
       part%n_in = .FALSE.
 
     END SUBROUTINE ionization
@@ -212,6 +212,116 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE symmetryAxis
 
+    !Secondary emission of electrons by impacting particle.
+    SUBROUTINE secondaryEmission(edge, part)
+      USE moduleSpecies
+      USE moduleRandom
+      USE moduleConstParam
+      IMPLICIT NONE
+
+      CLASS(meshEdge), INTENT(inout):: edge
+      CLASS(particle), INTENT(inout):: part
+      REAL(8):: vRel, eRel
+      REAL(8), DIMENSION(1:3):: rElectron, XiElectron!Position of new electrons (impacting particle)
+      REAL(8), DIMENSION(1:3):: rIncident !Vector from imapcting particle position to particle position
+      REAL(8):: theta !incident angle
+      REAL(8):: yield
+      REAL(8):: energy !incident energy corrected by threshold and 
+      INTEGER:: nNewElectrons
+      REAL(8), ALLOCATABLE:: weight(:)
+      INTEGER:: p
+      INTEGER:: cell
+      TYPE(particle), POINTER:: secondaryElectron
+
+      SELECT TYPE(bound => edge%boundary%bTypes(part%species%n)%obj)
+      TYPE IS(boundarySEE)
+        !Get relative velocity
+        vRel = NORM2(part%v)
+        !Convert to relative energy
+        eRel = part%species%mass*vRel**2*5.D-1
+
+        !If energy is abound threshold calculate secondary electrons
+        IF (eRel >= bound%energyThreshold) THEN
+
+          !position of impacting ion
+          rElectron  = edge%intersection(part%r)
+          XiElectron = mesh%cells(part%cell)%obj%phy2log(rElectron)
+
+          !Calculate incident angle
+          rIncident = part%r - rElectron
+          theta = ACOS(DOT_PRODUCT(rIncident, edge%normal) / (NORM2(rIncident) * NORM2(edge%normal)))
+
+          !Calculate incident energy
+          energy = (eRel - bound%energyThreshold) * PI2 / (PI2 + theta**2) + bound%energyThreshold
+
+          !Get number of secondary electrons particles
+          yield = part%weight*bound%yield%get(eRel) * (1.D0 * theta**2 / PI2) !Check equation!
+          
+          !Convert the number to macro-particles
+          nNewElectrons = FLOOR(yield / bound%electron%weight)
+
+          !If the weight of new macro-particles is below the yield, correct adding a particle
+          IF (REAL(nNewElectrons) * bound%electron%weight < yield) THEN
+            nNewElectrons = nNewElectrons + 1
+            ALLOCATE(weight(1:nNewElectrons))
+            weight(1:nNewElectrons-1) = bound%electron%weight
+            weight(nNewElectrons)     = yield - SUM(weight(1:nNewElectrons-1))
+
+          ELSE
+            ALLOCATE(weight(1:nNewElectrons))
+            weight(1:nNewElectrons) = bound%electron%weight
+          
+          END IF
+
+          !New cell of origin
+          IF (ASSOCIATED(edge%e1)) THEN
+            cell = edge%e1%n
+
+          ELSEIF (ASSOCIATED(edge%e2)) THEN
+            cell = edge%e2%n
+
+          END IF
+
+          !Create the new electron macro-particles
+          DO p = 1, nNewElectrons
+            !Create new macro-particle
+            ALLOCATE(secondaryElectron)
+
+            !Assign species to electron
+            secondaryElectron%species => bound%electron
+
+            !Assign position to particle
+            secondaryElectron%r  = rElectron
+            secondaryElectron%Xi = XiElectron
+
+            !Assign cell to electron
+            secondaryElectron%cell = cell
+
+            !Assign weight
+            secondaryElectron%weight = weight(p)
+
+            !Assume particle is inside the numerical domain
+            secondaryElectron%n_in = .TRUE.
+
+            !Assign velocity
+            secondaryElectron%v = 2.D0 * edge%normal + 1.D-1 * (/ randomMaxwellian(), randomMaxwellian(), randomMaxwellian() /)
+
+            !Add particle to list
+            CALL partSurfaces%setLock()
+            CALL partSurfaces%add(secondaryElectron)
+            CALL partSurfaces%unsetLock()
+
+          END DO
+
+        END IF
+
+        !Absorb impacting particle
+        CALL absorption(edge, part)
+
+      END SELECT
+
+    END SUBROUTINE secondaryEmission
+
     !Points the boundary function to specific type
     SUBROUTINE pointBoundaryFunction(edge, s)
       USE moduleErrors
@@ -239,6 +349,9 @@ MODULE moduleMeshBoundary
       TYPE IS(boundaryAxis)
         edge%fBoundary(s)%apply => symmetryAxis
 
+      TYPE IS(boundarySEE)
+        edge%fBoundary(s)%apply => secondaryEmission
+
       CLASS DEFAULT
         CALL criticalError("Boundary type not defined in this geometry", 'pointBoundaryFunction')
 
diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90
index 83c815c..b3c1d12 100644
--- a/src/modules/moduleBoundary.f90
+++ b/src/modules/moduleBoundary.f90
@@ -46,6 +46,16 @@ MODULE moduleBoundary
 
   END TYPE boundaryIonization
 
+  !Secondary electron emission (by ion impact)
+  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundarySEE
+    !Yield as a function of ion energy
+    TYPE(table1D):: yield
+    CLASS(speciesGeneric), POINTER:: electron !Electron species for secondary emission
+    REAL(8):: energyThreshold
+    CONTAINS
+
+  END TYPE boundarySEE
+
   !Symmetry axis
   TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
     CONTAINS
@@ -137,4 +147,36 @@ MODULE moduleBoundary
 
     END SUBROUTINE initIonization
 
+    SUBROUTINE initSEE(boundary, tableFile, speciesID)
+      USE moduleRefParam
+      USE moduleConstParam
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: tableFile
+      INTEGER:: speciesID
+      INTEGER:: i
+
+      ALLOCATE(boundarySEE:: boundary)
+      SELECT TYPE(boundary)
+      TYPE IS(boundarySEE)
+        CALL boundary%yield%init(tableFile)
+        CALL boundary%yield%convert(eV2J/(m_ref*v_ref**2), 1.D0)
+        boundary%electron => species(speciesID)%obj
+        !Search for the threshold energy in the table
+        DO i = 1, SIZE(boundary%yield%f)
+          IF (boundary%yield%f(i) > 0.D0) THEN
+            boundary%energyThreshold = boundary%yield%x(i)
+
+            EXIT
+
+          END IF
+
+        END DO
+
+      END SELECT
+
+    END SUBROUTINE initSEE
+
 END MODULE moduleBoundary
diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90
index ccca930..f01c20b 100644
--- a/src/modules/moduleCollisions.f90
+++ b/src/modules/moduleCollisions.f90
@@ -164,8 +164,8 @@ MODULE moduleCollisions
 
       self%amount = amount
 
-      mass_i = species(i)%obj%m
-      mass_j = species(j)%obj%m
+      mass_i = species(i)%obj%mass
+      mass_j = species(j)%obj%mass
 
       ALLOCATE(self%collisions(1:self%amount))
 
@@ -227,8 +227,8 @@ MODULE moduleCollisions
       REAL(8):: m_i, m_j
       REAL(8), DIMENSION(1:3):: vCM, vp
 
-      m_i = part_i%species%m
-      m_j = part_j%species%m
+      m_i = part_i%species%mass
+      m_j = part_j%species%mass
       !Applies the collision
       vCM = velocityCM(part_i%weight*m_i, part_i%v, part_j%weight*m_j, part_j%v)
       vp  = vRel*randomDirectionVHS()
@@ -283,7 +283,7 @@ MODULE moduleCollisions
         END SELECT
 
         !momentum change per ionization process
-        collision%deltaV = sqrt(collision%eThreshold / collision%electron%m)
+        collision%deltaV = sqrt(collision%eThreshold / collision%electron%mass)
 
       END SELECT
 
@@ -307,7 +307,7 @@ MODULE moduleCollisions
       REAL(8), DIMENSION(1:3):: vChange
       TYPE(particle), POINTER:: newElectron => NULL(), remainingNeutral => NULL()
 
-      rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
+      rMass = reducedMass(part_i%weight*part_i%species%mass, part_j%weight*part_j%species%mass)
       eRel = rMass*vRel**2
       !Relative energy must be higher than threshold
       IF (eRel > self%eThreshold) THEN
diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90
index 166b6b2..da03ad8 100644
--- a/src/modules/moduleCoulomb.f90
+++ b/src/modules/moduleCoulomb.f90
@@ -61,7 +61,7 @@ MODULE moduleCoulomb
       self%sp_i => species(i)%obj
       self%sp_j => species(j)%obj
 
-      self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m
+      self%one_plus_massRatio_ij = 1.D0 + self%sp_i%mass/self%sp_j%mass
 
       Z_i = 0.D0
       Z_j = 0.D0
@@ -87,11 +87,11 @@ MODULE moduleCoulomb
 
       scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3)
 
-      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent
-      self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent
+      self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%mass**2) * scaleFactor !Missing density because it's cell dependent
+      self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%mass**2) * scaleFactor !Missing density because it's cell dependent
 
-      self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent
-      self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent
+      self%l2_j = self%sp_j%mass / 2.D0 !Missing temperature because it's cell dependent
+      self%l2_i = self%sp_i%mass / 2.D0 !Missing temperature because it's cell dependent
 
     END SUBROUTINE initInteractionCoulomb
 
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 4e57083..b5babfb 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -179,16 +179,21 @@ MODULE moduleInject
       IMPLICIT NONE
 
       INTEGER:: i
+      INTEGER, DIMENSION(1:nInject):: nMin, nMax
       
       !$OMP SINGLE
       nPartInj = 0
       DO i = 1, nInject
         IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
+          nMin(i)  = nPartInj + 1 
           nPartInj = nPartInj + inject(i)%nParticles
+          nMax(i)  = nPartInj
 
         END IF
 
       END DO
+      PRINT *, nMin
+      PRINT *, nMax
 
       IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
       ALLOCATE(partInj(1:nPartInj))
@@ -196,7 +201,7 @@ MODULE moduleInject
 
       DO i=1, nInject
         IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
-          CALL inject(i)%addParticles()
+          CALL inject(i)%addParticles(nMin(i), nMax(i))
 
         END IF
       END DO
@@ -254,10 +259,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      DO WHILE (v <= 0.D0)
-        v = self%vTh*randomMaxwellian()
-
-      END DO
+      v = self%vTh*randomHalfMaxwellian()
 
     END FUNCTION randomVelHalfMaxwellian
 
@@ -282,8 +284,8 @@ MODULE moduleInject
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: nMin, nMax !Min and Max index in partInj array
       INTEGER:: randomX
-      INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
       INTEGER:: i
       INTEGER:: n, sp
       CLASS(meshEdge), POINTER:: randomEdge
@@ -291,17 +293,6 @@ MODULE moduleInject
 
       !Insert particles
       !$OMP SINGLE
-      nMin = 0
-      DO i = 1, self%id -1
-        IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN 
-          nMin = nMin + inject(i)%nParticles
-
-        END IF
-
-      END DO
-
-      nMin = nMin + 1
-      nMax = nMin + self%nParticles - 1
       !Assign weight to particle.
       partInj(nMin:nMax)%weight = self%species%weight
       !Particle is considered to be outside the domain
@@ -334,17 +325,19 @@ MODULE moduleInject
 
         direction = self%n
 
-        partInj(n)%v = 0.D0
-
+        !Sample initial velocity
         partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
                                                  self%v(2)%obj%randomVel(), &
                                                  self%v(3)%obj%randomVel() /)
 
-        !If velocity is not in the right direction, invert it
-        IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN
-          partInj(n)%v = - partInj(n)%v
+        !For each direction, velocities have to agree with the direction of injection
+        DO i = 1, 3
+          DO WHILE (partInj(n)%v(i)*direction(i) < 0)
+            partInj(n)%v(i) = self%vMod*direction(i) + self%v(i)%obj%randomVel()
 
-        END IF
+          END DO
+
+        END DO
 
         !Obtain natural coordinates of particle in cell
         partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90
index ab08f08..d039247 100644
--- a/src/modules/moduleSpecies.f90
+++ b/src/modules/moduleSpecies.f90
@@ -4,12 +4,54 @@ MODULE moduleSpecies
   USE OMP_LIB
   IMPLICIT NONE
   
+  !Basic type that defines a macro-particle
+  TYPE:: particle
+    REAL(8):: r(1:3)                         !Position
+    REAL(8):: v(1:3)                         !Velocity
+    CLASS(speciesGeneric), POINTER:: species !Pointer to particle's species
+    INTEGER:: cell                           !Index of element in which the particle is located TODO: Make these pointers
+    INTEGER:: cellColl                       !Index of element in which the particle is located in the Collision Mesh
+    REAL(8):: Xi(1:3)                        !Logical coordinates of particle in cell
+    LOGICAL:: n_in                           !Flag that indicates if a particle is in the domain
+    REAL(8):: weight=0.D0                    !weight of particle
+
+  END TYPE particle
+
+  !Wrapper to store the particles per species
+  TYPE:: particleArray
+    TYPE(particle), ALLOCATABLE, DIMENSION(:):: p
+
+  END TYPE particleArray
+
+  !Array of pointers for the particles to be pushed
+  TYPE:: particlePointer
+    TYPE(particle), POINTER:: p
+
+  END TYPE particlePointer
+
+  !Arrays that contain the particles
+  TYPE(particleArray),   ALLOCATABLE, TARGET, DIMENSION(:):: particles           !array of particles in the domain
+  TYPE(particle),        ALLOCATABLE, TARGET, DIMENSION(:):: particlesInjection  !array of inject particles
+  TYPE(particlePointer), ALLOCATABLE,         DIMENSION(:):: particlesToPush     !particles pushed in each iteration
+  !Integers to track number of particles in domain
+  INTEGER, ALLOCATABLE, DIMENSION(:):: nParticles
+  INTEGER:: nParticlesTotal
+  !Number of injected particles
+  INTEGER, ALLOCATABLE, DIMENSION(:):: nParticlesInject
+  INTEGER:: nPariclesInjectTotal
+
+  !Generic species type
   TYPE, ABSTRACT:: speciesGeneric
-    CHARACTER(:), ALLOCATABLE:: name
-    REAL(8):: m=0.D0, weight=0.D0, qm=0.D0
-    INTEGER:: n=0
+    INTEGER:: n=0 !Index of species
+    CHARACTER(:), ALLOCATABLE:: name !Name of species
+    !Mass, default weight of species and charge over mass
+    REAL(8):: mass=0.D0, weight=0.D0, qm=0.D0
+    INTEGER:: every !How many interations between advancing the species
+    LOGICAL:: pushSpecies !Boolean to indicate if the species is moved in the iteration
+
   END TYPE speciesGeneric
 
+  !Neutral species
   TYPE, EXTENDS(speciesGeneric):: speciesNeutral
     CLASS(speciesGeneric), POINTER:: ion => NULL()
     CONTAINS
@@ -17,6 +59,7 @@ MODULE moduleSpecies
 
   END TYPE speciesNeutral
 
+  !Charged species
   TYPE, EXTENDS(speciesGeneric):: speciesCharged
     REAL(8):: q=0.D0
     CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL()
@@ -26,34 +69,17 @@ MODULE moduleSpecies
 
   END TYPE speciesCharged
 
+  !Wrapper for species
   TYPE:: speciesCont
     CLASS(speciesGeneric), ALLOCATABLE:: obj
 
   END TYPE
 
+  !Number of species
   INTEGER:: nSpecies
+  !Array for species
   TYPE(speciesCont), ALLOCATABLE, TARGET:: species(:)
 
-  TYPE particle
-    REAL(8):: r(1:3) !Position
-    REAL(8):: v(1:3) !Velocity
-    CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle
-    INTEGER:: cell !Index of element in which the particle is located
-    INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh
-    REAL(8):: Xi(1:3) !Logical coordinates of particle in element e_p.
-    LOGICAL:: n_in !Flag that indicates if a particle is in the domain
-    REAL(8):: weight=0.D0 !weight of particle
-
-  END TYPE particle
-
-  !Number of old particles
-  INTEGER:: nPartOld
-  !Number of injected particles
-  INTEGER:: nPartInj
-  !Arrays that contain the particles
-  TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration
-  TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles
-
   CONTAINS
     FUNCTION speciesName2Index(speciesName) RESULT(sp)
       USE moduleErrors
diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90
index 18fbb7f..16ed1b4 100644
--- a/src/modules/output/moduleOutput.f90
+++ b/src/modules/output/moduleOutput.f90
@@ -148,8 +148,8 @@ MODULE moduleOutput
         formatValues%density     = rawValues%den*tempVol
         formatValues%velocity(:) = tempVel
         IF (tensorTrace(tensorTemp) > 0.D0) THEN
-         formatValues%pressure    = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0
-         formatValues%temperature = formatValues%pressure/(formatValues%density*kb)
+          formatValues%pressure    = speciesIn%mass*tensorTrace(tensorTemp)*tempVol/3.D0
+          formatValues%temperature = formatValues%pressure/(formatValues%density*kb)
 
         END IF
       END IF
@@ -187,7 +187,7 @@ MODULE moduleOutput
 
         OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
 
-        WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
+        WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nParticlesTotal, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField
 
         CLOSE(20)
 
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index e539fed..96e3f78 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -3,7 +3,7 @@ MODULE moduleSolver
 
   !Generic type for pusher of particles
   TYPE, PUBLIC:: pusherGeneric
-    PROCEDURE(push_interafece), POINTER, NOPASS:: pushParticle => NULL()
+    PROCEDURE(push_interface), POINTER, NOPASS:: pushParticle => NULL()
     !Boolean to indicate if the species is moved in the iteration
     LOGICAL:: pushSpecies
     !How many interations between advancing the species
@@ -29,13 +29,13 @@ MODULE moduleSolver
 
   INTERFACE
     !Push a particle
-    PURE SUBROUTINE push_interafece(part, tauIn)
+    PURE SUBROUTINE push_interface(part, tauIn)
       USE moduleSpecies
 
       TYPE(particle), INTENT(inout):: part
       REAL(8), INTENT(in):: tauIn
 
-    END SUBROUTINE push_interafece
+    END SUBROUTINE push_interface
 
     !Solve the electromagnetic field
     SUBROUTINE solveEM_interface()
@@ -123,7 +123,12 @@ MODULE moduleSolver
       END SELECT
 
       self%pushSpecies = .FALSE.
-      self%every = INT(tauSp/tau)
+      self%every = FLOOR(tauSp/tau)
+      !Correct value if not fulfilled
+      IF (tau*REAL(self%every) < tauSp) THEN
+        self%every = self%every + 1
+
+      END IF
 
     END SUBROUTINE initPusher
 
@@ -164,24 +169,75 @@ MODULE moduleSolver
       USE moduleMesh
       IMPLICIT NONE
 
-      INTEGER:: n
-      INTEGER:: sp
+      INTEGER:: p
+      INTEGER:: s, sp
+      INTEGER:: e
 
-      !$OMP DO
-      DO n = 1, nPartOld
-        !Select species type
-        sp = partOld(n)%species%n
+      !$OMP SINGLE
+      !Allocate the array of particles to push
+      nSpeciesToPush = COUNT(solver%pusher(:)%pushSpecies)
+      ALLOCATE(particlesToPush(1:nSpeciesToPush))
+      ALLOCATE(nPartOldToPush(1:nSpeciesToPush))
+      !Point the arrays to the location of the particles
+      sp = 0
+      DO s = 1, nSpecies
         !Checks if the species sp is update this iteration
         IF (solver%pusher(sp)%pushSpecies) THEN
-          !Push particle
-          CALL solver%pusher(sp)%pushParticle(partOld(n), tau(sp))
-          !Find cell in wich particle reside
-          CALL solver%updateParticleCell(partOld(n))
+          sp = sp + 1
+          particlesToPush(sp)%p = partOld(s)%p
+          nPartOldToPush(sp)    = nPartOld(s)
 
         END IF
 
       END DO
-      !$OMP END DO
+
+      !Delete list of particles in cells for collisions if the particle is pushed
+      IF (listInCells) THEN
+        DO e = 1, mesh%numCells
+          DO s = 1, nSpecies
+            IF (solver%pusher(s)%pushSpecies) THEN 
+              CALL mesh%cells(e)%obj%listPart_in(s)%erase()
+              mesh%cells(e)%obj%totalWeight(s) = 0.D0
+
+            END IF
+
+          END DO
+
+        END DO
+
+      END IF
+
+      !Erase the list of particles inside the cell in coll mesh if the particle is pushed
+      IF (doubleMesh) THEN
+        DO e = 1, meshColl%numCells
+          DO s = 1, nSpecies
+            IF (solver%pusher(s)%pushSpecies) THEN 
+              CALL meshColl%cells(e)%obj%listPart_in(s)%erase()
+              meshColl%cells(e)%obj%totalWeight(s) = 0.D0
+
+            END IF
+
+          END DO
+
+        END DO
+    
+      END IF
+      !$OMP END SINGLE
+
+      !Now, push particles
+      !$OMP DO
+      DO sp = 1, nSpeciesToPush
+        DO p = 1, nPartOldToPush(sp)
+          !Push particle
+          CALL solver%pusher(sp)%pushParticle(particlesToPush(sp)%p(p), tau(sp))
+          !Find cell in which particle reside
+          CALL solver%updateParticleCell(particlesToPush(sp)%p(p))
+
+        END DO
+
+      END DO
+      !$END OMP DO
+
 
     END SUBROUTINE doPushes
 
@@ -242,7 +298,10 @@ MODULE moduleSolver
       !$OMP SECTION
       nOldIn = 0
       IF (ALLOCATED(partOld)) THEN
-        nOldIn = COUNT(partOld%n_in)
+        DO s = 1, nSpecies
+          nOldIn = nOldin + COUNT(partOld(s)%p(:)%n_in)
+
+        END DO
 
       END IF
       !$OMP SECTION
@@ -319,40 +378,6 @@ MODULE moduleSolver
 
       END DO
 
-      !$OMP SECTION
-      !Erase the list of particles inside the cell if particles have been pushed
-      IF (listInCells) THEN
-        DO s = 1, nSpecies
-          DO e = 1, mesh%numCells
-            IF (solver%pusher(s)%pushSpecies) THEN 
-              CALL mesh%cells(e)%obj%listPart_in(s)%erase()
-              mesh%cells(e)%obj%totalWeight(s) = 0.D0
-
-            END IF
-
-          END DO
-
-        END DO
-
-      END IF
-
-      !$OMP SECTION
-      !Erase the list of particles inside the cell in coll mesh
-      IF (doubleMesh) THEN
-        DO s = 1, nSpecies
-          DO e = 1, meshColl%numCells
-            IF (solver%pusher(s)%pushSpecies) THEN 
-              CALL meshColl%cells(e)%obj%listPart_in(s)%erase()
-              meshColl%cells(e)%obj%totalWeight(s) = 0.D0
-
-            END IF
-
-          END DO
-
-        END DO
-      
-      END IF
-
       !$OMP END SECTIONS
 
       !$OMP SINGLE
@@ -367,14 +392,17 @@ MODULE moduleSolver
       USE moduleMesh
       IMPLICIT NONE
 
-      INTEGER:: n
+      INTEGER:: n, sp
       CLASS(meshCell), POINTER:: cell
 
       !Loops over the particles to scatter them
       !$OMP DO
-      DO n = 1, nPartOld
-        cell => mesh%cells(partOld(n)%cell)%obj
-        CALL cell%scatter(cell%nNodes, partOld(n))
+      DO sp = 1, nSpeciesToPush
+        DO n = 1, nPartOldToPush(sp)
+          cell => mesh%cells(particlesToPush(sp)%p(n)%cell)%obj
+          CALL cell%scatter(cell%nNodes, particlesToPush(sp)%p(n))
+
+        END DO
 
       END DO
       !$OMP END DO
@@ -544,8 +572,8 @@ MODULE moduleSolver
 
         END IF
 
-        IF (nPartOld > 0) THEN
-          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns"
+        IF (nPartOldTotal > 0) THEN
+          WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOldTotal), "ns"
 
         END IF
         WRITE(*,*)
diff --git a/src/modules/solver/pusher/modulePusher.f90 b/src/modules/solver/pusher/modulePusher.f90
index f75c733..4c7e372 100644
--- a/src/modules/solver/pusher/modulePusher.f90
+++ b/src/modules/solver/pusher/modulePusher.f90
@@ -14,7 +14,7 @@ MODULE modulePusher
     END SUBROUTINE pushCartNeutral
 
     PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
-      USE moduleSPecies
+      USE moduleSpecies
       USE moduleMesh
       IMPLICIT NONE