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 167 by chuckv, Fri Nov 8 14:21:52 2002 UTC vs.
Revision 246 by chuckv, Fri Jan 24 21:36:52 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 +  real( kind = dp ), allocatable,dimension(:,:),save :: q0
46  
47 +  public :: check
48 +  public :: save_nlist
49 +  public :: wrap
50 +  public :: getBox
51 +  public :: getRcut
52 +  public :: getRlist
53 +  public :: getNlocal
54 + !  public :: setRcut
55 +
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 +
68 +
69 +
70 +
71 +
72 +
73 +
74 +
75   contains
76  
77 +  subroutine setSimulation(nLRParticles,box,rlist,rcut)
78 +    integer, intent(in) :: nLRParticles
79 +    real(kind = dp ), intent(in), dimension(3) :: box
80 +    real(kind = dp ), intent(in) :: rlist
81 +    real(kind = dp ), intent(in) :: rcut
82  
83 <  subroutine set_simulation(box,rlist,rcut)
84 <    
83 >    if( setsim ) return  ! simulation is already initialized
84 >    setSim = .true.
85  
86 +    thisSim%nLRParticles = nLRParticles
87 +    thisSim%box          = box
88 +    thisSim%rlist        = rlist
89 +    thisSim%rcut         = rcut
90 +    thisSim%rcutsq       = rcut * rcut
91 +    thisSim%rcut6        = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
92  
93 <  end subroutine set_simulation
93 >  end subroutine setSimulation
94  
95 +  function getNparticles() result(nparticles)
96 +    integer :: nparticles
97 +    nparticles = thisSim%nLRparticles
98 +  end function getNparticles
99  
100  
101    subroutine change_box_size(new_box_size)
102      real(kind=dp), dimension(3) :: new_box_size
103  
104 <    box = new_box_size
104 >    thisSim%box = new_box_size
105  
106    end subroutine change_box_size
107  
108  
109 +  function getBox_3d() result(thisBox)
110 +    real( kind = dp ), dimension(3) :: thisBox
111 +    thisBox = thisSim%box
112 +  end function getBox_3d
113 +
114 +  function getBox_dim(dim) result(thisBox)
115 +    integer, intent(in) :: dim
116 +    real( kind = dp ) :: thisBox
117 +
118 +    thisBox = thisSim%box(dim)
119 +  end function getBox_dim
120 +    
121 +  subroutine check(q,update_nlist)
122 +    real( kind = dp ), dimension(:,:)  :: q
123 +    integer :: i
124 +    real( kind = DP ) :: dispmx
125 +    logical, intent(out) :: update_nlist
126 +    real( kind = DP ) :: dispmx_tmp
127 +    real( kind = dp ) :: skin_thickness
128 +    integer :: natoms
129 +
130 +    natoms = thisSim%nLRparticles
131 +    skin_thickness = thisSim%rlist
132 +    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 +  subroutine save_nlist(q)
160 +    real(kind = dp ), dimension(:,:), intent(in)  :: q
161 +    integer :: list_size
162 +    
163 +
164 +    list_size = size(q)
165 +
166 +    if (.not. allocated(q0)) then
167 +       allocate(q0(3,list_size))
168 +    else if( list_size > size(q0)) then
169 +       deallocate(q0)
170 +       allocate(q0(3,list_size))
171 +    endif
172 +
173 +    q0 = q
174 +
175 +  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 +       this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
188 +    else
189 +       this_wrap = r
190 +    endif
191 +    
192 +    return
193 +  end function wrap_1d
194 +
195 +  function wrap_3d(r) result(this_wrap)
196 +    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 +       this_wrap = r - thisSim%box*nint(r/thisSim%box)
203 +    else
204 +       this_wrap = r
205 +    endif
206 +  end function wrap_3d
207 +
208 +  
209 +
210 +  subroutine getRcut(thisrcut,rcut2,rcut6,status)
211 +    real( kind = dp ), intent(out) :: thisrcut
212 +    real( kind = dp ), intent(out), optional :: rcut2
213 +    real( kind = dp ), intent(out), optional :: rcut6
214 +    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 +    thisrcut = thisSim%rcut
224 +    if(present(rcut2)) rcut2 = thisSim%rcutsq
225 +    if(present(rcut2)) rcut6 = thisSim%rcut6
226 +
227 +  end subroutine getRcut
228 +  
229 +  
230 +  
231 +
232 +  subroutine getRlist(thisrlist,rlist2,status)
233 +    real( kind = dp ), intent(out) :: thisrlist
234 +    real( kind = dp ), intent(out), optional :: rlist2
235 +
236 +    integer, optional :: status
237 +
238 +    if (present(status)) status = 0
239 +
240 +    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   end module simulation

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines