--- trunk/mdtools/md_code/simulation_module.F90	2003/01/23 16:19:59	241
+++ trunk/mdtools/md_code/simulation_module.F90	2003/01/30 15:20:21	253
@@ -11,7 +11,7 @@ module simulation
 
   type, public :: simtype
      PRIVATE
-     SEQUENCE
+!     SEQUENCE
 !! Number of particles on this processor
      integer :: nLRparticles
 !! Periodic Box     
@@ -38,17 +38,23 @@ module simulation
   integer, allocatable, dimension(:) :: tag_column
 #endif
 
-
+!! WARNING: use_pbc hardcoded, fixme
+  logical :: use_pbc = .true.
   logical :: setSim = .false.
-  
-  real( kind = dp ), allocatable(:,:),save :: q0
 
+!! array for saving previous positions for neighbor lists.  
+  real( kind = dp ), allocatable,dimension(:,:),save :: q0
+
+
   public :: check
   public :: save_nlist
   public :: wrap
   public :: getBox
   public :: getRcut
-  public :: setRcut
+  public :: getRlist
+  public :: getNlocal
+  public :: setSimulation
+!  public :: setRcut
 
   interface wrap
      module procedure wrap_1d
@@ -76,17 +82,21 @@ contains
     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        = rcutsq * rcutsq * rcutsq
-
+    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)
@@ -103,7 +113,7 @@ contains
   end subroutine change_box_size
 
 
-  elemental function getBox_3d() result(thisBox)
+  function getBox_3d() result(thisBox)
     real( kind = dp ), dimension(3) :: thisBox
     thisBox = thisSim%box
   end function getBox_3d
@@ -115,20 +125,23 @@ contains
     thisBox = thisSim%box(dim)
   end function getBox_dim
     
-  subroutine check(update_nlist)
-  
+  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, nlocal
+    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 )
@@ -136,7 +149,8 @@ contains
     call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, & 
        mpi_max,mpi_comm_world,mpi_err)
 #else
-    do i = 1, natoms
+
+    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 )
@@ -159,7 +173,7 @@ contains
 
     if (.not. allocated(q0)) then
        allocate(q0(3,list_size))
-    else if( list_size > size(q0))
+    else if( list_size > size(q0)) then
        deallocate(q0)
        allocate(q0(3,list_size))
     endif
@@ -186,25 +200,25 @@ contains
     return
   end function wrap_1d
 
-  elemental function wrap_3d(r) result(this_wrap)
+  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(1:3) = r(1:3) - box(1:3)*nint(r(1:3)/box(1:3))
+       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 :: thisrcut6
+    real( kind = dp ), intent(out), optional :: rcut6
     integer, optional :: status
 
     if (present(status)) status = 0
@@ -214,13 +228,41 @@ contains
        return
     end if
     
-    thisrcut = rcut
-    if(present(rcut2)) rcut2 = rcutsq
-    if(present(rcut2)) rcut6 = rcut_6
+    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