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
    integer :: alloc_stat
    if( setsim ) return  ! simulation is already initialized
    setSim = .true.

    thisSim%nLRParticles = nLRParticles
    thisSim%box          = box
    thisSim%rlist        = rlist
    thisSIm%rlistsq      = 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


  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 - thisSim%rcut
    dispmx = 0.0E0_DP
    !! calculate the largest displacement of any atom in any direction
    
#ifdef MPI 
    dispmx_tmp = 0.0E0_DP
    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 )
    end do
    call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, & 
       mpi_max,mpi_comm_world,mpi_err)
#else

    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 )
    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(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
Revision:	264
Committed:	Tue Feb 4 20:16:08 2003 UTC (22 years, 9 months ago) by chuckv
File size:	6374 byte(s)
Log Message:	Fixed bug with pot energy.
#	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	252	integer :: alloc_stat
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	chuckv	253	thisSIm%rlistsq = rlist * rlist
93	chuckv	230	thisSim%rcut = rcut
94			thisSim%rcutsq = rcut * rcut
95	chuckv	246	thisSim%rcut6 = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
96	chuckv	252
97			if (.not. allocated(q0)) then
98			allocate(q0(3,nLRParticles),stat=alloc_stat)
99			endif
100	chuckv	230	end subroutine setSimulation
101	chuckv	194
102	chuckv	230	function getNparticles() result(nparticles)
103			integer :: nparticles
104			nparticles = thisSim%nLRparticles
105			end function getNparticles
106	chuckv	194
107
108	chuckv	167	subroutine change_box_size(new_box_size)
109			real(kind=dp), dimension(3) :: new_box_size
110
111	chuckv	230	thisSim%box = new_box_size
112	chuckv	167
113			end subroutine change_box_size
114
115	chuckv	222
116	chuckv	246	function getBox_3d() result(thisBox)
117	chuckv	222	real( kind = dp ), dimension(3) :: thisBox
118	chuckv	230	thisBox = thisSim%box
119	chuckv	222	end function getBox_3d
120
121			function getBox_dim(dim) result(thisBox)
122			integer, intent(in) :: dim
123			real( kind = dp ) :: thisBox
124
125	chuckv	230	thisBox = thisSim%box(dim)
126	chuckv	222	end function getBox_dim
127
128	chuckv	246	subroutine check(q,update_nlist)
129			real( kind = dp ), dimension(:,:) :: q
130	chuckv	222	integer :: i
131			real( kind = DP ) :: dispmx
132			logical, intent(out) :: update_nlist
133			real( kind = DP ) :: dispmx_tmp
134	chuckv	246	real( kind = dp ) :: skin_thickness
135			integer :: natoms
136
137			natoms = thisSim%nLRparticles
138	chuckv	264	skin_thickness = thisSim%rlist - thisSim%rcut
139	chuckv	222	dispmx = 0.0E0_DP
140			!! calculate the largest displacement of any atom in any direction
141
142			#ifdef MPI
143			dispmx_tmp = 0.0E0_DP
144	chuckv	252	do i = 1, thisSim%nLRparticles
145	chuckv	222	dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
146			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
147			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
148			end do
149			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
150			mpi_max,mpi_comm_world,mpi_err)
151			#else
152	chuckv	252
153			do i = 1, thisSim%nLRparticles
154	chuckv	222	dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
155			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
156			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
157			end do
158			#endif
159
160			!! a conservative test of list skin crossings
161			dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
162
163			update_nlist = (dispmx.gt.(skin_thickness))
164
165			end subroutine check
166
167	chuckv	230	subroutine save_nlist(q)
168			real(kind = dp ), dimension(:,:), intent(in) :: q
169			integer :: list_size
170
171	chuckv	222
172	chuckv	230	list_size = size(q)
173
174			if (.not. allocated(q0)) then
175			allocate(q0(3,list_size))
176	chuckv	246	else if( list_size > size(q0)) then
177	chuckv	230	deallocate(q0)
178			allocate(q0(3,list_size))
179			endif
180
181			q0 = q
182
183	chuckv	222	end subroutine save_nlist
184
185
186			function wrap_1d(r,dim) result(this_wrap)
187
188
189			real( kind = DP ) :: r
190			real( kind = DP ) :: this_wrap
191			integer :: dim
192
193			if (use_pbc) then
194			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
195	chuckv	230	this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
196	chuckv	222	else
197			this_wrap = r
198			endif
199
200			return
201			end function wrap_1d
202
203	chuckv	246	function wrap_3d(r) result(this_wrap)
204	chuckv	222	real( kind = dp ), dimension(3), intent(in) :: r
205			real( kind = dp ), dimension(3) :: this_wrap
206
207
208			if (use_pbc) then
209			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
210	chuckv	246	this_wrap = r - thisSim%box*nint(r/thisSim%box)
211	chuckv	222	else
212			this_wrap = r
213			endif
214			end function wrap_3d
215
216	chuckv	246
217
218	chuckv	230	subroutine getRcut(thisrcut,rcut2,rcut6,status)
219			real( kind = dp ), intent(out) :: thisrcut
220			real( kind = dp ), intent(out), optional :: rcut2
221	chuckv	246	real( kind = dp ), intent(out), optional :: rcut6
222	chuckv	230	integer, optional :: status
223
224			if (present(status)) status = 0
225
226			if (.not.setSim ) then
227			if (present(status)) status = -1
228			return
229			end if
230
231	chuckv	246	thisrcut = thisSim%rcut
232			if(present(rcut2)) rcut2 = thisSim%rcutsq
233	chuckv	252	if(present(rcut6)) rcut6 = thisSim%rcut6
234	chuckv	230
235			end subroutine getRcut
236	chuckv	222
237	chuckv	246
238
239
240			subroutine getRlist(thisrlist,rlist2,status)
241			real( kind = dp ), intent(out) :: thisrlist
242			real( kind = dp ), intent(out), optional :: rlist2
243	chuckv	240
244	chuckv	246	integer, optional :: status
245	chuckv	240
246	chuckv	246	if (present(status)) status = 0
247	chuckv	240
248	chuckv	246	if (.not.setSim ) then
249			if (present(status)) status = -1
250			return
251			end if
252
253			thisrlist = thisSim%rlist
254			if(present(rlist2)) rlist2 = thisSim%rlistsq
255
256	chuckv	253
257	chuckv	246	end subroutine getRlist
258
259
260
261			pure function getNlocal() result(nlocal)
262			integer :: nlocal
263			nlocal = thisSim%nLRparticles
264			end function getNlocal
265
266
267
268	chuckv	167	end module simulation