--- trunk/mdtools/md_code/simulation_module.F90	2002/12/04 21:19:38	194
+++ trunk/mdtools/md_code/simulation_module.F90	2003/02/24 21:25:15	281
@@ -1,74 +1,287 @@
 module simulation
   use definitions, ONLY :dp
-  use force_wrappers, ONLY : alloc_force_wrappers
+#ifdef IS_MPI
+  use mpiSimulation
+#endif
 
   implicit none
   PRIVATE
 
+#define __FORTRAN90
+#include "../headers/fsimulation.h"
 
-  
-  real ( kind = dp ), public, dimension(3) :: box
+  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
 
-  real ( kind = dp ), public :: rlist
-  real ( kind = dp ), public :: rcut
-  real ( kind = dp ), public :: rlistsq
-  real ( kind = dp ), public :: rcutsq
+!! WARNING: use_pbc hardcoded, fixme
+   logical :: setSim = .false.
 
-  integer,public, allocatable, dimension(:)   :: point
-  integer,public, allocatable, dimension(:)   :: list
+!! array for saving previous positions for neighbor lists.  
+  real( kind = dp ), allocatable,dimension(:,:),save :: q0
 
 
-!MPI dependent routines
+  public :: check
+  public :: save_nlist
+  public :: wrap
+  public :: getBox
+  public :: getRcut
+  public :: getRlist
+  public :: getNlocal
+  public :: setSimulation
+  public :: isEnsemble
+  public :: isPBC
+  public :: getStringLen
+  public :: returnMixingRules
 
-#ifdef IS_MPI
-! Universal routines: All types of force calculations will need these arrays
-! Arrays specific to a type of force calculation should be declared in that module.
-  real( kind = dp ), allocatable, dimension(:,:) :: qRow
-  real( kind = dp ), allocatable, dimension(:,:) :: qColumn
+!  public :: setRcut
 
-  real( kind = dp ), allocatable, dimension(:,:) :: fRow
-  real( kind = dp ), allocatable, dimension(:,:) :: fColumn
+  interface wrap
+     module procedure wrap_1d
+     module procedure wrap_3d
+  end interface
 
-  real( kind = dp ), allocatable, dimension(:,:) :: tRow
-  real( kind = dp ), allocatable, dimension(:,:) :: tColumn
+  interface getBox
+     module procedure getBox_3d
+     module procedure getBox_dim
+  end interface
 
-#endif
 
 
 
 
 
 
+
+
+
 contains
 
-#ifdef MPI
-! Allocated work arrays for MPI
-  subroutine allocate_mpi_arrays(nDimensions,numComponents)
-    integer, intent(in) :: nDimensions
-    integer, intent(in) :: numComponents
+  subroutine setSimulation(nLRParticles,box,rlist,rcut,ensemble,mixingRule,use_pbc)
+    integer, intent(in) :: nLRParticles 
+    real(kind = dp ), intent(in), dimension(3) :: box
+    real(kind = dp ), intent(in) :: rlist
+    real(kind = dp ), intent(in) :: rcut
+    character( len = stringLen), intent(in)  :: ensemble
+    character( len = stringLen), intent(in)  :: mixingRule
+    logical, intent(in) :: use_pbc
+    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
+    
+    thisSim%ensemble = ensemble
+    thisSim%mixingRule = mixingRule
+    thisSim%use_pbc = use_pbc
 
+    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
 
 
-  end subroutine allocate_mpi_arrays
+  subroutine change_box_size(new_box_size)
+    real(kind=dp), dimension(3) :: 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 - thisSim%rcut
+    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
+    
 
-  subroutine set_simulation(box,rlist,rcut)
+    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
 
-  end subroutine set_simulation
+    
+    if (this_sim%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
 
-  subroutine change_box_size(new_box_size)
-    real(kind=dp), dimension(3) :: new_box_size
+    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
 
-    box = new_box_size
+  end subroutine getRcut
+  
+  
+  
+ 
+  subroutine getRlist(thisrlist,rlist2,status)
+    real( kind = dp ), intent(out) :: thisrlist
+    real( kind = dp ), intent(out), optional :: rlist2
 
-  end subroutine change_box_size
+    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
+
+
+  function isEnsemble(this_ensemble) result(is_this_ensemble)
+    character(len = *) :: this_ensemble
+    logical :: is_this_enemble
+    is_this_ensemble = .false.
+    if (this_ensemble == thisSim%ensemble) is_this_ensemble = .true.
+  end function isEnsemble
+
+  function returnEnsemble() result(thisEnsemble)
+    character (len = len(thisSim%ensemble)) :: thisEnsemble
+    thisEnsemble = thisSim%ensemble
+  end function returnEnsemble
+
+  function returnMixingRules() result(thisMixingRule)
+    character (len = len(thisSim%ensemble)) :: thisMixingRule
+    thisMixingRule = thisSim%MixingRule
+  end function returnMixingRules
+
+  function isPBC() result(PBCset)
+    logical :: PBCset
+    PBCset = .false.
+    if (thisSim%use_pbc) PBCset = .true.
+  end function isPBC
+
+  pure function getStringLen() result (thislen)
+    integer :: thislen    
+    thislen = string_len
+  end function setStringLen
+
 end module simulation