diff --git a/lennard.f90 b/lennard.f90 index 3c10725..dee7dc6 100644 --- a/lennard.f90 +++ b/lennard.f90 @@ -687,9 +687,866 @@ function comp_fric(r,v,fric_term) result(f) f = (fric_term/(r(1)**2 + r(2)**2))*(v(1)*r(1)+v(2)*r(2))*[r(1),r(2)] end function comp_fric - + ! atualiza forças - subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) + subroutine comp_FT(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) + use linkedlist + use mod1 + use data + use mod0 + use mpi + + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + type(container), allocatable,dimension(:,:),intent(in) :: malha + real(dp), intent(in) :: GField(2), t + real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut, fric_term !fric_term = força de ficção + real(dp) :: x1(2),v1(2),x2(2),v2(2), rs1, rs2, coss, sine + integer :: i,j,ct = 0, dox(2), doy(2) !,ptr, ptrn + integer, intent(in) :: mesh(:),domx(2),domy(2), id + real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2 + type(list_t), pointer :: node, next_node + type(data_ptr) :: ptr,ptrn + integer ( kind = 4 ), intent(in) :: ids(8) + + !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 + if (domy(1) > 1) doy(1) = domy(1) - 1 + if (domx(2) < mesh(1)+2) dox(2) = domx(2) + 1 + if (domy(2) < mesh(2)+2) doy(2) = domy(2) + 1 + else + dox = domx + doy = domy + end if + + do i = doy(1),doy(2) ! i é linha + do j = dox(1),dox(2) + node => list_next(malha(i,j)%list) + ! if (associated(node)) then + ! ptr = transfer(list_get(node), ptr) + ! ! print*, "L 253", ptr%p%n + ! end if + do while (associated(node)) + ! print*, "Encontrada", ct, "celulas, id=", id + ptr = transfer(list_get(node), ptr) !particula selecionada + ptr%p%flag = .true. ! indica ao comp_x que a partícula precisa ser calculada + x1 = ptr%p%x + v1 = ptr%p%v + m1 = propriedade(ptr%p%grupo)%m + 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 + if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then + ptr%p%F = ptr%p%F + GField*m1 + end if + ! print '("x1 =", f10.6, " ", f10.6, " n ", i2, " cell ",i2," ",i2)',x1(1),x1(2), ptr%p%n,i,j + ! if (id == 0) read(*,*) + !calcular a força desta com todas as outras partículas + next_node => list_next(node) ! próxima partícula da célula + node => list_next(node) ! a ser computado com a particula selecionada + + ! NA PRÓPRIA CELULA + 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 + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a *epsil_b ) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + rs1 + rs2 + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a *epsil_b ) + else + rcut = r_cut*sigma_a + rs1 + rs2 + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 389 r", r, "id",id + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + !Células ao redor !i é linha, j é coluna + if (i /= mesh(2)+2) then ! se não for a última linha + if (j == mesh(1)+2) then ! se for a última coluna + node => list_next(malha(i+1,j)%list) !interagirá com a próxima linha apenas + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 666 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + if (j /= 1) then !se não for a primeira coluna + node => list_next(malha(i+1,j-1)%list) !interagirá com a próxima linha apenas + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 710 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + else + !interagirá com a próxima linha e coluna, e na diagonal + node => list_next(malha(i,j+1)%list) + + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 757 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + node => list_next(malha(i+1,j)%list) !interagirá com a próxima linha + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 798 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + node => list_next(malha(i+1,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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 841 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + if (j /= 1) then !se não for a primeira coluna + node => list_next(malha(i+1,j-1)%list) !interagirá com a próxima linha apenas + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 885 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + end if + + else ! 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 931 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + node => next_node + end do + end do + end do + !print*, 'Linha 143' + ! if (id == 0) read(*,*) + ! call MPI_barrier(MPI_COMM_WORLD, ierr) + end subroutine comp_FT + + ! atualiza forças + subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) + use linkedlist + use mod1 + use data + use mod0 + use mpi + + type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade + type(container), allocatable,dimension(:,:),intent(in) :: malha + real(dp), intent(in) :: GField(2), t + real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut, fric_term !fric_term = força de ficção + real(dp) :: x1(2),v1(2),x2(2),v2(2), rs1, rs2, coss, sine + integer :: i,j,ct = 0, dox(2), doy(2) !,ptr, ptrn + integer, intent(in) :: mesh(:),domx(2),domy(2), id + real(dp) :: Fi(2)=0,r, aux2(2),fR(2), fric_term1, fric_term2 + type(list_t), pointer :: node, next_node + type(data_ptr) :: ptr,ptrn + integer ( kind = 4 ), intent(in) :: ids(8) + + !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 + if (domy(1) > 1) doy(1) = domy(1) - 1 + if (domx(2) < mesh(1)+2) dox(2) = domx(2) + 1 + if (domy(2) < mesh(2)+2) doy(2) = domy(2) + 1 + else + dox = domx + doy = domy + end if + + do i = doy(1),doy(2) ! i é linha + do j = dox(1),dox(2) + node => list_next(malha(i,j)%list) + ! if (associated(node)) then + ! ptr = transfer(list_get(node), ptr) + ! ! print*, "L 253", ptr%p%n + ! end if + do while (associated(node)) + ! print*, "Encontrada", ct, "celulas, id=", id + ptr = transfer(list_get(node), ptr) !particula selecionada + ptr%p%flag = .true. ! indica ao comp_x que a partícula precisa ser calculada + x1 = ptr%p%x + v1 = ptr%p%v + m1 = propriedade(ptr%p%grupo)%m + 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 + if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then + ptr%p%F = ptr%p%F + GField*m1 + end if + ! print '("x1 =", f10.6, " ", f10.6, " n ", i2, " cell ",i2," ",i2)',x1(1),x1(2), ptr%p%n,i,j + ! if (id == 0) read(*,*) + !calcular a força desta com todas as outras partículas + next_node => list_next(node) ! próxima partícula da célula + node => list_next(node) ! a ser computado com a particula selecionada + + ! NA PRÓPRIA CELULA + 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 + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a *epsil_b ) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + rs1 + rs2 + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a *epsil_b ) + else + rcut = r_cut*sigma_a + rs1 + rs2 + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 389 r", r, "id",id + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + !Células ao redor !i é linha, j é coluna + if (i /= mesh(2)+2) then ! se não for a última linha + if (j == mesh(1)+2) then ! se for a última coluna + node => list_next(malha(i+1,j)%list) !interagirá com a próxima linha apenas + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 666 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + if (j /= 1) then !se não for a primeira coluna + node => list_next(malha(i+1,j-1)%list) !interagirá com a próxima linha apenas + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 710 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + else + !interagirá com a próxima linha e coluna, e na diagonal + node => list_next(malha(i,j+1)%list) + + do while (associated(node)) + ptrn = transfer(list_get(node), ptrn) !outra particula selecionada + sigma_b = propriedade(ptrn%p%grupo)%sigma + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 757 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + node => list_next(malha(i+1,j)%list) !interagirá com a próxima linha + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 798 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + node => list_next(malha(i+1,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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 841 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + + if (j /= 1) then !se não for a primeira coluna + node => list_next(malha(i+1,j-1)%list) !interagirá com a próxima linha apenas + 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 885 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + end if + + else ! 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) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else if (sigma_a < sigma_b) then + rcut = r_cut*sigma_b + (rs1 + rs2) + sigma = 0.5*(sigma_a + sigma_b) + epsil = sqrt(epsil_a*epsil_b) + else + rcut = r_cut*sigma_a + (rs1 + rs2) + 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) + coss = (x1(1)-x2(1))/r + sine = (x1(2)-x2(2))/r + r = r - rs1 - rs2 !raio + ! print*, "L 931 r", r, "id",id + 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] + 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) + end if + ptr%p%F = aux2 + ptr%p%F +fR + ptrn%p%F = -aux2 + ptrn%p%F - fR + end if + node => list_next(node) ! próxima partícula da célula + end do + end if + node => next_node + end do + end do + end do + !print*, 'Linha 143' + ! if (id == 0) read(*,*) + ! 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) use linkedlist use mod1 use data @@ -699,17 +1556,16 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) type(prop_grupo), allocatable,dimension(:),intent(in) :: propriedade type(container), allocatable,dimension(:,:),intent(in) :: malha real(dp), intent(in) :: GField(2), t - real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut, fric_term !fric_term = força de ficção + real(dp) :: sigma, epsil, sigma_a, epsil_a,sigma_b, epsil_b, rcut,r_cut real(dp) :: x1(2),v1(2),x2(2),v2(2), rs1, rs2, coss, sine integer :: i,j,ct = 0, dox(2), doy(2) !,ptr, ptrn 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) type(list_t), pointer :: node, next_node type(data_ptr) :: ptr,ptrn integer ( kind = 4 ), intent(in) :: ids(8) !Lennard Jones - fR = [0,0] dox = domx doy = domy if (sum(ids(1:4)) > -4) then @@ -738,15 +1594,11 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) v1 = ptr%p%v m1 = propriedade(ptr%p%grupo)%m 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 if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then ptr%p%F = ptr%p%F + GField*m1 end if - ! print '("x1 =", f10.6, " ", f10.6, " n ", i2, " cell ",i2," ",i2)',x1(1),x1(2), ptr%p%n,i,j - ! if (id == 0) read(*,*) !calcular a força desta com todas as outras partículas next_node => list_next(node) ! próxima partícula da célula node => list_next(node) ! a ser computado com a particula selecionada @@ -757,7 +1609,6 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) 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 @@ -785,13 +1636,8 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) ! 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) - end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F end if node => list_next(node) ! próxima partícula da célula end do @@ -805,7 +1651,6 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) 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) @@ -831,13 +1676,8 @@ 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] - 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) - end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F end if node => list_next(node) ! próxima partícula da célula end do @@ -875,13 +1715,8 @@ 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] - 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) - end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F end if node => list_next(node) ! próxima partícula da célula end do @@ -922,13 +1757,8 @@ 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] - 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) - end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F end if node => list_next(node) ! próxima partícula da célula end do @@ -939,7 +1769,6 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) 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) @@ -964,13 +1793,8 @@ 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] - 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) - end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F end if node => list_next(node) ! próxima partícula da célula end do @@ -981,7 +1805,6 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) 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) @@ -1007,13 +1830,8 @@ 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] - 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) - end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F end if node => list_next(node) ! próxima partícula da célula end do @@ -1025,7 +1843,6 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) 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) @@ -1051,13 +1868,8 @@ 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] - 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) - end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F end if node => list_next(node) ! próxima partícula da célula end do @@ -1071,7 +1883,6 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) 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) @@ -1097,13 +1908,8 @@ 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] - 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) - end if - ptr%p%F = aux2 + ptr%p%F +fR - ptrn%p%F = -aux2 + ptrn%p%F - fR + ptr%p%F = aux2 + ptr%p%F + ptrn%p%F = -aux2 + ptrn%p%F end if node => list_next(node) ! próxima partícula da célula end do @@ -1115,9 +1921,9 @@ subroutine comp_F(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) !print*, 'Linha 143' ! if (id == 0) read(*,*) ! call MPI_barrier(MPI_COMM_WORLD, ierr) - end subroutine comp_F + end subroutine comp_F_thermo - subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) + subroutine comp_FT_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id, t) use linkedlist use mod1 use data @@ -1204,9 +2010,7 @@ subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id 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 + ptr%p%F = aux2 + ptr%p%F ptrn%p%F = -aux2 + ptrn%p%F end if @@ -1492,7 +2296,7 @@ subroutine comp_F_thermo(GField, mesh,malha,propriedade,r_cut,domx,domy, ids, id !print*, 'Linha 143' ! if (id == 0) read(*,*) ! call MPI_barrier(MPI_COMM_WORLD, ierr) - end subroutine comp_F_thermo + end subroutine comp_FT_thermo ! atualiza posições subroutine comp_x(icell,jcell,malha,N,mesh,propriedade, dx_max,t,dt,ids,LT,domx,domy,wall,id, np) @@ -3512,9 +4316,9 @@ subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic if (propriedade(ptr%p%grupo)%ismolecule) then ! Termostato afeta ! passo preditor - ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xih*ptr%p%v(1) - ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xih*ptr%p%v(2) - Khot = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + Khot + ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xih*ptr%p%v(1)/2 + ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xih*ptr%p%v(2)/2 + Khot = propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + Khot numph = numph + 1 else ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) @@ -3537,7 +4341,7 @@ subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic ! passo preditor ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xic*ptr%p%v(1)/2 ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xic*ptr%p%v(2)/2 - Kcold = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + Kcold + Kcold = propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + Kcold numpc = numpc + 1 else ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) @@ -3552,58 +4356,7 @@ subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic end do end do - ! do i = hot_cells(3), hot_cells(4) !(domy(1)+bsh(3)), (domy(2)+bsh(4)) - ! do j = hot_cells(1), hot_cells(2) ! (domx(1)+bsh(1)), (domx(2)+bsh(2)) - ! !if (j >= hot_cells(1) .and. j <= hot_cells(2) .and. i >= hot_cells(3) .and. i <= hot_cells(4)) then - ! ! Termostato afeta QUENTE - ! node => list_next(malha(i,j)%list) - ! do while (associated(node)) - ! ptr = transfer(list_get(node), ptr) - ! if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then - ! if (propriedade(ptr%p%grupo)%ismolecule) then - ! ! Termostato afeta - ! ! passo preditor - ! ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xih*ptr%p%v(1) - ! ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xih*ptr%p%v(2) - ! Khot = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + Khot - ! numph = numph + 1 - ! else - ! ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - ! ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - - ! end if - ! end if - ! ptr%p%F = [0,0] - ! node => list_next(node) - ! end do - ! end do - ! end do - ! do i = cold_cells(3), cold_cells(4) !(domy(1)+bsh(3)), (domy(2)+bsh(4)) - ! do j = cold_cells(1), cold_cells(2) ! (domx(1)+bsh(1)), (domx(2)+bsh(2)) - ! ! Termostato afeta FRIO - ! ! else if (j >= cold_cells(1) .and. j <= cold_cells(2) .and. i >= cold_cells(3) .and. i <= cold_cells(4)) then - ! node => list_next(malha(i,j)%list) - ! do while (associated(node)) - ! ptr = transfer(list_get(node), ptr) - ! if (propriedade(ptr%p%grupo)%x_lockdelay <= t) then - ! if (propriedade(ptr%p%grupo)%ismolecule) then - ! ! Termostato afeta - ! ! passo preditor - ! ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xic*ptr%p%v(1)/2 - ! ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - dt*xic*ptr%p%v(2)/2 - ! Kcold = 0.5*propriedade(ptr%p%grupo)%m*(ptr%p%v(1)**2+ptr%p%v(2)**2) + Kcold - ! numpc = numpc + 1 - ! else - ! ptr%p%v(1) = ptr%p%v(1) + ptr%p%F(1)*dt/(2*propriedade(ptr%p%grupo)%m) - ! ptr%p%v(2) = ptr%p%v(2) + ptr%p%F(2)*dt/(2*propriedade(ptr%p%grupo)%m) - ! end if - ! end if - ! ptr%p%F = [0,0] - ! node => list_next(node) - ! end do - ! !end if - ! end do - ! end do + if (np > 1) then aux = [Khot, numph, Kcold, numpc] @@ -3616,10 +4369,9 @@ subroutine comp_v_thermo_pred(malha,mesh,dt,t,propriedade,domx,domy,Mc, xih, xic Kcold = auxres(3)+auxres(7)+auxres(11)+auxres(15) numpc = auxres(4)+auxres(8)+auxres(12)+auxres(16) end if + ! print*, "khot,numph", khot, numph, "kcold,numpc",kcold, numpc + ! print*, "Thot", Khot/numph, "Tcold", Kcold/numpc - - ! Khot = Khot/numph - ! Kcold = Kcold/numpc ! Este M aqui é o grau de acoplamento ! 2*N ao invés de 3*N pq estamos 2D diff --git a/potential_analysis.py b/potential_analysis.py index 99c2bec..30a7a04 100644 --- a/potential_analysis.py +++ b/potential_analysis.py @@ -159,14 +159,14 @@ def einstein_condutividade(delta_e,na): if yp == mesh[1]: yp = yp - 1 particle_map[xp][yp].append( pos_vel.loc[nn,'n'] ) - density_map[xp,yp,step-a] += 1 + density_map[xp,yp,step-stepini] += 1 # "print C" for i in range(region[0],region[1]): for j in range(region[2],region[3]): for nn in range(len(particle_map[i][j])): n1 = particle_map[i][j][nn] m = mass[pos_vel.loc[n1,'tipo']] - Kmap[i,j,step-a] += (pos_vel.loc[n1,'v_x']**2+pos_vel.loc[n1,'v_y']**2)*m/(2) + Kmap[i,j,step-stepini] += (pos_vel.loc[n1,'v_x']**2+pos_vel.loc[n1,'v_y']**2)*m/(2) # Vmap[i,j,step] += pos_vel.loc[n1,'u'] if pbar: bar.update(step-stepini) diff --git a/saida.f90 b/saida.f90 index 7b15948..dec7af6 100644 --- a/saida.f90 +++ b/saida.f90 @@ -84,7 +84,7 @@ subroutine vec2csv(v,n,d,prop,step,t,nimpre,start) timep = time else - if (step > 0) then + if (step > 2) then print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", f10.3, "s" )',prop//extensao//'.'//trim(passo),t,etc else print '("Salvo arquivo ", A, " t = ", f10.3, " ETC: ", "unknown" )',prop//extensao//'.'//trim(passo),t diff --git a/settings.ini b/settings.ini index 3b534f7..08b77cf 100755 --- a/settings.ini +++ b/settings.ini @@ -6,7 +6,7 @@ nimpre = 100 # quntidade de saidas N = 4100 # número de partículas Ntype = 2 # número de tipos de partícula - t_fim = 35000 # tempo de execução + t_fim = 40000 # tempo de execução nimpre_init = 0 # começa a simular a partir no tempo referente ao snapshot (nimpre). dt = 0.001 # passo de tempo dimX = 10000 # Dimensões da região de cálculo