ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/simulation_module.F90
Revision: 1217
Committed: Tue Jun 1 21:45:22 2004 UTC (20 years, 1 month ago) by gezelter
File size: 15116 byte(s)
Log Message:
Bug fix (fixes of skipList and neighbor list under MPI)

File Contents

# User Rev Content
1 mmeineke 377 !! Fortran interface to C entry plug.
2    
3     module simulation
4     use definitions
5     use neighborLists
6     use force_globals
7     use vector_class
8     use atype_module
9 gezelter 1150 use switcheroo
10 mmeineke 377 #ifdef IS_MPI
11     use mpiSimulation
12     #endif
13    
14     implicit none
15     PRIVATE
16    
17     #define __FORTRAN90
18     #include "fSimulation.h"
19 gezelter 1150 #include "fSwitchingFunction.h"
20 mmeineke 377
21 gezelter 747 type (simtype), public, save :: thisSim
22 mmeineke 377
23     logical, save :: simulation_setup_complete = .false.
24    
25 mmeineke 491 integer, public, save :: nLocal, nGlobal
26 tim 1198 integer, public, save :: nGroups, nGroupGlobal
27 mmeineke 377 integer, public, save :: nExcludes_Global = 0
28     integer, public, save :: nExcludes_Local = 0
29     integer, allocatable, dimension(:,:), public :: excludesLocal
30 gezelter 1183 integer, allocatable, dimension(:), public :: excludesGlobal
31     integer, allocatable, dimension(:), public :: molMembershipList
32 tim 1198 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 gezelter 1183 integer, allocatable, dimension(:), public :: nSkipsForAtom
39     integer, allocatable, dimension(:,:), public :: skipsForAtom
40 tim 1198 real(kind=dp), allocatable, dimension(:), public :: mfactRow
41     real(kind=dp), allocatable, dimension(:), public :: mfactCol
42     real(kind=dp), allocatable, dimension(:), public :: mfactLocal
43 mmeineke 377
44 mmeineke 569 real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
45 gezelter 570 logical, public, save :: boxIsOrthorhombic
46 mmeineke 377
47     public :: SimulationSetup
48     public :: getNlocal
49     public :: setBox
50     public :: getDielect
51     public :: SimUsesPBC
52     public :: SimUsesLJ
53 gezelter 941 public :: SimUsesCharges
54 mmeineke 377 public :: SimUsesDipoles
55     public :: SimUsesSticky
56     public :: SimUsesRF
57     public :: SimUsesGB
58     public :: SimUsesEAM
59     public :: SimRequiresPrepairCalc
60     public :: SimRequiresPostpairCalc
61     public :: SimUsesDirectionalAtoms
62    
63     contains
64    
65 mmeineke 491 subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
66 gezelter 483 CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
67 gezelter 1214 CmolMembership, Cmfact, CnGroups, CglobalGroupMembership, &
68 mmeineke 377 status)
69    
70     type (simtype) :: setThisSim
71 mmeineke 491 integer, intent(inout) :: CnGlobal, CnLocal
72     integer, dimension(CnLocal),intent(inout) :: c_idents
73 mmeineke 377
74     integer :: CnLocalExcludes
75     integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
76     integer :: CnGlobalExcludes
77     integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
78 mmeineke 491 integer, dimension(CnGlobal),intent(in) :: CmolMembership
79 mmeineke 377 !! Result status, success = 0, status = -1
80     integer, intent(out) :: status
81 gezelter 1183 integer :: i, j, me, thisStat, alloc_stat, myNode, id1, id2
82 tim 1198 integer :: ia
83 tim 1144
84     !! mass factors used for molecular cutoffs
85 gezelter 1150 real ( kind = dp ), dimension(CnLocal) :: Cmfact
86 tim 1198 integer, intent(in):: CnGroups
87 gezelter 1214 integer, dimension(CnGlobal), intent(in):: CglobalGroupMembership
88     integer :: maxSkipsForAtom, glPointer
89 tim 1144
90 mmeineke 377 #ifdef IS_MPI
91     integer, allocatable, dimension(:) :: c_idents_Row
92     integer, allocatable, dimension(:) :: c_idents_Col
93 gezelter 1214 integer :: nAtomsInRow, nGroupsInRow, aid
94     integer :: nAtomsInCol, nGroupsInCol, gid
95 mmeineke 377 #endif
96    
97     simulation_setup_complete = .false.
98     status = 0
99    
100     ! copy C struct into fortran type
101 mmeineke 491
102     nLocal = CnLocal
103     nGlobal = CnGlobal
104 tim 1198 nGroups = CnGroups
105 mmeineke 491
106 mmeineke 377 thisSim = setThisSim
107 mmeineke 620
108 mmeineke 377 nExcludes_Global = CnGlobalExcludes
109     nExcludes_Local = CnLocalExcludes
110    
111 mmeineke 491 call InitializeForceGlobals(nLocal, thisStat)
112 mmeineke 377 if (thisStat /= 0) then
113 chuckv 480 write(default_error,*) "SimSetup: InitializeForceGlobals error"
114 mmeineke 377 status = -1
115     return
116     endif
117    
118     call InitializeSimGlobals(thisStat)
119     if (thisStat /= 0) then
120 chuckv 480 write(default_error,*) "SimSetup: InitializeSimGlobals error"
121 mmeineke 377 status = -1
122     return
123     endif
124    
125     #ifdef IS_MPI
126     ! We can only set up forces if mpiSimulation has been setup.
127     if (.not. isMPISimSet()) then
128     write(default_error,*) "MPI is not set"
129     status = -1
130     return
131     endif
132 tim 1198 nAtomsInRow = getNatomsInRow(plan_atom_row)
133     nAtomsInCol = getNatomsInCol(plan_atom_col)
134     nGroupsInRow = getNgroupsInRow(plan_group_row)
135     nGroupsInCol = getNgroupsInCol(plan_group_col)
136 mmeineke 377 mynode = getMyNode()
137    
138 tim 1198 allocate(c_idents_Row(nAtomsInRow),stat=alloc_stat)
139 mmeineke 377 if (alloc_stat /= 0 ) then
140     status = -1
141     return
142     endif
143    
144 tim 1198 allocate(c_idents_Col(nAtomsInCol),stat=alloc_stat)
145 mmeineke 377 if (alloc_stat /= 0 ) then
146     status = -1
147     return
148     endif
149    
150 tim 1198 call gather(c_idents, c_idents_Row, plan_atom_row)
151     call gather(c_idents, c_idents_Col, plan_atom_col)
152 mmeineke 377
153 tim 1198 do i = 1, nAtomsInRow
154 mmeineke 377 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
155     atid_Row(i) = me
156     enddo
157    
158 tim 1198 do i = 1, nAtomsInCol
159 mmeineke 377 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
160     atid_Col(i) = me
161     enddo
162    
163     !! free temporary ident arrays
164     if (allocated(c_idents_Col)) then
165     deallocate(c_idents_Col)
166     end if
167     if (allocated(c_idents_Row)) then
168     deallocate(c_idents_Row)
169     endif
170 tim 1198
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 gezelter 1214 allocate(mfactLocal(nLocal),stat=alloc_stat)
205 tim 1198 if (alloc_stat /= 0 ) then
206     status = -1
207     return
208     endif
209 gezelter 1214
210     glPointer = 1
211 gezelter 1203
212 gezelter 1214 do i = 1, nGroupsInRow
213 gezelter 1203
214 gezelter 1214 gid = GroupRowToGlobal(i)
215     groupStartRow(i) = glPointer
216    
217 gezelter 1203 do j = 1, nAtomsInRow
218 gezelter 1214 aid = AtomRowToGlobal(j)
219     if (CglobalGroupMembership(aid) .eq. gid) then
220     groupListRow(glPointer) = j
221     glPointer = glPointer + 1
222 gezelter 1203 endif
223     enddo
224     enddo
225 gezelter 1214 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 gezelter 1203 do j = 1, nAtomsInCol
235 gezelter 1214 aid = AtomColToGlobal(j)
236     if (CglobalGroupMembership(aid) .eq. gid) then
237     groupListCol(glPointer) = j
238     glPointer = glPointer + 1
239 gezelter 1203 endif
240     enddo
241     enddo
242 gezelter 1214 groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1
243    
244     mfactLocal = Cmfact
245    
246 tim 1198 call gather(mfactLocal, mfactRow, plan_atom_row)
247     call gather(mfactLocal, mfactCol, plan_atom_col)
248 mmeineke 377
249 tim 1198 if (allocated(mfactLocal)) then
250     deallocate(mfactLocal)
251     end if
252     #else
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 gezelter 1214 allocate(mfactLocal(nLocal),stat=alloc_stat)
284     if (alloc_stat /= 0 ) then
285     status = -1
286     return
287     endif
288    
289     glPointer = 1
290 tim 1198 do i = 1, nGroups
291 gezelter 1214 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 gezelter 1199 groupStartRow(nGroups+1) = nLocal + 1
302     groupStartCol(nGroups+1) = nLocal + 1
303 gezelter 1214
304 tim 1198 do i = 1, nLocal
305     mfactRow(i) = Cmfact(i)
306     mfactCol(i) = Cmfact(i)
307     end do
308    
309 chuckv 648 #endif
310    
311 tim 1198
312 chuckv 648 ! We build the local atid's for both mpi and nonmpi
313 gezelter 490 do i = 1, nLocal
314 mmeineke 377
315     me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
316     atid(i) = me
317    
318     enddo
319    
320     do i = 1, nExcludes_Local
321     excludesLocal(1,i) = CexcludesLocal(1,i)
322     excludesLocal(2,i) = CexcludesLocal(2,i)
323     enddo
324 gezelter 1183
325 gezelter 1217 #ifdef IS_MPI
326     allocate(nSkipsForAtom(nAtomsInRow), stat=alloc_stat)
327     #else
328     allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
329     #endif
330     if (alloc_stat /= 0 ) then
331     thisStat = -1
332     write(*,*) 'Could not allocate nSkipsForAtom array'
333     return
334     endif
335    
336 gezelter 1183 maxSkipsForAtom = 0
337 tim 1206 #ifdef IS_MPI
338     do j = 1, nAtomsInRow
339     #else
340 gezelter 1183 do j = 1, nLocal
341 tim 1206 #endif
342 gezelter 1183 nSkipsForAtom(j) = 0
343     #ifdef IS_MPI
344 tim 1198 id1 = AtomRowToGlobal(j)
345 gezelter 1183 #else
346     id1 = j
347     #endif
348     do i = 1, nExcludes_Local
349     if (excludesLocal(1,i) .eq. id1 ) then
350     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
351    
352     if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
353     maxSkipsForAtom = nSkipsForAtom(j)
354     endif
355     endif
356     if (excludesLocal(2,i) .eq. id1 ) then
357     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
358    
359     if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
360     maxSkipsForAtom = nSkipsForAtom(j)
361     endif
362     endif
363     end do
364     enddo
365    
366 tim 1206 #ifdef IS_MPI
367     allocate(skipsForAtom(nAtomsInRow, maxSkipsForAtom), stat=alloc_stat)
368     #else
369 gezelter 1183 allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
370 tim 1206 #endif
371 gezelter 1183 if (alloc_stat /= 0 ) then
372     write(*,*) 'Could not allocate skipsForAtom array'
373     return
374     endif
375    
376 tim 1206 #ifdef IS_MPI
377     do j = 1, nAtomsInRow
378     #else
379 gezelter 1183 do j = 1, nLocal
380 tim 1206 #endif
381 gezelter 1183 nSkipsForAtom(j) = 0
382     #ifdef IS_MPI
383 tim 1198 id1 = AtomRowToGlobal(j)
384 gezelter 1183 #else
385     id1 = j
386     #endif
387     do i = 1, nExcludes_Local
388     if (excludesLocal(1,i) .eq. id1 ) then
389     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
390 tim 1206 ! exclude lists have global ID's so this line is
391     ! the same in MPI and non-MPI
392 gezelter 1183 id2 = excludesLocal(2,i)
393     skipsForAtom(j, nSkipsForAtom(j)) = id2
394     endif
395 gezelter 1217 if (excludesLocal(2, i) .eq. id1 ) then
396 gezelter 1183 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
397 tim 1206 ! exclude lists have global ID's so this line is
398     ! the same in MPI and non-MPI
399 gezelter 1183 id2 = excludesLocal(1,i)
400     skipsForAtom(j, nSkipsForAtom(j)) = id2
401     endif
402     end do
403     enddo
404 mmeineke 377
405     do i = 1, nExcludes_Global
406     excludesGlobal(i) = CexcludesGlobal(i)
407     enddo
408 mmeineke 435
409 gezelter 490 do i = 1, nGlobal
410 gezelter 483 molMemberShipList(i) = CmolMembership(i)
411 gezelter 1150 enddo
412    
413 mmeineke 377 if (status == 0) simulation_setup_complete = .true.
414    
415     end subroutine SimulationSetup
416    
417 mmeineke 569 subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
418 gezelter 570 real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
419 mmeineke 569 integer :: cBoxIsOrthorhombic
420 mmeineke 377 integer :: smallest, status, i
421 gezelter 570
422 mmeineke 569 Hmat = cHmat
423     HmatInv = cHmatInv
424 gezelter 570 if (cBoxIsOrthorhombic .eq. 0 ) then
425 mmeineke 569 boxIsOrthorhombic = .false.
426 gezelter 570 else
427     boxIsOrthorhombic = .true.
428 mmeineke 569 endif
429    
430 mmeineke 377 return
431 mmeineke 569 end subroutine setBox
432 mmeineke 377
433     function getDielect() result(dielect)
434     real( kind = dp ) :: dielect
435     dielect = thisSim%dielect
436     end function getDielect
437    
438     function SimUsesPBC() result(doesit)
439     logical :: doesit
440     doesit = thisSim%SIM_uses_PBC
441     end function SimUsesPBC
442    
443     function SimUsesLJ() result(doesit)
444     logical :: doesit
445     doesit = thisSim%SIM_uses_LJ
446     end function SimUsesLJ
447    
448     function SimUsesSticky() result(doesit)
449     logical :: doesit
450     doesit = thisSim%SIM_uses_sticky
451     end function SimUsesSticky
452    
453 gezelter 941 function SimUsesCharges() result(doesit)
454     logical :: doesit
455     doesit = thisSim%SIM_uses_charges
456     end function SimUsesCharges
457    
458 mmeineke 377 function SimUsesDipoles() result(doesit)
459     logical :: doesit
460     doesit = thisSim%SIM_uses_dipoles
461     end function SimUsesDipoles
462    
463     function SimUsesRF() result(doesit)
464     logical :: doesit
465     doesit = thisSim%SIM_uses_RF
466     end function SimUsesRF
467    
468     function SimUsesGB() result(doesit)
469     logical :: doesit
470     doesit = thisSim%SIM_uses_GB
471     end function SimUsesGB
472    
473     function SimUsesEAM() result(doesit)
474     logical :: doesit
475     doesit = thisSim%SIM_uses_EAM
476     end function SimUsesEAM
477    
478     function SimUsesDirectionalAtoms() result(doesit)
479     logical :: doesit
480     doesit = thisSim%SIM_uses_dipoles .or. thisSim%SIM_uses_sticky .or. &
481     thisSim%SIM_uses_GB .or. thisSim%SIM_uses_RF
482     end function SimUsesDirectionalAtoms
483    
484     function SimRequiresPrepairCalc() result(doesit)
485     logical :: doesit
486     doesit = thisSim%SIM_uses_EAM
487     end function SimRequiresPrepairCalc
488    
489     function SimRequiresPostpairCalc() result(doesit)
490     logical :: doesit
491     doesit = thisSim%SIM_uses_RF
492     end function SimRequiresPostpairCalc
493    
494     subroutine InitializeSimGlobals(thisStat)
495     integer, intent(out) :: thisStat
496     integer :: alloc_stat
497    
498     thisStat = 0
499    
500     call FreeSimGlobals()
501    
502     allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
503     if (alloc_stat /= 0 ) then
504     thisStat = -1
505     return
506     endif
507    
508     allocate(excludesGlobal(nExcludes_Global), stat=alloc_stat)
509     if (alloc_stat /= 0 ) then
510     thisStat = -1
511     return
512     endif
513 chuckv 482
514 mmeineke 491 allocate(molMembershipList(nGlobal), stat=alloc_stat)
515 chuckv 482 if (alloc_stat /= 0 ) then
516     thisStat = -1
517     return
518     endif
519 mmeineke 377
520     end subroutine InitializeSimGlobals
521    
522     subroutine FreeSimGlobals()
523    
524     !We free in the opposite order in which we allocate in.
525 gezelter 483
526 gezelter 1183 if (allocated(skipsForAtom)) deallocate(skipsForAtom)
527 gezelter 1217 if (allocated(nSkipsForAtom)) deallocate(nSkipsForAtom)
528 gezelter 1214 if (allocated(mfactLocal)) deallocate(mfactLocal)
529 tim 1198 if (allocated(mfactCol)) deallocate(mfactCol)
530     if (allocated(mfactRow)) deallocate(mfactRow)
531     if (allocated(groupListCol)) deallocate(groupListCol)
532     if (allocated(groupListRow)) deallocate(groupListRow)
533     if (allocated(groupStartCol)) deallocate(groupStartCol)
534     if (allocated(groupStartRow)) deallocate(groupStartRow)
535 gezelter 483 if (allocated(molMembershipList)) deallocate(molMembershipList)
536 mmeineke 377 if (allocated(excludesGlobal)) deallocate(excludesGlobal)
537     if (allocated(excludesLocal)) deallocate(excludesLocal)
538 gezelter 1150
539 mmeineke 377 end subroutine FreeSimGlobals
540 gezelter 1150
541 mmeineke 491 pure function getNlocal() result(n)
542     integer :: n
543     n = nLocal
544 mmeineke 377 end function getNlocal
545    
546 gezelter 1150
547 mmeineke 377 end module simulation