implementado rotação - Em testes

debug_lennard.sh

@@ -18,9 +18,9 @@ else
         mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -debug -O0 -shared-intel -nocheck -fno-omit-frame-pointer -fasynchronous-unwind-tables -fexceptions -o lennard.out
         echo "output: lennard.out"
     else
-        echo "compiling with intel mpiifort without options"
-        mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -debug -o lennard.out
-        echo "output: lennard"
+        echo "compiling with intel mpiifort without optimization"
+        mpiifort mod1.f90 linkedlist.f90 mod0.f90 saida.f90 data.f90 matprint.f90 m_config.f90 randnormal.f90 lennard.f90 -O0 -g -traceback -check all -o lennard.out
+        echo "output: lennard.out"
     fi
 
 fi

mod0.f90

@@ -15,6 +15,7 @@ module mod0
             real(dp)  :: x_lockdelay    ! só vai poder mudar de posição a partir de t = x_lockdelay
             real(dp) :: fric_term
             logical :: ismolecule
+            integer :: quant
         end type prop_grupo
 
             ! LSTR  lista de transferência pro MPI

pack.sh

@@ -3,7 +3,7 @@ read -e -p "Enter particle generator file: " VAR1
 
 if [ "$VAR1" = "" ]; then
     echo "No particle generator given."
-    tar -cvf lennard_files.tar pos_e_vel_inicial_usando_sim_anterior.py csv2vtk_particles.py lennard settings.ini verify_settings.py visualizar.py
+    tar -cvf lennard_files.tar potential_analysis.py pos_e_vel_inicial_usando_sim_anterior.py csv2vtk_particles.py lennard settings.ini verify_settings.py visualizar.py
 else
     tar -cvf lennard_files.tar pos_e_vel_inicial_usando_sim_anterior.py csv2vtk_particles.py lennard $VAR1 settings.ini verify_settings.py visualizar.py
 fi

rot_par.py

@@ -0,0 +1,80 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.ticker import MaxNLocator
+
+sigma_ref = 0.341
+
+delta = 0.005
+sig0 = 1
+ep0 = 1
+rs = 5 # raio solido
+rm = .245 / sigma_ref #raio do atomo
+# .245 é a distância entra atomos em estrutura hexagonal rombica
+
+x = y = np.arange(-(3*sig0+rs), 3*sig0+rs, delta)
+X, Y = np.meshgrid(x, y)
+
+r = np.sqrt(X**2 + Y**2)
+np.place(r,r<rs,np.nan)
+gamma = np.arctan(X/Y)
+
+
+A = (0.071/sigma_ref)/2 # no epsilon
+B = (0.071/sigma_ref)/2 # no sigma
+alpha = beta = 2*np.pi*rs/rm
+ph = 2*np.pi/alpha
+#beta = 6
+
+V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0*(1+B*np.sin(beta*gamma))/(r-rs) )**12 - (sig0*(1+B*np.sin(beta*gamma))/(r-rs))**6  )
+
+Fr = 24*ep0*(1+A*np.sin(alpha*gamma+ph))*(1/(r-rs)**2) *  ( (sig0*(1+B*np.sin(beta*gamma))/(r-rs))**6 * (1-2* (sig0*(1+B*np.sin(beta*gamma))/(r-rs))**6 )) *  (r-rs)
+
+Fg = (1/(r-rs))* ( (4*ep0*(A*np.sin(alpha*gamma+ph) + 1 ))* (sig0/(r-rs))**6*(B*np.sin(beta*gamma)+1)**5 * \
+     ( (sig0/(r-rs))**6 * 12*beta*B*np.cos(beta*gamma)*(B*np.sin(beta*gamma)+1)**6 - 6*beta*B*np.cos(beta*gamma)) + \
+    4*alpha*A*ep0*np.cos(alpha*gamma+ph)*(sig0/(r-rs))**6*(B*np.sin(beta*gamma)+1)**6 * \
+    ((sig0/(r-rs))**6 * (B*np.sin(beta*gamma)+1)**6 - 1) ) 
+
+
+# Fg = (1/r)* ( (4*ep0*(A*np.sin(alpha*gamma) + 1 ))* \
+#      ( (sig0/r)**12 * 12*beta*B*np.cos(beta*gamma)*(B*np.sin(beta*gamma)+1)**11 - (sig0/r)**6 * 6*beta*B*np.cos(beta*gamma)*(B*np.sin(beta*gamma)+1)**5) + \
+#     4*alpha*A*ep0*np.cos(alpha*gamma)*((sig0/r)**12 * (B*np.sin(beta*gamma)+1) - (sig0/r)**6*(B*np.sin(beta*gamma)+1)**6) ) 
+
+
+
+# pot 
+np.place(V,abs(V) > 20, np.nan)
+
+plt.figure()
+plt.plot(X[int(len(X)/2),:] ,V[int(len(X)/2),:])
+plt.title("potential")
+
+
+fig1, ax1 = plt.subplots()
+cs1 = ax1.contourf(X, Y, V,levels=20)
+cbar = fig1.colorbar(cs1)
+ax1.set_title("potential")
+
+#force 
+
+np.place(Fr,abs(Fr) > 5000, np.nan)
+
+# plt.figure()
+# plt.plot(X[int(len(X)/2),:] ,Fr[int(len(X)/2),:])
+# plt.title("force")
+
+
+cmap = plt.get_cmap('PiYG')
+fig2, ax2 = plt.subplots()
+cs2 = ax2.contourf(X, Y, Fr,levels=20, cmap=cmap)
+cbar = fig2.colorbar(cs2)
+ax2.set_title('radial force')
+
+
+np.place(Fg,abs(Fg) > 5000, np.nan)
+fig3, ax3 = plt.subplots()
+cs3 = ax3.contourf(X, Y, Fg, levels=20, cmap=cmap)
+cbar = fig3.colorbar(cs3)
+ax3.set_title('tangential force')
+
+
+plt.show()

saida.f90

@@ -49,7 +49,7 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start)
             end do          
         else if (d == 7) 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)
+                write(10,fmt7) v(i,1),v(i,2),v(i,3),v(i,4),v(i,5),v(i,6),v(i,7)
             end do          
         end if
 

settings.ini

@@ -3,18 +3,18 @@
 # Os outros parâmetros são adimensionais. Checar anotações.txt T' = T*k_b/epsilon
 # Argonio: epsilon/kb = 120K, sigma = 0.341nm 
 [global]
-    nimpre  = 100 # quntidade de saidas 
-    N = 4100  # número de partículas 
+    nimpre  = 200 # quntidade de saidas 
+    N = 2  # número de partículas 
     Ntype = 2 # número de tipos de partícula
-    t_fim = 40000  # tempo de execução 
+    t_fim = 20  # tempo de execução 
     nimpre_init = 0 # começa a simular a partir no tempo referente ao snapshot (nimpre). 
     dt = 0.001 # passo de tempo
-    dimX = 10000 # Dimensões da região de cálculo
-    dimY = 1250 #
-    mesh = 200 30 # 300 60 # malha(elementos)
+    dimX = 600 # Dimensões da região de cálculo
+    dimY = 600 #
+    mesh = 8 8 # 300 60 # malha(elementos)
     rcut = 3 # *sigma
-    wall = 'ppee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste
-    termostato_nose_hoover = .true. #liga termostato Nose Hoover. Este funciona durante toda execução.
+    wall = 'eeee' # condição Elastica ou Periodica nas paredes norte, sul, leste, oeste
+    termostato_nose_hoover = .false. #liga termostato Nose Hoover. Este funciona durante toda execução.
     termostato_vel_scaling = .false.  # liga termostato por velocity scaling
     Mc = 1
     Td = -1 #temperatura desejada global constante (-1 = não aplicar) kb = 1.38064852E-23, Td = (2/(3*kb))*E/epsilon(literatura)
@@ -24,40 +24,44 @@
     Td_hot = 2 # temperatura nas cold cells constante. Td_hot = -1 desliga termostato
     Td_cold = 0.8 # v = sqrt(2*T) temperatura nas hot cells constante. Td_cold = -1 desliga termostato
     force_lingradT = -1 # (NUMERO MULTIPLO DE 4 RECOMENDADO!) se > 0 irá dividir região no numero de subregiões entradas na direação X e será forçado um gradiente de temperatura linear. A região dos banhos será ignorada. A temperatura será controlado por um termostato velocity scaling
-    vd = 1.4 1.4 # velocidade distribuida com Maxwell Boltzmann = sqrt(Td') 
+    vd = 0 0 # velocidade distribuida com Maxwell Boltzmann = sqrt(Td') 
     NMPT = 1000 # 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 = .false. # 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. 
 # se ismolecule = .true. então a partícula é levada em conta no termostato 
+# pr = constantes [A, B, alpha, beta] do lennard jones rugoso que simulará um cristal de atomos. Se o valor for dado, a particula poderá girar. 
+#      o grupo com pr deve ser o último e deverá haver apenas um. Ver em rot_par.py como será o potencial
 [par_0]
-    x = 'moleculas8x1_sc.csv' 
-    v = 0 0 # velocidade global
+    x = 'mole.csv' 
+    v = 5 0 # velocidade global
     v_file = '%%v_file_0.csv' # velocidade de cada partícula
     m = 1
     nome = 'g1'
     sigma = 1 # Lennard Jones
     epsilon = 1 # Lennard Jones
-    quantidade = 4000
+    quantidade = 1
     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 não funciona quando o termostato nose hoover está ligado
     ismolecule = .true.			
 [par_1]
-    x = 'particulas8x1_sc.csv' 
-    v = 0 0 # velocidade global
+    x = 'parti.csv' 
+    v = 5 0 # velocidade global
     v_file = '%%v_file_1.csv' # velocidade de cada partícula
     m = 60
     nome = 'g2'
     sigma = 0.94 # Lennard Jones
     epsilon = 0.428 # Lennard Jones
-    quantidade = 100
+    quantidade = 1
     x_lockdelay = 1000 # só vai poder mudar de posição a partir de t = x_lockdelay
     rs = 5 # raio sólido. Posição de partículas na superfície de um sólido de raio rs
     fric_term = 0 # fricção artifical	
-    ismolecule = .false.						
+    ismolecule = .false.
+    pr = 0.104105 0.104105 43.72584 43.72584 0.071847 # constantes [A, B, alpha, beta, ph] do lennard jones rugoso que simulará um cristal de atomos 
+
 #[par_2]
 #    x = 'p_g.csv'
 #    v = 0 0