diff --git a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90 b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90
index 9b1a442..203de0c 100644
--- a/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90
+++ b/src/modules/mesh/inout/vtk/moduleMeshOutputVTK.f90
@@ -1,22 +1,44 @@
 MODULE moduleMeshOutputVTK
 
+  CHARACTER(LEN=6):: prefix = 'OUTPUT'
+
   CONTAINS
 
-    SUBROUTINE writeFileHeader(self, fileID)
+    PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName)
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CHARACTER(*), INTENT(in):: prefix, suffix, extension
+      INTEGER, INTENT(in), OPTIONAL:: t
+      CHARACTER (LEN=iterationDigits):: tString
+      CHARACTER(:), ALLOCATABLE:: fileName
+
+      IF (PRESENT(t)) THEN
+        WRITE(tString, iterationFormat) t
+        fileName = prefix //  '_' // tString // '_' // suffix // '.' // extension
+
+      ELSE
+        fileName = prefix //  '_' // suffix // '.' // extension
+
+      END IF
+
+    END FUNCTION formatFileName
+
+    SUBROUTINE writeHeader(nNodes, nCells, fileID)
       USE moduleMesh
       IMPLICIT NONE
 
-      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes, nCells
       INTEGER, INTENT(in):: fileID
 
       WRITE(fileID,"(A)") '<?xml version="1.0"?>'
       WRITE(fileID,"(2X, A)") '<VTKFile type="UnstructuredGrid">'
       WRITE(fileID,"(4X, A,ES20.6E3,A)") '<UnstructuredGrid>'
-      WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', self%numNodes, '" NumberOfCells="', self%numCells, '">'
+      WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', nNodes, '" NumberOfCells="', nCells, '">'
 
-    END SUBROUTINE writeFileHeader
+    END SUBROUTINE writeHeader
 
-    SUBROUTINE writeFileFooter(fileID)
+    SUBROUTINE writeFooter(fileID)
       IMPLICIT NONE
 
       INTEGER, INTENT(in):: fileID
@@ -25,7 +47,7 @@ MODULE moduleMeshOutputVTK
       WRITE(fileID,"(4X, A)") '</UnstructuredGrid>'
       WRITE(fileID,"(2X, A)") '</VTKFile>'
 
-    END SUBROUTINE writeFileFooter
+    END SUBROUTINE writeFooter
 
     FUNCTION getCellType(cell) RESULT(indexType)
       USE moduleMesh3DCart
@@ -74,12 +96,12 @@ MODULE moduleMeshOutputVTK
 
     END FUNCTION getCellType
 
-    SUBROUTINE writeFileMesh(self, fileID)
+    SUBROUTINE writeMesh(self, fileID)
       USE moduleMesh
       USE moduleRefParam
       IMPLICIT NONE
 
-      CLASS(meshParticles), INTENT(in):: self
+      CLASS(meshGeneric), INTENT(in):: self
       INTEGER, INTENT(in):: fileID
       CHARACTER(LEN=25):: nodeFormat
       CHARACTER(LEN=25):: cellFormat
@@ -123,104 +145,252 @@ MODULE moduleMeshOutputVTK
       WRITE(fileID, "(10X, A)") '</DataArray>'
       WRITE(fileID, "(8X, A)") '</Cells>'
 
-    END SUBROUTINE writeFileMesh
+    END SUBROUTINE writeMesh
 
-    SUBROUTINE writeNodeData(self, fileID, output)
+    SUBROUTINE writeSpeciesOutput(self, fileID, speciesIndex)
       USE moduleMesh
       USE moduleOutput
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
       INTEGER, INTENT(in):: fileID
+      INTEGER, INTENT(in):: speciesIndex
       TYPE(outputFormat):: output(1:self%numNodes)
       CHARACTER(LEN=25):: nodeFormat
       INTEGER:: n
 
-      WRITE(fileID,"(A)") '<PointData>'
-      WRITE(fileID,"(A)") '<DataArray type="Float64" Name="Density (m^-3)" NumberOfComponents="1">'
+      DO n = 1, self%numNodes
+        CALL calculateOutput(self%nodes(n)%obj%output(speciesIndex), output(n), self%nodes(n)%obj%v, species(speciesIndex)%obj)
+
+      END DO
+
+      WRITE(fileID,"(8X,A)") '<PointData>'
+      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Density (m^-3)" NumberOfComponents="1">'
       WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
       WRITE(fileID,nodeFormat) (output(n)%density, n  = 1, self%numNodes)
-      WRITE(fileID,"(A)") '</DataArray>'
-      WRITE(fileID,"(A)") '<DataArray type="Float64" Name="Velocity (m s^-1)" NumberOfComponents="3">'
+      WRITE(fileID,"(10X,A)") '</DataArray>'
+      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Velocity (m s^-1)" NumberOfComponents="3">'
       WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
       WRITE(fileID,nodeFormat) (output(n)%velocity(1:3), n  = 1, self%numNodes)
-      WRITE(fileID,"(A)") '</DataArray>'
-      WRITE(fileID,"(A)") '<DataArray type="Float64" Name="Pressure (Pa)" NumberOfComponents="1">'
+      WRITE(fileID,"(10X,A)") '</DataArray>'
+      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Pressure (Pa)" NumberOfComponents="1">'
       WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
       WRITE(fileID,nodeFormat) (output(n)%pressure, n  = 1, self%numNodes)
-      WRITE(fileID,"(A)") '</DataArray>'
-      WRITE(fileID,"(A)") '<DataArray type="Float64" Name="Temperature (K)" NumberOfComponents="1">'
+      WRITE(fileID,"(10X,A)") '</DataArray>'
+      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Temperature (K)" NumberOfComponents="1">'
       WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
       WRITE(fileID,nodeFormat) (output(n)%temperature, n  = 1, self%numNodes)
-      WRITE(fileID,"(A)") '</DataArray>'
-      WRITE(fileID,"(A)") '</PointData>'
+      WRITE(fileID,"(10X,A)") '</DataArray>'
+      WRITE(fileID,"(8X,A)") '</PointData>'
 
-    END SUBROUTINE writeNodeData
+    END SUBROUTINE writeSpeciesOutput
+
+    SUBROUTINE writeCollOutput(self,fileID)
+      USE moduleMesh
+      USE moduleCollisions
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: fileID
+      INTEGER:: k, c, n
+      CHARACTER(:), ALLOCATABLE:: title
+      CHARACTER (LEN=2):: cString
+      CHARACTER(LEN=25):: cellFormat
+
+      WRITE(fileID,"(8X,A)") '<CellData>'
+      DO k = 1, nCollPairs
+        DO c = 1, interactionMatrix(k)%amount
+          WRITE(cString, "(I2)") c
+          title = 'Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
+          WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="',title, '" NumberOfComponents="1">'
+          WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(I10))"
+          WRITE(fileID, cellFormat) (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells)
+          WRITE(fileID, "(10X, A)") '</DataArray>'
+
+        END DO
+      END DO
+      WRITE(fileID,"(8X,A)") '</CellData>'
+
+    END SUBROUTINE writeCollOutput
+
+    SUBROUTINE writeEM(self, fileID)
+      USE moduleMesh
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: fileID
+      CHARACTER(LEN=25):: nodeFormat
+      CHARACTER(LEN=25):: cellFormat
+      INTEGER:: n
+      REAL(8):: Xi(1:3)
+
+      Xi = (/ 0.D0, 0.D0, 0.D0 /)
+
+      !Points in nodes
+      WRITE(fileID,"(8X,A)") '<PointData>'
+      !Electric potential
+      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Potential (V)" NumberOfComponents="1">'
+      WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%phi*Volt_ref, n  = 1, self%numNodes)
+      WRITE(fileID,"(10X,A)") '</DataArray>'
+      !Magnetic Field
+      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Magnetic Field (T)" NumberOfComponents="3">'
+      WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
+      WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%B*B_ref, n  = 1, self%numNodes)
+      WRITE(fileID,"(10X,A)") '</DataArray>'
+      WRITE(fileID,"(8X,A)") '</PointData>'
+
+      !Cell Data
+      WRITE(fileID,"(8X,A)") '<CellData>'
+      !Electric field
+      WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="Electric Field (V m^-1)" NumberOfComponents="3">'
+      WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(3(ES20.6E3)))"
+      WRITE(fileID, cellFormat) (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
+      WRITE(fileID,"(10X,A)") '</DataArray>'
+      WRITE(fileID,"(8X,A)") '</CellData>'
+
+    END SUBROUTINE writeEM
+
+    SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
+      USE moduleCaseParam
+      USE moduleOutput
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      INTEGER:: fileID
+      INTEGER, INTENT(in):: t
+      CHARACTER(*):: fileNameStep, fileNameCollection
+
+      IF (t == tInitial) THEN
+        !Create collection file
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
+        OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection)
+        WRITE (fileID + 1, "(A)") '<VTKFile type="Collection">'
+        WRITE (fileID + 1, "(2X, A)") '<Collection>'
+        CLOSE(fileID + 1)
+
+      END IF
+
+      !Write iteration file in collection
+      OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection, ACCESS='APPEND')
+      WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
+
+      !Close collection file
+      IF (t == tFinal) THEN
+        WRITE (fileID + 1, "(2X, A)") '</Collection>'
+        WRITE (fileID + 1, "(A)") '</VTKFile>'
+
+      END IF
+      CLOSE(fileID + 1)
+
+    END SUBROUTINE writeCollection
     
     SUBROUTINE printOutputVTK(self,t)
       USE moduleMesh
-      USE moduleRefParam
       USE moduleSpecies
-      USE moduleOutput
-      USE moduleCaseParam
       IMPLICIT NONE
 
       CLASS(meshParticles), INTENT(in):: self
       INTEGER, INTENT(in):: t
       INTEGER:: n, i, fileID
       CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
-      CHARACTER (LEN=iterationDigits):: tstring
       TYPE(outputFormat):: output(1:self%numNodes)
 
       fileID = 60
 
       DO i = 1, nSpecies
-        WRITE(tstring, iterationFormat) t
-        fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.vtu'
-
+        fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t)
         WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
         OPEN (fileID, file = path // folder // '/' //  fileName)
 
-        fileNameCollection = 'OUTPUT_Collection_' // species(i)%obj%name // '.pvd'
-        IF (t == tInitial) THEN
-          !Create collection file
-          WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
-          OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection)
-          WRITE (fileID + 1, "(A)") '<VTKFile type="Collection">'
-          WRITE (fileID + 1, "(2X, A)") '<Collection>'
-          CLOSE(fileID + 1)
+        CALL writeHeader(self%numNodes, self%numCells, fileID)
 
-        END IF
+        CALL writeMesh(self, fileID)
 
-        !Write iteration file in collection
-        OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection, ACCESS='APPEND')
-        WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)"), '<DataSet timestep="', t*ti_ref,'" file="', fileName,'"/>'
+        CALL writeSpeciesOutput(self, fileID, i)
 
-        !Close collection file
-        IF (t == tFinal) THEN
-          WRITE (fileID + 1, "(2X, A)") '</Collection>'
-          WRITE (fileID + 1, "(A)") '</VTKFile>'
-
-        END IF
-        CLOSE(fileID + 1)
-
-        CALL writeFileHeader(self, fileID)
-
-        CALL writeFileMesh(self, fileID)
-
-        DO n = 1, self%numNodes
-          CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
-
-        END DO
-
-        CALL writeNodeData(self, fileID, output)
-
-        CALL writeFileFooter(fileID)
+        CALL writeFooter(fileID)
 
         CLOSE(fileID)
 
+        !Write collection file for time plotting
+        fileNameCollection = formatFileName(prefix, 'Collection_' // species(i)%obj%name, 'pvd')
+        CALL writeCollection(fileID, t, fileName, filenameCollection)
+
       END DO
 
     END SUBROUTINE printOutputVTK
 
+    SUBROUTINE printCollVTK(self,t)
+      USE moduleMesh
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+      INTEGER:: n, i, fileID
+      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
+      CHARACTER (LEN=iterationDigits):: tstring
+      TYPE(outputFormat):: output(1:self%numNodes)
+
+      fileID = 62
+
+      IF (collOutput) THEN
+        fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (fileID, file = path // folder // '/' //  fileName)
+
+        CALL writeHeader(self%numNodes, self%numCells, fileID)
+
+        CALL writeMesh(self, fileID)
+
+        CALL writeCollOutput(self, fileID)
+
+        CALL writeFooter(fileID)
+
+        CLOSE(fileID)
+
+        !Write collection file for time plotting
+        fileNameCollection = formatFileName(prefix, 'Collection_Collisions', 'pvd')
+        CALL writeCollection(fileID, t, fileName, filenameCollection)
+
+      END IF
+
+    END SUBROUTINE printCollVTK
+
+    SUBROUTINE printEMVTK(self, t)
+      USE moduleMesh
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in):: self
+      INTEGER, INTENT(in):: t
+      INTEGER:: fileID
+      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
+
+      fileID = 64
+
+      IF (emOutput) THEN
+        fileName = formatFileName(prefix, 'EMField', 'vtu', t)
+        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
+        OPEN (fileID, file = path // folder // '/' //  fileName)
+
+        CALL writeHeader(self%numNodes, self%numCells, fileID)
+
+        CALL writeMesh(self, fileID)
+
+        CALL writeEM(self, fileID)
+
+        CALL writeFooter(fileID)
+
+        CLOSE(fileID)
+
+        !Write collection file for time plotting
+        fileNameCollection = formatFileName(prefix, 'Collection_EMField', 'pvd')
+        CALL writeCollection(fileID, t, fileName, filenameCollection)
+
+      END IF
+
+    END SUBROUTINE printEMVTK
+
 END MODULE moduleMeshOutputVTK
diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90
index 813aeb3..def5ca4 100644
--- a/src/modules/solver/moduleSolver.f90
+++ b/src/modules/solver/moduleSolver.f90
@@ -528,9 +528,11 @@ MODULE moduleSolver
 
         CALL outputProbes(t)
         CALL mesh%printOutput(t)
-        CALL printOutputVTK(mesh, t)
+        CALL printOutputVTK(mesh, t) !TEMPORARY TO TEST VTK OUTPUT
         IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
+        IF (ASSOCIATED(meshForMCC)) CALL printCollVTK(meshForMCC,t) !TEMPORARY TO TEST VTK OUTPUT
         CALL mesh%printEM(t)
+        CALL printEMVTK(mesh, t) !TEMPORARY TO TEST VTK OUTPUT
         WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal
         WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
         IF (t == 0) THEN