--- trunk/OOPSE/libmdtools/simulation_module.F90	2003/03/21 17:42:12	378
+++ trunk/OOPSE/libmdtools/simulation_module.F90	2003/07/01 21:33:45	569
@@ -6,7 +6,6 @@ module simulation
   use force_globals
   use vector_class
   use atype_module
-  use lj
 #ifdef IS_MPI
   use mpiSimulation
 #endif
@@ -21,23 +20,22 @@ module simulation
 
   logical, save :: simulation_setup_complete = .false.
 
-  integer, public, save :: natoms
+  integer, public, save :: nLocal, nGlobal
   integer, public, save :: nExcludes_Global = 0
   integer, public, save :: nExcludes_Local = 0
   integer, allocatable, dimension(:,:), public :: excludesLocal
   integer, allocatable, dimension(:), public :: excludesGlobal
+  integer, allocatable, dimension(:), public :: molMembershipList
 
   real(kind=dp), save :: rcut2 = 0.0_DP
   real(kind=dp), save :: rcut6 = 0.0_DP
   real(kind=dp), save :: rlist2 = 0.0_DP
-  real(kind=dp), public, dimension(3), save :: box
-
+  real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
+  logical, save :: boxIsOrthorhombic
   
   public :: SimulationSetup
   public :: getNlocal
   public :: setBox
-  public :: setBox_3d
-  public :: getBox
   public :: setRcut
   public :: getRcut
   public :: getRlist
@@ -54,31 +52,23 @@ module simulation
   public :: SimRequiresPrepairCalc
   public :: SimRequiresPostpairCalc
   public :: SimUsesDirectionalAtoms
-
-  interface getBox
-     module procedure getBox_3d
-     module procedure getBox_1d
-  end interface
   
-  interface setBox
-     module procedure setBox_3d
-     module procedure setBox_1d
-  end interface
-  
 contains
   
-  subroutine SimulationSetup(setThisSim, nComponents, c_idents, &
+  subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
        CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
+       CmolMembership, &
        status)    
 
     type (simtype) :: setThisSim
-    integer, intent(inout) :: nComponents
-    integer, dimension(nComponents),intent(inout) :: c_idents
+    integer, intent(inout) :: CnGlobal, CnLocal
+    integer, dimension(CnLocal),intent(inout) :: c_idents
 
     integer :: CnLocalExcludes
     integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
     integer :: CnGlobalExcludes
     integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
+    integer, dimension(CnGlobal),intent(in) :: CmolMembership 
     !!  Result status, success = 0, status = -1
     integer, intent(out) :: status
     integer :: i, me, thisStat, alloc_stat, myNode
@@ -93,8 +83,11 @@ contains
     status = 0
 
     ! copy C struct into fortran type
+
+    nLocal = CnLocal
+    nGlobal = CnGlobal
+
     thisSim = setThisSim
-    natoms = nComponents
     rcut2 = thisSim%rcut * thisSim%rcut
     rcut6 = rcut2 * rcut2 * rcut2
     rlist2 = thisSim%rlist * thisSim%rlist
@@ -103,14 +96,16 @@ contains
     nExcludes_Global = CnGlobalExcludes
     nExcludes_Local = CnLocalExcludes
 
-    call InitializeForceGlobals(natoms, thisStat)
+    call InitializeForceGlobals(nLocal, thisStat)
     if (thisStat /= 0) then
+       write(default_error,*) "SimSetup: InitializeForceGlobals error"
        status = -1
        return
     endif
 
     call InitializeSimGlobals(thisStat)
     if (thisStat /= 0) then
+       write(default_error,*) "SimSetup: InitializeSimGlobals error"
        status = -1
        return
     endif
@@ -160,7 +155,7 @@ contains
     endif
     
 #else
-    do i = 1, nComponents
+    do i = 1, nLocal
        
        me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
        atid(i) = me
@@ -168,13 +163,8 @@ contains
     enddo
 #endif
 
-    !! Create neighbor lists
-    call expandNeighborList(nComponents, thisStat)
-    if (thisStat /= 0) then
-       status = -1
-       return
-    endif
 
+
     do i = 1, nExcludes_Local
        excludesLocal(1,i) = CexcludesLocal(1,i)
        excludesLocal(2,i) = CexcludesLocal(2,i)
@@ -183,64 +173,40 @@ contains
     do i = 1, nExcludes_Global
        excludesGlobal(i) = CexcludesGlobal(i)
     enddo
-    
+
+    do i = 1, nGlobal
+       molMemberShipList(i) = CmolMembership(i)
+     enddo
+
     if (status == 0) simulation_setup_complete = .true.
     
   end subroutine SimulationSetup
   
-  subroutine setBox_3d(new_box_size)
-    real(kind=dp), dimension(3) :: new_box_size
+  subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
+    real(kind=dp), dimension(3,3) :: cHmat, cHamtInv
+    integer :: cBoxIsOrthorhombic
     integer :: smallest, status, i
 
-    thisSim%box = new_box_size
-    box = thisSim%box
+    Hmat = cHmat
+    HmatInv = cHmatInv
+    if(cBoxisOrthorhombic .eq. 0 )then
+       boxIsOrthorhombic = .false.
+    else boxIsOrthorhombic = .true.
+    endif
+    
 
-    smallest = 1
-    do i = 2, 3
-       if (new_box_size(i) .lt. new_box_size(smallest)) smallest = i
-    end do
-    if (thisSim%rcut .gt. 0.5_dp * new_box_size(smallest)) &
-         call setRcut(0.5_dp * new_box_size(smallest), status)
     return     
-  end subroutine setBox_3d
+  end subroutine setBox
 
-  subroutine setBox_1d(dim, new_box_size)
-    integer :: dim, status
-    real(kind=dp) :: new_box_size
-    thisSim%box(dim) = new_box_size
-    box(dim) = thisSim%box(dim)
-    if (thisSim%rcut .gt. 0.5_dp * new_box_size) &
-         call setRcut(0.5_dp * new_box_size, status)
-  end subroutine setBox_1d
-
   subroutine setRcut(new_rcut, status) 
     real(kind = dp) :: new_rcut
     integer :: myStatus, status
     thisSim%rcut = new_rcut
     rcut2 = thisSim%rcut * thisSim%rcut
     rcut6 = rcut2 * rcut2 * rcut2
-    myStatus = 0
-    call LJ_new_rcut(new_rcut, myStatus)
-    if (myStatus .ne. 0) then
-       write(default_error, *) 'LJ module refused our rcut!'
-       status = -1
-       return
-    endif
     status = 0
     return
   end subroutine setRcut
-
-  function getBox_3d() result(thisBox)
-    real( kind = dp ), dimension(3) :: thisBox
-    thisBox = thisSim%box
-  end function getBox_3d
-  
-  function getBox_1d(dim) result(thisBox)
-    integer, intent(in) :: dim
-    real( kind = dp ) :: thisBox
-    
-    thisBox = thisSim%box(dim)
-  end function getBox_1d
      
   subroutine getRcut(thisrcut,rc2,rc6,status)
     real( kind = dp ), intent(out) :: thisrcut
@@ -280,6 +246,7 @@ contains
   function getRrf() result(rrf)
     real( kind = dp ) :: rrf
     rrf = thisSim%rrf
+    write(*,*) 'getRrf = ', rrf, thisSim%rrf
   end function getRrf
   
   function getRt() result(rt)
@@ -362,21 +329,28 @@ contains
        thisStat = -1
        return
     endif
+
+    allocate(molMembershipList(nGlobal), stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       thisStat = -1
+       return
+    endif
     
   end subroutine InitializeSimGlobals
   
   subroutine FreeSimGlobals()
     
     !We free in the opposite order in which we allocate in.
-    
+
+    if (allocated(molMembershipList)) deallocate(molMembershipList)     
     if (allocated(excludesGlobal)) deallocate(excludesGlobal)
     if (allocated(excludesLocal)) deallocate(excludesLocal)
 
   end subroutine FreeSimGlobals
 
-  pure function getNlocal() result(nlocal)
-    integer :: nlocal
-    nlocal = natoms
+  pure function getNlocal() result(n)
+    integer :: n
+    n = nLocal
   end function getNlocal