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 :: 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
  

  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:	284
Committed:	Tue Feb 25 21:30:09 2003 UTC (21 years, 4 months ago) by chuckv
File size:	5505 byte(s)
Log Message:	Added neighbor list module that auto expands...
#	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 :: wrap
30			public :: getBox
31	chuckv	230	public :: getRcut
32	chuckv	246	public :: getRlist
33			public :: getNlocal
34	chuckv	249	public :: setSimulation
35	chuckv	281	public :: isEnsemble
36			public :: isPBC
37			public :: getStringLen
38			public :: returnMixingRules
39
40	chuckv	246	! public :: setRcut
41	chuckv	222
42			interface wrap
43			module procedure wrap_1d
44			module procedure wrap_3d
45			end interface
46
47			interface getBox
48			module procedure getBox_3d
49			module procedure getBox_dim
50			end interface
51
52
53	chuckv	185
54
55
56	chuckv	194
57
58
59
60
61	chuckv	167	contains
62
63	chuckv	281	subroutine setSimulation(nLRParticles,box,rlist,rcut,ensemble,mixingRule,use_pbc)
64	chuckv	230	integer, intent(in) :: nLRParticles
65	mmeineke	241	real(kind = dp ), intent(in), dimension(3) :: box
66	chuckv	230	real(kind = dp ), intent(in) :: rlist
67			real(kind = dp ), intent(in) :: rcut
68	chuckv	281	character( len = stringLen), intent(in) :: ensemble
69			character( len = stringLen), intent(in) :: mixingRule
70			logical, intent(in) :: use_pbc
71	chuckv	252	integer :: alloc_stat
72	chuckv	230	if( setsim ) return ! simulation is already initialized
73			setSim = .true.
74	chuckv	194
75	chuckv	230	thisSim%nLRParticles = nLRParticles
76			thisSim%box = box
77			thisSim%rlist = rlist
78	chuckv	253	thisSIm%rlistsq = rlist * rlist
79	chuckv	230	thisSim%rcut = rcut
80			thisSim%rcutsq = rcut * rcut
81	chuckv	246	thisSim%rcut6 = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
82	chuckv	252
83	chuckv	281	thisSim%ensemble = ensemble
84			thisSim%mixingRule = mixingRule
85			thisSim%use_pbc = use_pbc
86
87	chuckv	252	if (.not. allocated(q0)) then
88			allocate(q0(3,nLRParticles),stat=alloc_stat)
89			endif
90	chuckv	230	end subroutine setSimulation
91	chuckv	194
92	chuckv	230	function getNparticles() result(nparticles)
93			integer :: nparticles
94			nparticles = thisSim%nLRparticles
95			end function getNparticles
96	chuckv	194
97
98	chuckv	167	subroutine change_box_size(new_box_size)
99			real(kind=dp), dimension(3) :: new_box_size
100
101	chuckv	230	thisSim%box = new_box_size
102	chuckv	167
103			end subroutine change_box_size
104
105	chuckv	222
106	chuckv	246	function getBox_3d() result(thisBox)
107	chuckv	222	real( kind = dp ), dimension(3) :: thisBox
108	chuckv	230	thisBox = thisSim%box
109	chuckv	222	end function getBox_3d
110
111			function getBox_dim(dim) result(thisBox)
112			integer, intent(in) :: dim
113			real( kind = dp ) :: thisBox
114
115	chuckv	230	thisBox = thisSim%box(dim)
116	chuckv	222	end function getBox_dim
117
118
119			function wrap_1d(r,dim) result(this_wrap)
120
121
122			real( kind = DP ) :: r
123			real( kind = DP ) :: this_wrap
124			integer :: dim
125
126			if (use_pbc) then
127			! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
128	chuckv	230	this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
129	chuckv	222	else
130			this_wrap = r
131			endif
132
133			return
134			end function wrap_1d
135
136	chuckv	246	function wrap_3d(r) result(this_wrap)
137	chuckv	222	real( kind = dp ), dimension(3), intent(in) :: r
138			real( kind = dp ), dimension(3) :: this_wrap
139
140
141	chuckv	281	if (this_sim%use_pbc) then
142	chuckv	222	! this_wrap = r - box(dim)dsign(1.0E0_DP,r)int(abs(r/box(dim)) + 0.5E0_DP)
143	chuckv	246	this_wrap = r - thisSim%box*nint(r/thisSim%box)
144	chuckv	222	else
145			this_wrap = r
146			endif
147			end function wrap_3d
148
149	chuckv	246
150
151	chuckv	230	subroutine getRcut(thisrcut,rcut2,rcut6,status)
152			real( kind = dp ), intent(out) :: thisrcut
153			real( kind = dp ), intent(out), optional :: rcut2
154	chuckv	246	real( kind = dp ), intent(out), optional :: rcut6
155	chuckv	230	integer, optional :: status
156
157			if (present(status)) status = 0
158
159			if (.not.setSim ) then
160			if (present(status)) status = -1
161			return
162			end if
163
164	chuckv	246	thisrcut = thisSim%rcut
165			if(present(rcut2)) rcut2 = thisSim%rcutsq
166	chuckv	252	if(present(rcut6)) rcut6 = thisSim%rcut6
167	chuckv	230
168			end subroutine getRcut
169	chuckv	222
170	chuckv	246
171
172
173			subroutine getRlist(thisrlist,rlist2,status)
174			real( kind = dp ), intent(out) :: thisrlist
175			real( kind = dp ), intent(out), optional :: rlist2
176	chuckv	240
177	chuckv	246	integer, optional :: status
178	chuckv	240
179	chuckv	246	if (present(status)) status = 0
180	chuckv	240
181	chuckv	246	if (.not.setSim ) then
182			if (present(status)) status = -1
183			return
184			end if
185
186			thisrlist = thisSim%rlist
187			if(present(rlist2)) rlist2 = thisSim%rlistsq
188
189	chuckv	253
190	chuckv	246	end subroutine getRlist
191
192
193
194			pure function getNlocal() result(nlocal)
195			integer :: nlocal
196			nlocal = thisSim%nLRparticles
197			end function getNlocal
198
199
200	chuckv	281	function isEnsemble(this_ensemble) result(is_this_ensemble)
201			character(len = *) :: this_ensemble
202			logical :: is_this_enemble
203			is_this_ensemble = .false.
204			if (this_ensemble == thisSim%ensemble) is_this_ensemble = .true.
205			end function isEnsemble
206	chuckv	246
207	chuckv	281	function returnEnsemble() result(thisEnsemble)
208			character (len = len(thisSim%ensemble)) :: thisEnsemble
209			thisEnsemble = thisSim%ensemble
210			end function returnEnsemble
211
212			function returnMixingRules() result(thisMixingRule)
213			character (len = len(thisSim%ensemble)) :: thisMixingRule
214			thisMixingRule = thisSim%MixingRule
215			end function returnMixingRules
216
217			function isPBC() result(PBCset)
218			logical :: PBCset
219			PBCset = .false.
220			if (thisSim%use_pbc) PBCset = .true.
221			end function isPBC
222
223			pure function getStringLen() result (thislen)
224			integer :: thislen
225			thislen = string_len
226			end function setStringLen
227
228	chuckv	167	end module simulation