--- trunk/mdtools/md_code/simulation_module.F90	2003/01/27 18:28:11	247
+++ trunk/mdtools/md_code/simulation_module.F90	2003/01/30 15:20:21	253
@@ -53,6 +53,7 @@ module simulation
   public :: getRcut
   public :: getRlist
   public :: getNlocal
+  public :: setSimulation
 !  public :: setRcut
 
   interface wrap
@@ -81,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        = 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)
@@ -136,7 +141,7 @@ contains
     
 #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 )
@@ -144,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 )
@@ -224,7 +230,7 @@ contains
     
     thisrcut = thisSim%rcut
     if(present(rcut2)) rcut2 = thisSim%rcutsq
-    if(present(rcut2)) rcut6 = thisSim%rcut6
+    if(present(rcut6)) rcut6 = thisSim%rcut6
 
   end subroutine getRcut
   
@@ -247,6 +253,7 @@ contains
     thisrlist = thisSim%rlist
     if(present(rlist2)) rlist2 = thisSim%rlistsq
 
+
   end subroutine getRlist