module modlist

  use global

  implicit none


  

contains



  !........................................!
  !     constructor for class list         !
  !........................................!
  subroutine init_list()
    
    implicit none
    integer :: i
    integer :: aver
    
    print*,param%box,param%rcut
        
    listvar%rn = (param%box/int(param%box/param%rcut))
    listvar%ncellT = nint(param%vol/(listvar%rn**3))
    listvar%ncellD = nint(listvar%ncellT**(1.0d0/3.0d0))
    listvar%alloc_size = int(ceiling(listvar%ncellT*1.2))
    
    print*,'--------------- INTIALIZING LISTS-----------------------------------'
    print*,'           rn                :',listvar%rn
    print*,'           ncellT            :',listvar%ncellT
    print*,'           ncellD            :',listvar%ncellD
    print*,'           MAKE SURE NCELLD CUBED IS NCELLT'
    print*,'--------------- END INTIALIZING LISTS-------------------------------'
    

    !--- try these arrays
    ! attributes align:64::min_x,min_y,min_z,cnum
    allocate(min_x(0:13*listvar%alloc_size),min_y(0:13*listvar%alloc_size),min_z(0:13*listvar%alloc_size))
    allocate(cnum(0:13*listvar%alloc_size))





    allocate(cl(0:listvar%alloc_size))
    allocate(cv_total(13,0:listvar%alloc_size))
    allocate(start(0:listvar%alloc_size),end(0:listvar%alloc_size))


    !--- allocate temp variables used for cell sort
    allocate(temp_num(0:listvar%alloc_size))
   

    aver = cell_max()
    aver = int(ceiling(1.4*aver))

    print*,'AVERAGRE',aver
    do i =0, listvar%alloc_size

       allocate(cl(i)%loc_r(aver))
       cl(i)%size = aver
    enddo
    


    !--- allocate atom cell construct
    allocate(atom(param%np))


    
  end subroutine init_list

  


  !......................................................................!
  !     determines the neigh neighbors for total energy calculation      !
  !......................................................................!  
  subroutine  neighcell_total()
    
    implicit none
    integer :: offsets(3,13)
    integer :: in,i,count
    integer :: ic,itel,icz,icy,icx
    integer :: ix,iy,iz
    integer :: iccx,iccy,iccz
    integer :: x,y,z

    !--- build offsets array
    offsets(1,1)=1; offsets(2,1)=0;  offsets(3,1) = 0
    offsets(1,2)=1; offsets(2,2)=1;  offsets(3,2) = 0
    offsets(1,3)=1; offsets(2,3)=-1; offsets(3,3) = 0
    offsets(1,4)=0; offsets(2,4)=-1; offsets(3,4) = 0
    
    offsets(1,5)=1; offsets(2,5)=0;  offsets(3,5) = 1
    offsets(1,6)=1; offsets(2,6)=1;  offsets(3,6) = 1
    offsets(1,7)=1; offsets(2,7)=-1; offsets(3,7) = 1
    offsets(1,8)=0; offsets(2,8)=-1; offsets(3,8) = 1
    
    offsets(1,9)=1;  offsets(2,9)=0;  offsets(3,9) = -1
    offsets(1,10)=1; offsets(2,10)=1; offsets(3,10) = -1
    offsets(1,11)=1; offsets(2,11)=-1;offsets(3,11) = -1
    offsets(1,12)=0; offsets(2,12)=-1;offsets(3,12) = -1
    
    offsets(1,13) = 0; offsets(2,13)=0; offsets(3,13) = 1
    

    count = 0
    do ic = 0,listvar%ncellT-1
             
       
       !--- obtain x,y,z coordinates of cell of interest: ic
       !--- this is doing using integer arithmetic 
       icz = ic/(listvar%ncellD*listvar%ncellD)
       icy = (ic - icz*listvar%ncellD*listvar%ncellD)/(listvar%ncellD)
       icx = ic -icy*listvar%ncellD - icz*listvar%ncellD*listvar%ncellD
       
       
       
       do i = 1,13
          
          
          ix = offsets(1,i); iy = offsets(2,i); iz = offsets(3,i)
          
          iccz = icz+iz
          
          if(iccz.lt.0) then
             z = iccz + listvar%ncellD                              
          else if (iccz.ge.listvar%ncellD) then
             z = iccz - listvar%ncellD                              
          else
             z = iccz
          endif
          
          
          iccy = icy+iy
          
          if(iccy.lt.0) then
             y = iccy + listvar%ncellD                                            
          else if (iccy.ge.listvar%ncellD) then
             y = iccy - listvar%ncellD                                                                                 
          else
             y = iccy
          endif
          
          
          iccx = icx+ix
          if(iccx.lt.0) then
             x = iccx + listvar%ncellD             
          else if (iccx.ge.listvar%ncellD) then
             x = iccx - listvar%ncellD                      
          else
             x = iccx    
          endif
          
          in = (x+y*listvar%ncellD + z*listvar%ncellD*listvar%ncellD)   
          
          cnum(count) = in

         
          cv_total(i,ic)%cnum = in
          
          !--- store min im in cv array
          if(iccx.lt.0) then  
             min_x(count) = param%box                 
          else if (iccx.ge.listvar%ncellD) then
             min_x(count) = -param%box             
          else
             min_x(count) = 0.0d0        
          endif
          
          if(iccy.lt.0) then  
             min_y(count) = param%box                 
          else if (iccy.ge.listvar%ncellD) then
             min_y(count) = -param%box             
          else
             min_y(count) = 0.0d0        
          endif
          
          if(iccz.lt.0) then  
             min_z(count) = param%box                 
          else if (iccz.ge.listvar%ncellD) then
             min_z(count) = -param%box             
          else
             min_z(count) = 0.0d0        
          endif

          count = count + 1
       enddo
     
    enddo
    
  end subroutine neighcell_total
  

  !........................................!
  ! compute cell interaction indicies      !
  ! compute miniomage for each cell pair   !
  ! compute cell coordinates for each cell
  !........................................!
  subroutine compute_cell_neighbors()
    implicit none
       
    call neighcell_total()


  end subroutine compute_cell_neighbors
          


  !........................................!
  !     determines cell number             !
  !........................................!
  integer function cell(x,y,z)
    
    implicit none
    double precision :: x,y,z
    integer :: ix,iy,iz
    
    if(x.lt.param%box)then    
       ix = int(x/listvar%rn)
    else
       ix = listvar%ncellD-1
    endif

    if(y.lt.param%box)then    
       iy = int(y/listvar%rn)
    else
       iy = listvar%ncellD-1
    endif
    
    if(z.lt.param%box)then    
       iz = int(z/listvar%rn)
    else
       iz = listvar%ncellD-1
    endif

   
    cell = ix+iy*listvar%ncellD+iz*listvar%ncellD*listvar%ncellD
    
  end function cell
  
  
 


  !.....................................................................!
  ! finds cell with largest number of atoms for initial array allocation
  !.....................................................................!
  integer function cell_max()
    implicit none
    integer :: temp_cell(param%np)
    integer :: i,icell

    temp_num(:) = 0
  
    cell_max = 0
    do i = 1 ,param%np

       icell = cell( x(i),y(i),z(i))

       temp_cell(i) = icell
       temp_num(icell) = temp_num(icell) + 1     

       if(temp_num(icell).gt.cell_max)then
          cell_max = temp_num(icell)
       endif

    enddo

  end function cell_max



  !.....................................................................................!
  ! want to assign each atom to appropriate cell
  ! generate list of atoms in each cell for GONET algorithm
  ! do not want to deallocate and reallocate her
  !......................................................................................!
  subroutine total_cell_sort_init(step)
    implicit none
    integer :: step
    integer :: i,icell
    integer :: new_size

    temp_num(:) = 0
      
    do i = 1 ,param%np

       icell = cell(x(i),y(i),z(i))

       temp_num(icell) = temp_num(icell) + 1     

       !--- keep particles current cell
       atom(i)%cell = icell 
       atom(i)%loc  = temp_num(icell)       
    enddo



    do i=0, listvar%ncellT - 1

       !--- number of molecules in cell 
       cl(i)%num = temp_num(i)
       
       !--- while sorting here, check to see if reallocation is needed
       if(cl(i)%num.gt.cl(i)%size) then
          
          new_size = int(ceiling(1.4*cl(i)%num))
          

          deallocate(cl(i)%loc_r)
          allocate(cl(i)%loc_r(new_size))
          cl(i)%size = new_size
                   
       endif
    enddo

   

    do i= 1,param%np
       !--- get cell for particle
       icell = atom(i)%cell
        
       !--- place particles in cell
       cl(icell)%loc_r(atom(i)%loc)%x = x(i)   
       cl(icell)%loc_r(atom(i)%loc)%y = y(i)     
       cl(icell)%loc_r(atom(i)%loc)%z = z(i)     
  
    enddo


  end subroutine total_cell_sort_init

  subroutine build_r()
    integer :: i,j,counter
    integer :: counter2
    counter = 0 
    counter2 = 0
    do i = 0,listvar%ncellT-1  
 
       start(i) = counter + 1

       do j = 1 ,cl(i)%num
          
          counter = counter + 1
          
          x(counter)    = cl(i)%loc_r(j)%x
          y(counter)    = cl(i)%loc_r(j)%y
          z(counter)    = cl(i)%loc_r(j)%z

          counter2 = counter2+1
          flat(counter2)   = x(counter)

          counter2 = counter2+1
          flat(counter2) = y(counter)

          counter2 = counter2+1
          flat(counter2) = z(counter)

       enddo
       end(i) = counter
    enddo
  end subroutine build_r

  !----------------------------------------------------------------------------------
  !*** subroutine vlist()
  
  subroutine verlet_build(step)
    
    implicit none
    double precision, save :: time
    double precision :: dx,dy,dz
    double precision :: x1,y1,z1,x2,y2,z2
    double precision :: minx, miny, minz
    double precision :: dr,dri,dr2,dr2i,dr6i,dr12i,energy
    integer :: i,j,k,l,step
    integer :: cs,ce,cnews,cnewe
    integer :: T1, T2, clock_rate, clock_max
    integer :: tid,count,delta
    integer :: c1s,c1e,c2s,c2e,cell_neigh
  

    allocate(numN(param%np))
    allocate(verletList(500,param%np))

    numN(:) = 0
    
    call system_clock(T1,clock_rate,clock_max)

    do k=0,listvar%ncellT-1
       

       cs   = start(k)
       ce   = end(k)
      
       
       do i= cs, ce-1  
          !---load particle coodinates
          x1 = x(i); y1 = y(i); z1 =z(i)
        
          do j= i+1,ce
             
             !---load particle coodinates
             x2 = x(j); y2 = y(j); z2 = z(j)
            
             !--- calculate displacements
             dx = x2-x1 
             dy = y2-y1
             dz = z2-z1
             
             dr2 = dx*dx+dy*dy+dz*dz
             
             if(dr2.lt.param%rcut2) then
                
                numN(i) = numN(i) + 1
                verletList(numN(i),i) = j
             
             endif
             
          enddo
       enddo
    enddo
    
    do i = 0,listvar%ncellT-1
       
       
       c1s   = start(i)
       c1e   = end(i)
       
       !--- loop over cell neighbors
       do k =  i*13,13*i+12    
          
          !---determine this cell
          cell_neigh = cnum(k)
          
          !---minimum image
          minx = min_x(k)
          miny = min_y(k)
          minz = min_z(k)
          
          
          c2s = start(cell_neigh)
          c2e = end(cell_neigh)  
          
          
          !--- loop over particles in each cell         
          do j = c1s,c1e   
             x1 = x(j); y1 = y(j); z1 = z(j)
             
             do l = c2s,c2e
                x2 = x(l); y2 = y(l); z2 = z(l)
                
                !--- obtain displacements
                !--- apply minimum image here
                dx = x2-x1-minx
                dy = y2-y1-miny
                dz = z2-z1-minz
                
                dr2 = dx*dx+dy*dy+dz*dz
                
                tid = 0.0d0
                if(dr2.lt.param%rcut2) then
                   
                   numN(j) = numN(j) + 1
                   verletList(numN(j),j) = l
                   
                                      
                endif
             enddo
          enddo
       enddo
    enddo

    call system_clock(T2,clock_rate,clock_max)

    print*,'TOTAL TIME NO VECTORIZATION:',real(T2-T1)/real(clock_rate)
    

  end subroutine verlet_build

     
end module modlist
