Append an array into an array

Help
1

I need to keep track of all pairs of particles that have the distance less then a certain value, i want the pairs to be allocated in an array but i don’t know how to append, for example, if the distance between the particles 1 and 8 satisfies the condition, my array id_pairs_collide recive the pair (1,8) so id_pairs_collide(1) = (1,8) and keep doing it for all pair of particles so my array shoul be something like id_pairs_collide = ((1,2), (2,7),(6,11)) for example.
But the only way i managed to make it was in a 1D array, so my it is something like id_pairs_collide = (1,8,2,7,6,11). The code is a simullation of particles in a gas, here is what i have for now:

program collision
    !se rf_mod, only: rf
    implicit none
    integer :: n_particles = 16, i, j
    real :: radius=0.06
    real(8), allocatable :: r(:,:), v(:,:),x(:),y(:), x_pairs(:,:), y_pairs(:,:), ids(:), dx_pairs(:), dy_pairs(:), cos, d_pairs(:)
    integer, allocatable ::  id_pairs(:,:), id_pairs_collide(:)
    
    !GENARATING THE PARTICLES
    allocate(r(2,n_particles)) !position of the particles
    allocate(v(2,n_particles)) !vellocities of the particles
    allocate(ids(n_particles)) 

    do i=1,n_particles !Vector "numerating" each particle
        ids(i) = i
    end do
    
    call random_number(r)
    !do i=1, 2 !randomly positioning each particle
        !do j=1, n_particles
            r(i,j) = rf()
        !end do
    !end do
    
    !open(40, file = 'pos.txt', status='unknown')        
        do j =1, n_particles
                if(r(1,j) > 0.5) then
                    write(40,1)r(:,j)!, 'red'
                    1 format(2f20.6)
                else
                    write(40,2)r(:,j)!, 'green'
                    2 format(2f20.6)
                end if
        end do
    !close(40)
    !call system ('gnuplot -p grafico.plt')
    
    !initial vellocities
    !if the particle is in the left corner, v=100.
    !if it's in the right, v=-100
    v(2,:) = 0
    do i=1, n_particles
        if(r(1,i) < 0.5) then
            v(1,i) = 100
        else 
            v(1,i) = -100
        end if
    end do 
    
    !DISTANCE BETWEEN ALL THE PAIR OF PARTICLES
    
    !all pair of particles
    id_pairs = combination(ids)

    allocate(x(n_particles))
    allocate(y(n_particles))
    x = r(1,:) !x coordinates
    y = r(2,:) !y coordinates
    x_pairs = combination(x) 
    y_pairs = combination(y)
    
    dx_pairs = diff(x_pairs) !distance between x coordinates of every particle
    dy_pairs = diff(y_pairs) !distance between y coordinates of every particle

    d_pairs = sqrt(dx_pairs**2 + dy_pairs**2) !total distance between all pairs of particles
    
    do i=1, size(d_pairs)
        if(d_pairs(i) < 2*radius) then 
            print *,d_pairs(i), id_pairs(:,i)
            id_pairs_collide = [id_pairs_collide, id_pairs(:,i)]
        end if
    enddo

    
    !Now, we need to calculate the velocities of each particle, iterate, compare distances and change the vellocities


    !open(40, file = 'comb.txt', status='unknown')        
            do i=1,size(d_pairs)
                    write(40,*)d_pairs(i)
            enddo
    !close(40)
    !open(40, file = 'comb2.txt', status='unknown')        
    !do i=1,size(x_pairs)/2
            write(40,*)x_pairs(:,i)
    !enddo

contains
function combination(arg) result(comb)
    implicit none
    real(8), dimension(:,:), allocatable :: comb
    real(8), intent(in), dimension(:) :: arg
    integer :: M, N, counter, i, j

    N=size(arg)
    M=int(N*(N-1)/2)
    allocate(comb(2,M))
    counter=1
    do i=1,N
        do j=i+1,N
            comb(1,counter)=arg(i)
            comb(2,counter)=arg(j)
            counter = counter + 1
        enddo
    enddo
    
end function combination

function diff(v)
    implicit none
    real(8), intent(in), dimension(:,:) :: v
    real(8), allocatable:: diff(:), v1(:), v2(:)
    integer :: n 

    n = size(v)/2
    allocate(v1(n))
    allocate(v2(n))
    v1 = v(1,:)
    v2 = v(2,:)
    diff = abs(v1 - v2)
end function diff

end program
1 Like
2

Will making it objective-oriented be useful?
Like

 type :: pair
    integer :: i ! the label of particle i
    integer :: j ! the label of particle j
    real :: r ! the distance btw particle i and j
 end type

Then, say you have n particles, so you have n(n-1)/2 different pairs, you create an array of type pair with n(n-1)/2 elements.

type(pair) :: pair_array(  n*(n-1)/2  )

So each elements means a pair of of particle ij, and the index of the array element is the label of this pair.
You can calculate the each pair’s distance r, say the Kth pair, the r is stored as

 pair_array(k)%r

and you know the corresponding label of the two particles are

pair_array(k)%i
pair_array(k)%j

Sorry if there are some sloppy grammar errors, but roughly that is that, you can do whatever you want. I am just saying that perhaps objective-oriented thinking will be useful in this case.

PS.
Correct me if I am wrong, I see

do i=1,N
    do j=i+1,N
...

in your code, perhaps the first line should be

do i=1,N-1
3

I don’t have an answer to your specific question but I have general suggestions. Before writing more code, you would want to read about using either Cell Lists or a Voronoi library to accelerate the neighbor detection. Computer Simulation of Liquids by Allen and Tildesley is a good resource. They talk about neighbor lists and other patterns that you might find beneficial for particle simulations.

2 Likes
4

You can grow an allocatable array of derived types, for example

module m
implicit none
type :: point
   integer :: i, j
end type point
end module m
!
program main
use m, only: point
implicit none
type(point), allocatable :: v(:)
v = [point(1,2)]
v = [v,point(2,7)]
v = [v,point(6,11)]
print "(2(1x,i0))", v
end program main
! output:
!  1 2
!  2 7
!  6 11

although it will be faster to use move_alloc to grow the allocatable array.

1 Like
5

If the objective is just to append an array to another array, move_alloc could be useful (as already mentioned by @Beliavsky). On the other hand, if the objective is to evaluate the dynamics of neighboring particles (pairs of particles), a model linked-lists cell could be useful (as already mentioned by @Abhilash).

This code is a fragment of likend-cell list model in a stochastic simulation. Although the code is a bit old (2015) it still works :slight_smile: . See melquiades/m_precells.f90 at main · aslozada/melquiades · GitHub


x%m_celli(1) = x%m_cell(1)/x%m_edge(1)
  x%m_celli(2) = x%m_cell(2)/x%m_edge(2)
  x%m_celli(3) = x%m_cell(3)/x%m_edge(3)

  write(*,'("---------------------------------------------------------")')
  write(*,'("          Linked-Cells method specfications")')
  write(ival,'(i8)') y%m_ncellc
  write(*,'(" Total number of cells          : ",a8)') adjustl(ival)
  write(*,'(" Size of subcells               : ",3f8.1," [Å]")') x%m_edge(:)/x%m_cell(:)
  write(*,'("---------------------------------------------------------")')

 else
   x%m_ndiv(1) = int(x%m_edge(1))
   x%m_ndiv(2) = int(x%m_edge(2))
   x%m_ndiv(3) = int(x%m_edge(3))

   do while((x%m_edge(1)/x%m_ndiv(1)) < y%m_cutoff .and. (x%m_edge(2)/x%m_ndiv(2)) < y%m_cutoff&
           & .and. (x%m_edge(3)/x%m_ndiv(3)) < y%m_cutoff)

   x%m_ndiv(1) = x%m_ndiv(1) - 1
   x%m_ndiv(2) = x%m_ndiv(2) - 1
   x%m_ndiv(3) = x%m_ndiv(3) - 1

   x%m_cell(1) = x%m_ndiv(1)
   x%m_cell(2) = x%m_ndiv(2)
   x%m_cell(3) = x%m_ndiv(3)

   end do
2 Likes