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

Comparing trunk/OOPSE/libmdtools/simulation_module.F90 (file contents):
Revision 460 by chuckv, Fri Apr 4 22:22:30 2003 UTC vs.
Revision 1206 by tim, Thu May 27 19:51:18 2004 UTC

# Line 6 | Line 6 | module simulation
6    use force_globals
7    use vector_class
8    use atype_module
9 +  use switcheroo
10   #ifdef IS_MPI
11    use mpiSimulation
12   #endif
# Line 15 | Line 16 | module simulation
16  
17   #define __FORTRAN90
18   #include "fSimulation.h"
19 + #include "fSwitchingFunction.h"
20  
21 <  type (simtype), public :: thisSim
21 >  type (simtype), public, save :: thisSim
22  
23    logical, save :: simulation_setup_complete = .false.
24  
25 <  integer, public, save :: natoms
25 >  integer, public, save :: nLocal, nGlobal
26 >  integer, public, save :: nGroups, nGroupGlobal
27    integer, public, save :: nExcludes_Global = 0
28    integer, public, save :: nExcludes_Local = 0
29    integer, allocatable, dimension(:,:), public :: excludesLocal
30 <  integer, allocatable, dimension(:), public :: excludesGlobal
30 >  integer, allocatable, dimension(:),   public :: excludesGlobal
31 >  integer, allocatable, dimension(:),   public :: molMembershipList
32 >  integer, allocatable, dimension(:),   public :: groupListRow
33 >  integer, allocatable, dimension(:),   public :: groupStartRow
34 >  integer, allocatable, dimension(:),   public :: groupListCol
35 >  integer, allocatable, dimension(:),   public :: groupStartCol
36 >  integer, allocatable, dimension(:),   public :: groupListLocal
37 >  integer, allocatable, dimension(:),   public :: groupStartLocal
38 >  integer, allocatable, dimension(:),   public :: nSkipsForAtom
39 >  integer, allocatable, dimension(:,:), public :: skipsForAtom
40 >  real(kind=dp), allocatable, dimension(:), public :: mfactRow
41 >  real(kind=dp), allocatable, dimension(:), public :: mfactCol
42 >  real(kind=dp), allocatable, dimension(:), public :: mfactLocal
43  
44 <  real(kind=dp), save :: rcut2 = 0.0_DP
45 <  real(kind=dp), save :: rcut6 = 0.0_DP
31 <  real(kind=dp), save :: rlist2 = 0.0_DP
32 <  real(kind=dp), public, dimension(3), save :: box
33 <
44 >  real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
45 >  logical, public, save :: boxIsOrthorhombic
46    
47    public :: SimulationSetup
48    public :: getNlocal
49    public :: setBox
38 –  public :: setBox_3d
39 –  public :: getBox
40 –  public :: setRcut
41 –  public :: getRcut
42 –  public :: getRlist
43 –  public :: getRrf
44 –  public :: getRt
50    public :: getDielect
51    public :: SimUsesPBC
52    public :: SimUsesLJ
53 +  public :: SimUsesCharges
54    public :: SimUsesDipoles
55    public :: SimUsesSticky
56    public :: SimUsesRF
# Line 53 | Line 59 | module simulation
59    public :: SimRequiresPrepairCalc
60    public :: SimRequiresPostpairCalc
61    public :: SimUsesDirectionalAtoms
56 –
57 –  interface getBox
58 –     module procedure getBox_3d
59 –     module procedure getBox_1d
60 –  end interface
62    
62 –  interface setBox
63 –     module procedure setBox_3d
64 –     module procedure setBox_1d
65 –  end interface
66 –  
63   contains
64    
65 <  subroutine SimulationSetup(setThisSim, nComponents, c_idents, &
65 >  subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
66         CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
67 +       CmolMembership, Cmfact, CnGroups, CgroupList, CgroupStart, &
68         status)    
69  
70      type (simtype) :: setThisSim
71 <    integer, intent(inout) :: nComponents
72 <    integer, dimension(nComponents),intent(inout) :: c_idents
71 >    integer, intent(inout) :: CnGlobal, CnLocal
72 >    integer, dimension(CnLocal),intent(inout) :: c_idents
73  
74      integer :: CnLocalExcludes
75      integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
76      integer :: CnGlobalExcludes
77      integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
78 +    integer, dimension(CnGlobal),intent(in) :: CmolMembership
79      !!  Result status, success = 0, status = -1
80      integer, intent(out) :: status
81 <    integer :: i, me, thisStat, alloc_stat, myNode
81 >    integer :: i, j, me, thisStat, alloc_stat, myNode, id1, id2
82 >    integer :: ia
83 >
84 >    !! mass factors used for molecular cutoffs
85 >    real ( kind = dp ), dimension(CnLocal) :: Cmfact
86 >    integer, intent(in):: CnGroups
87 >    integer, dimension(CnLocal),intent(in) :: CgroupList
88 >    integer, dimension(CnGroups),intent(in) :: CgroupStart
89 >    integer :: maxSkipsForAtom
90 >
91   #ifdef IS_MPI
92      integer, allocatable, dimension(:) :: c_idents_Row
93      integer, allocatable, dimension(:) :: c_idents_Col
94 <    integer :: nrow
95 <    integer :: ncol
94 >    integer :: nAtomsInRow, nGroupsInRow
95 >    integer :: nAtomsInCol, nGroupsInCol
96   #endif  
97  
98      simulation_setup_complete = .false.
99      status = 0
100  
101      ! copy C struct into fortran type
102 +
103 +    nLocal = CnLocal
104 +    nGlobal = CnGlobal
105 +    nGroups = CnGroups
106 +
107      thisSim = setThisSim
96 –    natoms = nComponents
97 –    rcut2 = thisSim%rcut * thisSim%rcut
98 –    rcut6 = rcut2 * rcut2 * rcut2
99 –    rlist2 = thisSim%rlist * thisSim%rlist
100 –    box = thisSim%box
108  
109      nExcludes_Global = CnGlobalExcludes
110      nExcludes_Local = CnLocalExcludes
111  
112 <    call InitializeForceGlobals(natoms, thisStat)
112 >    call InitializeForceGlobals(nLocal, thisStat)
113      if (thisStat /= 0) then
114 +       write(default_error,*) "SimSetup: InitializeForceGlobals error"
115         status = -1
116         return
117      endif
118  
119      call InitializeSimGlobals(thisStat)
120      if (thisStat /= 0) then
121 +       write(default_error,*) "SimSetup: InitializeSimGlobals error"
122         status = -1
123         return
124      endif
# Line 121 | Line 130 | contains
130         status = -1
131         return
132      endif
133 <    nrow = getNrow(plan_row)
134 <    ncol = getNcol(plan_col)
133 >    nAtomsInRow = getNatomsInRow(plan_atom_row)
134 >    nAtomsInCol = getNatomsInCol(plan_atom_col)
135 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
136 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
137      mynode = getMyNode()
138      
139 <    allocate(c_idents_Row(nrow),stat=alloc_stat)
139 >    allocate(c_idents_Row(nAtomsInRow),stat=alloc_stat)
140      if (alloc_stat /= 0 ) then
141         status = -1
142         return
143      endif
144  
145 <    allocate(c_idents_Col(ncol),stat=alloc_stat)
145 >    allocate(c_idents_Col(nAtomsInCol),stat=alloc_stat)
146      if (alloc_stat /= 0 ) then
147         status = -1
148         return
149      endif
150  
151 <    call gather(c_idents, c_idents_Row, plan_row)
152 <    call gather(c_idents, c_idents_Col, plan_col)
151 >    call gather(c_idents, c_idents_Row, plan_atom_row)
152 >    call gather(c_idents, c_idents_Col, plan_atom_col)
153  
154 <    do i = 1, nrow
154 >    do i = 1, nAtomsInRow
155         me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
156         atid_Row(i) = me
157      enddo
158  
159 <    do i = 1, ncol
159 >    do i = 1, nAtomsInCol
160         me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
161         atid_Col(i) = me
162      enddo
# Line 157 | Line 168 | contains
168      if (allocated(c_idents_Row)) then
169         deallocate(c_idents_Row)
170      endif
171 +  
172 + #endif
173 +
174 + #ifdef IS_MPI
175 +    allocate(groupStartRow(nGroupsInRow+1),stat=alloc_stat)
176 +    if (alloc_stat /= 0 ) then
177 +       status = -1
178 +       return
179 +    endif
180 +    allocate(groupStartCol(nGroupsInCol+1),stat=alloc_stat)
181 +    if (alloc_stat /= 0 ) then
182 +       status = -1
183 +       return
184 +    endif
185 +    allocate(groupListRow(nAtomsInRow),stat=alloc_stat)
186 +    if (alloc_stat /= 0 ) then
187 +       status = -1
188 +       return
189 +    endif
190 +    allocate(groupListCol(nAtomsInCol),stat=alloc_stat)
191 +    if (alloc_stat /= 0 ) then
192 +       status = -1
193 +       return
194 +    endif
195 +    allocate(mfactRow(nAtomsInRow),stat=alloc_stat)
196 +    if (alloc_stat /= 0 ) then
197 +       status = -1
198 +       return
199 +    endif
200 +    allocate(mfactCol(nAtomsInCol),stat=alloc_stat)
201 +    if (alloc_stat /= 0 ) then
202 +       status = -1
203 +       return
204 +    endif
205 +    allocate(groupStartLocal(nGroups), stat=alloc_stat)
206 +    if (alloc_stat /= 0 ) then
207 +       status = -1
208 +       return
209 +    endif
210 +    allocate(groupListLocal(nLocal), stat=alloc_stat)
211 +    if (alloc_stat /= 0 ) then
212 +       status = -1
213 +       return
214 +    endif    
215 +    allocate(mfactLocal(nLocal), stat=alloc_stat)
216 +    if (alloc_stat /= 0 ) then
217 +       status = -1
218 +       return
219 +    endif    
220 +            
221 +    groupStartLocal = CgroupStart
222 +    groupListLocal = CgroupList
223 +    mfactLocal = Cmfact        
224 +            
225 +    call gather(groupStartLocal, groupStartRow, plan_group_row)
226 +    call gather(groupStartLocal, groupStartCol, plan_group_col)  
227 +    groupStartRow(nGroupsInRow+1) = nAtomsInRow + 1
228 +    groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1
229 +    call gather(groupListLocal,  groupListRow,  plan_atom_row)
230 +    call gather(groupListLocal,  groupListCol,  plan_atom_col)
231 +
232 +    ! C passes us groupList as globalID numbers for the atoms
233 +    ! we need to remap these onto the row and column ids on this
234 +    ! processor.  This is a linear search, but is only done once
235 +    ! (we hope)
236 +
237 +    do i = 1, nAtomsInRow
238 +       do j = 1, nAtomsInRow
239 +          if (AtomRowToGlobal(j) .eq. groupListRow(i)) then
240 +             groupListRow(i) = j
241 +          endif
242 +       enddo
243 +    enddo
244 +    do i = 1, nAtomsInCol
245 +       do j = 1, nAtomsInCol
246 +          if (AtomColToGlobal(j) .eq. groupListCol(i)) then
247 +             groupListCol(i) = j
248 +          endif
249 +       enddo
250 +    enddo
251 +    call gather(mfactLocal,      mfactRow,      plan_atom_row)
252 +    call gather(mfactLocal,      mfactCol,      plan_atom_col)
253      
254 +    if (allocated(mfactLocal)) then
255 +       deallocate(mfactLocal)
256 +    end if
257 +    if (allocated(groupListLocal)) then
258 +       deallocate(groupListLocal)
259 +    endif
260 +    if (allocated(groupStartLocal)) then
261 +       deallocate(groupStartLocal)
262 +    endif
263   #else
264 <    do i = 1, nComponents
264 >    allocate(groupStartRow(nGroups+1),stat=alloc_stat)
265 >    if (alloc_stat /= 0 ) then
266 >       status = -1
267 >       return
268 >    endif
269 >    allocate(groupStartCol(nGroups+1),stat=alloc_stat)
270 >    if (alloc_stat /= 0 ) then
271 >       status = -1
272 >       return
273 >    endif
274 >    allocate(groupListRow(nLocal),stat=alloc_stat)
275 >    if (alloc_stat /= 0 ) then
276 >       status = -1
277 >       return
278 >    endif
279 >    allocate(groupListCol(nLocal),stat=alloc_stat)
280 >    if (alloc_stat /= 0 ) then
281 >       status = -1
282 >       return
283 >    endif
284 >    allocate(mfactRow(nLocal),stat=alloc_stat)
285 >    if (alloc_stat /= 0 ) then
286 >       status = -1
287 >       return
288 >    endif
289 >    allocate(mfactCol(nLocal),stat=alloc_stat)
290 >    if (alloc_stat /= 0 ) then
291 >       status = -1
292 >       return
293 >    endif
294 >    do i = 1, nGroups
295 >       groupStartRow(i) = CgroupStart(i)
296 >       groupStartCol(i) = CgroupStart(i)
297 >    end do
298 >    groupStartRow(nGroups+1) = nLocal + 1
299 >    groupStartCol(nGroups+1) = nLocal + 1
300 >    do i = 1, nLocal
301 >       groupListRow(i) = CgroupList(i)
302 >       groupListCol(i) = CgroupList(i)
303 >       mfactRow(i) = Cmfact(i)
304 >       mfactCol(i) = Cmfact(i)
305 >    end do
306 >    
307 > #endif
308 >
309 >
310 > ! We build the local atid's for both mpi and nonmpi
311 >    do i = 1, nLocal
312        
313         me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
314         atid(i) = me
315    
316      enddo
317 +
318 +    do i = 1, nExcludes_Local
319 +       excludesLocal(1,i) = CexcludesLocal(1,i)
320 +       excludesLocal(2,i) = CexcludesLocal(2,i)
321 +    enddo
322 +
323 +    maxSkipsForAtom = 0
324 + #ifdef IS_MPI
325 +    do j = 1, nAtomsInRow
326 + #else
327 +    do j = 1, nLocal
328   #endif
329 +       nSkipsForAtom(j) = 0
330 + #ifdef IS_MPI
331 +       id1 = AtomRowToGlobal(j)
332 + #else
333 +       id1 = j
334 + #endif
335 +       do i = 1, nExcludes_Local
336 +          if (excludesLocal(1,i) .eq. id1 ) then
337 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
338  
339 <    !! Create neighbor lists
340 <    call expandNeighborList(nComponents, thisStat)
341 <    if (thisStat /= 0) then
342 <       status = -1
339 >             if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
340 >                maxSkipsForAtom = nSkipsForAtom(j)
341 >             endif
342 >          endif
343 >          if (excludesLocal(2,i) .eq. id1 ) then
344 >             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
345 >
346 >             if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
347 >                maxSkipsForAtom = nSkipsForAtom(j)
348 >             endif
349 >          endif
350 >       end do
351 >    enddo
352 >
353 > #ifdef IS_MPI
354 >    allocate(skipsForAtom(nAtomsInRow, maxSkipsForAtom), stat=alloc_stat)
355 > #else
356 >    allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
357 > #endif
358 >    if (alloc_stat /= 0 ) then
359 >       write(*,*) 'Could not allocate skipsForAtom array'
360         return
361      endif
362  
363 <
364 <    do i = 1, nExcludes_Local
365 <       excludesLocal(1,i) = CexcludesLocal(1,i)
366 <       excludesLocal(2,i) = CexcludesLocal(2,i)
363 > #ifdef IS_MPI
364 >    do j = 1, nAtomsInRow
365 > #else
366 >    do j = 1, nLocal
367 > #endif
368 >       nSkipsForAtom(j) = 0
369 > #ifdef IS_MPI
370 >       id1 = AtomRowToGlobal(j)
371 > #else
372 >       id1 = j
373 > #endif
374 >       do i = 1, nExcludes_Local
375 >          if (excludesLocal(1,i) .eq. id1 ) then
376 >             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
377 >             ! exclude lists have global ID's so this line is
378 >             ! the same in MPI and non-MPI
379 >             id2 = excludesLocal(2,i)
380 >             skipsForAtom(j, nSkipsForAtom(j)) = id2
381 >          endif
382 >          if (excludesLocal(2, i) .eq. id2 ) then
383 >             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
384 >             ! exclude lists have global ID's so this line is
385 >             ! the same in MPI and non-MPI
386 >             id2 = excludesLocal(1,i)
387 >             skipsForAtom(j, nSkipsForAtom(j)) = id2
388 >          endif
389 >       end do
390      enddo
391      
392      do i = 1, nExcludes_Global
393         excludesGlobal(i) = CexcludesGlobal(i)
394      enddo
395  
396 +    do i = 1, nGlobal
397 +       molMemberShipList(i) = CmolMembership(i)
398 +    enddo
399 +    
400      if (status == 0) simulation_setup_complete = .true.
401      
402    end subroutine SimulationSetup
403    
404 <  subroutine setBox_3d(new_box_size)
405 <    real(kind=dp), dimension(3) :: new_box_size
404 >  subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
405 >    real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
406 >    integer :: cBoxIsOrthorhombic
407      integer :: smallest, status, i
194 –
195 –    thisSim%box = new_box_size
196 –    box = thisSim%box
197 –
198 –    return    
199 –  end subroutine setBox_3d
200 –
201 –  subroutine setBox_1d(dim, new_box_size)
202 –    integer :: dim, status
203 –    real(kind=dp) :: new_box_size
204 –    thisSim%box(dim) = new_box_size
205 –    box(dim) = thisSim%box(dim)
206 –  end subroutine setBox_1d
207 –
208 –  subroutine setRcut(new_rcut, status)
209 –    real(kind = dp) :: new_rcut
210 –    integer :: myStatus, status
211 –    thisSim%rcut = new_rcut
212 –    rcut2 = thisSim%rcut * thisSim%rcut
213 –    rcut6 = rcut2 * rcut2 * rcut2
214 –    status = 0
215 –    return
216 –  end subroutine setRcut
217 –
218 –  function getBox_3d() result(thisBox)
219 –    real( kind = dp ), dimension(3) :: thisBox
220 –    thisBox = thisSim%box
221 –  end function getBox_3d
222 –  
223 –  function getBox_1d(dim) result(thisBox)
224 –    integer, intent(in) :: dim
225 –    real( kind = dp ) :: thisBox
408      
409 <    thisBox = thisSim%box(dim)
410 <  end function getBox_1d
411 <    
412 <  subroutine getRcut(thisrcut,rc2,rc6,status)
413 <    real( kind = dp ), intent(out) :: thisrcut
414 <    real( kind = dp ), intent(out), optional :: rc2
415 <    real( kind = dp ), intent(out), optional :: rc6
234 <    integer, optional :: status
235 <
236 <    if (present(status)) status = 0
409 >    Hmat = cHmat
410 >    HmatInv = cHmatInv
411 >    if (cBoxIsOrthorhombic .eq. 0 ) then
412 >       boxIsOrthorhombic = .false.
413 >    else
414 >       boxIsOrthorhombic = .true.
415 >    endif
416      
417 <    if (.not.simulation_setup_complete ) then
418 <       if (present(status)) status = -1
240 <       return
241 <    end if
242 <    
243 <    thisrcut = thisSim%rcut
244 <    if(present(rc2)) rc2 = rcut2
245 <    if(present(rc6)) rc6 = rcut6
246 <  end subroutine getRcut
247 <  
248 <  subroutine getRlist(thisrlist,rl2,status)
249 <    real( kind = dp ), intent(out) :: thisrlist
250 <    real( kind = dp ), intent(out), optional :: rl2
417 >    return    
418 >  end subroutine setBox
419  
252 –    integer, optional :: status
253 –
254 –    if (present(status)) status = 0
255 –
256 –    if (.not.simulation_setup_complete ) then
257 –       if (present(status)) status = -1
258 –       return
259 –    end if
260 –    
261 –    thisrlist = thisSim%rlist
262 –    if(present(rl2)) rl2 = rlist2
263 –  end subroutine getRlist
264 –
265 –  function getRrf() result(rrf)
266 –    real( kind = dp ) :: rrf
267 –    rrf = thisSim%rrf
268 –    write(*,*) 'getRrf = ', rrf, thisSim%rrf
269 –  end function getRrf
270 –  
271 –  function getRt() result(rt)
272 –    real( kind = dp ) :: rt
273 –    rt = thisSim%rt
274 –  end function getRt
275 –
420    function getDielect() result(dielect)
421      real( kind = dp ) :: dielect
422      dielect = thisSim%dielect
# Line 293 | Line 437 | contains
437      doesit = thisSim%SIM_uses_sticky
438    end function SimUsesSticky
439  
440 +  function SimUsesCharges() result(doesit)
441 +    logical :: doesit
442 +    doesit = thisSim%SIM_uses_charges
443 +  end function SimUsesCharges
444 +
445    function SimUsesDipoles() result(doesit)
446      logical :: doesit
447      doesit = thisSim%SIM_uses_dipoles
# Line 336 | Line 485 | contains
485      thisStat = 0
486      
487      call FreeSimGlobals()    
488 +
489 +    allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
490 +    if (alloc_stat /= 0 ) then
491 +       thisStat = -1
492 +       return
493 +    endif
494      
495      allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
496      if (alloc_stat /= 0 ) then
# Line 348 | Line 503 | contains
503         thisStat = -1
504         return
505      endif
506 +
507 +    allocate(molMembershipList(nGlobal), stat=alloc_stat)
508 +    if (alloc_stat /= 0 ) then
509 +       thisStat = -1
510 +       return
511 +    endif
512      
513    end subroutine InitializeSimGlobals
514    
515    subroutine FreeSimGlobals()
516      
517      !We free in the opposite order in which we allocate in.
518 <    
518 >
519 >    if (allocated(skipsForAtom)) deallocate(skipsForAtom)
520 >    if (allocated(mfactCol)) deallocate(mfactCol)
521 >    if (allocated(mfactRow)) deallocate(mfactRow)
522 >    if (allocated(groupListCol)) deallocate(groupListCol)    
523 >    if (allocated(groupListRow)) deallocate(groupListRow)    
524 >    if (allocated(groupStartCol)) deallocate(groupStartCol)
525 >    if (allocated(groupStartRow)) deallocate(groupStartRow)    
526 >    if (allocated(molMembershipList)) deallocate(molMembershipList)    
527      if (allocated(excludesGlobal)) deallocate(excludesGlobal)
528      if (allocated(excludesLocal)) deallocate(excludesLocal)
529 <
529 >    
530    end subroutine FreeSimGlobals
531 <
532 <  pure function getNlocal() result(nlocal)
533 <    integer :: nlocal
534 <    nlocal = natoms
531 >  
532 >  pure function getNlocal() result(n)
533 >    integer :: n
534 >    n = nLocal
535    end function getNlocal
367 –
536    
537 +  
538   end module simulation

Diff Legend

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