diff --git a/comp_FT.f90 b/comp_FT.f90 index 70ede8d..a39656a 100644 --- a/comp_FT.f90 +++ b/comp_FT.f90 @@ -1,4 +1,4 @@ -subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,domy, ids, id, t,dt,partlst) +subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,domy, ids, id, wall, t,dt,partlst) use linkedlist use mod1 use data @@ -20,6 +20,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, type(data_ptr) :: ptr,ptrn integer, intent(out), dimension(:) :: partlst integer ( kind = 4 ), intent(in) :: ids(8) + character(4), intent(in) :: wall partlst = 0 A = pr(1) @@ -59,7 +60,19 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, dox = domx doy = domy end if - + + ! em caso de parede elástica, não precisa ver as celulas fantasma + if (wall(1:2) == 'ee') then + if (doy(1) == 1) doy(1) = 2 + if (doy(2) == mesh(2)+2) doy(2) = mesh(2)+1 + end if + + if (wall(3:4) == 'ee') then + if (dox(1) == 1) dox(1) = 2 + if (dox(2) == mesh(1)+2) dox(2) = mesh(1)+1 + end if + ! + do i = doy(1),doy(2) ! i é linha do j = dox(1),dox(2) node => list_next(malha(i,j)%list) @@ -117,6 +130,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, sigma = sigma_a epsil = epsil_a end if + x2 = ptrn%p%x v2 = ptrn%p%v r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) @@ -133,20 +147,33 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if 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) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) & + + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if Tor(-pa +ptr%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] -fR ! print*, "r1 >0" ! print*, "ptr%p%F" , ptr%p%F, "particula" ! print*, "ptrn%p%F" , ptrn%p%F @@ -159,14 +186,26 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, else if (rs2 > 0 .and. rs1 == 0) then gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) & + + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if Tor(-pa +ptrn%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) - ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] - ptrn%p%F = -aux2 + ptrn%p%F + ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] +fR + ptrn%p%F = -aux2 + ptrn%p%F -fR ! print*, "r2 > 0" ! print*, "ptr%p%F" , ptr%p%F ! print*, "ptrn%p%F" , ptrn%p%F, "particula" @@ -184,11 +223,11 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptrn%p%F = -aux2 + ptrn%p%F -fR end if end if node => list_next(node) ! próxima partícula da célula @@ -219,6 +258,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, epsil = epsil_a end if + x2 = ptrn%p%x v2 = ptrn%p%v r = sqrt((x1(1)-x2(1))**2 + (x1(2)-x2(2))**2) @@ -235,32 +275,58 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%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) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) *& + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if Tor(-pa +ptr%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] -fR ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if Tor(-pa +ptrn%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) - ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] - ptrn%p%F = -aux2 + ptrn%p%F + ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] +fR + ptrn%p%F = -aux2 + ptrn%p%F -fR partlst(ptrn%p%n - pa) = 1 else ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos @@ -272,7 +338,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR @@ -322,30 +388,56 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if - ptr%p%F = aux2 + ptr%p%F +fR + ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR - else if (rs1 > 0 .and. rs2 == 0) then + else if (rs1 > 0 .and. rs2 == 0) then gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss,& + -(x1(2)-x2(2)+B*sin(beta * gamma))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptr%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] -fR ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss,& + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptrn%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) - ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] - ptrn%p%F = -aux2 + ptrn%p%F + ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] +fR + ptrn%p%F = -aux2 + ptrn%p%F -fR partlst(ptrn%p%n - pa) = 1 else ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos @@ -357,7 +449,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR @@ -409,30 +501,56 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR - else if (rs1 > 0 .and. rs2 == 0) then + else if (rs1 > 0 .and. rs2 == 0) then gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptr%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] -fR ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss,& + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptrn%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) - ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] - ptrn%p%F = -aux2 + ptrn%p%F + ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] +fR + ptrn%p%F = -aux2 + ptrn%p%F -fR partlst(ptrn%p%n - pa) = 1 else ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos @@ -444,7 +562,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR @@ -491,32 +609,58 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR - else if (rs1 > 0 .and. rs2 == 0) then + else if (rs1 > 0 .and. rs2 == 0) then gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptr%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] -fR ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptrn%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) - ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] - ptrn%p%F = -aux2 + ptrn%p%F + ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] +fR + ptrn%p%F = -aux2 + ptrn%p%F -fR partlst(ptrn%p%n - pa) = 1 else ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos @@ -528,7 +672,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR @@ -576,31 +720,54 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if 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) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptr%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] -fR ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptrn%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) - ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] - ptrn%p%F = -aux2 + ptrn%p%F + ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] +fR + ptrn%p%F = -aux2 + ptrn%p%F -fR partlst(ptrn%p%n - pa) = 1 else ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos @@ -612,7 +779,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR @@ -660,33 +827,59 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then - fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR - else if (rs1 > 0 .and. rs2 == 0) then + else if (rs1 > 0 .and. rs2 == 0) then gamma = theta(-pa +ptr%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptr%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] -fR ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptrn%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) - ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] - ptrn%p%F = -aux2 + ptrn%p%F + ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] +fR + ptrn%p%F = -aux2 + ptrn%p%F -fR partlst(ptrn%p%n - pa) = 1 else ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos @@ -698,7 +891,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR @@ -710,7 +903,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, end if end if - else ! se for a última lina, só interage com a celua ao lado + else if (j /= mesh(1) + 2) then ! se for a última lina, só interage com a celua ao lado node => list_next(malha(i,j+1)%list) !interagirá com a próxima linha e coluna do while (associated(node)) ptrn = transfer(list_get(node), ptrn) !outra particula selecionada @@ -730,7 +923,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, else rcut = r_cut*sigma_a + (rs1 + rs2) sigma = sigma_a - epsil = epsil__a + epsil = epsil_a end if x2 = ptrn%p%x @@ -749,37 +942,58 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if 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) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) & + + alpha * A * epsil * cos(alpha * gamma) * & + (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) - ! print*, "r1 > 0 gamma", gamma - ! print*, "aux2, aux3",aux2, aux3 - ! print*, "x1 =", x1 - ! print*, "x2 =", x2 - ! read(*,*) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if + Tor(-pa +ptr%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptr%p%n) - ptr%p%F = aux2 + ptr%p%F - ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss] -fR ! partlst(ptr%p%n - pa) = 1 else if (rs2 > 0 .and. rs1 == 0) then gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine) - aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + aux2 = - ( (6 * epsil * sigma**6 * (A * sin(alpha * gamma) + 1) * & + ((B * sin(beta * gamma) + r)**6 - 2 * sigma**6))/(B * sin(beta * gamma) + r)**13 ) * & + [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine] + + aux3 = (1/(B * sin(beta * gamma) + r)) * (epsil * (A * sin(alpha * gamma) + 1)* & + ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - & + (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + & + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 & + - sigma**6/(B * sin(beta * gamma) + r)**6)) - aux3 = (1/r) * (epsil * (A * sin(alpha * gamma) + 1)* ((6 * beta * B * sigma**6 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**7 - (12 * beta * B * sigma**12 * cos(beta * gamma))/(B * sin(beta * gamma) + r)**13) + alpha * A * epsil * cos(alpha * gamma) * (sigma**12/(B * sin(beta * gamma) + r)**12 - sigma**6/(B * sin(beta * gamma) + r)**6)) + fric_term = (fric_term1+fric_term2)/2 + if (fric_term > 0) then + fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2+B*sin(beta * gamma))*coss, & + -(x1(2)-x2(2))+(rs1+rs2+B*sin(beta * gamma))*sine], & + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) + end if Tor(-pa +ptrn%p%n) = -aux3*(r+rs1+rs2) + Tor(-pa +ptrn%p%n) - ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] - ptrn%p%F = -aux2 + ptrn%p%F + ptr%p%F = aux2 + ptr%p%F + aux3*[sine,-coss] +fR + ptrn%p%F = -aux2 + ptrn%p%F -fR partlst(ptrn%p%n - pa) = 1 else ! Não vou estudar aqui a interação entre duas partículas, elas vão se repelir como átomos @@ -791,7 +1005,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx, fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR @@ -946,7 +1160,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR ! ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -956,7 +1170,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -979,7 +1193,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1005,7 +1219,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR @@ -1056,7 +1270,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR ! ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1066,7 +1280,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1))*(sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1))*(sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)* & ! (B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & @@ -1081,7 +1295,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1100,7 +1314,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR @@ -1150,7 +1364,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR ! ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1160,7 +1374,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* & ! (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)* & ! (B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & @@ -1176,7 +1390,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* & ! (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 & ! - 6*beta*B*cos(beta*gamma)) + & @@ -1197,7 +1411,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR @@ -1249,7 +1463,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR ! ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1259,7 +1473,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1273,7 +1487,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)*( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))*( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1292,7 +1506,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR @@ -1339,7 +1553,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR ! ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1350,7 +1564,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1364,7 +1578,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1383,7 +1597,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR @@ -1431,7 +1645,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR ! ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1441,7 +1655,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1455,7 +1669,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1474,7 +1688,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR @@ -1523,7 +1737,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR ! ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1533,7 +1747,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1547,7 +1761,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1566,7 +1780,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR @@ -1617,7 +1831,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR ! ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1627,7 +1841,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1647,7 +1861,7 @@ end subroutine comp_FT ! 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] -! aux3 = (1/r)* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & +! aux3 = (1/(B * sin(beta * gamma) + r))* ( (4*epsil*(A*sin(alpha*gamma+ph) + 1 ))* (sigma/r)**6*(B*sin(beta*gamma)+1)**5 * & ! ((sigma/r)**6 * 12*beta*B*cos(beta*gamma)*(B*sin(beta*gamma)+1)**6 - 6*beta*B*cos(beta*gamma)) + & ! 4*alpha*A*epsil*cos(alpha*gamma+ph)*(sigma/r)**6*(B*sin(beta*gamma)+1)**6 * & ! ((sigma/r)**6 * (B*sin(beta*gamma)+1)**6 - 1)) @@ -1666,7 +1880,7 @@ end subroutine comp_FT ! fric_term = (fric_term1+fric_term2)/2 ! if (fric_term > 0) then ! fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & -! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) +! [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) ! end if ! ptr%p%F = aux2 + ptr%p%F +fR diff --git a/distribution_particle.py b/distribution_particle.py index 388dff9..d8ef9e9 100644 --- a/distribution_particle.py +++ b/distribution_particle.py @@ -1,25 +1,45 @@ import numpy as np -import vtk from vtk.util.numpy_support import vtk_to_numpy import configparser - +import pandas as pd +import os.path +from glob import glob +import matplotlib.pyplot as plt +try: + import vtk +except: + print("Vtk module not found! If using anaconda try:\nconda install -c anaconda vtk\n OR\n pip install vtk") +## -------- Funções --------- ## 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 +#---------------------------------------------------------------# + +# escolhe a pasta + +dirname = os.getcwd() #os.path.dirname(os.path.abspath(__file__)) +dirlist = glob(dirname + "/*/") +print("Choose a folder there the results are contained:\nNo | Folder") +for a in range(len(dirlist)): + print("{} | {}\n".format(a,dirlist[a])) +a = int(input("Enter the number of the folder\n")) +res_dir = dirlist[a] +res_name = res_dir.split('/')[-2] # nome da pasta com o resultado + # LÊ A CONFIGURAÇÃO config = configparser.ConfigParser() -config.read('settings.ini') -print('Reading settings.ini') +config.read(res_dir + 'settings.txt') +print('Reading settings.txt') 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]) +dimx = float(config['global']['dimX'].split()[0]) +dimy = float(config['global']['dimY'].split()[0]) quant = [] @@ -29,56 +49,83 @@ def trapz(x,y): 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) +pgroup = input("Enter the particle group (>= 0, < ntype):[{}] ".format(ntype-1)) +if (pgroup == ''): + pgroup = str(ntype-1) + print("pgroup = {}\n".format(pgroup)) +stepini = input("Enter the initial step:[0] ") +if (stepini == ''): + stepini = 0 + print("stepini = 0\n") +else: + stepini = int(stepini) +stepfim = input("Enter the final step :[{}] ".format(nimpre)) +if (stepfim == ''): + stepfim = nimpre + print("stepfim = {}".format(nimpre)) +else: + stepfim = int(stepfim) + +aux1 = input("Number of files to use (will skip some if less than given above):[{}] ".format(nimpre/2)) +if aux1 == '': + aux1 = nimpre/2 + print("jump = {}\n".format(int(nimpre/2))) +else: + aux1 = int(aux1) +jump = int((stepfim - stepini)/aux1) + +x = np.zeros((stepfim - stepini, quant[int(pgroup)], 2)) # LÊ OS VTK -grupo = "./result/grupo" + pgroup + "_" +grupo = res_dir + "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 +proximo = stepini +for file in range(stepini,stepfim): + if file == proximo: + reader.SetFileName(grupo+str(file)+'.vtu') + 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 + proximo += jump + #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() +aux1 = input("Enter the number of subdivisions and the direction (x or y).\nIf direction not given, will be x.\n [10 x]").split() +if len(aux1) == 0: + aux1 = ['10'] + print("There will be 10 subdivisions.\n") direc = 0 if len(aux1) > 1: - if aux[0] == 'y': + if aux1[0] == 'y': direc = 1 -ndensiy = np.zeros((nfiles,aux1)) +ndensiy = np.zeros((nfiles,int(aux1[0]))) 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 + ndensiy[file, int(x[file,i,direc]//dhist) ] += 1 -n2density = (np.sum(ndensity,axis=0)/nfiles)/(quant[int(pgroup)]) +ndensity2 = (np.sum(ndensiy,axis=0)/nfiles)/(quant[int(pgroup)]) xx = np.linspace(0,dimx, int(aux1[0])) -plt.plot(n2density) +plt.plot(ndensity2) du = quant[int(pgroup)]/dimx #distribuição uniforme -DS = trapz(x,n2density*np.log(n2density/du)) - - - +DS = trapz(xx,ndensity2*np.log(ndensity2/du)) +print("Relative distribution relative to the Uniform: {}\n".format(DS)) +plt.show() #N = reader.GetNumberOfPoints() # numero de pontos @@ -91,3 +138,36 @@ def trapz(x,y): +##-- SALVA dados em arquivos csv --## +if os.path.exists('Distributions.csv'): + csv_input = pd.read_csv('Distributions.csv') + headers_list = list(csv_input) + + laux1 = True # variavel auxiliar lógica + while laux1: + for name in headers_list: + if name == res_name: + print("A file with name {} already exists.") + res_name = input('Enter a new name or ctrl+C to abort:\n') + laux1 = True + break + else: + laux1 = False + + csv_input[res_name] = ndensity2 +else: + csv_input = pd.DataFrame(ndensity2,columns=[res_name]) + + +csv_input.to_csv('Distributions.csv', index=False) + +#--# + + +f=open('Distrib_entropies.csv', "a+") +f.write("{}, {}\n".format(res_name,DS)) + + + + + diff --git a/lennard.f90 b/lennard.f90 index cb19598..3d192f2 100644 --- a/lennard.f90 +++ b/lennard.f90 @@ -9,7 +9,7 @@ function magtorque(angle, hfield, t) result(Tor) real(dp) :: t, Tor if (t > hfield(5)) then - Tor = hfield(1)*(cos(angle)*hfield(3) - sin(angle)*hfield(2))*sin(t*hfield(4)) + Tor = hfield(1)*(cos(angle)*hfield(3) - sin(angle)*hfield(2))*cos(t*hfield(4)) else Tor = 0 end if @@ -67,12 +67,12 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade type(container), allocatable,dimension(:,:),intent(in) :: malha 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) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut 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, dox(2), doy(2) integer, intent(in) :: mesh(:),domx(2),domy(2), id - real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2 + real(dp) :: Fi(2)=0,r, aux2(2),fR(2) type(list_t), pointer :: node, next_node type(data_ptr) :: ptr,ptrn integer ( kind = 4 ), intent(in) :: ids(8) @@ -108,7 +108,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) m1 = propriedade(ptr%p%grupo)%m p1 = [v1(1), v1(2)]*m1 rs1 = propriedade(ptr%p%grupo)%rs !raio sólido - fric_term1 = propriedade(ptr%p%grupo)%fric_term sigma_a = propriedade(ptr%p%grupo)%sigma ! rcut = r_cut*sigma epsil_a = propriedade(ptr%p%grupo)%epsilon @@ -178,7 +177,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) sigma_b = propriedade(ptrn%p%grupo)%sigma epsil_b = propriedade(ptrn%p%grupo)%epsilon rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term ! Lorenz-Betherlot rule for mixing epsilon sigma if (sigma_a > sigma_b) then rcut = r_cut*sigma_a + (rs1 + rs2) @@ -221,7 +219,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) sigma_b = propriedade(ptrn%p%grupo)%sigma epsil_b = propriedade(ptrn%p%grupo)%epsilon rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term ! Lorenz-Betherlot rule for mixing epsilon sigma if (sigma_a > sigma_b) then rcut = r_cut*sigma_a + (rs1 + rs2) @@ -267,7 +264,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) sigma_b = propriedade(ptrn%p%grupo)%sigma epsil_b = propriedade(ptrn%p%grupo)%epsilon rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term ! Lorenz-Betherlot rule for mixing epsilon sigma if (sigma_a > sigma_b) then rcut = r_cut*sigma_a + (rs1 + rs2) @@ -311,7 +307,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) sigma_b = propriedade(ptrn%p%grupo)%sigma epsil_b = propriedade(ptrn%p%grupo)%epsilon rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term ! Lorenz-Betherlot rule for mixing epsilon sigma if (sigma_a > sigma_b) then rcut = r_cut*sigma_a + (rs1 + rs2) @@ -353,7 +348,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) sigma_b = propriedade(ptrn%p%grupo)%sigma epsil_b = propriedade(ptrn%p%grupo)%epsilon rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term ! Lorenz-Betherlot rule for mixing epsilon sigma if (sigma_a > sigma_b) then rcut = r_cut*sigma_a + (rs1 + rs2) @@ -397,7 +391,6 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) sigma_b = propriedade(ptrn%p%grupo)%sigma epsil_b = propriedade(ptrn%p%grupo)%epsilon rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term ! Lorenz-Betherlot rule for mixing epsilon sigma if (sigma_a > sigma_b) then rcut = r_cut*sigma_a + (rs1 + rs2) @@ -436,14 +429,13 @@ subroutine comp_pot(mesh,malha,propriedade,r_cut,domx,domy, ids, id, t, nRfu) end if end if - else ! se for a última lina, só interage com a celua ao lado + else if (j /= mesh(1) + 2) then ! se for a última lina, só interage com a celua ao lado node => list_next(malha(i,j+1)%list) !interagirá com a próxima linha e coluna do while (associated(node)) ptrn = transfer(list_get(node), ptrn) !outra particula selecionada sigma_b = propriedade(ptrn%p%grupo)%sigma epsil_b = propriedade(ptrn%p%grupo)%epsilon rs2 = propriedade(ptrn%p%grupo)%rs - fric_term2 = propriedade(ptrn%p%grupo)%fric_term ! Lorenz-Betherlot rule for mixing epsilon sigma if (sigma_a > sigma_b) then rcut = r_cut*sigma_a + (rs1 + rs2) @@ -696,12 +688,13 @@ subroutine corr_K(malha,domx,domy,cold_cells,hot_cells,Td_hot,Td_cold,propriedad end subroutine corr_K !força de ficção - function comp_fric(r,v,fric_term) result(f) + function comp_fric(r,v,fric_term,sigma) result(f) use mod1 - real(dp), intent(in) :: r(2),v(2),fric_term - real(dp) :: f(2) + real(dp), intent(in) :: r(2),v(2),fric_term,sigma + real(dp) :: f(2),aux - f = (fric_term/(r(1)**2 + r(2)**2))*(v(1)*r(1)+v(2)*r(2))*[r(1),r(2)] + aux = (0.56*sigma**2*fric_term/(r(1)**2 + r(2)**2)**2) + if (aux > 0) f = aux*(v(1)*r(1)+v(2)*r(2))*[r(1),r(2)] end function comp_fric ! atualiza forças @@ -709,7 +702,7 @@ end function comp_fric include 'comp_FT.f90' ! atualiza forças - subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) + subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id,wall, t) use linkedlist use mod1 use data @@ -727,11 +720,13 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) type(list_t), pointer :: node, next_node type(data_ptr) :: ptr,ptrn integer ( kind = 4 ), intent(in) :: ids(8) + character(4), intent(in) :: wall !Lennard Jones fR = [0,0] dox = domx doy = domy + if (sum(ids(1:4)) > -4) then !caso paralelo if (domx(1) > 1) dox(1) = domx(1) - 1 @@ -742,7 +737,18 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) dox = domx doy = domy end if - + ! em caso de parede elástica, não precisa ver as celulas fantasma + if (wall(1:2) == 'ee') then + if (doy(1) == 1) doy(1) = 2 + if (doy(2) == mesh(2)+2) doy(2) = mesh(2)+1 + end if + + if (wall(3:4) == 'ee') then + if (dox(1) == 1) dox(1) = 2 + if (dox(2) == mesh(1)+2) dox(2) = mesh(1)+1 + end if + + do i = doy(1),doy(2) ! i é linha do j = dox(1),dox(2) node => list_next(malha(i,j)%list) @@ -802,13 +808,11 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) 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] - ! print*, "L 395 r", r, "id",id - !print '("(x1(1)-x2(1)) ", f7.3 ," (rs1+rs2)*coss ", f7.3, " (x1(2)-x2(2)) ", f7.3, " (rs1+rs2)*sine ", f7.3 )', (x1(1)-x2(1)), & - !(rs1+rs2)*coss, (x1(2)-x2(2)), (rs1+rs2)*sine + fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -854,7 +858,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -898,7 +902,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -945,7 +949,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -987,7 +991,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1030,7 +1034,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1074,7 +1078,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1084,7 +1088,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) end if end if - else ! se for a última lina, só interage com a celua ao lado + else if (j /= mesh(1) + 2) then ! se for a última lina, só interage com a celua ao lado node => list_next(malha(i,j+1)%list) !interagirá com a próxima linha e coluna do while (associated(node)) ptrn = transfer(list_get(node), ptrn) !outra particula selecionada @@ -1104,7 +1108,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) else rcut = r_cut*sigma_a + (rs1 + rs2) sigma = sigma_a - epsil = epsil__a + epsil = epsil_a end if x2 = ptrn%p%x @@ -1120,7 +1124,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) fric_term = (fric_term1+fric_term2)/2 if (fric_term > 0) then fR = comp_fric([-(x1(1)-x2(1))+ (rs1+rs2)*coss,-(x1(2)-x2(2))+(rs1+rs2)*sine], & - [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term) + [-(v1(1)-v2(1)),-(v1(2)-v2(2))],fric_term,sigma) end if ptr%p%F = aux2 + ptr%p%F +fR ptrn%p%F = -aux2 + ptrn%p%F - fR @@ -1137,7 +1141,7 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) ! call MPI_barrier(MPI_COMM_WORLD, ierr) end subroutine comp_F - subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) + subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, wall,t) use linkedlist use mod1 use data @@ -1155,12 +1159,13 @@ subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id type(list_t), pointer :: node, next_node type(data_ptr) :: ptr,ptrn integer ( kind = 4 ), intent(in) :: ids(8) + character(4), intent(in) :: wall !Lennard Jones dox = domx doy = domy if (sum(ids(1:4)) > -4) then - !caso paralelo + ! print*, "caso paralelo" if (domx(1) > 1) dox(1) = domx(1) - 1 if (domy(1) > 1) doy(1) = domy(1) - 1 if (domx(2) < mesh(1)+2) dox(2) = domx(2) + 1 @@ -1169,7 +1174,18 @@ subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id dox = domx doy = domy end if - + + ! em caso de parede elástica, não precisa ver as celulas fantasma + if (wall(1:2) == 'ee') then + if (doy(1) == 1) doy(1) = 2 + if (doy(2) == mesh(2)+2) doy(2) = mesh(2)+1 + end if + + if (wall(3:4) == 'ee') then + if (dox(1) == 1) dox(1) = 2 + if (dox(2) == mesh(1)+2) dox(2) = mesh(1)+1 + end if + do i = doy(1),doy(2) ! i é linha do j = dox(1),dox(2) node => list_next(malha(i,j)%list) @@ -1467,7 +1483,7 @@ subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id end if end if - else ! se for a última lina, só interage com a celua ao lado + else if (j /= mesh(1) + 2) then ! se for a última lina, só interage com a celua ao lado node => list_next(malha(i,j+1)%list) !interagirá com a próxima linha e coluna do while (associated(node)) ptrn = transfer(list_get(node), ptrn) !outra particula selecionada @@ -1551,7 +1567,6 @@ subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx, east = wall(3:3) west = wall(4:4) ! if (id == 0) read(*,*) - ! Esvazia as celulas emprestadas if (domy(1) > 1) then i = domy(1)-1 @@ -4902,6 +4917,7 @@ subroutine walls(icell,jcell,mesh,malha,domx,domy,wall,subx,suby,np,id,mic) end if end do end do + i = 2 ! copiada pra celula fantasma do lado oposto do j = domx(1),domx(2) previous_node => malha(i,j)%list @@ -4987,9 +5003,9 @@ subroutine walls(icell,jcell,mesh,malha,domx,domy,wall,subx,suby,np,id,mic) ! print*, "n", ptr%p%x allocate(ptrn%p) ptrn%p = ptr%p - ptrn%p%x(1) = ptr%p%x(1) + jcell(mesh(1)+1) - ptr%p%flag = .true. - call list_insert(malha(i,mesh(1)+2)%list, data=transfer(ptr,list_data)) + ptrn%p%x(1) = ptrn%p%x(1) + jcell(mesh(1)+1) + ptrn%p%flag = .true. + call list_insert(malha(i,mesh(1)+2)%list, data=transfer(ptrn,list_data)) node => list_next(node) end do end do @@ -5058,7 +5074,7 @@ subroutine walls(icell,jcell,mesh,malha,domx,domy,wall,subx,suby,np,id,mic) ptrn%p = ptr%p ptrn%p%x(1) = ptr%p%x(1) - jcell(mesh(1)+1) ptrn%p%flag = .true. - call list_insert(malha(i,1)%list, data=transfer(ptr,list_data)) + call list_insert(malha(i,1)%list, data=transfer(ptrn,list_data)) node => list_next(node) end do end do @@ -5460,7 +5476,8 @@ program main end if i = 0 - if (Td == 0) then + if (Td == 0 .and. (termostato_nose_hoover .or. termostato_vel_scaling)) then + print*, "TD = 0. The desired temperature will be defined by the given velocities." 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 @@ -5687,13 +5704,24 @@ program main theta = PI*theta partlst = 0 omega = 0 + + if (Td > 0) then + print*, "Global temperature specified. Td = ", Td + else if (Td_cold > 0) then + print*, "Regional temperatures specified:", Td_cold, Td_hot + else + print*, "No temperature specified!" + end if ! Calcula a rotação da partícula. Velocity scaling. do while (t_fim > t) ! COMPF - ! if (id == 0) read(*,*) - ! call MPI_barrier(MPI_COMM_WORLD, ierr) - ! print*, "L 6544" - call comp_FT(GField,hfield,theta,Tor,pr, mesh,malha,propriedade,rcut,domx,domy,ids,id,t,dt,partlst) !altera Força + ! if (id == 0) then + ! ! print*, t + ! read(*,*) + ! end if + call MPI_barrier(MPI_COMM_WORLD, ierr) + ! print*, "L 6544", id + call comp_FT(GField,hfield,theta,Tor,pr, mesh,malha,propriedade,rcut,domx,domy,ids,id,wall,t,dt,partlst) !altera Força ! print*, "tor", tor ! print*, "omega",omega ! print*, "theta",theta @@ -5701,14 +5729,14 @@ program main ! IDS são os ids das regiões vizinhas, vetor(4) int posições, ! DOMX e DOMY são vetores(2) int com o domínio (quat. de celulas) ! que o processo vai cuidar (subdivisões) - ! print*, "L 6408" + ! print*, "L 6408", id if (propriedade(gruporot)%x_lockdelay > t) omega = 0 theta = theta*partlst omega = omega*partlst ! COMP X call comp_xT(icell,jcell,malha,N,mesh,propriedade,Tor,theta,omega, dx_max,t,dt,ids,LT,domx,domy,wall,id, np) ! altera posição - ! print*, "L 6460" + ! print*, "L 6460", id ! print*, "tor", tor ! print*, "omega",omega ! print*, "theta",theta @@ -5756,9 +5784,7 @@ program main ! aux5(num_rot+1:2*num_rot) = omega ! " " ", omega ! aux5(2*num_rot+1:3*num_rot) = Tor ! " " ", Tor end if - end if - ! print*, "L 6609", np - if (np > 1) then + call MPI_BCAST(aux5, 3*num_rot, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) if (id /= 0) then theta = aux5(1:num_rot) @@ -5767,11 +5793,21 @@ program main end if end if + ! do aux4 = 0, 3 + ! if (id == aux4) then + ! print*, "ID #### THETA", id + ! print*, theta + ! print*, "########" + ! end if + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + ! end do + + ! print*, "L 6619" call walls(icell,jcell,mesh,malha,domx,domy,wall,subx,suby,np,id,mic) ! altera posição e malha ! print*, "L 6621" ! COMP F - call comp_FT(GField,hfield,theta,Tor,pr, mesh,malha,propriedade,rcut,domx,domy,ids,id,t,dt,partlst) !altera força + call comp_FT(GField,hfield,theta,Tor,pr, mesh,malha,propriedade,rcut,domx,domy,ids,id,wall,t,dt,partlst) !altera força ! print*, "L 6470",tor ! print*, "L 6625" ! print*, "tor", tor @@ -5811,8 +5847,7 @@ program main if (termostato_vel_scaling .and. i == interv_Td(j2)) then if (Td > 0) then call corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,t) - end if - if (Td_hot > 0 .and. Td_cold > 0) then + else if (Td_cold > 0) then call corr_K(malha,domx,domy,cold_cells,hot_cells,Td_hot,Td_cold,propriedade, np,id,t) end if ! j2 = j2 + 1 @@ -5997,7 +6032,7 @@ program main !print*, "L 1990" ! COMPF - call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t) !altera Força + call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,wall,t) !altera Força ! IDS são os ids das regiões vizinhas, vetor(4) int posições, ! DOMX e DOMY são vetores(2) int com o domínio (quat. de celulas) ! que o processo vai cuidar (subdivisões) @@ -6014,7 +6049,7 @@ program main call comp_v_thermo_pred(malha,mesh,dt,t,np,propriedade,domx,domy, Mc, xih, xic, Td_hot, Td_cold, hot_cells, cold_cells) ! COMP F - call comp_F_thermo(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t) !altera força + call comp_F_thermo(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,wall,t) !altera força ! COMP V call comp_v_thermo_corr(malha,mesh,dt,t,propriedade,domx,domy, xih, xic, hot_cells, cold_cells) @@ -6177,10 +6212,18 @@ program main else ! TERMOSTATO SCALING OU SEM TEMOSTATO if (id==0) print*, "CASE: SCALING OR NO TERMOSTAT" + if (Td > 0) then + print*, "Global temperature specified. Td = ", Td + else if (Td_cold > 0) then + print*, "Regional temperatures specified:", Td_cold, Td_hot + else + print*, "No temperature specified!" + end if + do while (t_fim > t) ! COMPF ! print*, "L 6151" - call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t) !altera Força + call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,wall,t) !altera Força ! IDS são os ids das regiões vizinhas, vetor(4) int posições, ! DOMX e DOMY são vetores(2) int com o domínio (quat. de celulas) ! que o processo vai cuidar (subdivisões) @@ -6195,12 +6238,10 @@ program main call walls(icell,jcell,mesh,malha,domx,domy,wall,subx,suby,np,id,mic) ! altera posição e malha ! COMP F - call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t) !altera força - ! print*, "L 1717" + call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,wall,t) !altera força ! COMP V call comp_v(malha,mesh,dt,t,propriedade,domx,domy) !altera velocidade - ! print*, "1721" t = t + dt cont2 = cont2+1 @@ -6225,8 +6266,7 @@ program main if (termostato_vel_scaling .and. i == interv_Td(j2)) then if (Td > 0) then call corr_Kglobal(malha,domx,domy,Td,propriedade, np,id,t) - end if - if (Td_hot > 0 .and. Td_cold > 0) then + else if (Td_cold > 0) then call corr_K(malha,domx,domy,cold_cells,hot_cells,Td_hot,Td_cold,propriedade, np,id,t) end if ! j2 = j2 + 1 diff --git a/part_gen_pequeno.py b/part_gen_pequeno.py index f1c3ae8..80a2e17 100644 --- a/part_gen_pequeno.py +++ b/part_gen_pequeno.py @@ -24,14 +24,14 @@ import matplotlib.pyplot as plt from random import randint -ff = 8 # formato da região ff = x/y -Lx = 2000 -arquivo2 = 'particulas8x1.csv' -arquivo1 = 'moleculas8x1.csv' - -ratio = 0.8 -N1 = 5000 -N2 = 130 +ff = 1 # formato da região ff = x/y +Lx = 200 +arquivo2 = 'parp.csv' +arquivo1 = 'molp.csv' + +ratio = 0.5 +N1 = 1000 +N2 = 30 L1 = Lx L2 = Lx/ff diff --git a/pos_e_vel_inicial_usando_sim_anterior.py b/pos_e_vel_inicial_usando_sim_anterior.py index 51c523e..673279f 100644 --- a/pos_e_vel_inicial_usando_sim_anterior.py +++ b/pos_e_vel_inicial_usando_sim_anterior.py @@ -17,6 +17,7 @@ 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')