diff --git a/distribution_particle.py b/distribution_particle.py new file mode 100644 index 0000000..388dff9 --- /dev/null +++ b/distribution_particle.py @@ -0,0 +1,93 @@ +import numpy as np +import vtk +from vtk.util.numpy_support import vtk_to_numpy +import configparser + +def trapz(x,y): + res = 0 + for i in range(1,len(x)): + res += (y[i-1]+y[i])*(x[i-1]+x[i])/2 + return res + +# LÊ A CONFIGURAÇÃO + +config = configparser.ConfigParser() +config.read('settings.ini') +print('Reading settings.ini') +N = int(config['global']['N'].split()[0]) +nimpre = int(config['global']['nimpre'].split()[0]) +ntype = int(config['global']['Ntype'].split()[0]) +nimpre_init = int(config['global']['nimpre_init'].split()[0]) +dimx = int(config['global']['dimX'].split()[0]) +dimy = int(config['global']['dimY'].split()[0]) + +quant = [] + + +for i in range(ntype): + quant.append(int(config['par_'+str(i)]['quantidade'].split()[0])) + +print("Initial step: {}. Final step: {}. Number of particle groups (ntype): {}\n".format(nimpre_init,nimpre,ntype)) + +pgroup = input("Enter the particle group (>= 0, < ntype): ") +stepini = input("Enter the initial step: ") +stepfim = input("Enter the final step: ") + +x = np.zeros((stepfim - stepini, quant[int(pgroup)], 2) + +# LÊ OS VTK + +grupo = "./result/grupo" + pgroup + "_" + +reader = vtk.vtkXMLUnstructuredGridReader() +nfiles = 0 +for file in range(stepini,stepfim) + reader.SetFileName(grupo+str(file)) + reader.Update() + data = reader.GetOutput() + + points = data.GetPoints() + x[nfiles,:,:] = vtk_to_numpy(points.GetData())[:,0:2]# coluna 1 = x, coluna 2 = y + nfiles += 1 + #u = vtk_to_numpy(data.GetPointData().GetArray(0)) # colunas: Vx, Vy, Vz + +# Todas as posições carregadas + +aux1 = input("Enter the number of subdivisions and the direction (x or y).\nIf direction not given, will be x.\n")).split() +direc = 0 +if len(aux1) > 1: + if aux[0] == 'y': + direc = 1 + + +ndensiy = np.zeros((nfiles,aux1)) +dhist = dimx/int(aux1[0]) # distancia no histograma + +for file in range(nfiles): + for i in range(quant[int(pgroup)]): + ndensiy[file, (x[file,i,direc]//dhist) ] += 1 + + +n2density = (np.sum(ndensity,axis=0)/nfiles)/(quant[int(pgroup)]) +xx = np.linspace(0,dimx, int(aux1[0])) +plt.plot(n2density) + +du = quant[int(pgroup)]/dimx #distribuição uniforme + + +DS = trapz(x,n2density*np.log(n2density/du)) + + + + + +#N = reader.GetNumberOfPoints() # numero de pontos +#npa = reader.GetNumberOfPointArrays() # numero de vetores de ponto + +# n_arrays = reader.GetNumberOfPointArrays() +# print("Available data:\n") +# for i in range(n_arrays): +# print(reader.GetPointArrayName(i)) + + + diff --git a/lennard.f90 b/lennard.f90 index 53155b1..9009eae 100644 --- a/lennard.f90 +++ b/lennard.f90 @@ -845,7 +845,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR else if (rs1 > 0 .and. rs2 == 0) then - gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & @@ -869,7 +869,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then - gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 * & (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -956,7 +956,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - fR else if (rs1 > 0 .and. rs2 == 0) then - gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * & (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -971,7 +971,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then - gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 & * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1050,7 +1050,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - fR else if (rs1 > 0 .and. rs2 == 0) then - gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 & * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1066,7 +1066,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then - gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1149,7 +1149,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - fR else if (rs1 > 0 .and. rs2 == 0) then - gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) & * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1163,7 +1163,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then - gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * & (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1240,7 +1240,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, else if (rs1 > 0 .and. rs2 == 0) then - gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & ((sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1254,7 +1254,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then - gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * & (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1331,7 +1331,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - fR else if (rs1 > 0 .and. rs2 == 0) then - gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * & (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1345,7 +1345,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then - gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1423,7 +1423,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - fR else if (rs1 > 0 .and. rs2 == 0) then - gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 & * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1437,7 +1437,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then - gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 * & (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -1516,7 +1516,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR else if (rs1 > 0 .and. rs2 == 0) then - gamma = theta(-pa +ptr%p%n) + atan(coss/sine) + gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * & ( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & @@ -1537,7 +1537,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then - gamma = theta(-pa +ptrn%p%n) + atan(coss/sine) + gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 * & (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * & [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] @@ -6516,6 +6516,7 @@ program main end if !! CASE OF RUGGED LENNARD JONES THAT CAN MAKE THE PARTICLE ROTATE if (abs(pr(4)+pr(3)) > 0) then ! PARTICLE ROTATES + print*, "CASE: Particle Rotates" if (id == 0) then deallocate(x,v) allocate(x(N,3),v(N,3)) @@ -6826,6 +6827,7 @@ program main end do else if (termostato_nose_hoover) then + print*, "CASE: NOSE HOOVER" ! TERMOSTATO NOSE HOOVER do while (t_fim > t) !print*, "L 1990" @@ -7010,6 +7012,7 @@ program main else ! TERMOSTATO SCALING OU SEM TEMOSTATO + print*, "CASE: SCALING OR NO TERMOSTAT" do while (t_fim > t) ! COMPF call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t) !altera Força diff --git a/rot_par.py b/rot_par.py index 35708bb..37f43cd 100644 --- a/rot_par.py +++ b/rot_par.py @@ -1,5 +1,5 @@ import numpy as np -from numpy import sin, cos +from numpy import sin, cos, heaviside import matplotlib.pyplot as plt @@ -9,14 +9,19 @@ def abssin(x): maxF = 40 sigma_ref = 0.341 -delta = 0.005 + +delta = 0.01 sig0 = 1 -ep0 = 4 +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(0, 3*sig0+rs, delta) +# caso com polimero +Np = 6 +deltad = 2 + +x = y = np.arange(-(3*sig0+rs), 3*sig0+rs, delta) X, Y = np.meshgrid(x, y) r = np.sqrt(X**2 + Y**2) @@ -24,9 +29,9 @@ def abssin(x): gamma = np.arctan2(Y,X) -A = 0.0001 #ep0/2 #(0.071/sigma_ref)/4 # no epsilon -B = 0.0001 #(0.071/sigma_ref)/8 # no sigma -alpha = beta = (2*np.pi*rs/rm) +A = 0 #ep0/8 #(0.071/sigma_ref)/4 # no epsilon +B = 0 #(0.071/sigma_ref)/8 # no sigma +alpha = beta = 6 # (2*np.pi*rs/rm)/2 ph = (2*np.pi/alpha) print("A {}, B {}, alpha {}, fase {}, periodo {}\n".format(A, B, alpha ,ph,2*np.pi/alpha)) @@ -35,17 +40,20 @@ def abssin(x): gmolec = np.linspace(0,2*np.pi/4,round((2*np.pi*rs/rm)/4)) #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 ) -#V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs) )**12 - (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs))**6 ) +# LJ rugoso +#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 ) +# com polimero +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 ) + (Np/(2*deltad))*((deltad*2 + sig0) + (r-rs) *(np.log((deltad*2 + sig0)/(r-rs) ) -1))*np.sin(beta*gamma)**2*heaviside(r-(rs+sig0) , 1)*(1-heaviside(r-(deltad*2+rs+sig0) , 1)) +#V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs) )**12 - (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs))**6 ) #V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0/(r-rs*(1+B*np.sin(beta*gamma))) )**12 - (sig0/(r-rs*(1+B*np.sin(beta*gamma))))**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) ) +#rugoso Fr = 24*ep0*(1+A*sin(alpha*gamma+ph))*(1/(r-rs)**2) * ((sig0*(1+B*sin(beta*gamma))/(r-rs))**6 * \ (1-2* (sig0*(1+B*sin(beta*gamma))/(r-rs))**6 )) * (r-rs) @@ -54,6 +62,10 @@ def abssin(x): 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) ) +# polimero +Fr = 24*ep0*(1/(r-rs)**2) * ((sig0/(r-rs))**6 *(1-2* (sig0/(r-rs))**6 )) * (r-rs) + (Np/(2*deltad))*((np.log((deltad*2 + sig0)/(r-rs) ) -1) + (deltad*2+sig0))*np.sin(beta*gamma)**2*heaviside(r-(rs+sig0) , 1)*(1-heaviside(r-(deltad*2+rs+sig0) , 1)) +Fg = 2*beta*cos(beta*gamma)*sin(beta*gamma)*(Np/(2*deltad))*((deltad*2 + sig0) + (r-rs) *(np.log((deltad*2 + sig0)/(r-rs) ) -1))*np.sin(beta*gamma)**2*heaviside(r-(rs+sig0) , 1)*(1-heaviside(r-(deltad*2+rs+sig0) , 1)) +### # Fr = 24*ep0*(1+A*np.sin(alpha*gamma+ph))*(1/(r-rs)**2) * ( (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs))**6 * (1-2* (sig0*(1+B*abs(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*abs(np.sin(beta*gamma))+1)**5 * \ @@ -66,9 +78,9 @@ def abssin(x): np.place(Fr,abs(Fr) > maxF, np.nan) np.place(Fg,abs(Fg) > maxF, np.nan) -plt.figure() -plt.plot(X[int(len(X)/2),:] ,V[int(len(X)/2),:]) -plt.title("potential") +# plt.figure() +# plt.plot(X[int(len(X)/2),:] ,V[int(len(X)/2),:]) +# plt.title("potential") fig1, ax1 = plt.subplots()