--- trunk/mdtools/md_code/simulation_module.F90	2002/11/22 20:39:33	185
+++ trunk/mdtools/md_code/simulation_module.F90	2003/01/30 15:20:21	253
@@ -1,35 +1,268 @@
 module simulation
   use definitions, ONLY :dp
+#ifdef IS_MPI
+  use mpiSimulation
+#endif
+
   implicit none
   PRIVATE
 
 
 
-  real ( kind = dp ), dimension(3) :: box
+  type, public :: simtype
+     PRIVATE
+!     SEQUENCE
+!! Number of particles on this processor
+     integer :: nLRparticles
+!! Periodic Box     
+     real ( kind = dp ), dimension(3) :: box
+!! List Cutoff     
+     real ( kind = dp ) :: rlist = 0.0_dp
+!! Radial cutoff
+     real ( kind = dp ) :: rcut  = 0.0_dp
+!! List cutoff squared
+     real ( kind = dp ) :: rlistsq = 0.0_dp
+!! Radial Cutoff squared
+     real ( kind = dp ) :: rcutsq  = 0.0_dp
+!! Radial Cutoff^6
+     real ( kind = dp ) :: rcut6  = 0.0_dp
 
+  end type simtype
 
+  type (simtype), public :: thisSim
+!! Tag for MPI calculations  
+  integer, allocatable, dimension(:) :: tag
 
+#ifdef IS_MPI
+  integer, allocatable, dimension(:) :: tag_row
+  integer, allocatable, dimension(:) :: tag_column
+#endif
 
+!! WARNING: use_pbc hardcoded, fixme
+  logical :: use_pbc = .true.
+  logical :: setSim = .false.
 
+!! array for saving previous positions for neighbor lists.  
+  real( kind = dp ), allocatable,dimension(:,:),save :: q0
 
 
-contains
+  public :: check
+  public :: save_nlist
+  public :: wrap
+  public :: getBox
+  public :: getRcut
+  public :: getRlist
+  public :: getNlocal
+  public :: setSimulation
+!  public :: setRcut
 
+  interface wrap
+     module procedure wrap_1d
+     module procedure wrap_3d
+  end interface
 
-  subroutine set_simulation(box,rlist,rcut)
-    
+  interface getBox
+     module procedure getBox_3d
+     module procedure getBox_dim
+  end interface
 
 
-  end subroutine set_simulation
 
 
 
+
+
+
+
+
+contains
+
+  subroutine setSimulation(nLRParticles,box,rlist,rcut)
+    integer, intent(in) :: nLRParticles 
+    real(kind = dp ), intent(in), dimension(3) :: box
+    real(kind = dp ), intent(in) :: rlist
+    real(kind = dp ), intent(in) :: rcut
+    integer :: alloc_stat
+    if( setsim ) return  ! simulation is already initialized
+    setSim = .true.
+
+    thisSim%nLRParticles = nLRParticles
+    thisSim%box          = box
+    thisSim%rlist        = rlist
+    thisSIm%rlistsq      = rlist * rlist
+    thisSim%rcut         = rcut
+    thisSim%rcutsq       = rcut * rcut
+    thisSim%rcut6        = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
+    
+    if (.not. allocated(q0)) then
+       allocate(q0(3,nLRParticles),stat=alloc_stat)
+    endif
+  end subroutine setSimulation
+
+  function getNparticles() result(nparticles)
+    integer :: nparticles
+    nparticles = thisSim%nLRparticles
+  end function getNparticles
+
+
   subroutine change_box_size(new_box_size)
     real(kind=dp), dimension(3) :: new_box_size
 
-    box = new_box_size
+    thisSim%box = new_box_size
 
   end subroutine change_box_size
 
 
+  function getBox_3d() result(thisBox)
+    real( kind = dp ), dimension(3) :: thisBox
+    thisBox = thisSim%box
+  end function getBox_3d
+
+  function getBox_dim(dim) result(thisBox)
+    integer, intent(in) :: dim
+    real( kind = dp ) :: thisBox
+
+    thisBox = thisSim%box(dim)
+  end function getBox_dim
+    
+  subroutine check(q,update_nlist)
+    real( kind = dp ), dimension(:,:)  :: q 
+    integer :: i
+    real( kind = DP ) :: dispmx 
+    logical, intent(out) :: update_nlist
+    real( kind = DP ) :: dispmx_tmp
+    real( kind = dp ) :: skin_thickness
+    integer :: natoms
+
+    natoms = thisSim%nLRparticles
+    skin_thickness = thisSim%rlist
+    dispmx = 0.0E0_DP
+    !! calculate the largest displacement of any atom in any direction
+    
+#ifdef MPI 
+    dispmx_tmp = 0.0E0_DP
+    do i = 1, thisSim%nLRparticles
+       dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
+       dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
+       dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
+    end do
+    call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, & 
+       mpi_max,mpi_comm_world,mpi_err)
+#else
+
+    do i = 1, thisSim%nLRparticles
+       dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
+       dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
+       dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
+    end do
+#endif
+  
+    !! a conservative test of list skin crossings
+    dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
+    
+    update_nlist = (dispmx.gt.(skin_thickness))
+    
+  end subroutine check
+  
+  subroutine save_nlist(q)
+    real(kind = dp ), dimension(:,:), intent(in)  :: q
+    integer :: list_size
+    
+
+    list_size = size(q)
+
+    if (.not. allocated(q0)) then
+       allocate(q0(3,list_size))
+    else if( list_size > size(q0)) then
+       deallocate(q0)
+       allocate(q0(3,list_size))
+    endif
+
+    q0 = q
+
+  end subroutine save_nlist
+  
+
+  function wrap_1d(r,dim) result(this_wrap)
+    
+    
+    real( kind = DP ) :: r
+    real( kind = DP ) :: this_wrap
+    integer           :: dim
+    
+    if (use_pbc) then
+       !     this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP) 
+       this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
+    else
+       this_wrap = r
+    endif
+    
+    return
+  end function wrap_1d
+
+  function wrap_3d(r) result(this_wrap)
+    real( kind = dp ), dimension(3), intent(in) :: r
+    real( kind = dp ), dimension(3) :: this_wrap
+
+    
+    if (use_pbc) then
+       !     this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP) 
+       this_wrap = r - thisSim%box*nint(r/thisSim%box)
+    else
+       this_wrap = r
+    endif
+  end function wrap_3d
+
+  
+ 
+  subroutine getRcut(thisrcut,rcut2,rcut6,status)
+    real( kind = dp ), intent(out) :: thisrcut
+    real( kind = dp ), intent(out), optional :: rcut2
+    real( kind = dp ), intent(out), optional :: rcut6
+    integer, optional :: status
+
+    if (present(status)) status = 0
+
+    if (.not.setSim ) then
+       if (present(status)) status = -1
+       return
+    end if
+    
+    thisrcut = thisSim%rcut
+    if(present(rcut2)) rcut2 = thisSim%rcutsq
+    if(present(rcut6)) rcut6 = thisSim%rcut6
+
+  end subroutine getRcut
+  
+  
+  
+ 
+  subroutine getRlist(thisrlist,rlist2,status)
+    real( kind = dp ), intent(out) :: thisrlist
+    real( kind = dp ), intent(out), optional :: rlist2
+
+    integer, optional :: status
+
+    if (present(status)) status = 0
+
+    if (.not.setSim ) then
+       if (present(status)) status = -1
+       return
+    end if
+    
+    thisrlist = thisSim%rlist
+    if(present(rlist2)) rlist2 = thisSim%rlistsq
+
+
+  end subroutine getRlist
+  
+  
+
+ pure function getNlocal() result(nlocal)
+    integer :: nlocal
+    nlocal = thisSim%nLRparticles
+  end function getNlocal
+
+
+
 end module simulation