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.
  
  real( kind = dp ), allocatable,dimension(:,:),save :: q0

  public :: check
  public :: save_nlist
  public :: wrap
  public :: getBox
  public :: getRcut
  public :: getRlist
  public :: getNlocal
!  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:	246
Committed:	Fri Jan 24 21:36:52 2003 UTC (22 years, 9 months ago) by chuckv
File size:	6077 byte(s)
Log Message:	lj_FF.F90 and simulation_module now compile. Logical bug still exists in lj_FF.F90 lj_atypes_list.
#	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
45	chuckv	246	real( kind = dp ), allocatable,dimension(:,:),save :: q0
46	chuckv	167
47	chuckv	222	public :: check
48			public :: save_nlist
49			public :: wrap
50			public :: getBox
51	chuckv	230	public :: getRcut
52	chuckv	246	public :: getRlist
53			public :: getNlocal
54			! public :: setRcut
55	chuckv	222
56			interface wrap
57			module procedure wrap_1d
58			module procedure wrap_3d
59			end interface
60
61			interface getBox
62			module procedure getBox_3d
63			module procedure getBox_dim
64			end interface
65
66
67	chuckv	185
68
69
70	chuckv	194
71
72
73
74
75	chuckv	167	contains
76
77	chuckv	230	subroutine setSimulation(nLRParticles,box,rlist,rcut)
78			integer, intent(in) :: nLRParticles
79	mmeineke	241	real(kind = dp ), intent(in), dimension(3) :: box
80	chuckv	230	real(kind = dp ), intent(in) :: rlist
81			real(kind = dp ), intent(in) :: rcut
82	chuckv	167
83	chuckv	230	if( setsim ) return ! simulation is already initialized
84			setSim = .true.
85	chuckv	194
86	chuckv	230	thisSim%nLRParticles = nLRParticles
87			thisSim%box = box
88			thisSim%rlist = rlist
89			thisSim%rcut = rcut
90			thisSim%rcutsq = rcut * rcut
91	chuckv	246	thisSim%rcut6 = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
92	chuckv	194
93	chuckv	230	end subroutine setSimulation
94	chuckv	194
95	chuckv	230	function getNparticles() result(nparticles)
96			integer :: nparticles
97			nparticles = thisSim%nLRparticles
98			end function getNparticles
99	chuckv	194
100
101	chuckv	167	subroutine change_box_size(new_box_size)
102			real(kind=dp), dimension(3) :: new_box_size
103
104	chuckv	230	thisSim%box = new_box_size
105	chuckv	167
106			end subroutine change_box_size
107
108	chuckv	222
109	chuckv	246	function getBox_3d() result(thisBox)
110	chuckv	222	real( kind = dp ), dimension(3) :: thisBox
111	chuckv	230	thisBox = thisSim%box
112	chuckv	222	end function getBox_3d
113
114			function getBox_dim(dim) result(thisBox)
115			integer, intent(in) :: dim
116			real( kind = dp ) :: thisBox
117
118	chuckv	230	thisBox = thisSim%box(dim)
119	chuckv	222	end function getBox_dim
120
121	chuckv	246	subroutine check(q,update_nlist)
122			real( kind = dp ), dimension(:,:) :: q
123	chuckv	222	integer :: i
124			real( kind = DP ) :: dispmx
125			logical, intent(out) :: update_nlist
126			real( kind = DP ) :: dispmx_tmp
127	chuckv	246	real( kind = dp ) :: skin_thickness
128			integer :: natoms
129
130			natoms = thisSim%nLRparticles
131			skin_thickness = thisSim%rlist
132	chuckv	222	dispmx = 0.0E0_DP
133			!! calculate the largest displacement of any atom in any direction
134
135			#ifdef MPI
136			dispmx_tmp = 0.0E0_DP
137			do i = 1, nlocal
138			dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
139			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
140			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
141			end do
142			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
143			mpi_max,mpi_comm_world,mpi_err)
144			#else
145			do i = 1, natoms
146			dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
147			dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
148			dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
149			end do
150			#endif
151
152			!! a conservative test of list skin crossings
153			dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
154
155			update_nlist = (dispmx.gt.(skin_thickness))
156
157			end subroutine check
158
159	chuckv	230	subroutine save_nlist(q)
160			real(kind = dp ), dimension(:,:), intent(in) :: q
161			integer :: list_size
162
163	chuckv	222
164	chuckv	230	list_size = size(q)
165
166			if (.not. allocated(q0)) then
167			allocate(q0(3,list_size))
168	chuckv	246	else if( list_size > size(q0)) then
169	chuckv	230	deallocate(q0)
170			allocate(q0(3,list_size))
171			endif
172
173			q0 = q
174
175	chuckv	222	end subroutine save_nlist
176
177
178			function wrap_1d(r,dim) result(this_wrap)
179
180
181			real( kind = DP ) :: r
182			real( kind = DP ) :: this_wrap
183			integer :: dim
184
185			if (use_pbc) then
186			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
187	chuckv	230	this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
188	chuckv	222	else
189			this_wrap = r
190			endif
191
192			return
193			end function wrap_1d
194
195	chuckv	246	function wrap_3d(r) result(this_wrap)
196	chuckv	222	real( kind = dp ), dimension(3), intent(in) :: r
197			real( kind = dp ), dimension(3) :: this_wrap
198
199
200			if (use_pbc) then
201			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
202	chuckv	246	this_wrap = r - thisSim%box*nint(r/thisSim%box)
203	chuckv	222	else
204			this_wrap = r
205			endif
206			end function wrap_3d
207
208	chuckv	246
209
210	chuckv	230	subroutine getRcut(thisrcut,rcut2,rcut6,status)
211			real( kind = dp ), intent(out) :: thisrcut
212			real( kind = dp ), intent(out), optional :: rcut2
213	chuckv	246	real( kind = dp ), intent(out), optional :: rcut6
214	chuckv	230	integer, optional :: status
215
216			if (present(status)) status = 0
217
218			if (.not.setSim ) then
219			if (present(status)) status = -1
220			return
221			end if
222
223	chuckv	246	thisrcut = thisSim%rcut
224			if(present(rcut2)) rcut2 = thisSim%rcutsq
225			if(present(rcut2)) rcut6 = thisSim%rcut6
226	chuckv	230
227			end subroutine getRcut
228	chuckv	222
229	chuckv	246
230
231
232			subroutine getRlist(thisrlist,rlist2,status)
233			real( kind = dp ), intent(out) :: thisrlist
234			real( kind = dp ), intent(out), optional :: rlist2
235	chuckv	240
236	chuckv	246	integer, optional :: status
237	chuckv	240
238	chuckv	246	if (present(status)) status = 0
239	chuckv	240
240	chuckv	246	if (.not.setSim ) then
241			if (present(status)) status = -1
242			return
243			end if
244
245			thisrlist = thisSim%rlist
246			if(present(rlist2)) rlist2 = thisSim%rlistsq
247
248			end subroutine getRlist
249
250
251
252			pure function getNlocal() result(nlocal)
253			integer :: nlocal
254			nlocal = thisSim%nLRparticles
255			end function getNlocal
256
257
258
259	chuckv	167	end module simulation