mdtools/md_code/simulation_module.F90

module simulation
  use definitions, ONLY :dp
#ifdef IS_MPI
  use mpiSimulation
#endif

  implicit none
  PRIVATE


  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

  end type simtype

  type (simtype), public :: thisSim
!! Tag for MPI calculations  
  integer, allocatable, dimension(:) :: tag

#ifdef IS_MPI
  integer, allocatable, dimension(:) :: tag_row
  integer, allocatable, dimension(:) :: tag_column
#endif

!! 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
     module procedure wrap_3d
  end interface

  interface getBox
     module procedure getBox_3d
     module procedure getBox_dim
  end interface


contains

  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

    if( setsim ) return  ! simulation is already initialized
    setSim = .true.

    thisSim%nLRParticles = nLRParticles
    thisSim%box          = box
    thisSim%rlist        = rlist
    thisSim%rcut         = rcut
    thisSim%rcutsq       = rcut * rcut
    thisSim%rcut6        = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq

  end subroutine setSimulation

  function getNparticles() result(nparticles)
    integer :: nparticles
    nparticles = thisSim%nLRparticles
  end function getNparticles


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

    thisSim%box = new_box_size

  end subroutine change_box_size


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

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

    thisBox = thisSim%box(dim)
  end function getBox_dim
    
  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
       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(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
  

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

  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 = 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(rcut2)) 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
Revision:	249
Committed:	Mon Jan 27 21:28:19 2003 UTC (21 years, 5 months ago) by chuckv
File size:	6163 byte(s)
Log Message:	For some unknown reason the Single processor builds. Has not been tested!
#	User	Rev	Content
1	chuckv	167	module simulation
2			use definitions, ONLY :dp
3	chuckv	230	#ifdef IS_MPI
4			use mpiSimulation
5			#endif
6	chuckv	194
7	chuckv	167	implicit none
8			PRIVATE
9
10
11
12	chuckv	230	type, public :: simtype
13			PRIVATE
14	chuckv	246	! SEQUENCE
15	chuckv	230	!! Number of particles on this processor
16			integer :: nLRparticles
17			!! Periodic Box
18			real ( kind = dp ), dimension(3) :: box
19			!! List Cutoff
20			real ( kind = dp ) :: rlist = 0.0_dp
21			!! Radial cutoff
22			real ( kind = dp ) :: rcut = 0.0_dp
23			!! List cutoff squared
24			real ( kind = dp ) :: rlistsq = 0.0_dp
25			!! Radial Cutoff squared
26			real ( kind = dp ) :: rcutsq = 0.0_dp
27			!! Radial Cutoff^6
28			real ( kind = dp ) :: rcut6 = 0.0_dp
29	chuckv	167
30	chuckv	230	end type simtype
31	chuckv	222
32	chuckv	230	type (simtype), public :: thisSim
33			!! Tag for MPI calculations
34			integer, allocatable, dimension(:) :: tag
35	chuckv	222
36	chuckv	230	#ifdef IS_MPI
37			integer, allocatable, dimension(:) :: tag_row
38			integer, allocatable, dimension(:) :: tag_column
39			#endif
40	chuckv	167
41	chuckv	246	!! WARNING: use_pbc hardcoded, fixme
42			logical :: use_pbc = .true.
43	chuckv	230	logical :: setSim = .false.
44	chuckv	247
45			!! array for saving previous positions for neighbor lists.
46	chuckv	246	real( kind = dp ), allocatable,dimension(:,:),save :: q0
47	chuckv	167
48	chuckv	247
49	chuckv	222	public :: check
50			public :: save_nlist
51			public :: wrap
52			public :: getBox
53	chuckv	230	public :: getRcut
54	chuckv	246	public :: getRlist
55			public :: getNlocal
56	chuckv	249	public :: setSimulation
57	chuckv	246	! public :: setRcut
58	chuckv	222
59			interface wrap
60			module procedure wrap_1d
61			module procedure wrap_3d
62			end interface
63
64			interface getBox
65			module procedure getBox_3d
66			module procedure getBox_dim
67			end interface
68
69
70	chuckv	185
71
72
73	chuckv	194
74
75
76
77
78	chuckv	167	contains
79
80	chuckv	230	subroutine setSimulation(nLRParticles,box,rlist,rcut)
81			integer, intent(in) :: nLRParticles
82	mmeineke	241	real(kind = dp ), intent(in), dimension(3) :: box
83	chuckv	230	real(kind = dp ), intent(in) :: rlist
84			real(kind = dp ), intent(in) :: rcut
85	chuckv	167
86	chuckv	230	if( setsim ) return ! simulation is already initialized
87			setSim = .true.
88	chuckv	194
89	chuckv	230	thisSim%nLRParticles = nLRParticles
90			thisSim%box = box
91			thisSim%rlist = rlist
92			thisSim%rcut = rcut
93			thisSim%rcutsq = rcut * rcut
94	chuckv	246	thisSim%rcut6 = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
95	chuckv	194
96	chuckv	230	end subroutine setSimulation
97	chuckv	194
98	chuckv	230	function getNparticles() result(nparticles)
99			integer :: nparticles
100			nparticles = thisSim%nLRparticles
101			end function getNparticles
102	chuckv	194
103
104	chuckv	167	subroutine change_box_size(new_box_size)
105			real(kind=dp), dimension(3) :: new_box_size
106
107	chuckv	230	thisSim%box = new_box_size
108	chuckv	167
109			end subroutine change_box_size
110
111	chuckv	222
112	chuckv	246	function getBox_3d() result(thisBox)
113	chuckv	222	real( kind = dp ), dimension(3) :: thisBox
114	chuckv	230	thisBox = thisSim%box
115	chuckv	222	end function getBox_3d
116
117			function getBox_dim(dim) result(thisBox)
118			integer, intent(in) :: dim
119			real( kind = dp ) :: thisBox
120
121	chuckv	230	thisBox = thisSim%box(dim)
122	chuckv	222	end function getBox_dim
123
124	chuckv	246	subroutine check(q,update_nlist)
125			real( kind = dp ), dimension(:,:) :: q
126	chuckv	222	integer :: i
127			real( kind = DP ) :: dispmx
128			logical, intent(out) :: update_nlist
129			real( kind = DP ) :: dispmx_tmp
130	chuckv	246	real( kind = dp ) :: skin_thickness
131			integer :: natoms
132
133			natoms = thisSim%nLRparticles
134			skin_thickness = thisSim%rlist
135	chuckv	222	dispmx = 0.0E0_DP
136			!! calculate the largest displacement of any atom in any direction
137
138			#ifdef MPI
139			dispmx_tmp = 0.0E0_DP
140			do i = 1, nlocal
141			dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
142			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
143			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
144			end do
145			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
146			mpi_max,mpi_comm_world,mpi_err)
147			#else
148			do i = 1, natoms
149			dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
150			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
151			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
152			end do
153			#endif
154
155			!! a conservative test of list skin crossings
156			dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
157
158			update_nlist = (dispmx.gt.(skin_thickness))
159
160			end subroutine check
161
162	chuckv	230	subroutine save_nlist(q)
163			real(kind = dp ), dimension(:,:), intent(in) :: q
164			integer :: list_size
165
166	chuckv	222
167	chuckv	230	list_size = size(q)
168
169			if (.not. allocated(q0)) then
170			allocate(q0(3,list_size))
171	chuckv	246	else if( list_size > size(q0)) then
172	chuckv	230	deallocate(q0)
173			allocate(q0(3,list_size))
174			endif
175
176			q0 = q
177
178	chuckv	222	end subroutine save_nlist
179
180
181			function wrap_1d(r,dim) result(this_wrap)
182
183
184			real( kind = DP ) :: r
185			real( kind = DP ) :: this_wrap
186			integer :: dim
187
188			if (use_pbc) then
189			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
190	chuckv	230	this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
191	chuckv	222	else
192			this_wrap = r
193			endif
194
195			return
196			end function wrap_1d
197
198	chuckv	246	function wrap_3d(r) result(this_wrap)
199	chuckv	222	real( kind = dp ), dimension(3), intent(in) :: r
200			real( kind = dp ), dimension(3) :: this_wrap
201
202
203			if (use_pbc) then
204			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
205	chuckv	246	this_wrap = r - thisSim%box*nint(r/thisSim%box)
206	chuckv	222	else
207			this_wrap = r
208			endif
209			end function wrap_3d
210
211	chuckv	246
212
213	chuckv	230	subroutine getRcut(thisrcut,rcut2,rcut6,status)
214			real( kind = dp ), intent(out) :: thisrcut
215			real( kind = dp ), intent(out), optional :: rcut2
216	chuckv	246	real( kind = dp ), intent(out), optional :: rcut6
217	chuckv	230	integer, optional :: status
218
219			if (present(status)) status = 0
220
221			if (.not.setSim ) then
222			if (present(status)) status = -1
223			return
224			end if
225
226	chuckv	246	thisrcut = thisSim%rcut
227			if(present(rcut2)) rcut2 = thisSim%rcutsq
228			if(present(rcut2)) rcut6 = thisSim%rcut6
229	chuckv	230
230			end subroutine getRcut
231	chuckv	222
232	chuckv	246
233
234
235			subroutine getRlist(thisrlist,rlist2,status)
236			real( kind = dp ), intent(out) :: thisrlist
237			real( kind = dp ), intent(out), optional :: rlist2
238	chuckv	240
239	chuckv	246	integer, optional :: status
240	chuckv	240
241	chuckv	246	if (present(status)) status = 0
242	chuckv	240
243	chuckv	246	if (.not.setSim ) then
244			if (present(status)) status = -1
245			return
246			end if
247
248			thisrlist = thisSim%rlist
249			if(present(rlist2)) rlist2 = thisSim%rlistsq
250
251			end subroutine getRlist
252
253
254
255			pure function getNlocal() result(nlocal)
256			integer :: nlocal
257			nlocal = thisSim%nLRparticles
258			end function getNlocal
259
260
261
262	chuckv	167	end module simulation