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 480 by chuckv, Tue Apr 8 17:16:22 2003 UTC vs.
Revision 1214 by gezelter, Tue Jun 1 18:42:58 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, CglobalGroupMembership, &
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(CnGlobal), intent(in):: CglobalGroupMembership
88 >    integer :: maxSkipsForAtom, glPointer
89 >
90   #ifdef IS_MPI
91      integer, allocatable, dimension(:) :: c_idents_Row
92      integer, allocatable, dimension(:) :: c_idents_Col
93 <    integer :: nrow
94 <    integer :: ncol
93 >    integer :: nAtomsInRow, nGroupsInRow, aid
94 >    integer :: nAtomsInCol, nGroupsInCol, gid
95   #endif  
96  
97      simulation_setup_complete = .false.
98      status = 0
99  
100      ! copy C struct into fortran type
101 +
102 +    nLocal = CnLocal
103 +    nGlobal = CnGlobal
104 +    nGroups = CnGroups
105 +
106      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
107  
108      nExcludes_Global = CnGlobalExcludes
109      nExcludes_Local = CnLocalExcludes
110  
111 <    call InitializeForceGlobals(natoms, thisStat)
111 >    call InitializeForceGlobals(nLocal, thisStat)
112      if (thisStat /= 0) then
113         write(default_error,*) "SimSetup: InitializeForceGlobals error"
114         status = -1
# Line 123 | Line 129 | contains
129         status = -1
130         return
131      endif
132 <    nrow = getNrow(plan_row)
133 <    ncol = getNcol(plan_col)
132 >    nAtomsInRow = getNatomsInRow(plan_atom_row)
133 >    nAtomsInCol = getNatomsInCol(plan_atom_col)
134 >    nGroupsInRow = getNgroupsInRow(plan_group_row)
135 >    nGroupsInCol = getNgroupsInCol(plan_group_col)
136      mynode = getMyNode()
137      
138 <    allocate(c_idents_Row(nrow),stat=alloc_stat)
138 >    allocate(c_idents_Row(nAtomsInRow),stat=alloc_stat)
139      if (alloc_stat /= 0 ) then
140         status = -1
141         return
142      endif
143  
144 <    allocate(c_idents_Col(ncol),stat=alloc_stat)
144 >    allocate(c_idents_Col(nAtomsInCol),stat=alloc_stat)
145      if (alloc_stat /= 0 ) then
146         status = -1
147         return
148      endif
149  
150 <    call gather(c_idents, c_idents_Row, plan_row)
151 <    call gather(c_idents, c_idents_Col, plan_col)
150 >    call gather(c_idents, c_idents_Row, plan_atom_row)
151 >    call gather(c_idents, c_idents_Col, plan_atom_col)
152  
153 <    do i = 1, nrow
153 >    do i = 1, nAtomsInRow
154         me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
155         atid_Row(i) = me
156      enddo
157  
158 <    do i = 1, ncol
158 >    do i = 1, nAtomsInCol
159         me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
160         atid_Col(i) = me
161      enddo
# Line 159 | Line 167 | contains
167      if (allocated(c_idents_Row)) then
168         deallocate(c_idents_Row)
169      endif
170 +  
171 + #endif
172 +
173 + #ifdef IS_MPI
174 +    allocate(groupStartRow(nGroupsInRow+1),stat=alloc_stat)
175 +    if (alloc_stat /= 0 ) then
176 +       status = -1
177 +       return
178 +    endif
179 +    allocate(groupStartCol(nGroupsInCol+1),stat=alloc_stat)
180 +    if (alloc_stat /= 0 ) then
181 +       status = -1
182 +       return
183 +    endif
184 +    allocate(groupListRow(nAtomsInRow),stat=alloc_stat)
185 +    if (alloc_stat /= 0 ) then
186 +       status = -1
187 +       return
188 +    endif
189 +    allocate(groupListCol(nAtomsInCol),stat=alloc_stat)
190 +    if (alloc_stat /= 0 ) then
191 +       status = -1
192 +       return
193 +    endif
194 +    allocate(mfactRow(nAtomsInRow),stat=alloc_stat)
195 +    if (alloc_stat /= 0 ) then
196 +       status = -1
197 +       return
198 +    endif
199 +    allocate(mfactCol(nAtomsInCol),stat=alloc_stat)
200 +    if (alloc_stat /= 0 ) then
201 +       status = -1
202 +       return
203 +    endif
204 +    allocate(mfactLocal(nLocal),stat=alloc_stat)
205 +    if (alloc_stat /= 0 ) then
206 +       status = -1
207 +       return
208 +    endif
209      
210 +    glPointer = 1
211 +
212 +    do i = 1, nGroupsInRow
213 +
214 +       gid = GroupRowToGlobal(i)
215 +       groupStartRow(i) = glPointer      
216 +
217 +       do j = 1, nAtomsInRow
218 +          aid = AtomRowToGlobal(j)
219 +          if (CglobalGroupMembership(aid) .eq. gid) then
220 +             groupListRow(glPointer) = j
221 +             glPointer = glPointer + 1
222 +          endif
223 +       enddo
224 +    enddo
225 +    groupStartRow(nGroupsInRow+1) = nAtomsInRow + 1
226 +
227 +    glPointer = 1
228 +
229 +    do i = 1, nGroupsInCol
230 +
231 +       gid = GroupColToGlobal(i)
232 +       groupStartCol(i) = glPointer      
233 +
234 +       do j = 1, nAtomsInCol
235 +          aid = AtomColToGlobal(j)
236 +          if (CglobalGroupMembership(aid) .eq. gid) then
237 +             groupListCol(glPointer) = j
238 +             glPointer = glPointer + 1
239 +          endif
240 +       enddo
241 +    enddo
242 +    groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1
243 +
244 +    mfactLocal = Cmfact        
245 +
246 +    call gather(mfactLocal,      mfactRow,      plan_atom_row)
247 +    call gather(mfactLocal,      mfactCol,      plan_atom_col)
248 +    
249 +    if (allocated(mfactLocal)) then
250 +       deallocate(mfactLocal)
251 +    end if
252   #else
253 <    do i = 1, nComponents
253 >    allocate(groupStartRow(nGroups+1),stat=alloc_stat)
254 >    if (alloc_stat /= 0 ) then
255 >       status = -1
256 >       return
257 >    endif
258 >    allocate(groupStartCol(nGroups+1),stat=alloc_stat)
259 >    if (alloc_stat /= 0 ) then
260 >       status = -1
261 >       return
262 >    endif
263 >    allocate(groupListRow(nLocal),stat=alloc_stat)
264 >    if (alloc_stat /= 0 ) then
265 >       status = -1
266 >       return
267 >    endif
268 >    allocate(groupListCol(nLocal),stat=alloc_stat)
269 >    if (alloc_stat /= 0 ) then
270 >       status = -1
271 >       return
272 >    endif
273 >    allocate(mfactRow(nLocal),stat=alloc_stat)
274 >    if (alloc_stat /= 0 ) then
275 >       status = -1
276 >       return
277 >    endif
278 >    allocate(mfactCol(nLocal),stat=alloc_stat)
279 >    if (alloc_stat /= 0 ) then
280 >       status = -1
281 >       return
282 >    endif
283 >    allocate(mfactLocal(nLocal),stat=alloc_stat)
284 >    if (alloc_stat /= 0 ) then
285 >       status = -1
286 >       return
287 >    endif
288 >
289 >    glPointer = 1
290 >    do i = 1, nGroups
291 >       groupStartRow(i) = glPointer      
292 >       groupStartCol(i) = glPointer
293 >       do j = 1, nLocal
294 >          if (CglobalGroupMembership(j) .eq. i) then
295 >             groupListRow(glPointer) = j
296 >             groupListCol(glPointer) = j
297 >             glPointer = glPointer + 1
298 >          endif
299 >       enddo
300 >    enddo
301 >    groupStartRow(nGroups+1) = nLocal + 1
302 >    groupStartCol(nGroups+1) = nLocal + 1
303 >
304 >    do i = 1, nLocal
305 >       mfactRow(i) = Cmfact(i)
306 >       mfactCol(i) = Cmfact(i)
307 >    end do
308 >    
309 > #endif
310 >
311 >
312 > ! We build the local atid's for both mpi and nonmpi
313 >    do i = 1, nLocal
314        
315         me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
316         atid(i) = me
317    
318      enddo
170 – #endif
319  
172 –
173 –
320      do i = 1, nExcludes_Local
321         excludesLocal(1,i) = CexcludesLocal(1,i)
322         excludesLocal(2,i) = CexcludesLocal(2,i)
323      enddo
324 +
325 +    maxSkipsForAtom = 0
326 + #ifdef IS_MPI
327 +    do j = 1, nAtomsInRow
328 + #else
329 +    do j = 1, nLocal
330 + #endif
331 +       nSkipsForAtom(j) = 0
332 + #ifdef IS_MPI
333 +       id1 = AtomRowToGlobal(j)
334 + #else
335 +       id1 = j
336 + #endif
337 +       do i = 1, nExcludes_Local
338 +          if (excludesLocal(1,i) .eq. id1 ) then
339 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
340 +
341 +             if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
342 +                maxSkipsForAtom = nSkipsForAtom(j)
343 +             endif
344 +          endif
345 +          if (excludesLocal(2,i) .eq. id1 ) then
346 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
347 +
348 +             if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
349 +                maxSkipsForAtom = nSkipsForAtom(j)
350 +             endif
351 +          endif
352 +       end do
353 +    enddo
354 +
355 + #ifdef IS_MPI
356 +    allocate(skipsForAtom(nAtomsInRow, maxSkipsForAtom), stat=alloc_stat)
357 + #else
358 +    allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
359 + #endif
360 +    if (alloc_stat /= 0 ) then
361 +       write(*,*) 'Could not allocate skipsForAtom array'
362 +       return
363 +    endif
364 +
365 + #ifdef IS_MPI
366 +    do j = 1, nAtomsInRow
367 + #else
368 +    do j = 1, nLocal
369 + #endif
370 +       nSkipsForAtom(j) = 0
371 + #ifdef IS_MPI
372 +       id1 = AtomRowToGlobal(j)
373 + #else
374 +       id1 = j
375 + #endif
376 +       do i = 1, nExcludes_Local
377 +          if (excludesLocal(1,i) .eq. id1 ) then
378 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
379 +             ! exclude lists have global ID's so this line is
380 +             ! the same in MPI and non-MPI
381 +             id2 = excludesLocal(2,i)
382 +             skipsForAtom(j, nSkipsForAtom(j)) = id2
383 +          endif
384 +          if (excludesLocal(2, i) .eq. id2 ) then
385 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
386 +             ! exclude lists have global ID's so this line is
387 +             ! the same in MPI and non-MPI
388 +             id2 = excludesLocal(1,i)
389 +             skipsForAtom(j, nSkipsForAtom(j)) = id2
390 +          endif
391 +       end do
392 +    enddo
393      
394      do i = 1, nExcludes_Global
395         excludesGlobal(i) = CexcludesGlobal(i)
396      enddo
397  
398 +    do i = 1, nGlobal
399 +       molMemberShipList(i) = CmolMembership(i)
400 +    enddo
401 +    
402      if (status == 0) simulation_setup_complete = .true.
403      
404    end subroutine SimulationSetup
405    
406 <  subroutine setBox_3d(new_box_size)
407 <    real(kind=dp), dimension(3) :: new_box_size
406 >  subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
407 >    real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
408 >    integer :: cBoxIsOrthorhombic
409      integer :: smallest, status, i
190 –
191 –    thisSim%box = new_box_size
192 –    box = thisSim%box
193 –
194 –    return    
195 –  end subroutine setBox_3d
196 –
197 –  subroutine setBox_1d(dim, new_box_size)
198 –    integer :: dim, status
199 –    real(kind=dp) :: new_box_size
200 –    thisSim%box(dim) = new_box_size
201 –    box(dim) = thisSim%box(dim)
202 –  end subroutine setBox_1d
203 –
204 –  subroutine setRcut(new_rcut, status)
205 –    real(kind = dp) :: new_rcut
206 –    integer :: myStatus, status
207 –    thisSim%rcut = new_rcut
208 –    rcut2 = thisSim%rcut * thisSim%rcut
209 –    rcut6 = rcut2 * rcut2 * rcut2
210 –    status = 0
211 –    return
212 –  end subroutine setRcut
213 –
214 –  function getBox_3d() result(thisBox)
215 –    real( kind = dp ), dimension(3) :: thisBox
216 –    thisBox = thisSim%box
217 –  end function getBox_3d
218 –  
219 –  function getBox_1d(dim) result(thisBox)
220 –    integer, intent(in) :: dim
221 –    real( kind = dp ) :: thisBox
410      
411 <    thisBox = thisSim%box(dim)
412 <  end function getBox_1d
413 <    
414 <  subroutine getRcut(thisrcut,rc2,rc6,status)
415 <    real( kind = dp ), intent(out) :: thisrcut
416 <    real( kind = dp ), intent(out), optional :: rc2
417 <    real( kind = dp ), intent(out), optional :: rc6
230 <    integer, optional :: status
231 <
232 <    if (present(status)) status = 0
411 >    Hmat = cHmat
412 >    HmatInv = cHmatInv
413 >    if (cBoxIsOrthorhombic .eq. 0 ) then
414 >       boxIsOrthorhombic = .false.
415 >    else
416 >       boxIsOrthorhombic = .true.
417 >    endif
418      
419 <    if (.not.simulation_setup_complete ) then
420 <       if (present(status)) status = -1
236 <       return
237 <    end if
238 <    
239 <    thisrcut = thisSim%rcut
240 <    if(present(rc2)) rc2 = rcut2
241 <    if(present(rc6)) rc6 = rcut6
242 <  end subroutine getRcut
243 <  
244 <  subroutine getRlist(thisrlist,rl2,status)
245 <    real( kind = dp ), intent(out) :: thisrlist
246 <    real( kind = dp ), intent(out), optional :: rl2
419 >    return    
420 >  end subroutine setBox
421  
248 –    integer, optional :: status
249 –
250 –    if (present(status)) status = 0
251 –
252 –    if (.not.simulation_setup_complete ) then
253 –       if (present(status)) status = -1
254 –       return
255 –    end if
256 –    
257 –    thisrlist = thisSim%rlist
258 –    if(present(rl2)) rl2 = rlist2
259 –  end subroutine getRlist
260 –
261 –  function getRrf() result(rrf)
262 –    real( kind = dp ) :: rrf
263 –    rrf = thisSim%rrf
264 –    write(*,*) 'getRrf = ', rrf, thisSim%rrf
265 –  end function getRrf
266 –  
267 –  function getRt() result(rt)
268 –    real( kind = dp ) :: rt
269 –    rt = thisSim%rt
270 –  end function getRt
271 –
422    function getDielect() result(dielect)
423      real( kind = dp ) :: dielect
424      dielect = thisSim%dielect
# Line 289 | Line 439 | contains
439      doesit = thisSim%SIM_uses_sticky
440    end function SimUsesSticky
441  
442 +  function SimUsesCharges() result(doesit)
443 +    logical :: doesit
444 +    doesit = thisSim%SIM_uses_charges
445 +  end function SimUsesCharges
446 +
447    function SimUsesDipoles() result(doesit)
448      logical :: doesit
449      doesit = thisSim%SIM_uses_dipoles
# Line 332 | Line 487 | contains
487      thisStat = 0
488      
489      call FreeSimGlobals()    
490 +
491 +    allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
492 +    if (alloc_stat /= 0 ) then
493 +       thisStat = -1
494 +       return
495 +    endif
496      
497      allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
498      if (alloc_stat /= 0 ) then
# Line 344 | Line 505 | contains
505         thisStat = -1
506         return
507      endif
508 +
509 +    allocate(molMembershipList(nGlobal), stat=alloc_stat)
510 +    if (alloc_stat /= 0 ) then
511 +       thisStat = -1
512 +       return
513 +    endif
514      
515    end subroutine InitializeSimGlobals
516    
517    subroutine FreeSimGlobals()
518      
519      !We free in the opposite order in which we allocate in.
520 <    
520 >
521 >    if (allocated(skipsForAtom)) deallocate(skipsForAtom)
522 >    if (allocated(mfactLocal)) deallocate(mfactLocal)
523 >    if (allocated(mfactCol)) deallocate(mfactCol)
524 >    if (allocated(mfactRow)) deallocate(mfactRow)
525 >    if (allocated(groupListCol)) deallocate(groupListCol)    
526 >    if (allocated(groupListRow)) deallocate(groupListRow)    
527 >    if (allocated(groupStartCol)) deallocate(groupStartCol)
528 >    if (allocated(groupStartRow)) deallocate(groupStartRow)    
529 >    if (allocated(molMembershipList)) deallocate(molMembershipList)    
530      if (allocated(excludesGlobal)) deallocate(excludesGlobal)
531      if (allocated(excludesLocal)) deallocate(excludesLocal)
532 <
532 >    
533    end subroutine FreeSimGlobals
534 <
535 <  pure function getNlocal() result(nlocal)
536 <    integer :: nlocal
537 <    nlocal = natoms
534 >  
535 >  pure function getNlocal() result(n)
536 >    integer :: n
537 >    n = nLocal
538    end function getNlocal
363 –
539    
540 +  
541   end module simulation

Diff Legend

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