ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/simulation_module.F90
Revision: 1203
Committed: Thu May 27 18:59:17 2004 UTC (20 years, 1 month ago) by gezelter
File size: 14978 byte(s)
Log Message:
Cutoff group changes 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 tim 1198 CmolMembership, Cmfact, CnGroups, CgroupList, CgroupStart, &
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 1150 integer, dimension(CnLocal),intent(in) :: CgroupList
88 tim 1198 integer, dimension(CnGroups),intent(in) :: CgroupStart
89 gezelter 1183 integer :: maxSkipsForAtom
90 tim 1144
91 mmeineke 377 #ifdef IS_MPI
92     integer, allocatable, dimension(:) :: c_idents_Row
93     integer, allocatable, dimension(:) :: c_idents_Col
94 tim 1198 integer :: nAtomsInRow, nGroupsInRow
95     integer :: nAtomsInCol, nGroupsInCol
96 mmeineke 377 #endif
97    
98     simulation_setup_complete = .false.
99     status = 0
100    
101     ! copy C struct into fortran type
102 mmeineke 491
103     nLocal = CnLocal
104     nGlobal = CnGlobal
105 tim 1198 nGroups = CnGroups
106 mmeineke 491
107 mmeineke 377 thisSim = setThisSim
108 mmeineke 620
109 mmeineke 377 nExcludes_Global = CnGlobalExcludes
110     nExcludes_Local = CnLocalExcludes
111    
112 mmeineke 491 call InitializeForceGlobals(nLocal, thisStat)
113 mmeineke 377 if (thisStat /= 0) then
114 chuckv 480 write(default_error,*) "SimSetup: InitializeForceGlobals error"
115 mmeineke 377 status = -1
116     return
117     endif
118    
119     call InitializeSimGlobals(thisStat)
120     if (thisStat /= 0) then
121 chuckv 480 write(default_error,*) "SimSetup: InitializeSimGlobals error"
122 mmeineke 377 status = -1
123     return
124     endif
125    
126     #ifdef IS_MPI
127     ! We can only set up forces if mpiSimulation has been setup.
128     if (.not. isMPISimSet()) then
129     write(default_error,*) "MPI is not set"
130     status = -1
131     return
132     endif
133 tim 1198 nAtomsInRow = getNatomsInRow(plan_atom_row)
134     nAtomsInCol = getNatomsInCol(plan_atom_col)
135     nGroupsInRow = getNgroupsInRow(plan_group_row)
136     nGroupsInCol = getNgroupsInCol(plan_group_col)
137 mmeineke 377 mynode = getMyNode()
138    
139 tim 1198 allocate(c_idents_Row(nAtomsInRow),stat=alloc_stat)
140 mmeineke 377 if (alloc_stat /= 0 ) then
141     status = -1
142     return
143     endif
144    
145 tim 1198 allocate(c_idents_Col(nAtomsInCol),stat=alloc_stat)
146 mmeineke 377 if (alloc_stat /= 0 ) then
147     status = -1
148     return
149     endif
150    
151 tim 1198 call gather(c_idents, c_idents_Row, plan_atom_row)
152     call gather(c_idents, c_idents_Col, plan_atom_col)
153 mmeineke 377
154 tim 1198 do i = 1, nAtomsInRow
155 mmeineke 377 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
156     atid_Row(i) = me
157     enddo
158    
159 tim 1198 do i = 1, nAtomsInCol
160 mmeineke 377 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
161     atid_Col(i) = me
162     enddo
163    
164     !! free temporary ident arrays
165     if (allocated(c_idents_Col)) then
166     deallocate(c_idents_Col)
167     end if
168     if (allocated(c_idents_Row)) then
169     deallocate(c_idents_Row)
170     endif
171 tim 1198
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 gezelter 1203 call gather(groupStartLocal, groupStartCol, plan_group_col)
227 gezelter 1199 groupStartRow(nGroupsInRow+1) = nAtomsInRow + 1
228     groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1
229 tim 1198 call gather(groupListLocal, groupListRow, plan_atom_row)
230     call gather(groupListLocal, groupListCol, plan_atom_col)
231 gezelter 1203
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 tim 1198 call gather(mfactLocal, mfactRow, plan_atom_row)
252     call gather(mfactLocal, mfactCol, plan_atom_col)
253 mmeineke 377
254 tim 1198 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     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 gezelter 1199 groupStartRow(nGroups+1) = nLocal + 1
299     groupStartCol(nGroups+1) = nLocal + 1
300 tim 1198 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 chuckv 648 #endif
308    
309 tim 1198
310 chuckv 648 ! We build the local atid's for both mpi and nonmpi
311 gezelter 490 do i = 1, nLocal
312 mmeineke 377
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 gezelter 1183
323     maxSkipsForAtom = 0
324     do j = 1, nLocal
325     nSkipsForAtom(j) = 0
326     #ifdef IS_MPI
327 tim 1198 id1 = AtomRowToGlobal(j)
328 gezelter 1183 #else
329     id1 = j
330     #endif
331     do i = 1, nExcludes_Local
332     if (excludesLocal(1,i) .eq. id1 ) then
333     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
334    
335     if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
336     maxSkipsForAtom = nSkipsForAtom(j)
337     endif
338     endif
339     if (excludesLocal(2,i) .eq. id1 ) then
340     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
341    
342     if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
343     maxSkipsForAtom = nSkipsForAtom(j)
344     endif
345     endif
346     end do
347     enddo
348    
349     allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
350     if (alloc_stat /= 0 ) then
351     write(*,*) 'Could not allocate skipsForAtom array'
352     return
353     endif
354    
355     do j = 1, nLocal
356     nSkipsForAtom(j) = 0
357     #ifdef IS_MPI
358 tim 1198 id1 = AtomRowToGlobal(j)
359 gezelter 1183 #else
360     id1 = j
361     #endif
362     do i = 1, nExcludes_Local
363     if (excludesLocal(1,i) .eq. id1 ) then
364     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
365     #ifdef IS_MPI
366 tim 1198 id2 = AtomColToGlobal(excludesLocal(2,i))
367 gezelter 1183 #else
368     id2 = excludesLocal(2,i)
369     #endif
370     skipsForAtom(j, nSkipsForAtom(j)) = id2
371     endif
372     if (excludesLocal(2, i) .eq. id2 ) then
373     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
374     #ifdef IS_MPI
375 tim 1198 id2 = AtomColToGlobal(excludesLocal(1,i))
376 gezelter 1183 #else
377     id2 = excludesLocal(1,i)
378     #endif
379     skipsForAtom(j, nSkipsForAtom(j)) = id2
380     endif
381     end do
382     enddo
383 mmeineke 377
384     do i = 1, nExcludes_Global
385     excludesGlobal(i) = CexcludesGlobal(i)
386     enddo
387 mmeineke 435
388 gezelter 490 do i = 1, nGlobal
389 gezelter 483 molMemberShipList(i) = CmolMembership(i)
390 gezelter 1150 enddo
391    
392 mmeineke 377 if (status == 0) simulation_setup_complete = .true.
393    
394     end subroutine SimulationSetup
395    
396 mmeineke 569 subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
397 gezelter 570 real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
398 mmeineke 569 integer :: cBoxIsOrthorhombic
399 mmeineke 377 integer :: smallest, status, i
400 gezelter 570
401 mmeineke 569 Hmat = cHmat
402     HmatInv = cHmatInv
403 gezelter 570 if (cBoxIsOrthorhombic .eq. 0 ) then
404 mmeineke 569 boxIsOrthorhombic = .false.
405 gezelter 570 else
406     boxIsOrthorhombic = .true.
407 mmeineke 569 endif
408    
409 mmeineke 377 return
410 mmeineke 569 end subroutine setBox
411 mmeineke 377
412     function getDielect() result(dielect)
413     real( kind = dp ) :: dielect
414     dielect = thisSim%dielect
415     end function getDielect
416    
417     function SimUsesPBC() result(doesit)
418     logical :: doesit
419     doesit = thisSim%SIM_uses_PBC
420     end function SimUsesPBC
421    
422     function SimUsesLJ() result(doesit)
423     logical :: doesit
424     doesit = thisSim%SIM_uses_LJ
425     end function SimUsesLJ
426    
427     function SimUsesSticky() result(doesit)
428     logical :: doesit
429     doesit = thisSim%SIM_uses_sticky
430     end function SimUsesSticky
431    
432 gezelter 941 function SimUsesCharges() result(doesit)
433     logical :: doesit
434     doesit = thisSim%SIM_uses_charges
435     end function SimUsesCharges
436    
437 mmeineke 377 function SimUsesDipoles() result(doesit)
438     logical :: doesit
439     doesit = thisSim%SIM_uses_dipoles
440     end function SimUsesDipoles
441    
442     function SimUsesRF() result(doesit)
443     logical :: doesit
444     doesit = thisSim%SIM_uses_RF
445     end function SimUsesRF
446    
447     function SimUsesGB() result(doesit)
448     logical :: doesit
449     doesit = thisSim%SIM_uses_GB
450     end function SimUsesGB
451    
452     function SimUsesEAM() result(doesit)
453     logical :: doesit
454     doesit = thisSim%SIM_uses_EAM
455     end function SimUsesEAM
456    
457     function SimUsesDirectionalAtoms() result(doesit)
458     logical :: doesit
459     doesit = thisSim%SIM_uses_dipoles .or. thisSim%SIM_uses_sticky .or. &
460     thisSim%SIM_uses_GB .or. thisSim%SIM_uses_RF
461     end function SimUsesDirectionalAtoms
462    
463     function SimRequiresPrepairCalc() result(doesit)
464     logical :: doesit
465     doesit = thisSim%SIM_uses_EAM
466     end function SimRequiresPrepairCalc
467    
468     function SimRequiresPostpairCalc() result(doesit)
469     logical :: doesit
470     doesit = thisSim%SIM_uses_RF
471     end function SimRequiresPostpairCalc
472    
473     subroutine InitializeSimGlobals(thisStat)
474     integer, intent(out) :: thisStat
475     integer :: alloc_stat
476    
477     thisStat = 0
478    
479     call FreeSimGlobals()
480 gezelter 1183
481     allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
482     if (alloc_stat /= 0 ) then
483     thisStat = -1
484     return
485     endif
486 mmeineke 377
487     allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
488     if (alloc_stat /= 0 ) then
489     thisStat = -1
490     return
491     endif
492    
493     allocate(excludesGlobal(nExcludes_Global), stat=alloc_stat)
494     if (alloc_stat /= 0 ) then
495     thisStat = -1
496     return
497     endif
498 chuckv 482
499 mmeineke 491 allocate(molMembershipList(nGlobal), stat=alloc_stat)
500 chuckv 482 if (alloc_stat /= 0 ) then
501     thisStat = -1
502     return
503     endif
504 mmeineke 377
505     end subroutine InitializeSimGlobals
506    
507     subroutine FreeSimGlobals()
508    
509     !We free in the opposite order in which we allocate in.
510 gezelter 483
511 gezelter 1183 if (allocated(skipsForAtom)) deallocate(skipsForAtom)
512 tim 1198 if (allocated(mfactCol)) deallocate(mfactCol)
513     if (allocated(mfactRow)) deallocate(mfactRow)
514     if (allocated(groupListCol)) deallocate(groupListCol)
515     if (allocated(groupListRow)) deallocate(groupListRow)
516     if (allocated(groupStartCol)) deallocate(groupStartCol)
517     if (allocated(groupStartRow)) deallocate(groupStartRow)
518 gezelter 483 if (allocated(molMembershipList)) deallocate(molMembershipList)
519 mmeineke 377 if (allocated(excludesGlobal)) deallocate(excludesGlobal)
520     if (allocated(excludesLocal)) deallocate(excludesLocal)
521 gezelter 1150
522 mmeineke 377 end subroutine FreeSimGlobals
523 gezelter 1150
524 mmeineke 491 pure function getNlocal() result(n)
525     integer :: n
526     n = nLocal
527 mmeineke 377 end function getNlocal
528    
529 gezelter 1150
530 mmeineke 377 end module simulation