adimensional e gospindo dados print_TC

anotações.txt

@@ -0,0 +1,37 @@
+Equações adimensionais para MD
+m' = adimensional
+m* = valor dimensional de referência
+[] vetor
+
+m' = m/m* (massa)
+sigma' = sigma/sigma* (distância)
+epsilon' = epsion/epsilon* (energia = força*[distância]) No caso do LJ, epsilon*[distância]/distância^2 = força
+[xi'] = [xi]/sigma* (distância)
+Dimensões:
+
+comprimento    sigma (LJ)
+energia        epsilon (LJ)
+massa          m
+Velocidade     sqrt(epsilon/m)
+força          epsilon/sigma 
+pressão        epsilon/sigma^3
+temperatura    epsilon*k_b
+
+Adimensionais
+
+[rij] = [rij]/sigma*
+T' = T*k_b/epsilon
+E' = E/epsilon*
+V' = V/epsilon* (energia potencial) 
+t' = t/alpha com alpha = sqrt( (m^2*sigma^2/epsilon*) )
+
+O que dá pra fazer aqui para admensionalizar. 
+Escolher uma das particulas como referência, esta terá valores de
+sigma, m, epsilon = 1 
+A velocidade das partículas determinada por T ou E será o que fará as particulas
+se aglutinarem mais ou menos. 
+
+
+
+Ek =  T'* epsilon 
+T' = Ek/epsilon (numero grande) = Ek'

csv2vtk_particles.py

@@ -46,7 +46,7 @@
                 aux = True
                 i = i+1       
 
-print('The results will be saved at {}'.format(folder))
+print('\nThe results will be saved at {}\n'.format(folder))
 
 #Ler a quantidade de arquivos de settings.
 

lennard.f90

@@ -18,7 +18,6 @@ subroutine MaxwellBoltzmann(malha,mesh,factor)
 
         do i = 2,mesh(2)+1
             do j = 2,mesh(1)+1
-               ! print *, 'posição', i, ',', j
                 node => list_next(malha(i,j)%list)
                 do while (associated(node))
                     ptr = transfer(list_get(node), ptr)
@@ -535,19 +534,14 @@ subroutine corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,tt)
         integer ( kind = 4 ) :: np, id
         real(dp), intent(in) :: tt
         !calcula temperatura atual
-        T = (2/(3*kb))*comp_Kglobal(malha,domx,domy,propriedade,np,id,tt)
+        ! T = (2/(Nf*kb))*Ekin com Nf = número de graus de liberadade 
+        T = (2/(2*kb))*comp_Kglobal(malha,domx,domy,propriedade,np,id,tt)
 
-        ! print*, "T =", T, "Td =", Td
-        ! if (id  == 0) read(*,*)    
-        ! call MPI_barrier(MPI_COMM_WORLD, ierr)
         beta = sqrt(Td/T) ! aqui o T é o Beta
         do i = domy(1),domy(2)
             do j = domx(1),domx(2)
                ! print *, 'posição', i, ',', j
                 node => list_next(malha(i,j)%list)
-                !  print*, i,j
-                !  print*, 'is associated?', associated(node)
-                !  print*, associated(node)
                 do while (associated(node))
                     ptr = transfer(list_get(node), ptr)
                     if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then
@@ -628,9 +622,9 @@ subroutine corr_K(malha,domx,domy,cold_cells,hot_cells,Td_hot,Td_cold,propriedad
         integer ( kind = 4 ) :: np, id
         real(dp), intent(in) :: tt
         !calcula temperatura atual
-
-        T_hot = (2/(3*kb))*comp_K(malha,domx,domy,hot_cells,propriedade,np,id,tt)
-        T_cold = (2/(3*kb))*comp_K(malha,domx,domy,cold_cells,propriedade,np,id,tt)
+        ! T = (2/(Nf*kb))*Ekin com Nf = número de graus de liberadade 
+        T_hot =  comp_K(malha,domx,domy,hot_cells,propriedade,np,id,tt)
+        T_cold = comp_K(malha,domx,domy,cold_cells,propriedade,np,id,tt)
         ! print*, "T =", T, "Td =", Td
         ! if (id  == 0) read(*,*)    
         ! call MPI_barrier(MPI_COMM_WORLD, ierr)
@@ -2803,7 +2797,7 @@ program main
     i = 0
 
     if (Td == 0) then
-        Td = (2/(3*N*kb))*sum(0.5*propriedade(ptr%p%grupo)%m*(v(:,1)**2+v(:,2)**2))
+        Td = (2/(2*N*kb))*sum(0.5*propriedade(ptr%p%grupo)%m*(v(:,1)**2+v(:,2)**2))
     else if (Td < 0 .and. (Td_hot*Td_hot) <= 0) then
         interv_Td = -1
     end if
@@ -2841,7 +2835,8 @@ program main
     deallocate(grupo)
 
     ! adiciona às particulas velocidade inicial de acordo com distribuição maxwell boltzmann
-    if (vd(1) /= 0 .and. vd(2) /= 0 .and. vd(1) /= 0) call MaxwellBoltzmann(malha,mesh,sqrt(vd))
+    if (vd(1) /= 0 .and. vd(2) /= 0) call MaxwellBoltzmann(malha,mesh, vd)
+    if (vd(1) == 0 .and. vd(2) == 0 .and. Td > 0) call MaxwellBoltzmann(malha,mesh,sqrt([Td,Td]))
 
     ! adicionamos uma celula para corrigir o fato de estarmos usando fantasmas
     cold_cells = cold_cells + 1

pos_e_vel_inicial_usando_sim_anterior.py

@@ -73,7 +73,7 @@
 #         velocities.append(line)
 
 # verificar se não é desejavel fazer backup de arquivos de posição e velocidade
-a == 'n'
+a = 'n'
 i = 0
 for f in x_files + v_files:
     if os.path.isfile(f) and a == 'n':

settings.ini

@@ -1,38 +1,40 @@
 # PARAMETROS GLOBAIS
+# Convenção: O primeiro grupo de partículas (part_0) deve ser o de referência (m = epsilon = sigma = 1)
+# Os outros parâmetros são adimensionais. Checar anotações.txt T' = T*k_b/epsilon
 [global]
     nimpre  = 50 # quntidade de saidas 
-    N = 108 #9 # número de partículas 
+    N = 1089 # número de partículas 
     Ntype = 1 # número de tipos de partícula
-    t_fim = 5  # tempo de execução 
+    t_fim = 10  # 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
     dimY = 410 #
     mesh = 82 82 # malha(elementos)
     rcut = 5 # *sigma
     wall = 'eeee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste
-    Td = -1 # 2.897189213660259e+24 #temperatura desejada global constante (-1 = não aplicar) kb = 1.38064852E-23, Td = (2/(3*kb))*E
+    Td = 3 # 2.897189213660259e+24 #temperatura desejada global constante (-1 = não aplicar) kb = 1.38064852E-23, Td = (2/(3*kb))*E/epsilon(literatura)
     temp_Td = 0 1000 100 # Tempo de velocity scaling [inicial, final, iterações_por_scaling] 
     cold_cells = 1 5 1 52 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
     hot_cells = 145 150 1 52 # [cellx_ini, cellx_fin, celly_ini, celly_fin]
     Td_hot = -1 #100e+23  # temperatura nas cold cells constante. Td_hot = -1 desliga termostato
     Td_cold = -1 #10e+23 # temperatura nas hot cells constante. Td_cold = -1 desliga termostato
-    vd = 2 2 # velocidade distribuida com Maxwell Boltzmann
+    vd = 0 0 # velocidade distribuida com Maxwell Boltzmann = sqrt(Td') 
     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 = 1 # imprimir dados para calcular coeficiente de transporte
+    print_TC = 0 # 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. 
 [par_0]
     x = 'bin_peq1.csv' 
     v = 0 0 # velocidade global
-    v_file = '%vgrupo0.csv' # velocidade de cada partícula
+    v_file = '%v_file_0.csv' # velocidade de cada partícula
     m = 1
     nome = 'g1'
-    sigma = 2 # Lennard Jones
+    sigma = 1 # Lennard Jones
     epsilon = 1 # Lennard Jones
-    quantidade = 108 #9
+    quantidade = 1089
     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