--- trunk/OOPSE/libmdtools/simulation_module.F90 2003/03/24 21:55:34 394 +++ trunk/OOPSE/libmdtools/simulation_module.F90 2004/05/27 19:51:18 1206 @@ -6,7 +6,7 @@ module simulation use force_globals use vector_class use atype_module - use lj + use switcheroo #ifdef IS_MPI use mpiSimulation #endif @@ -16,36 +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 :: natoms + 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 :: 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,63 +59,66 @@ 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, Cmfact, CnGroups, CgroupList, CgroupStart, & 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 + 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(CnLocal),intent(in) :: CgroupList + integer, dimension(CnGroups),intent(in) :: CgroupStart + integer :: maxSkipsForAtom + #ifdef IS_MPI integer, allocatable, dimension(:) :: c_idents_Row integer, allocatable, dimension(:) :: c_idents_Col - integer :: nrow - integer :: ncol + integer :: nAtomsInRow, nGroupsInRow + integer :: nAtomsInCol, nGroupsInCol #endif simulation_setup_complete = .false. status = 0 ! copy C struct into fortran type + + nLocal = CnLocal + nGlobal = CnGlobal + nGroups = CnGroups + thisSim = setThisSim - natoms = nComponents - rcut2 = thisSim%rcut * thisSim%rcut - rcut6 = rcut2 * rcut2 * rcut2 - rlist2 = thisSim%rlist * thisSim%rlist - box = thisSim%box 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 @@ -122,31 +130,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 @@ -158,137 +168,255 @@ 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(groupStartLocal(nGroups), stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif + allocate(groupListLocal(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 + + groupStartLocal = CgroupStart + groupListLocal = CgroupList + mfactLocal = Cmfact + + call gather(groupStartLocal, groupStartRow, plan_group_row) + call gather(groupStartLocal, groupStartCol, plan_group_col) + groupStartRow(nGroupsInRow+1) = nAtomsInRow + 1 + groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1 + call gather(groupListLocal, groupListRow, plan_atom_row) + call gather(groupListLocal, groupListCol, plan_atom_col) + + ! C passes us groupList as globalID numbers for the atoms + ! we need to remap these onto the row and column ids on this + ! processor. This is a linear search, but is only done once + ! (we hope) + + do i = 1, nAtomsInRow + do j = 1, nAtomsInRow + if (AtomRowToGlobal(j) .eq. groupListRow(i)) then + groupListRow(i) = j + endif + enddo + enddo + do i = 1, nAtomsInCol + do j = 1, nAtomsInCol + if (AtomColToGlobal(j) .eq. groupListCol(i)) then + groupListCol(i) = j + endif + enddo + enddo + call gather(mfactLocal, mfactRow, plan_atom_row) + call gather(mfactLocal, mfactCol, plan_atom_col) + if (allocated(mfactLocal)) then + deallocate(mfactLocal) + end if + if (allocated(groupListLocal)) then + deallocate(groupListLocal) + endif + if (allocated(groupStartLocal)) then + deallocate(groupStartLocal) + endif #else - do i = 1, nComponents + 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 + do i = 1, nGroups + groupStartRow(i) = CgroupStart(i) + groupStartCol(i) = CgroupStart(i) + end do + groupStartRow(nGroups+1) = nLocal + 1 + groupStartCol(nGroups+1) = nLocal + 1 + do i = 1, nLocal + groupListRow(i) = CgroupList(i) + groupListCol(i) = CgroupList(i) + 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 + + do i = 1, nExcludes_Local + excludesLocal(1,i) = CexcludesLocal(1,i) + excludesLocal(2,i) = CexcludesLocal(2,i) + enddo + + 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 - !! Create neighbor lists - call expandNeighborList(nComponents, thisStat) - if (thisStat /= 0) then - status = -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 - - do i = 1, nExcludes_Local - excludesLocal(1,i) = CexcludesLocal(1,i) - excludesLocal(2,i) = CexcludesLocal(2,i) +#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. id2 ) 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) 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, cHmatInv + integer :: cBoxIsOrthorhombic integer :: smallest, status, i - - thisSim%box = new_box_size - box = thisSim%box - - 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 - - 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 + + Hmat = cHmat + HmatInv = cHmatInv + if (cBoxIsOrthorhombic .eq. 0 ) then + boxIsOrthorhombic = .false. + else + boxIsOrthorhombic = .true. 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 - real( kind = dp ), intent(out), optional :: rc2 - real( kind = dp ), intent(out), optional :: rc6 - integer, optional :: status + return + end subroutine setBox - if (present(status)) status = 0 - - 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 - - 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 @@ -309,6 +437,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 @@ -352,6 +485,12 @@ contains thisStat = 0 call FreeSimGlobals() + + allocate(nSkipsForAtom(nLocal), stat=alloc_stat) + if (alloc_stat /= 0 ) then + thisStat = -1 + return + endif allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat) if (alloc_stat /= 0 ) then @@ -364,22 +503,36 @@ 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(skipsForAtom)) deallocate(skipsForAtom) + 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(nlocal) - integer :: nlocal - nlocal = natoms + + pure function getNlocal() result(n) + integer :: n + n = nLocal end function getNlocal - + end module simulation