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 194 by chuckv, Wed Dec 4 21:19:38 2002 UTC vs.
Revision 249 by chuckv, Mon Jan 27 21:28:19 2003 UTC

# Line 1 | Line 1
1   module simulation
2    use definitions, ONLY :dp
3 <  use force_wrappers, ONLY : alloc_force_wrappers
3 > #ifdef IS_MPI
4 >  use mpiSimulation
5 > #endif
6  
7    implicit none
8    PRIVATE
9  
10  
9 –  
10 –  real ( kind = dp ), public, dimension(3) :: box
11  
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 <  real ( kind = dp ), public :: rlist
14 <  real ( kind = dp ), public :: rcut
15 <  real ( kind = dp ), public :: rlistsq
16 <  real ( kind = dp ), public :: rcutsq
30 >  end type simtype
31  
32 <  integer,public, allocatable, dimension(:)   :: point
33 <  integer,public, allocatable, dimension(:)   :: list
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 < !MPI dependent routines
41 > !! WARNING: use_pbc hardcoded, fixme
42 >  logical :: use_pbc = .true.
43 >  logical :: setSim = .false.
44  
45 < #ifdef IS_MPI
46 < ! Universal routines: All types of force calculations will need these arrays
26 < ! Arrays specific to a type of force calculation should be declared in that module.
27 <  real( kind = dp ), allocatable, dimension(:,:) :: qRow
28 <  real( kind = dp ), allocatable, dimension(:,:) :: qColumn
45 > !! array for saving previous positions for neighbor lists.  
46 >  real( kind = dp ), allocatable,dimension(:,:),save :: q0
47  
30 –  real( kind = dp ), allocatable, dimension(:,:) :: fRow
31 –  real( kind = dp ), allocatable, dimension(:,:) :: fColumn
48  
49 <  real( kind = dp ), allocatable, dimension(:,:) :: tRow
50 <  real( kind = dp ), allocatable, dimension(:,:) :: tColumn
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 < #endif
59 >  interface wrap
60 >     module procedure wrap_1d
61 >     module procedure wrap_3d
62 >  end interface
63  
64 +  interface getBox
65 +     module procedure getBox_3d
66 +     module procedure getBox_dim
67 +  end interface
68  
69  
70  
71  
72  
43 – contains
73  
45 – #ifdef MPI
46 – ! Allocated work arrays for MPI
47 –  subroutine allocate_mpi_arrays(nDimensions,numComponents)
48 –    integer, intent(in) :: nDimensions
49 –    integer, intent(in) :: numComponents
74  
75  
76  
77  
78 + contains
79  
80 <  end subroutine allocate_mpi_arrays
81 < #endif
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  
86 <  subroutine set_simulation(box,rlist,rcut)
87 <    
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%rcut         = rcut
93 +    thisSim%rcutsq       = rcut * rcut
94 +    thisSim%rcut6        = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
95  
96 <  end subroutine set_simulation
96 >  end subroutine setSimulation
97  
98 +  function getNparticles() result(nparticles)
99 +    integer :: nparticles
100 +    nparticles = thisSim%nLRparticles
101 +  end function getNparticles
102  
103  
104    subroutine change_box_size(new_box_size)
105      real(kind=dp), dimension(3) :: new_box_size
106  
107 <    box = new_box_size
107 >    thisSim%box = new_box_size
108  
109    end subroutine change_box_size
110  
111  
112 +  function getBox_3d() result(thisBox)
113 +    real( kind = dp ), dimension(3) :: thisBox
114 +    thisBox = thisSim%box
115 +  end function getBox_3d
116 +
117 +  function getBox_dim(dim) result(thisBox)
118 +    integer, intent(in) :: dim
119 +    real( kind = dp ) :: thisBox
120 +
121 +    thisBox = thisSim%box(dim)
122 +  end function getBox_dim
123 +    
124 +  subroutine check(q,update_nlist)
125 +    real( kind = dp ), dimension(:,:)  :: q
126 +    integer :: i
127 +    real( kind = DP ) :: dispmx
128 +    logical, intent(out) :: update_nlist
129 +    real( kind = DP ) :: dispmx_tmp
130 +    real( kind = dp ) :: skin_thickness
131 +    integer :: natoms
132 +
133 +    natoms = thisSim%nLRparticles
134 +    skin_thickness = thisSim%rlist
135 +    dispmx = 0.0E0_DP
136 +    !! calculate the largest displacement of any atom in any direction
137 +    
138 + #ifdef MPI
139 +    dispmx_tmp = 0.0E0_DP
140 +    do i = 1, nlocal
141 +       dispmx_tmp = max( abs ( q(1,i) - q0(1,i) ), dispmx )
142 +       dispmx_tmp = max( abs ( q(2,i) - q0(2,i) ), dispmx )
143 +       dispmx_tmp = max( abs ( q(3,i) - q0(3,i) ), dispmx )
144 +    end do
145 +    call mpi_allreduce(dispmx_tmp,dispmx,1,mpi_double_precision, &
146 +       mpi_max,mpi_comm_world,mpi_err)
147 + #else
148 +    do i = 1, natoms
149 +       dispmx = max( abs ( q(1,i) - q0(1,i) ), dispmx )
150 +       dispmx = max( abs ( q(2,i) - q0(2,i) ), dispmx )
151 +       dispmx = max( abs ( q(3,i) - q0(3,i) ), dispmx )
152 +    end do
153 + #endif
154 +  
155 +    !! a conservative test of list skin crossings
156 +    dispmx = 2.0E0_DP * sqrt (3.0E0_DP * dispmx * dispmx)
157 +    
158 +    update_nlist = (dispmx.gt.(skin_thickness))
159 +    
160 +  end subroutine check
161 +  
162 +  subroutine save_nlist(q)
163 +    real(kind = dp ), dimension(:,:), intent(in)  :: q
164 +    integer :: list_size
165 +    
166 +
167 +    list_size = size(q)
168 +
169 +    if (.not. allocated(q0)) then
170 +       allocate(q0(3,list_size))
171 +    else if( list_size > size(q0)) then
172 +       deallocate(q0)
173 +       allocate(q0(3,list_size))
174 +    endif
175 +
176 +    q0 = q
177 +
178 +  end subroutine save_nlist
179 +  
180 +
181 +  function wrap_1d(r,dim) result(this_wrap)
182 +    
183 +    
184 +    real( kind = DP ) :: r
185 +    real( kind = DP ) :: this_wrap
186 +    integer           :: dim
187 +    
188 +    if (use_pbc) then
189 +       !     this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP)
190 +       this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
191 +    else
192 +       this_wrap = r
193 +    endif
194 +    
195 +    return
196 +  end function wrap_1d
197 +
198 +  function wrap_3d(r) result(this_wrap)
199 +    real( kind = dp ), dimension(3), intent(in) :: r
200 +    real( kind = dp ), dimension(3) :: this_wrap
201 +
202 +    
203 +    if (use_pbc) then
204 +       !     this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP)
205 +       this_wrap = r - thisSim%box*nint(r/thisSim%box)
206 +    else
207 +       this_wrap = r
208 +    endif
209 +  end function wrap_3d
210 +
211 +  
212 +
213 +  subroutine getRcut(thisrcut,rcut2,rcut6,status)
214 +    real( kind = dp ), intent(out) :: thisrcut
215 +    real( kind = dp ), intent(out), optional :: rcut2
216 +    real( kind = dp ), intent(out), optional :: rcut6
217 +    integer, optional :: status
218 +
219 +    if (present(status)) status = 0
220 +
221 +    if (.not.setSim ) then
222 +       if (present(status)) status = -1
223 +       return
224 +    end if
225 +    
226 +    thisrcut = thisSim%rcut
227 +    if(present(rcut2)) rcut2 = thisSim%rcutsq
228 +    if(present(rcut2)) rcut6 = thisSim%rcut6
229 +
230 +  end subroutine getRcut
231 +  
232 +  
233 +  
234 +
235 +  subroutine getRlist(thisrlist,rlist2,status)
236 +    real( kind = dp ), intent(out) :: thisrlist
237 +    real( kind = dp ), intent(out), optional :: rlist2
238 +
239 +    integer, optional :: status
240 +
241 +    if (present(status)) status = 0
242 +
243 +    if (.not.setSim ) then
244 +       if (present(status)) status = -1
245 +       return
246 +    end if
247 +    
248 +    thisrlist = thisSim%rlist
249 +    if(present(rlist2)) rlist2 = thisSim%rlistsq
250 +
251 +  end subroutine getRlist
252 +  
253 +  
254 +
255 + pure function getNlocal() result(nlocal)
256 +    integer :: nlocal
257 +    nlocal = thisSim%nLRparticles
258 +  end function getNlocal
259 +
260 +
261 +
262   end module simulation

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines