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, 1 month ago) by tim
File size: 15231 byte(s)
Log Message:
Bug fix for SkipList

File Contents

# Content
1 !! 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 use switcheroo
10 #ifdef IS_MPI
11 use mpiSimulation
12 #endif
13
14 implicit none
15 PRIVATE
16
17 #define __FORTRAN90
18 #include "fSimulation.h"
19 #include "fSwitchingFunction.h"
20
21 type (simtype), public, save :: thisSim
22
23 logical, save :: simulation_setup_complete = .false.
24
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
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), public, dimension(3,3), save :: Hmat, HmatInv
45 logical, public, save :: boxIsOrthorhombic
46
47 public :: SimulationSetup
48 public :: getNlocal
49 public :: setBox
50 public :: getDielect
51 public :: SimUsesPBC
52 public :: SimUsesLJ
53 public :: SimUsesCharges
54 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 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) :: 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, 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 :: 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
108
109 nExcludes_Global = CnGlobalExcludes
110 nExcludes_Local = CnLocalExcludes
111
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
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 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(nAtomsInRow),stat=alloc_stat)
140 if (alloc_stat /= 0 ) then
141 status = -1
142 return
143 endif
144
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_atom_row)
152 call gather(c_idents, c_idents_Col, plan_atom_col)
153
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, nAtomsInCol
160 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
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 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 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 #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(cHmat, cHmatInv, cBoxIsOrthorhombic)
405 real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
406 integer :: cBoxIsOrthorhombic
407 integer :: smallest, status, i
408
409 Hmat = cHmat
410 HmatInv = cHmatInv
411 if (cBoxIsOrthorhombic .eq. 0 ) then
412 boxIsOrthorhombic = .false.
413 else
414 boxIsOrthorhombic = .true.
415 endif
416
417 return
418 end subroutine setBox
419
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 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
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
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
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
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
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
530 end subroutine FreeSimGlobals
531
532 pure function getNlocal() result(n)
533 integer :: n
534 n = nLocal
535 end function getNlocal
536
537
538 end module simulation