--- trunk/OOPSE/libmdtools/simulation_module.F90	2003/04/11 18:46:37	491
+++ trunk/OOPSE/libmdtools/simulation_module.F90	2004/06/01 21:45:22	1217
@@ -6,6 +6,7 @@ module simulation
   use force_globals
   use vector_class
   use atype_module
+  use switcheroo
 #ifdef IS_MPI
   use mpiSimulation
 #endif
@@ -15,37 +16,41 @@ module simulation
 
 #define __FORTRAN90
 #include "fSimulation.h"
+#include "fSwitchingFunction.h"
 
-  type (simtype), public :: thisSim
+  type (simtype), public, save :: thisSim
 
   logical, save :: simulation_setup_complete = .false.
 
   integer, public, save :: nLocal, nGlobal
+  integer, public, save :: nGroups, nGroupGlobal
   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
+  integer, allocatable, dimension(:),   public :: excludesGlobal
+  integer, allocatable, dimension(:),   public :: molMembershipList
+  integer, allocatable, dimension(:),   public :: groupListRow
+  integer, allocatable, dimension(:),   public :: groupStartRow
+  integer, allocatable, dimension(:),   public :: groupListCol
+  integer, allocatable, dimension(:),   public :: groupStartCol
+  integer, allocatable, dimension(:),   public :: groupListLocal
+  integer, allocatable, dimension(:),   public :: groupStartLocal
+  integer, allocatable, dimension(:),   public :: nSkipsForAtom
+  integer, allocatable, dimension(:,:), public :: skipsForAtom
+  real(kind=dp), allocatable, dimension(:), public :: mfactRow
+  real(kind=dp), allocatable, dimension(:), public :: mfactCol
+  real(kind=dp), allocatable, dimension(:), public :: mfactLocal
 
-  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, public, save :: boxIsOrthorhombic
   
   public :: SimulationSetup
   public :: getNlocal
   public :: setBox
-  public :: setBox_3d
-  public :: getBox
-  public :: setRcut
-  public :: getRcut
-  public :: getRlist
-  public :: getRrf
-  public :: getRt
   public :: getDielect
   public :: SimUsesPBC
   public :: SimUsesLJ
+  public :: SimUsesCharges
   public :: SimUsesDipoles
   public :: SimUsesSticky
   public :: SimUsesRF
@@ -54,22 +59,12 @@ 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, CnGlobal, CnLocal, c_idents, &
        CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
-       CmolMembership, &
+       CmolMembership, Cmfact, CnGroups, CglobalGroupMembership, &
        status)    
 
     type (simtype) :: setThisSim
@@ -83,12 +78,20 @@ contains
     integer, dimension(CnGlobal),intent(in) :: CmolMembership 
     !!  Result status, success = 0, status = -1
     integer, intent(out) :: status
-    integer :: i, me, thisStat, alloc_stat, myNode
+    integer :: i, j, me, thisStat, alloc_stat, myNode, id1, id2
+    integer :: ia
+
+    !! mass factors used for molecular cutoffs
+    real ( kind = dp ), dimension(CnLocal) :: Cmfact
+    integer, intent(in):: CnGroups
+    integer, dimension(CnGlobal), intent(in):: CglobalGroupMembership
+    integer :: maxSkipsForAtom, glPointer
+
 #ifdef IS_MPI
     integer, allocatable, dimension(:) :: c_idents_Row
     integer, allocatable, dimension(:) :: c_idents_Col
-    integer :: nrow
-    integer :: ncol
+    integer :: nAtomsInRow, nGroupsInRow, aid
+    integer :: nAtomsInCol, nGroupsInCol, gid
 #endif  
 
     simulation_setup_complete = .false.
@@ -98,12 +101,9 @@ contains
 
     nLocal = CnLocal
     nGlobal = CnGlobal
+    nGroups = CnGroups
 
     thisSim = setThisSim
-    rcut2 = thisSim%rcut * thisSim%rcut
-    rcut6 = rcut2 * rcut2 * rcut2
-    rlist2 = thisSim%rlist * thisSim%rlist
-    box = thisSim%box
 
     nExcludes_Global = CnGlobalExcludes
     nExcludes_Local = CnLocalExcludes
@@ -129,31 +129,33 @@ contains
        status = -1
        return
     endif
-    nrow = getNrow(plan_row)
-    ncol = getNcol(plan_col)
+    nAtomsInRow = getNatomsInRow(plan_atom_row)
+    nAtomsInCol = getNatomsInCol(plan_atom_col)
+    nGroupsInRow = getNgroupsInRow(plan_group_row)
+    nGroupsInCol = getNgroupsInCol(plan_group_col)
     mynode = getMyNode()
     
-    allocate(c_idents_Row(nrow),stat=alloc_stat)
+    allocate(c_idents_Row(nAtomsInRow),stat=alloc_stat)
     if (alloc_stat /= 0 ) then
        status = -1
        return
     endif
 
-    allocate(c_idents_Col(ncol),stat=alloc_stat)
+    allocate(c_idents_Col(nAtomsInCol),stat=alloc_stat)
     if (alloc_stat /= 0 ) then
        status = -1
        return
     endif
 
-    call gather(c_idents, c_idents_Row, plan_row)
-    call gather(c_idents, c_idents_Col, plan_col)
+    call gather(c_idents, c_idents_Row, plan_atom_row)
+    call gather(c_idents, c_idents_Col, plan_atom_col)
 
-    do i = 1, nrow
+    do i = 1, nAtomsInRow
        me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
        atid_Row(i) = me
     enddo
 
-    do i = 1, ncol
+    do i = 1, nAtomsInCol
        me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
        atid_Col(i) = me
     enddo
@@ -165,22 +167,240 @@ contains
     if (allocated(c_idents_Row)) then
        deallocate(c_idents_Row)
     endif
+   
+#endif
+
+#ifdef IS_MPI
+    allocate(groupStartRow(nGroupsInRow+1),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(groupStartCol(nGroupsInCol+1),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(groupListRow(nAtomsInRow),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(groupListCol(nAtomsInCol),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(mfactRow(nAtomsInRow),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(mfactCol(nAtomsInCol),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(mfactLocal(nLocal),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
     
+    glPointer = 1
+
+    do i = 1, nGroupsInRow 
+
+       gid = GroupRowToGlobal(i)
+       groupStartRow(i) = glPointer       
+
+       do j = 1, nAtomsInRow
+          aid = AtomRowToGlobal(j)
+          if (CglobalGroupMembership(aid) .eq. gid) then
+             groupListRow(glPointer) = j
+             glPointer = glPointer + 1
+          endif
+       enddo
+    enddo
+    groupStartRow(nGroupsInRow+1) = nAtomsInRow + 1
+
+    glPointer = 1
+
+    do i = 1, nGroupsInCol
+
+       gid = GroupColToGlobal(i)
+       groupStartCol(i) = glPointer       
+
+       do j = 1, nAtomsInCol
+          aid = AtomColToGlobal(j)
+          if (CglobalGroupMembership(aid) .eq. gid) then
+             groupListCol(glPointer) = j
+             glPointer = glPointer + 1
+          endif
+       enddo
+    enddo
+    groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1
+
+    mfactLocal = Cmfact        
+
+    call gather(mfactLocal,      mfactRow,      plan_atom_row)
+    call gather(mfactLocal,      mfactCol,      plan_atom_col)
+    
+    if (allocated(mfactLocal)) then
+       deallocate(mfactLocal)
+    end if
 #else
+    allocate(groupStartRow(nGroups+1),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(groupStartCol(nGroups+1),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(groupListRow(nLocal),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(groupListCol(nLocal),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(mfactRow(nLocal),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(mfactCol(nLocal),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+    allocate(mfactLocal(nLocal),stat=alloc_stat)
+    if (alloc_stat /= 0 ) then
+       status = -1
+       return
+    endif
+
+    glPointer = 1
+    do i = 1, nGroups
+       groupStartRow(i) = glPointer       
+       groupStartCol(i) = glPointer
+       do j = 1, nLocal
+          if (CglobalGroupMembership(j) .eq. i) then
+             groupListRow(glPointer) = j
+             groupListCol(glPointer) = j
+             glPointer = glPointer + 1
+          endif
+       enddo
+    enddo
+    groupStartRow(nGroups+1) = nLocal + 1
+    groupStartCol(nGroups+1) = nLocal + 1
+
     do i = 1, nLocal
+       mfactRow(i) = Cmfact(i)
+       mfactCol(i) = Cmfact(i)
+    end do
+    
+#endif
+
+
+! We build the local atid's for both mpi and nonmpi
+    do i = 1, nLocal
        
        me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
        atid(i) = me
   
     enddo
-#endif
 
-
-
     do i = 1, nExcludes_Local
        excludesLocal(1,i) = CexcludesLocal(1,i)
        excludesLocal(2,i) = CexcludesLocal(2,i)
     enddo
+
+#ifdef IS_MPI
+    allocate(nSkipsForAtom(nAtomsInRow), stat=alloc_stat)
+#else
+    allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
+#endif
+    if (alloc_stat /= 0 ) then
+       thisStat = -1
+       write(*,*) 'Could not allocate nSkipsForAtom array'
+       return
+    endif
+
+    maxSkipsForAtom = 0
+#ifdef IS_MPI
+    do j = 1, nAtomsInRow
+#else
+    do j = 1, nLocal
+#endif
+       nSkipsForAtom(j) = 0
+#ifdef IS_MPI
+       id1 = AtomRowToGlobal(j)
+#else 
+       id1 = j
+#endif
+       do i = 1, nExcludes_Local
+          if (excludesLocal(1,i) .eq. id1 ) then
+             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
+
+             if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
+                maxSkipsForAtom = nSkipsForAtom(j)
+             endif
+          endif
+          if (excludesLocal(2,i) .eq. id1 ) then
+             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
+
+             if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
+                maxSkipsForAtom = nSkipsForAtom(j)
+             endif
+          endif
+       end do
+    enddo
+
+#ifdef IS_MPI
+    allocate(skipsForAtom(nAtomsInRow, maxSkipsForAtom), stat=alloc_stat)
+#else
+    allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
+#endif
+    if (alloc_stat /= 0 ) then
+       write(*,*) 'Could not allocate skipsForAtom array'
+       return
+    endif
+
+#ifdef IS_MPI
+    do j = 1, nAtomsInRow
+#else
+    do j = 1, nLocal
+#endif
+       nSkipsForAtom(j) = 0
+#ifdef IS_MPI
+       id1 = AtomRowToGlobal(j)
+#else 
+       id1 = j
+#endif
+       do i = 1, nExcludes_Local
+          if (excludesLocal(1,i) .eq. id1 ) then
+             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
+             ! exclude lists have global ID's so this line is 
+             ! the same in MPI and non-MPI
+             id2 = excludesLocal(2,i)
+             skipsForAtom(j, nSkipsForAtom(j)) = id2
+          endif
+          if (excludesLocal(2, i) .eq. id1 ) then
+             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
+             ! exclude lists have global ID's so this line is 
+             ! the same in MPI and non-MPI
+             id2 = excludesLocal(1,i)
+             skipsForAtom(j, nSkipsForAtom(j)) = id2
+          endif
+       end do
+    enddo
     
     do i = 1, nExcludes_Global
        excludesGlobal(i) = CexcludesGlobal(i)
@@ -188,97 +408,28 @@ contains
 
     do i = 1, nGlobal
        molMemberShipList(i) = CmolMembership(i)
-     enddo
-
+    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, cHmatInv
+    integer :: cBoxIsOrthorhombic
     integer :: smallest, status, i
-
-    thisSim%box = new_box_size
-    box = thisSim%box
-
-    return     
-  end subroutine setBox_3d
-
-  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)
-  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
-    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
-    real( kind = dp ), intent(out), optional :: rc2
-    real( kind = dp ), intent(out), optional :: rc6
-    integer, optional :: status
-
-    if (present(status)) status = 0
+    Hmat = cHmat
+    HmatInv = cHmatInv
+    if (cBoxIsOrthorhombic .eq. 0 ) then
+       boxIsOrthorhombic = .false.
+    else 
+       boxIsOrthorhombic = .true.
+    endif
     
-    if (.not.simulation_setup_complete ) then
-       if (present(status)) status = -1
-       return
-    end if
-    
-    thisrcut = thisSim%rcut
-    if(present(rc2)) rc2 = rcut2
-    if(present(rc6)) rc6 = rcut6
-  end subroutine getRcut
-  
-  subroutine getRlist(thisrlist,rl2,status)
-    real( kind = dp ), intent(out) :: thisrlist
-    real( kind = dp ), intent(out), optional :: rl2
+    return     
+  end subroutine setBox
 
-    integer, optional :: status
-
-    if (present(status)) status = 0
-
-    if (.not.simulation_setup_complete ) then
-       if (present(status)) status = -1
-       return
-    end if
-    
-    thisrlist = thisSim%rlist
-    if(present(rl2)) rl2 = rlist2
-  end subroutine getRlist
-
-  function getRrf() result(rrf)
-    real( kind = dp ) :: rrf
-    rrf = thisSim%rrf
-    write(*,*) 'getRrf = ', rrf, thisSim%rrf
-  end function getRrf
-  
-  function getRt() result(rt)
-    real( kind = dp ) :: rt
-    rt = thisSim%rt
-  end function getRt
-
   function getDielect() result(dielect)
     real( kind = dp ) :: dielect
     dielect = thisSim%dielect
@@ -299,6 +450,11 @@ contains
     doesit = thisSim%SIM_uses_sticky
   end function SimUsesSticky
 
+  function SimUsesCharges() result(doesit)
+    logical :: doesit
+    doesit = thisSim%SIM_uses_charges
+  end function SimUsesCharges
+
   function SimUsesDipoles() result(doesit)
     logical :: doesit
     doesit = thisSim%SIM_uses_dipoles
@@ -367,16 +523,25 @@ contains
     
     !We free in the opposite order in which we allocate in.
 
+    if (allocated(skipsForAtom)) deallocate(skipsForAtom)
+    if (allocated(nSkipsForAtom)) deallocate(nSkipsForAtom)
+    if (allocated(mfactLocal)) deallocate(mfactLocal)
+    if (allocated(mfactCol)) deallocate(mfactCol)
+    if (allocated(mfactRow)) deallocate(mfactRow)
+    if (allocated(groupListCol)) deallocate(groupListCol)     
+    if (allocated(groupListRow)) deallocate(groupListRow)    
+    if (allocated(groupStartCol)) deallocate(groupStartCol)
+    if (allocated(groupStartRow)) deallocate(groupStartRow)     
     if (allocated(molMembershipList)) deallocate(molMembershipList)     
     if (allocated(excludesGlobal)) deallocate(excludesGlobal)
     if (allocated(excludesLocal)) deallocate(excludesLocal)
-
+    
   end subroutine FreeSimGlobals
-
+  
   pure function getNlocal() result(n)
     integer :: n
     n = nLocal
   end function getNlocal
-
   
+  
 end module simulation