ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/simulation_module.F90
(Generate patch)

Comparing trunk/mdtools/md_code/simulation_module.F90 (file contents):
Revision 185 by chuckv, Fri Nov 22 20:39:33 2002 UTC vs.
Revision 253 by chuckv, Thu Jan 30 15:20:21 2003 UTC

# Line 1 | Line 1
1   module simulation
2    use definitions, ONLY :dp
3 + #ifdef IS_MPI
4 +  use mpiSimulation
5 + #endif
6 +
7    implicit none
8    PRIVATE
9  
10  
11  
12 <  real ( kind = dp ), dimension(3) :: box
12 >  type, public :: simtype
13 >     PRIVATE
14 > !     SEQUENCE
15 > !! 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  
30 +  end type simtype
31  
32 +  type (simtype), public :: thisSim
33 + !! Tag for MPI calculations  
34 +  integer, allocatable, dimension(:) :: tag
35  
36 + #ifdef IS_MPI
37 +  integer, allocatable, dimension(:) :: tag_row
38 +  integer, allocatable, dimension(:) :: tag_column
39 + #endif
40  
41 + !! WARNING: use_pbc hardcoded, fixme
42 +  logical :: use_pbc = .true.
43 +  logical :: setSim = .false.
44  
45 + !! array for saving previous positions for neighbor lists.  
46 +  real( kind = dp ), allocatable,dimension(:,:),save :: q0
47  
48  
49 < contains
49 >  public :: check
50 >  public :: save_nlist
51 >  public :: wrap
52 >  public :: getBox
53 >  public :: getRcut
54 >  public :: getRlist
55 >  public :: getNlocal
56 >  public :: setSimulation
57 > !  public :: setRcut
58  
59 +  interface wrap
60 +     module procedure wrap_1d
61 +     module procedure wrap_3d
62 +  end interface
63  
64 <  subroutine set_simulation(box,rlist,rcut)
65 <    
64 >  interface getBox
65 >     module procedure getBox_3d
66 >     module procedure getBox_dim
67 >  end interface
68  
69  
23  end subroutine set_simulation
70  
71  
72  
73 +
74 +
75 +
76 +
77 +
78 + contains
79 +
80 +  subroutine setSimulation(nLRParticles,box,rlist,rcut)
81 +    integer, intent(in) :: nLRParticles
82 +    real(kind = dp ), intent(in), dimension(3) :: box
83 +    real(kind = dp ), intent(in) :: rlist
84 +    real(kind = dp ), intent(in) :: rcut
85 +    integer :: alloc_stat
86 +    if( setsim ) return  ! simulation is already initialized
87 +    setSim = .true.
88 +
89 +    thisSim%nLRParticles = nLRParticles
90 +    thisSim%box          = box
91 +    thisSim%rlist        = rlist
92 +    thisSIm%rlistsq      = rlist * rlist
93 +    thisSim%rcut         = rcut
94 +    thisSim%rcutsq       = rcut * rcut
95 +    thisSim%rcut6        = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
96 +    
97 +    if (.not. allocated(q0)) then
98 +       allocate(q0(3,nLRParticles),stat=alloc_stat)
99 +    endif
100 +  end subroutine setSimulation
101 +
102 +  function getNparticles() result(nparticles)
103 +    integer :: nparticles
104 +    nparticles = thisSim%nLRparticles
105 +  end function getNparticles
106 +
107 +
108    subroutine change_box_size(new_box_size)
109      real(kind=dp), dimension(3) :: new_box_size
110  
111 <    box = new_box_size
111 >    thisSim%box = new_box_size
112  
113    end subroutine change_box_size
114  
115  
116 +  function getBox_3d() result(thisBox)
117 +    real( kind = dp ), dimension(3) :: thisBox
118 +    thisBox = thisSim%box
119 +  end function getBox_3d
120 +
121 +  function getBox_dim(dim) result(thisBox)
122 +    integer, intent(in) :: dim
123 +    real( kind = dp ) :: thisBox
124 +
125 +    thisBox = thisSim%box(dim)
126 +  end function getBox_dim
127 +    
128 +  subroutine check(q,update_nlist)
129 +    real( kind = dp ), dimension(:,:)  :: q
130 +    integer :: i
131 +    real( kind = DP ) :: dispmx
132 +    logical, intent(out) :: update_nlist
133 +    real( kind = DP ) :: dispmx_tmp
134 +    real( kind = dp ) :: skin_thickness
135 +    integer :: natoms
136 +
137 +    natoms = thisSim%nLRparticles
138 +    skin_thickness = thisSim%rlist
139 +    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 +    do i = 1, thisSim%nLRparticles
145 +       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 +
153 +    do i = 1, thisSim%nLRparticles
154 +       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 +  subroutine save_nlist(q)
168 +    real(kind = dp ), dimension(:,:), intent(in)  :: q
169 +    integer :: list_size
170 +    
171 +
172 +    list_size = size(q)
173 +
174 +    if (.not. allocated(q0)) then
175 +       allocate(q0(3,list_size))
176 +    else if( list_size > size(q0)) then
177 +       deallocate(q0)
178 +       allocate(q0(3,list_size))
179 +    endif
180 +
181 +    q0 = q
182 +
183 +  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 +       this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
196 +    else
197 +       this_wrap = r
198 +    endif
199 +    
200 +    return
201 +  end function wrap_1d
202 +
203 +  function wrap_3d(r) result(this_wrap)
204 +    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 +       this_wrap = r - thisSim%box*nint(r/thisSim%box)
211 +    else
212 +       this_wrap = r
213 +    endif
214 +  end function wrap_3d
215 +
216 +  
217 +
218 +  subroutine getRcut(thisrcut,rcut2,rcut6,status)
219 +    real( kind = dp ), intent(out) :: thisrcut
220 +    real( kind = dp ), intent(out), optional :: rcut2
221 +    real( kind = dp ), intent(out), optional :: rcut6
222 +    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 +    thisrcut = thisSim%rcut
232 +    if(present(rcut2)) rcut2 = thisSim%rcutsq
233 +    if(present(rcut6)) rcut6 = thisSim%rcut6
234 +
235 +  end subroutine getRcut
236 +  
237 +  
238 +  
239 +
240 +  subroutine getRlist(thisrlist,rlist2,status)
241 +    real( kind = dp ), intent(out) :: thisrlist
242 +    real( kind = dp ), intent(out), optional :: rlist2
243 +
244 +    integer, optional :: status
245 +
246 +    if (present(status)) status = 0
247 +
248 +    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 +
257 +  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   end module simulation

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines