diff --git a/csv2vtk_particles.py b/csv2vtk_particles.py
index f3d8966..4f94941 100755
--- a/csv2vtk_particles.py
+++ b/csv2vtk_particles.py
@@ -13,6 +13,7 @@
 import numpy as np
 import csv, os, shutil
 import configparser
+from zipfile import ZipFile
 try:
     import progressbar
     pbar = True 
@@ -28,13 +29,14 @@
 folder = 'result'
 i = 1
 
+
 while aux:
     try:
         os.mkdir(folder)
         aux = False
     except FileExistsError:
         print('Folder {} exists'.format(folder))
-        opt = input('Overwrite? [y/n]')
+        opt = input('Overwrite? [y/n] ')
         if opt == 'y' or 'opt' == 'Y':
             aux =False 
         else:
@@ -54,6 +56,7 @@
 N = int(config['global']['N'].split()[0])
 nimpre =  int(config['global']['nimpre'].split()[0])
 ntype = int(config['global']['Ntype'].split()[0])
+print_TC = int(config['global']['print_TC'].split()[0])
 quant = []
 rs = [] # raio sólido
 sigma = []
@@ -82,11 +85,6 @@
         rsol[j+k] = rs[i]
     k = sum(quant[0:i+1])
         
-
-# a = input("a")
-# for i in range(len(rsol)):
-#     print("{} {}".format(i,rsol[i]))
-
 # Ler os arquivos e colocá-los em vetores
 
 x = np.zeros(N)
@@ -103,9 +101,12 @@
 
 nID = np.linspace(1,N,N)
 
+zip_positions = ZipFile(via+'/'+folder+'/positions.zip','w')
+zip_velocities = ZipFile(via+'/'+folder+'/velocities.zip','w')
+
 print('Converting...')
 if pbar:
-    bar = progressbar.ProgressBar(max_value=nimpre)
+    bar = progressbar.ProgressBar(max_value=nimpre+1)
 for fnum in range(0,nimpre+1):
     with open('temp/position.csv.'+str(fnum),encoding='utf-8') as file_locus:
         csv_lector = csv.reader(file_locus,delimiter = ',')
@@ -114,8 +115,8 @@
             x[i] = linea[0]
             y[i] = linea[1]
             i = i+1
-    shutil.move(via+'/temp/position.csv.'+str(fnum),via+'/'+folder+'/position.csv.'+str(fnum)) 
-    
+    zip_positions.write(via+'/temp/position.csv.'+str(fnum), 'position.csv.'+str(fnum))
+    # shutil.move(via+'/temp/position.csv.'+str(fnum),via+'/'+folder+'/position.csv.'+str(fnum)) 
     with open('temp/velocity.csv.'+str(fnum),encoding='utf-8') as file_velocitas:
         csv_lector = csv.reader(file_velocitas,delimiter = ',')
         i = 0
@@ -124,7 +125,8 @@
             vy[i] = linea[1]
             
             i = i+1
-    shutil.move(via+'/temp/velocity.csv.'+str(fnum),via+'/'+folder+'/velocity.csv.'+str(fnum))
+    zip_velocities.write(via+'/temp/velocity.csv.'+str(fnum),'velocity.csv.'+str(fnum))
+    # shutil.move(via+'/temp/velocity.csv.'+str(fnum),via+'/'+folder+'/velocity.csv.'+str(fnum))
 
     fin = 0
     grupo = 0
@@ -144,8 +146,21 @@
 
     if pbar:
         bar.update(fnum)
+
+
+zip_positions.close()
+zip_velocities.close()
+
 shutil.rmtree('temp')
-        
+
+if print_TC == 1:
+    a = os.listdir('temp2')
+    zip_rFup = ZipFile(via+'/'+folder+'/rFuP.zip','w')
+    for f in a:
+        zip_rFup.write(via+'/temp2/'+f,f)
+        # shutil.move(via+'/temp2/'+f,via+'/'+folder + '/' + f)
+    zip_rFup.close()
+    shutil.rmtree('temp2')
         
 with open("settings.txt","a") as settingstxt:
     settingstxt.write("[out_files]\n")
@@ -168,6 +183,7 @@
     with open('.gitignore','a+') as file:
         file.write(folder+'\n')
 
+ 
     #with open('cell.csv.'+str(fnum),encoding='utf-8') as file_velocitas:
         #csv_lector = csv.reader(file_velocitas,delimiter = ',')
         #i = 0
diff --git a/lennard.f90 b/lennard.f90
index 2e1fa01..75d1864 100755
--- a/lennard.f90
+++ b/lennard.f90
@@ -33,7 +33,7 @@ subroutine MaxwellBoltzmann(malha,mesh,factor)
     end subroutine MaxwellBoltzmann
     
     !computa o potencial total 
-    subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nRfu) 
+    subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) 
         !Computa e imprime o potencial e forças em cada partícula, intermoleculares
 
         use linkedlist
@@ -44,19 +44,18 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
         
         type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade
         type(container), allocatable,dimension(:,:),intent(in) :: malha
-        real(dp), intent(in) :: t
-        real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut, fric_term !fric_term = força de ficção
+        real(dp), intent(in) :: t,r_cut
+        real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut, fric_term !fric_term = força de ficção
         real(dp) :: x1(2),v1(2),x2(2),p1(2), rs1, rs2, coss, sine, u   
         real(dp), intent(inout), allocatable, dimension(:) :: nRfu
         integer :: i,j,ct = 0, n1, n2 !,ptr, ptrn
-        integer, intent(in) :: mesh(:),domx(2),domy(2)
+        integer, intent(in) :: mesh(:),domx(2),domy(2), id
         real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2, dox(2), doy(2)
         type(list_t), pointer :: node, next_node
         type(data_ptr) :: ptr,ptrn
         integer ( kind = 4 ), intent(in) :: ids(8)
-        integer, intent(out) :: aux1 
+        
  
-        aux1 = 1
         !Lennard Jones
         dox = domx 
         doy = domy 
@@ -76,7 +75,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                 node => list_next(malha(i,j)%list)
                 if (associated(node)) then
                     ptr = transfer(list_get(node), ptr)
-                    ! print*, "L 253", ptr%p%n
                 end if
                 do while (associated(node))
                     ! print*, "Encontrada", ct, "celulas, id=", id
@@ -98,10 +96,15 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                     !calcular a força desta com todas as outras partículas
                     next_node => list_next(node) ! próxima partícula da célula
                     node => list_next(node) ! a ser computado com a particula selecionada
-             
-                    ! NA PRÓPRIA CELULA 
+
+                    if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                        nRfu(n1*6-5) = real(ptr%p%n,kind(0.d0))
+                    end if
+                    nRfu(n1*6-1:n1*6) = [p1(1), p1(2)] 
                     do while (associated(node))
+                        ! NA PRÓPRIA CELULA 
                         ptrn = transfer(list_get(node), ptrn) !outra particula selecionada
+                        ! print*, "L 108"
                         sigma_b = propriedade(ptrn%p%grupo)%sigma
                         epsil_b = propriedade(ptrn%p%grupo)%epsilon 
                         rs2 = propriedade(ptrn%p%grupo)%rs 
@@ -126,17 +129,21 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                         sine = (x1(2)-x2(2))/r 
                         r = r - rs1 - rs2 !raio
                         ! print*, "L 389 r", r, "id",id
+                        ! print*, "L 129"
                         if (r <= rcut) then
                             aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & 
                             [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]  
                             u = 4*epsil*((sigma/r)**12 - (sigma/r)**6)  
-                            nRfu(n1:n1+5) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u, p1(1), p1(2)] &
-                            + nRfu(n1:n1+5)
-                            nRfu(n2:n2+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
-                            + nRfu(n2:n2+3)
-
+                            nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                            + nRfu(n1*6-4:n1*6-2)
+                            nRfu(n2*6-4:n2*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] & 
+                            + nRfu(n2*6-4:n2*6-2)
+                            if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                                nRfu(n2*6-5) = real(n2,kind(0.d0))
+                            end if
+                            ! print*, "nrfu2", nRfu(n1*6-5:n1*6)
+                            !read(*,*)
                         end if
-                        aux1 = aux1 + 5
                         node => list_next(node) ! próxima partícula da célula
                     end do
 
@@ -165,6 +172,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                     epsil = epsil_a
                                 end if 
                                 x2 = ptrn%p%x
+                                n2 = ptrn%p%n
                                 r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2)
                                 coss = (x1(1)-x2(1))/r 
                                 sine = (x1(2)-x2(2))/r 
@@ -173,10 +181,13 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                 if (r <= rcut) then
                                     aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & 
                                         [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
-                                    nRfu(n1:n1+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
-                                        + nRfu(n1:n1+3)
-                                    nRfu(n2:n2+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
-                                        + nRfu(n2:n2+3)
+                                    nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                                        + nRfu(n1*6-4:n1*6-2)
+                                    nRfu(n2*6-4:n2*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] & 
+                                        + nRfu(n2*6-4:n2*6-2)
+                                    if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                                        nRfu(n2*6-5) = real(n2,kind(0.d0))
+                                    end if
                                 end if
                                 node => list_next(node) ! próxima partícula da célula
                             end do
@@ -205,6 +216,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                     end if 
 
                                     x2 = ptrn%p%x
+                                    n2 = ptrn%p%n
                                     r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2)
                                     coss = (x1(1)-x2(1))/r 
                                     sine = (x1(2)-x2(2))/r 
@@ -214,10 +226,13 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                         aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & 
                                             [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
 
-                                        nRfu(n1:n1+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
-                                        + nRfu(n1:n1+3)
-                                        nRfu(n2:n2+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
-                                        + nRfu(n2:n2+3)
+                                        nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                                        + nRfu(n1*6-4:n1*6-2)
+                                        nRfu(n2*6-4:n2*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] & 
+                                        + nRfu(n2*6-4:n2*6-2)
+                                        if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                                            nRfu(n2*6-5) = real(n2,kind(0.d0))
+                                        end if
                                     end if
                                     node => list_next(node) ! próxima partícula da célula
                                 end do                            
@@ -225,7 +240,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                         else
                              !interagirá com a próxima linha e coluna, e na diagonal
                             node => list_next(malha(i,j+1)%list) 
-                            
                             do while (associated(node))
                                 ptrn = transfer(list_get(node), ptrn) !outra particula selecionada
                                 sigma_b = propriedade(ptrn%p%grupo)%sigma
@@ -248,6 +262,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                 end if 
 
                                 x2 = ptrn%p%x
+                                n2 = ptrn%p%n
                                 r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2)
                                 coss = (x1(1)-x2(1))/r 
                                 sine = (x1(2)-x2(2))/r 
@@ -257,10 +272,13 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                     aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & 
                                         [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] 
 
-                                    nRfu(n1:n1+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
-                                    + nRfu(n1:n1+3)
-                                    nRfu(n2:n2+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
-                                    + nRfu(n2:n2+3)
+                                    nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                                    + nRfu(n1*6-4:n1*6-2)
+                                    nRfu(n2*6-4:n2*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] & 
+                                    + nRfu(n2*6-4:n2*6-2)
+                                    if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                                        nRfu(n2*6-5) = real(n2,kind(0.d0))
+                                    end if
                                 end if
                                 node => list_next(node) ! próxima partícula da célula                                
                             end do
@@ -287,6 +305,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                     epsil = epsil_a
                                 end if 
                                 x2 = ptrn%p%x
+                                n2 = ptrn%p%n
                                 r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2)
                                 coss = (x1(1)-x2(1))/r 
                                 sine = (x1(2)-x2(2))/r 
@@ -295,10 +314,13 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                 if (r <= rcut) then
                                     aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & 
                                         [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] 
-                                    nRfu(n1:n1+5) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
-                                        + nRfu(n1:n1+3)
-                                    nRfu(n2:n2+5) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
-                                        + nRfu(n2:n2+3)
+                                    nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                                        + nRfu(n1*6-4:n1*6-2)
+                                    nRfu(n2*6-4:n2*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] & 
+                                        + nRfu(n2*6-4:n2*6-2)
+                                    if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                                        nRfu(n2*6-5) = real(n2,kind(0.d0))
+                                    end if
                                 end if
                                 node => list_next(node) ! próxima partícula da célula                                
                             end do
@@ -326,6 +348,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                 end if 
 
                                 x2 = ptrn%p%x
+                                n2 = ptrn%p%n
                                 r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2)
                                 coss = (x1(1)-x2(1))/r 
                                 sine = (x1(2)-x2(2))/r 
@@ -334,10 +357,13 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                 if (r <= rcut) then
                                     aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & 
                                         [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] 
-                                    nRfu(n1:n1+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
-                                        + nRfu(n1:n1+3)
-                                    nRfu(n2:n2+5) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
-                                        + nRfu(n2:n2+3)
+                                    nRfu(n1*6-4:n1*6-2) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                                        + nRfu(n1*6-4:n1*6-2)
+                                    nRfu(n2*6-4:n2*6-2) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
+                                        + nRfu(n2*6-4:n2*6-2)
+                                    if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                                        nRfu(n2*6-5) = real(n2,kind(0.d0))
+                                    end if
                                 end if
                                 node => list_next(node) ! próxima partícula da célula                                
                             end do
@@ -366,6 +392,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                     end if 
 
                                     x2 = ptrn%p%x
+                                    n2 = ptrn%p%n
                                     r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2)
                                     coss = (x1(1)-x2(1))/r 
                                     sine = (x1(2)-x2(2))/r 
@@ -374,10 +401,13 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                                     if (r <= rcut) then
                                         aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & 
                                             [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] 
-                                        nRfu(n1:n1+5) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
-                                            + nRfu(n1:n1+3)
-                                        nRfu(n2:n2+5) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] & 
-                                            + nRfu(n2:n2+3)
+                                        nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                                            + nRfu(n1*6-4:n1*6-2)
+                                        nRfu(n2*6-4:n2*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] & 
+                                            + nRfu(n2*6-4:n2*6-2)
+                                        if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                                            nRfu(n2*6-5) = real(n2,kind(0.d0))
+                                        end if
                                     end if
                                     node => list_next(node) ! próxima partícula da célula
                                 end do                            
@@ -408,6 +438,7 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                             end if 
 
                             x2 = ptrn%p%x
+                            n2 = ptrn%p%n
                             r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2)
                             coss = (x1(1)-x2(1))/r 
                             sine = (x1(2)-x2(2))/r 
@@ -416,10 +447,13 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
                             if (r <= rcut) then
                                 aux2 = -(1/r**2)*(sigma/r)**6*(1-2*(sigma/r)**6)*24*epsil* & 
                                     [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] 
-                                nRfu(n1:n1+3) = [real(ptr%p%n,kind(0.d0)), aux2(2)*r*coss, aux2(1)*r*sine, u] &
-                                    + nRfu(n1:n1+3)
-                                nRfu(n2:n2+3) = [real(ptr%p%n,kind(0.d0)), -aux2(2)*r*coss, -aux2(1)*r*sine, u] & 
-                                    + nRfu(n2:n2+3)
+                                nRfu(n1*6-4:n1*6-2) = [aux2(2)*r*coss, aux2(1)*r*sine, u] &
+                                    + nRfu(n1*6-4:n1*6-2)
+                                nRfu(n2*6-4:n2*6-2) = [-aux2(2)*r*coss, -aux2(1)*r*sine, u] & 
+                                    + nRfu(n2*6-4:n2*6-2)
+                                if (i /= doy(1) .and. i /= doy(2) .and. j /= dox(1) .and. j /= dox(2)) then
+                                    nRfu(n2*6-5) = real(n2,kind(0.d0))
+                                end if
                             end if
                             node => list_next(node) ! próxima partícula da célula                                
                         end do
@@ -429,6 +463,8 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, aux1, nR
             end do
         end do
         aux1 = aux1 - 1
+        ! print*, "ID, AUX1", id, aux1
+        ! print*, "FIM"
     end subroutine comp_pot    
  
     function comp_Kglobal(malha,domx,domy,propriedade,np,id,t) result(K)
@@ -651,7 +687,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t)
         real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut, fric_term !fric_term = força de ficção
         real(dp) :: x1(2),v1(2),x2(2),v2(2), rs1, rs2, coss, sine   
         integer :: i,j,ct = 0 !,ptr, ptrn
-        integer, intent(in) :: mesh(:),domx(2),domy(2)
+        integer, intent(in) :: mesh(:),domx(2),domy(2), id
         real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2, dox(2), doy(2)
         type(list_t), pointer :: node, next_node
         type(data_ptr) :: ptr,ptrn
@@ -675,10 +711,10 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t)
         do i = doy(1),doy(2) ! i é linha
             do j = dox(1),dox(2)
                 node => list_next(malha(i,j)%list)
-                if (associated(node)) then
-                    ptr = transfer(list_get(node), ptr)
-                    ! print*, "L 253", ptr%p%n
-                end if
+                ! if (associated(node)) then
+                !     ptr = transfer(list_get(node), ptr)
+                !     ! print*, "L 253", ptr%p%n
+                ! end if
                 do while (associated(node))
                     ! print*, "Encontrada", ct, "celulas, id=", id
                     ptr = transfer(list_get(node), ptr) !particula selecionada
@@ -1067,7 +1103,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t)
     end subroutine comp_F
     
     ! atualiza posições
-    subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,domy,wall,mic,id, np)
+    subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,domy,wall,id, np)
         use linkedlist
         use mod1
         use data
@@ -1078,7 +1114,7 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,
         character(1) :: north, south, east, west
         type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade
         type(container), allocatable,dimension(:,:),intent(in) :: malha
-        integer, intent(inout) :: mic(:)
+        ! integer, intent(inout) :: mic(:)
         integer :: i,j,k, cell(2), status(MPI_STATUS_SIZE), count, destD,dx1,dx2,dy1,dy2
         integer, intent(in) :: N,mesh(2),domx(2),domy(2)
         type(list_t), pointer :: node, previous_node
@@ -1693,7 +1729,7 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,
                 previous_node => malha(i,j)%list
                 node => list_next(malha(i,j)%list)
     
-                if (associated(node)) ptr = transfer(list_get(node), ptr)
+                ! if (associated(node)) ptr = transfer(list_get(node), ptr)
                 ! if (associated(node))  print*, "transferencia diagonal", i,j, "ID", id
                 do while(associated(node))
                     ptr = transfer(list_get(node), ptr)
@@ -2531,7 +2567,7 @@ program main
     real(dp), dimension(:), allocatable :: icell,jcell, nxv, nxv_send, nRfu, nRfu_send !dimensões das celulas e vetor de resultado pra imprimir
     real(dp) :: t=0,t_fim,dt,printstep, sigma, epsil, rcut,aux2,start = 0,finish,Td,kb = 1.38064852E-23,vd(2)
     real(dp) :: GField(2), temp_Td(3), dimX, dimY, dx_max, Td_hot, Td_cold
-    integer,allocatable :: interv(:), interv_Td(:), grupo(:), rcounts(:), displs(:)) !, mic(:,:), mic_rcv(:) ! mic são vetores para contar quantas vezes atravessou a borda periodica
+    integer,allocatable :: interv(:), interv_Td(:), grupo(:), rcounts(:), displs(:) !, mic(:,:), mic_rcv(:) ! mic são vetores para contar quantas vezes atravessou a borda periodica
     type(container), allocatable,dimension(:,:) :: malha
     type(prop_grupo),allocatable,dimension(:) :: propriedade
     type(list_t), pointer :: node
@@ -2728,7 +2764,7 @@ program main
             do j = 1,quant
                 call CFG_get(my_cfg, particle//"%v", v(cont,:))    
                 grupo(cont) = i+1
-                read(20,*) x(cont,:) !,x(cont,2),x(cont,3)
+                read(20,*) x(cont,:) 
                 cont = cont+1
             end do
         end if
@@ -2921,10 +2957,10 @@ program main
     call clean_mesh(malha, mesh, domx, domy,id,.false.)
     ! print*, "bbb", id
   
-    if (id == 0) then
-        print*,'Press Return to continue'
-       read(*,*)
-    end if
+    ! if (id == 0) then
+    !     print*,'Press Return to continue'
+    !    read(*,*)
+    ! end if
     call MPI_barrier(MPI_COMM_WORLD, ierr)
     
     ! id = 0 é o processo raiz 
@@ -2941,6 +2977,11 @@ program main
         i = interv(j)
     end if
 
+    if (print_TC == 1) then
+        allocate(nRfu(N*6),nRfu_send(N*6),rFUp(N,6))
+        call system('mkdir temp2')
+    end if
+
     do while (t_fim > t)
         ! print*, "L 1990"
         call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t)  !altera Força
@@ -3024,33 +3065,28 @@ program main
             end if
             ! if (id == 0) read(*,*)
             ! call MPI_barrier(MPI_COMM_WORLD, ierr) 
-            j = j+1
+            
             deallocate(nxv,nxv_send)
 
             if (print_TC == 1) then 
-                allocate(nRfu(N*6),nRfu_send(N*6),rFUp(N,5))
-                nRfu = nan
-                ! call MPI_barrier(MPI_COMM_WORLD, ierr) 
-                ! print*, 'MANDANDO PRINTAR', id, i
-                ! CONTINUAR DAQUI
-                call comp_pot(mesh,malha,propriedade,rcut,domx,domy, ids, id, t, aux1, nRfu)
-                ! ! ! print*, "L nxv_send", nxv_send  
+                nRfu = 0.0_dp
+
+                call comp_pot(mesh,malha,propriedade,rcut,domx,domy, ids, id, t, nRfu) 
+
                 !! print*, "LINKED 2 VEC, id", id, "aux1", aux1
                 if (np > 1) then ! paralelo
                     aux3 = 1
                     do ii = 1,N
-                        if (.not. isnan(nRfu(ii*6-5))) then
+                        if (nRfu(ii*6-5) /= 0) then
                             nRfu_send(aux3:aux3+5) = nRfu(ii*6-5:ii*6)
                             aux3 = aux3 + 6
                         end if
                     end do
+                    aux1 = aux3 -1
 
                     ! call MPI_GATHER(sbuf, scount, MPI_integer, rbuf, rcount, MPI_integer, root, MPI_COMM_WORLD, ierr)
                     call MPI_GATHER(aux1,     1,    MPI_integer, rcounts, 1,   MPI_integer,  0,   MPI_COMM_WORLD, ierr)
     
-                    ! print*, "Aux1 =", aux1, "id =", id
-                    ! if (id == 0) print*, "rcounts =", rcounts
-                    ! !print*, 'L IMPRE >', id
                     displs = 0
     
                     if (id == 0) then
@@ -3059,37 +3095,49 @@ program main
                         end do
                         ! print*, "displs", displs
                     end if
-                    ! print*,'E_tot0 = ',(comp_pot(mesh,malha,propriedade,rcut) +comp_K(malha,mesh,propriedade))
+                    ! if (id == 0) then 
+                    !     print*, "L 3099"
+                    !     read(*,*)
+                    ! end if
+                    ! call MPI_barrier(MPI_COMM_WORLD, ierr) 
+
                     ! MPI_GATHERV    (sbuf,   scount,  stype, rbuf, rcounts,  displs,   rtype,  root,  comm,  ierr)
                     call MPI_GATHERV(nRfu_send, aux1, MPI_DOUBLE_PRECISION, nRfu, rcounts, &
                      displs, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
                 end if 
-                
+                ! print*, "AA"          
+
                 if (id == 0) then
                     do ii = 0, N-1
-                        rFUp(int(nRfu(ii*6+1)),:) = [nRfu(ii*6+2),nRfu(ii*6+3),nRfu(ii*6+4),nRfu(ii*6+5),nRfu(ii*6+6)]
+                        ! print*, [nRfu(ii*6+2),nRfu(ii*6+3),nRfu(ii*6+4),nRfu(ii*6+5),nRfu(ii*6+6)]
+                        rFUp(int(nRfu(ii*6+1)),:) = [nRfu(ii*6+1),nRfu(ii*6+2),nRfu(ii*6+3),nRfu(ii*6+4),nRfu(ii*6+5),nRfu(ii*6+6)]
                     end do
-                    
-                    call vec2csv(rFUp,N,5,'rF_u_P',j,t,nimpre,start)
+                    ! print*, "AAA"
+                    call vec2csv(rFUp,N,6,'rF_u_P',j,t,nimpre,start)
                     
                 end if 
-
-                deallocate(nRfu,nRfu_send,rFUp)
+                ! print*, "B"
+                ! read(*,*)
+                ! deallocate(nRfu,nRfu_send,rFUp)
+                ! print*, "C"
             end if 
-
             allocate(LT%lstrdb_N(NMPT*6),LT%lstrdb_S(NMPT*6),LT%lstrdb_E(NMPT*6),LT%lstrdb_W(NMPT*6),LT%lstrdb_D(NMPT*6), &
-                 LT%lstrint_N(NMPT*4),LT%lstrint_S(NMPT*4),LT%lstrint_E(NMPT*4),LT%lstrint_W(NMPT*4), LT%lstrint_D(NMPT*4))
+            LT%lstrint_N(NMPT*4),LT%lstrint_S(NMPT*4),LT%lstrint_E(NMPT*4),LT%lstrint_W(NMPT*4), LT%lstrint_D(NMPT*4))
+            j = j+1
         end if
         i = i+1
         ! print*, "L 1754 >", id
         ! if (id  == 0) read(*,*)    
         ! call MPI_barrier(MPI_COMM_WORLD, ierr)
+        ! print*, t
     end do
 
-    ! print*, 'L 2112', id 
     call clean_mesh(malha, mesh, domx, domy,id,.true.)
     deallocate(LT%lstrdb_N, LT%lstrdb_S, LT%lstrdb_E, LT%lstrdb_W, LT%lstrdb_D, &
             LT%lstrint_N, LT%lstrint_S, LT%lstrint_E, LT%lstrint_W, LT%lstrint_D)
+    ! print*, 'L 3104'
+    ! if (print_TC == 1) deallocate(nRfu,nRfu_send,rFUp)
+    
     if (id == 0)  then
         call cpu_time(finish)
         print '("Time = ",f10.3," seconds.")',(finish-start)
diff --git a/makefile b/makefile
index e61791e..6428a16 100644
--- a/makefile
+++ b/makefile
@@ -2,10 +2,10 @@
 # Defining variables
 objects = m_config.o linkedlist.o mod1.o saida.o lennard.o data.o mod0.o matprint.o randnormal.o 
 f90comp = mpiifort
-switch = -O3
+switch = -O3 -assume protect_parens,minus0 -prec-div -prec-sqrt 
 # Makefile
 execname: $(objects)
-	$(f90comp) -o lennard.out $(switch) $(objects)
+	$(f90comp) -o lennard $(switch) $(objects)
 randnormal.mod: randnormal.o randnormal.f90
 	$(f90comp) -c $(switch) randnormal.f90
 mod1.mod: mod1.o mod1.f90
diff --git a/part_gen_pequeno.py b/part_gen_pequeno.py
index 8a06492..6ec27c2 100644
--- a/part_gen_pequeno.py
+++ b/part_gen_pequeno.py
@@ -86,6 +86,8 @@
         p1[cont,:] = [i*spac+refe[0],j*spac+refe[1]]
         cont = cont+1   
 
+# np.random.shuffle(p1)        
+
 with open(arquivo1,'w') as file:
     for i in range(N):
         file.write('{},{}\n'.format(p1[i,0],p1[i,1]))
diff --git a/pos_e_vel_inicial_usando_sim_anterior.py b/pos_e_vel_inicial_usando_sim_anterior.py
index d2fe296..44fe718 100644
--- a/pos_e_vel_inicial_usando_sim_anterior.py
+++ b/pos_e_vel_inicial_usando_sim_anterior.py
@@ -7,6 +7,7 @@
 from glob import glob
 import os, shutil
 import configparser
+from zipfile import ZipFile
 
 dirname = os.getcwd() #os.path.dirname(os.path.abspath(__file__))
 dirlist = glob(dirname + "/*/")
@@ -16,8 +17,11 @@
 a = int(input("Enter the number of the folder\n"))
 res_dir = dirlist[a]
 
+zip_positions = ZipFile(res_dir+'/positions.zip','r')
+zip_velocities = ZipFile(res_dir+'/velocities.zip','r')
+
 config = configparser.ConfigParser()
-list_files = os.listdir(res_dir)
+len_list_files =  len(zip_positions.namelist()+zip_velocities.namelist())
 
 if os.path.isfile(res_dir + 'settings.txt'):
     config.read(res_dir + 'settings.txt')
@@ -26,7 +30,7 @@
     print('The still not converted to vtk. Probable incomplete simulation.\n'+\
           'settings.ini will be used.\n')
     config.read(dirname + '/settings.ini') 
-    nimpre = len(list_files)/2    
+    nimpre = len(len_list_files)/2    
 
     
 N = int(config['global']['N'].split()[0])
@@ -39,27 +43,42 @@
     quant.append(int(config['par_'+str(i)]['quantidade'].split()[0]))
     x_files.append(config['par_'+str(i)]['x'].split()[0])
     x_files[-1]  = x_files[-1][1:len(x_files[-1])-1] 
-    v_files.append(config['par_'+str(i)]['v_file'].split()[0])
-    v_files[-1]  = v_files[-1][1:len(v_files[-1])-1]
+    try:
+        v_files.append(config['par_'+str(i)]['v_file'].split()[0])
+        v_files[-1]  = v_files[-1][1:len(v_files[-1])-1]
+    except:
+        print("no velocity file used")
+        v_files.append('v_file_'+str(i))
+        pass 
     
 step = input('Extract step no. (max {}) '.format(nimpre-1))
 
 positions = []
-with open(res_dir + 'position.csv.'+step) as file:
+with zip_positions.open('position.csv.'+step) as file:
     for line in file:
         positions.append(line)
 velocities = []
-with open(res_dir + 'velocity.csv.'+step) as file:
+with zip_velocities.open('velocity.csv.'+step) as file:
     for line in file:
         velocities.append(line)
 
+
+# positions = []
+# with open(res_dir + 'position.csv.'+step) as file:
+#     for line in file:
+#         positions.append(line)
+# velocities = []
+# with open(res_dir + 'velocity.csv.'+step) as file:
+#     for line in file:
+#         velocities.append(line)
+
 # verificar se não é desejavel fazer backup de arquivos de posição e velocidade
 a == 'n'
 i = 0
 for f in x_files + v_files:
     if os.path.isfile(f) and a == 'n':
-        a = input('Old positions files exist. Backup? y/n')
-    if a != 'n' or a != 'N':
+        a = input('Old positions files exist. Backup? y/n ')
+    if a != 'n' or a != 'N': 
         try:
             if i == 0:
                 bkp = 'bkp'
@@ -80,14 +99,14 @@
             
 k = 0            
 for i in range(ntype):
-    with open(x_files[i],'w') as file:
+    with open(x_files[i],'wb') as file:
         for j in range(quant[i]):
             file.write(positions[k])
             k += 1
     
 k = 0        
 for i in range(ntype):
-    with open(v_files[i],'w') as file:
+    with open(v_files[i],'wb') as file:
         for j in range(quant[i]):
             file.write(velocities[k])
             k += 1        
diff --git a/saida.f90 b/saida.f90
index 61ec96f..4cd75e0 100644
--- a/saida.f90
+++ b/saida.f90
@@ -12,9 +12,15 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start)
         character(4) :: extensao = '.csv', passo
         real(dp) :: time
         real(dp), save :: timep, etc, dtimepp
+        character(LEN=*),parameter :: fmt5 = '(f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )'
+        character(LEN=*),parameter :: fmt6 = '(f5.0, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )'
 
         write(passo,'(i0)') step
-        open(10,file='temp/'//prop//extensao//'.'//trim(passo),status="replace")
+        if (d <= 3) then
+            open(10,file='temp/'//prop//extensao//'.'//trim(passo),status="replace")
+        else
+            open(10,file='temp2/'//prop//extensao//'.'//trim(passo),status="replace")
+        end if
         if (d == 1) then
             do i = 1,n
                 write(10,*) v(i,1)
@@ -33,9 +39,14 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start)
             end do
         else if (d == 5) then 
             do i = 1,n
-                write(10,*) v(i,1),',',v(i,2),',',v(i,3),',',v(i,4),',',v(i,5)
+                write(10,fmt5) v(i,1),v(i,2),v(i,3),v(i,4),v(i,5)
+            end do          
+        else if (d == 6) then 
+            do i = 1,n
+                write(10,fmt6) v(i,1),v(i,2),v(i,3),v(i,4),v(i,5),v(i,6)
             end do          
         end if
+        
         close(10)
 
         if (prop == "position") then 
diff --git a/settings.ini b/settings.ini
index e6ab3ff..ed3335a 100755
--- a/settings.ini
+++ b/settings.ini
@@ -1,9 +1,9 @@
 # PARAMETROS GLOBAIS
 [global]
-    nimpre  = 5000 # quntidade de saidas 
-    N = 1089 # número de partículas 
+    nimpre  = 50 # quntidade de saidas 
+    N = 108 #9 # número de partículas 
     Ntype = 1 # número de tipos de partícula
-    t_fim = 50  # tempo de execução 
+    t_fim = 5  # tempo de execução 
     nimpre_init = 0 # começa a simular a partir no tempo referente ao snapshot (nimpre). 
     dt = 0.0001 # passo de tempo
     dimX = 410 # Dimensões da região de cálculo
@@ -20,7 +20,7 @@
     vd = 2 2 # velocidade distribuida com Maxwell Boltzmann
     NMPT = 100 # Número máximo de partículas que poderá mudar de process. Se não souber estimar, usar -1 ou o valor de N.
     GField = 0 0 # Campo de força gravitacional que afeta todas as partículas
-    print_TC = 0 # imprimir dados para calcular coeficiente de transporte
+    print_TC = 1 # imprimir dados para calcular coeficiente de transporte
 # PARTICLE PROPERTIES
 # v é velocidade global, entrar dois número correspondentes ao vetor vx vy que serão aplicadas em todas as partículas
 # se v_file é dada, então ela será usada para dar a velocidade às partículas. Deve ser arquivo .csv. Deixar "%" se não usar. 
@@ -32,7 +32,7 @@
     nome = 'g1'
     sigma = 2 # Lennard Jones
     epsilon = 1 # Lennard Jones
-    quantidade = 1089
+    quantidade = 108 #9
     x_lockdelay = 0 # só vai poder mudar de posição a partir de t = x_lockdelay
     rs = 0 # raio sólido. Posição de partículas na superfície de um sólido de raio rs
     fric_term = 0 # fricção artifical				
diff --git a/teste.f90 b/teste.f90
index 4d90842..524c51a 100644
--- a/teste.f90
+++ b/teste.f90
@@ -1,11 +1,23 @@
 program teste
-    use mod1
+    ! use mod1
    
-    call system('mkdir temp')
-    open(10,file='temp/test.txt',status="replace")
-    write(10,*) 'lalalala'
-    print*,'a'
+    ! call system('mkdir temp')
+    ! open(10,file='temp/test.txt',status="replace")
+    ! write(10,*) 'lalalala'
+    ! print*,'a'
     
+    implicit none
+    real :: a,b,c,d
+    character(LEN=*), parameter :: fmt = '(f16.0, ", ",f32.16, ", ",f32.16, ", ",f32.16, ", ",f32.16 )'
+    character(170) :: out
+
+    a = 12.0
+    b = 1.2
+    c = sqrt(2.0)
+    d = -sqrt(3.0)
+    write(out,fmt) a,b,c,d,d
+    out = trim(out)
+    write(*,*) out