mdtools/md_code/simulation_module.F90

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

  implicit none
  PRIVATE

#define __FORTRAN90
#include "../headers/fsimulation.h"

  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 :: 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 :: isEnsemble
  public :: isPBC
  public :: getStringLen
  public :: returnMixingRules

!  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,ensemble,mixingRule,use_pbc)
    integer, intent(in) :: nLRParticles 
    real(kind = dp ), intent(in), dimension(3) :: box
    real(kind = dp ), intent(in) :: rlist
    real(kind = dp ), intent(in) :: rcut
    character( len = stringLen), intent(in)  :: ensemble
    character( len = stringLen), intent(in)  :: mixingRule
    logical, intent(in) :: use_pbc
    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
    
    thisSim%ensemble = ensemble
    thisSim%mixingRule = mixingRule
    thisSim%use_pbc = use_pbc

    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 (this_sim%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


  function isEnsemble(this_ensemble) result(is_this_ensemble)
    character(len = *) :: this_ensemble
    logical :: is_this_enemble
    is_this_ensemble = .false.
    if (this_ensemble == thisSim%ensemble) is_this_ensemble = .true.
  end function isEnsemble

  function returnEnsemble() result(thisEnsemble)
    character (len = len(thisSim%ensemble)) :: thisEnsemble
    thisEnsemble = thisSim%ensemble
  end function returnEnsemble

  function returnMixingRules() result(thisMixingRule)
    character (len = len(thisSim%ensemble)) :: thisMixingRule
    thisMixingRule = thisSim%MixingRule
  end function returnMixingRules

  function isPBC() result(PBCset)
    logical :: PBCset
    PBCset = .false.
    if (thisSim%use_pbc) PBCset = .true.
  end function isPBC

  pure function getStringLen() result (thislen)
    integer :: thislen    
    thislen = string_len
  end function setStringLen

end module simulation
Revision:	281
Committed:	Mon Feb 24 21:25:15 2003 UTC (21 years, 4 months ago) by chuckv
File size:	7148 byte(s)
Log Message:	Changes they are a comming....
#	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	chuckv	281	#define __FORTRAN90
11			#include "../headers/fsimulation.h"
12	chuckv	167
13	chuckv	230	type (simtype), public :: thisSim
14			!! Tag for MPI calculations
15			integer, allocatable, dimension(:) :: tag
16	chuckv	222
17	chuckv	230	#ifdef IS_MPI
18			integer, allocatable, dimension(:) :: tag_row
19			integer, allocatable, dimension(:) :: tag_column
20			#endif
21	chuckv	167
22	chuckv	246	!! WARNING: use_pbc hardcoded, fixme
23	chuckv	281	logical :: setSim = .false.
24	chuckv	247
25			!! array for saving previous positions for neighbor lists.
26	chuckv	246	real( kind = dp ), allocatable,dimension(:,:),save :: q0
27	chuckv	167
28	chuckv	247
29	chuckv	222	public :: check
30			public :: save_nlist
31			public :: wrap
32			public :: getBox
33	chuckv	230	public :: getRcut
34	chuckv	246	public :: getRlist
35			public :: getNlocal
36	chuckv	249	public :: setSimulation
37	chuckv	281	public :: isEnsemble
38			public :: isPBC
39			public :: getStringLen
40			public :: returnMixingRules
41
42	chuckv	246	! public :: setRcut
43	chuckv	222
44			interface wrap
45			module procedure wrap_1d
46			module procedure wrap_3d
47			end interface
48
49			interface getBox
50			module procedure getBox_3d
51			module procedure getBox_dim
52			end interface
53
54
55	chuckv	185
56
57
58	chuckv	194
59
60
61
62
63	chuckv	167	contains
64
65	chuckv	281	subroutine setSimulation(nLRParticles,box,rlist,rcut,ensemble,mixingRule,use_pbc)
66	chuckv	230	integer, intent(in) :: nLRParticles
67	mmeineke	241	real(kind = dp ), intent(in), dimension(3) :: box
68	chuckv	230	real(kind = dp ), intent(in) :: rlist
69			real(kind = dp ), intent(in) :: rcut
70	chuckv	281	character( len = stringLen), intent(in) :: ensemble
71			character( len = stringLen), intent(in) :: mixingRule
72			logical, intent(in) :: use_pbc
73	chuckv	252	integer :: alloc_stat
74	chuckv	230	if( setsim ) return ! simulation is already initialized
75			setSim = .true.
76	chuckv	194
77	chuckv	230	thisSim%nLRParticles = nLRParticles
78			thisSim%box = box
79			thisSim%rlist = rlist
80	chuckv	253	thisSIm%rlistsq = rlist * rlist
81	chuckv	230	thisSim%rcut = rcut
82			thisSim%rcutsq = rcut * rcut
83	chuckv	246	thisSim%rcut6 = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
84	chuckv	252
85	chuckv	281	thisSim%ensemble = ensemble
86			thisSim%mixingRule = mixingRule
87			thisSim%use_pbc = use_pbc
88
89	chuckv	252	if (.not. allocated(q0)) then
90			allocate(q0(3,nLRParticles),stat=alloc_stat)
91			endif
92	chuckv	230	end subroutine setSimulation
93	chuckv	194
94	chuckv	230	function getNparticles() result(nparticles)
95			integer :: nparticles
96			nparticles = thisSim%nLRparticles
97			end function getNparticles
98	chuckv	194
99
100	chuckv	167	subroutine change_box_size(new_box_size)
101			real(kind=dp), dimension(3) :: new_box_size
102
103	chuckv	230	thisSim%box = new_box_size
104	chuckv	167
105			end subroutine change_box_size
106
107	chuckv	222
108	chuckv	246	function getBox_3d() result(thisBox)
109	chuckv	222	real( kind = dp ), dimension(3) :: thisBox
110	chuckv	230	thisBox = thisSim%box
111	chuckv	222	end function getBox_3d
112
113			function getBox_dim(dim) result(thisBox)
114			integer, intent(in) :: dim
115			real( kind = dp ) :: thisBox
116
117	chuckv	230	thisBox = thisSim%box(dim)
118	chuckv	222	end function getBox_dim
119
120	chuckv	246	subroutine check(q,update_nlist)
121			real( kind = dp ), dimension(:,:) :: q
122	chuckv	222	integer :: i
123			real( kind = DP ) :: dispmx
124			logical, intent(out) :: update_nlist
125			real( kind = DP ) :: dispmx_tmp
126	chuckv	246	real( kind = dp ) :: skin_thickness
127			integer :: natoms
128
129			natoms = thisSim%nLRparticles
130	chuckv	264	skin_thickness = thisSim%rlist - thisSim%rcut
131	chuckv	222	dispmx = 0.0E0_DP
132			!! calculate the largest displacement of any atom in any direction
133
134			#ifdef MPI
135			dispmx_tmp = 0.0E0_DP
136	chuckv	252	do i = 1, thisSim%nLRparticles
137	chuckv	222	dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
138			dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
139			dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
140			end do
141			call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
142			mpi_max,mpi_comm_world,mpi_err)
143			#else
144	chuckv	252
145			do i = 1, thisSim%nLRparticles
146	chuckv	222	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	chuckv	281	if (this_sim%use_pbc) then
201	chuckv	222	! 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	chuckv	252	if(present(rcut6)) 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	chuckv	253
249	chuckv	246	end subroutine getRlist
250
251
252
253			pure function getNlocal() result(nlocal)
254			integer :: nlocal
255			nlocal = thisSim%nLRparticles
256			end function getNlocal
257
258
259	chuckv	281	function isEnsemble(this_ensemble) result(is_this_ensemble)
260			character(len = *) :: this_ensemble
261			logical :: is_this_enemble
262			is_this_ensemble = .false.
263			if (this_ensemble == thisSim%ensemble) is_this_ensemble = .true.
264			end function isEnsemble
265	chuckv	246
266	chuckv	281	function returnEnsemble() result(thisEnsemble)
267			character (len = len(thisSim%ensemble)) :: thisEnsemble
268			thisEnsemble = thisSim%ensemble
269			end function returnEnsemble
270
271			function returnMixingRules() result(thisMixingRule)
272			character (len = len(thisSim%ensemble)) :: thisMixingRule
273			thisMixingRule = thisSim%MixingRule
274			end function returnMixingRules
275
276			function isPBC() result(PBCset)
277			logical :: PBCset
278			PBCset = .false.
279			if (thisSim%use_pbc) PBCset = .true.
280			end function isPBC
281
282			pure function getStringLen() result (thislen)
283			integer :: thislen
284			thislen = string_len
285			end function setStringLen
286
287	chuckv	167	end module simulation