ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/simulation_module.F90
Revision: 1206
Committed: Thu May 27 19:51:18 2004 UTC (20 years, 3 months ago) by tim
File size: 15231 byte(s)
Log Message:
Bug fix for SkipList

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 tim 1206 #ifdef IS_MPI
325     do j = 1, nAtomsInRow
326     #else
327 gezelter 1183 do j = 1, nLocal
328 tim 1206 #endif
329 gezelter 1183 nSkipsForAtom(j) = 0
330     #ifdef IS_MPI
331 tim 1198 id1 = AtomRowToGlobal(j)
332 gezelter 1183 #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     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 tim 1206 #ifdef IS_MPI
354     allocate(skipsForAtom(nAtomsInRow, maxSkipsForAtom), stat=alloc_stat)
355     #else
356 gezelter 1183 allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
357 tim 1206 #endif
358 gezelter 1183 if (alloc_stat /= 0 ) then
359     write(*,*) 'Could not allocate skipsForAtom array'
360     return
361     endif
362    
363 tim 1206 #ifdef IS_MPI
364     do j = 1, nAtomsInRow
365     #else
366 gezelter 1183 do j = 1, nLocal
367 tim 1206 #endif
368 gezelter 1183 nSkipsForAtom(j) = 0
369     #ifdef IS_MPI
370 tim 1198 id1 = AtomRowToGlobal(j)
371 gezelter 1183 #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 tim 1206 ! exclude lists have global ID's so this line is
378     ! the same in MPI and non-MPI
379 gezelter 1183 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 tim 1206 ! exclude lists have global ID's so this line is
385     ! the same in MPI and non-MPI
386 gezelter 1183 id2 = excludesLocal(1,i)
387     skipsForAtom(j, nSkipsForAtom(j)) = id2
388     endif
389     end do
390     enddo
391 mmeineke 377
392     do i = 1, nExcludes_Global
393     excludesGlobal(i) = CexcludesGlobal(i)
394     enddo
395 mmeineke 435
396 gezelter 490 do i = 1, nGlobal
397 gezelter 483 molMemberShipList(i) = CmolMembership(i)
398 gezelter 1150 enddo
399    
400 mmeineke 377 if (status == 0) simulation_setup_complete = .true.
401    
402     end subroutine SimulationSetup
403    
404 mmeineke 569 subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
405 gezelter 570 real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
406 mmeineke 569 integer :: cBoxIsOrthorhombic
407 mmeineke 377 integer :: smallest, status, i
408 gezelter 570
409 mmeineke 569 Hmat = cHmat
410     HmatInv = cHmatInv
411 gezelter 570 if (cBoxIsOrthorhombic .eq. 0 ) then
412 mmeineke 569 boxIsOrthorhombic = .false.
413 gezelter 570 else
414     boxIsOrthorhombic = .true.
415 mmeineke 569 endif
416    
417 mmeineke 377 return
418 mmeineke 569 end subroutine setBox
419 mmeineke 377
420     function getDielect() result(dielect)
421     real( kind = dp ) :: dielect
422     dielect = thisSim%dielect
423     end function getDielect
424    
425     function SimUsesPBC() result(doesit)
426     logical :: doesit
427     doesit = thisSim%SIM_uses_PBC
428     end function SimUsesPBC
429    
430     function SimUsesLJ() result(doesit)
431     logical :: doesit
432     doesit = thisSim%SIM_uses_LJ
433     end function SimUsesLJ
434    
435     function SimUsesSticky() result(doesit)
436     logical :: doesit
437     doesit = thisSim%SIM_uses_sticky
438     end function SimUsesSticky
439    
440 gezelter 941 function SimUsesCharges() result(doesit)
441     logical :: doesit
442     doesit = thisSim%SIM_uses_charges
443     end function SimUsesCharges
444    
445 mmeineke 377 function SimUsesDipoles() result(doesit)
446     logical :: doesit
447     doesit = thisSim%SIM_uses_dipoles
448     end function SimUsesDipoles
449    
450     function SimUsesRF() result(doesit)
451     logical :: doesit
452     doesit = thisSim%SIM_uses_RF
453     end function SimUsesRF
454    
455     function SimUsesGB() result(doesit)
456     logical :: doesit
457     doesit = thisSim%SIM_uses_GB
458     end function SimUsesGB
459    
460     function SimUsesEAM() result(doesit)
461     logical :: doesit
462     doesit = thisSim%SIM_uses_EAM
463     end function SimUsesEAM
464    
465     function SimUsesDirectionalAtoms() result(doesit)
466     logical :: doesit
467     doesit = thisSim%SIM_uses_dipoles .or. thisSim%SIM_uses_sticky .or. &
468     thisSim%SIM_uses_GB .or. thisSim%SIM_uses_RF
469     end function SimUsesDirectionalAtoms
470    
471     function SimRequiresPrepairCalc() result(doesit)
472     logical :: doesit
473     doesit = thisSim%SIM_uses_EAM
474     end function SimRequiresPrepairCalc
475    
476     function SimRequiresPostpairCalc() result(doesit)
477     logical :: doesit
478     doesit = thisSim%SIM_uses_RF
479     end function SimRequiresPostpairCalc
480    
481     subroutine InitializeSimGlobals(thisStat)
482     integer, intent(out) :: thisStat
483     integer :: alloc_stat
484    
485     thisStat = 0
486    
487     call FreeSimGlobals()
488 gezelter 1183
489     allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
490     if (alloc_stat /= 0 ) then
491     thisStat = -1
492     return
493     endif
494 mmeineke 377
495     allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
496     if (alloc_stat /= 0 ) then
497     thisStat = -1
498     return
499     endif
500    
501     allocate(excludesGlobal(nExcludes_Global), stat=alloc_stat)
502     if (alloc_stat /= 0 ) then
503     thisStat = -1
504     return
505     endif
506 chuckv 482
507 mmeineke 491 allocate(molMembershipList(nGlobal), stat=alloc_stat)
508 chuckv 482 if (alloc_stat /= 0 ) then
509     thisStat = -1
510     return
511     endif
512 mmeineke 377
513     end subroutine InitializeSimGlobals
514    
515     subroutine FreeSimGlobals()
516    
517     !We free in the opposite order in which we allocate in.
518 gezelter 483
519 gezelter 1183 if (allocated(skipsForAtom)) deallocate(skipsForAtom)
520 tim 1198 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 gezelter 483 if (allocated(molMembershipList)) deallocate(molMembershipList)
527 mmeineke 377 if (allocated(excludesGlobal)) deallocate(excludesGlobal)
528     if (allocated(excludesLocal)) deallocate(excludesLocal)
529 gezelter 1150
530 mmeineke 377 end subroutine FreeSimGlobals
531 gezelter 1150
532 mmeineke 491 pure function getNlocal() result(n)
533     integer :: n
534     n = nLocal
535 mmeineke 377 end function getNlocal
536    
537 gezelter 1150
538 mmeineke 377 end module simulation