--- trunk/mdtools/md_code/simulation_module.F90 2003/01/02 21:45:45 222 +++ trunk/mdtools/md_code/simulation_module.F90 2003/01/28 22:16:55 252 @@ -1,33 +1,60 @@ module simulation use definitions, ONLY :dp - use force_wrappers, ONLY : alloc_force_wrappers +#ifdef IS_MPI + use mpiSimulation +#endif implicit none PRIVATE - integer :: nLR + type, public :: simtype + PRIVATE +! SEQUENCE +!! Number of particles on this processor + integer :: nLRparticles +!! Periodic Box + real ( kind = dp ), dimension(3) :: box +!! List Cutoff + real ( kind = dp ) :: rlist = 0.0_dp +!! Radial cutoff + real ( kind = dp ) :: rcut = 0.0_dp +!! List cutoff squared + real ( kind = dp ) :: rlistsq = 0.0_dp +!! Radial Cutoff squared + real ( kind = dp ) :: rcutsq = 0.0_dp +!! Radial Cutoff^6 + real ( kind = dp ) :: rcut6 = 0.0_dp - real ( kind = dp ), dimension(3) :: box + end type simtype + type (simtype), public :: thisSim +!! Tag for MPI calculations + integer, allocatable, dimension(:) :: tag - real ( kind = dp ), public :: rlist - real ( kind = dp ), public :: rcut - real ( kind = dp ), public :: rlistsq - real ( kind = dp ), public :: rcutsq +#ifdef IS_MPI + integer, allocatable, dimension(:) :: tag_row + integer, allocatable, dimension(:) :: tag_column +#endif - integer,public, allocatable, dimension(:) :: point - integer,public, allocatable, dimension(:) :: list +!! WARNING: use_pbc hardcoded, fixme + logical :: use_pbc = .true. + logical :: setSim = .false. +!! array for saving previous positions for neighbor lists. + real( kind = dp ), allocatable,dimension(:,:),save :: q0 - - public :: check public :: save_nlist public :: wrap public :: getBox + public :: getRcut + public :: getRlist + public :: getNlocal + public :: setSimulation +! public :: setRcut interface wrap module procedure wrap_1d @@ -40,84 +67,80 @@ module simulation end interface -!MPI dependent routines -#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 - real( kind = dp ), allocatable, dimension(:,:) :: fRow - real( kind = dp ), allocatable, dimension(:,:) :: fColumn - real( kind = dp ), allocatable, dimension(:,:) :: tRow - real( kind = dp ), allocatable, dimension(:,:) :: tColumn -#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) + integer, intent(in) :: nLRParticles + 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. - - - - - end subroutine allocate_mpi_arrays -#endif - - subroutine set_simulation(box,rlist,rcut) + thisSim%nLRParticles = nLRParticles + thisSim%box = box + thisSim%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) + integer :: nparticles + nparticles = thisSim%nLRparticles + end function getNparticles - end subroutine set_simulation - - subroutine change_box_size(new_box_size) real(kind=dp), dimension(3) :: new_box_size - box = new_box_size + thisSim%box = new_box_size end subroutine change_box_size - elemental function getBox_3d() result(thisBox) + function getBox_3d() result(thisBox) real( kind = dp ), dimension(3) :: thisBox - thisBox = box + thisBox = thisSim%box end function getBox_3d function getBox_dim(dim) result(thisBox) integer, intent(in) :: dim real( kind = dp ) :: thisBox - thisBox = box(dim) + thisBox = thisSim%box(dim) end function getBox_dim - subroutine check(update_nlist) - + 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 dispmx = 0.0E0_DP - !! calculate the largest displacement of any atom in any direction #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 ) @@ -125,7 +148,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 ) @@ -139,18 +163,22 @@ contains end subroutine check - subroutine save_nlist() - integer :: i -#ifdef MPI - do i = 1, nlocal -#else - do i = 1, natoms -#endif - q0(1,i) = q(1,i) - q0(2,i) = q(2,i) - q0(3,i) = q(3,i) - end do + subroutine save_nlist(q) + real(kind = dp ), dimension(:,:), intent(in) :: q + integer :: list_size + + 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 @@ -163,7 +191,7 @@ contains 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 - box(dim)*nint(r/box(dim)) + this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim)) else this_wrap = r endif @@ -171,18 +199,68 @@ contains return end function wrap_1d - elemental function wrap_3d(r) result(this_wrap) + function wrap_3d(r) result(this_wrap) real( kind = dp ), dimension(3), intent(in) :: r real( kind = dp ), dimension(3) :: this_wrap if (use_pbc) then ! this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP) - this_wrap(1:3) = r(1:3) - box(1:3)*nint(r(1:3)/box(1:3)) + 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 + + 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 + + end subroutine getRcut + + + + + subroutine getRlist(thisrlist,rlist2,status) + real( kind = dp ), intent(out) :: thisrlist + real( kind = dp ), intent(out), optional :: rlist2 + + 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 + + + end module simulation