mdtools/md_code/simulation_module.F90

module simulation
  use definitions, ONLY :dp
  use force_wrappers, ONLY : alloc_force_wrappers

  implicit none
  PRIVATE


  integer :: nLR


  real ( kind = dp ), dimension(3) :: box


  real ( kind = dp ), public :: rlist
  real ( kind = dp ), public :: rcut
  real ( kind = dp ), public :: rlistsq
  real ( kind = dp ), public :: rcutsq

  integer,public, allocatable, dimension(:)   :: point
  integer,public, allocatable, dimension(:)   :: list


  public :: check
  public :: save_nlist
  public :: wrap
  public :: getBox

  interface wrap
     module procedure wrap_1d
     module procedure wrap_3d
  end interface

  interface getBox
     module procedure getBox_3d
     module procedure getBox_dim
  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


  end subroutine allocate_mpi_arrays
#endif

  subroutine set_simulation(box,rlist,rcut)
    

  end subroutine set_simulation


  subroutine change_box_size(new_box_size)
    real(kind=dp), dimension(3) :: new_box_size

    box = new_box_size

  end subroutine change_box_size


  elemental function getBox_3d() result(thisBox)
    real( kind = dp ), dimension(3) :: thisBox
    thisBox = box
  end function getBox_3d

  function getBox_dim(dim) result(thisBox)
    integer, intent(in) :: dim
    real( kind = dp ) :: thisBox

    thisBox = box(dim)
  end function getBox_dim
    
  subroutine check(update_nlist)
  
    integer :: i
    real( kind = DP ) :: dispmx 
    logical, intent(out) :: update_nlist
    real( kind = DP ) :: dispmx_tmp
    
    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
       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, natoms
       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()
    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

  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 - box(dim)*nint(r/box(dim))
    else
       this_wrap = r
    endif
    
    return
  end function wrap_1d

  elemental 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))
    else
       this_wrap = r
    endif
  end function wrap_3d

  
end module simulation
Revision:	222
Committed:	Thu Jan 2 21:45:45 2003 UTC (21 years, 6 months ago) by chuckv
File size:	4168 byte(s)
Log Message:	More codeing on fortran force backend.
#	User	Rev	Content
1	chuckv	167	module simulation
2			use definitions, ONLY :dp
3	chuckv	194	use force_wrappers, ONLY : alloc_force_wrappers
4
5	chuckv	167	implicit none
6			PRIVATE
7
8
9	chuckv	222	integer :: nLR
10	chuckv	167
11
12	chuckv	222	real ( kind = dp ), dimension(3) :: box
13
14
15	chuckv	194	real ( kind = dp ), public :: rlist
16			real ( kind = dp ), public :: rcut
17			real ( kind = dp ), public :: rlistsq
18			real ( kind = dp ), public :: rcutsq
19	chuckv	167
20	chuckv	194	integer,public, allocatable, dimension(:) :: point
21			integer,public, allocatable, dimension(:) :: list
22	chuckv	167
23
24	chuckv	222
25
26
27			public :: check
28			public :: save_nlist
29			public :: wrap
30			public :: getBox
31
32			interface wrap
33			module procedure wrap_1d
34			module procedure wrap_3d
35			end interface
36
37			interface getBox
38			module procedure getBox_3d
39			module procedure getBox_dim
40			end interface
41
42
43	chuckv	194	!MPI dependent routines
44	chuckv	185
45	chuckv	194	#ifdef IS_MPI
46			! Universal routines: All types of force calculations will need these arrays
47			! Arrays specific to a type of force calculation should be declared in that module.
48			real( kind = dp ), allocatable, dimension(:,:) :: qRow
49			real( kind = dp ), allocatable, dimension(:,:) :: qColumn
50	chuckv	185
51	chuckv	194	real( kind = dp ), allocatable, dimension(:,:) :: fRow
52			real( kind = dp ), allocatable, dimension(:,:) :: fColumn
53	chuckv	185
54	chuckv	194	real( kind = dp ), allocatable, dimension(:,:) :: tRow
55			real( kind = dp ), allocatable, dimension(:,:) :: tColumn
56
57			#endif
58
59
60
61
62
63
64	chuckv	167	contains
65
66	chuckv	194	#ifdef MPI
67			! Allocated work arrays for MPI
68			subroutine allocate_mpi_arrays(nDimensions,numComponents)
69			integer, intent(in) :: nDimensions
70			integer, intent(in) :: numComponents
71	chuckv	167
72	chuckv	194
73
74
75
76			end subroutine allocate_mpi_arrays
77			#endif
78
79	chuckv	167	subroutine set_simulation(box,rlist,rcut)
80
81
82
83			end subroutine set_simulation
84
85
86
87			subroutine change_box_size(new_box_size)
88			real(kind=dp), dimension(3) :: new_box_size
89
90			box = new_box_size
91
92			end subroutine change_box_size
93
94	chuckv	222
95			elemental function getBox_3d() result(thisBox)
96			real( kind = dp ), dimension(3) :: thisBox
97			thisBox = box
98			end function getBox_3d
99
100			function getBox_dim(dim) result(thisBox)
101			integer, intent(in) :: dim
102			real( kind = dp ) :: thisBox
103
104			thisBox = box(dim)
105			end function getBox_dim
106
107			subroutine check(update_nlist)
108	chuckv	215
109	chuckv	222	integer :: i
110			real( kind = DP ) :: dispmx
111			logical, intent(out) :: update_nlist
112			real( kind = DP ) :: dispmx_tmp
113
114			dispmx = 0.0E0_DP
115
116			!! calculate the largest displacement of any atom in any direction
117
118			#ifdef MPI
119			dispmx_tmp = 0.0E0_DP
120			do i = 1, nlocal
121			dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
122			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
123			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
124			end do
125			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
126			mpi_max,mpi_comm_world,mpi_err)
127			#else
128			do i = 1, natoms
129			dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
130			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
131			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
132			end do
133			#endif
134
135			!! a conservative test of list skin crossings
136			dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
137
138			update_nlist = (dispmx.gt.(skin_thickness))
139
140			end subroutine check
141
142			subroutine save_nlist()
143			integer :: i
144			#ifdef MPI
145			do i = 1, nlocal
146			#else
147			do i = 1, natoms
148			#endif
149			q0(1,i) = q(1,i)
150			q0(2,i) = q(2,i)
151			q0(3,i) = q(3,i)
152			end do
153
154			end subroutine save_nlist
155
156
157			function wrap_1d(r,dim) result(this_wrap)
158
159
160			real( kind = DP ) :: r
161			real( kind = DP ) :: this_wrap
162			integer :: dim
163
164			if (use_pbc) then
165			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
166			this_wrap = r - box(dim)*nint(r/box(dim))
167			else
168			this_wrap = r
169			endif
170
171			return
172			end function wrap_1d
173
174			elemental function wrap_3d(r) result(this_wrap)
175			real( kind = dp ), dimension(3), intent(in) :: r
176			real( kind = dp ), dimension(3) :: this_wrap
177
178
179			if (use_pbc) then
180			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
181			this_wrap(1:3) = r(1:3) - box(1:3)*nint(r(1:3)/box(1:3))
182			else
183			this_wrap = r
184			endif
185			end function wrap_3d
186
187
188	chuckv	167	end module simulation