ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/simulation.F90
Revision: 2329
Committed: Mon Sep 26 16:42:53 2005 UTC (18 years, 9 months ago) by chuckv
File size: 21066 byte(s)
Log Message:
MPI fix for SimHasAtype in simulation module. We needed a seperate receive buffer.

File Contents

# Content
1 !!
2 !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 !!
4 !! The University of Notre Dame grants you ("Licensee") a
5 !! non-exclusive, royalty free, license to use, modify and
6 !! redistribute this software in source and binary code form, provided
7 !! that the following conditions are met:
8 !!
9 !! 1. Acknowledgement of the program authors must be made in any
10 !! publication of scientific results based in part on use of the
11 !! program. An acceptable form of acknowledgement is citation of
12 !! the article in which the program was described (Matthew
13 !! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 !! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 !! Parallel Simulation Engine for Molecular Dynamics,"
16 !! J. Comput. Chem. 26, pp. 252-271 (2005))
17 !!
18 !! 2. Redistributions of source code must retain the above copyright
19 !! notice, this list of conditions and the following disclaimer.
20 !!
21 !! 3. Redistributions in binary form must reproduce the above copyright
22 !! notice, this list of conditions and the following disclaimer in the
23 !! documentation and/or other materials provided with the
24 !! distribution.
25 !!
26 !! This software is provided "AS IS," without a warranty of any
27 !! kind. All express or implied conditions, representations and
28 !! warranties, including any implied warranty of merchantability,
29 !! fitness for a particular purpose or non-infringement, are hereby
30 !! excluded. The University of Notre Dame and its licensors shall not
31 !! be liable for any damages suffered by licensee as a result of
32 !! using, modifying or distributing the software or its
33 !! derivatives. In no event will the University of Notre Dame or its
34 !! licensors be liable for any lost revenue, profit or data, or for
35 !! direct, indirect, special, consequential, incidental or punitive
36 !! damages, however caused and regardless of the theory of liability,
37 !! arising out of the use of or inability to use software, even if the
38 !! University of Notre Dame has been advised of the possibility of
39 !! such damages.
40 !!
41
42 !! Fortran interface to C entry plug.
43
44 module simulation
45 use definitions
46 use neighborLists
47 use force_globals
48 use vector_class
49 use atype_module
50 use switcheroo
51 #ifdef IS_MPI
52 use mpiSimulation
53 #endif
54
55 implicit none
56 PRIVATE
57
58 #define __FORTRAN90
59 #include "brains/fSimulation.h"
60 #include "UseTheForce/fSwitchingFunction.h"
61 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
62
63 type (simtype), public, save :: thisSim
64
65 logical, save :: simulation_setup_complete = .false.
66
67 integer, public, save :: nLocal, nGlobal
68 integer, public, save :: nGroups, nGroupGlobal
69 integer, public, save :: nExcludes_Global = 0
70 integer, public, save :: nExcludes_Local = 0
71 integer, allocatable, dimension(:,:), public :: excludesLocal
72 integer, allocatable, dimension(:), public :: excludesGlobal
73 integer, allocatable, dimension(:), public :: molMembershipList
74 integer, allocatable, dimension(:), public :: groupListRow
75 integer, allocatable, dimension(:), public :: groupStartRow
76 integer, allocatable, dimension(:), public :: groupListCol
77 integer, allocatable, dimension(:), public :: groupStartCol
78 integer, allocatable, dimension(:), public :: groupListLocal
79 integer, allocatable, dimension(:), public :: groupStartLocal
80 integer, allocatable, dimension(:), public :: nSkipsForAtom
81 integer, allocatable, dimension(:,:), public :: skipsForAtom
82 real(kind=dp), allocatable, dimension(:), public :: mfactRow
83 real(kind=dp), allocatable, dimension(:), public :: mfactCol
84 real(kind=dp), allocatable, dimension(:), public :: mfactLocal
85
86 logical, allocatable, dimension(:) :: simHasAtypeMap
87 #ifdef IS_MPI
88 logical, allocatable, dimension(:) :: simHasAtypeMapTemp
89 #endif
90
91 real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
92 logical, public, save :: boxIsOrthorhombic
93
94 public :: SimulationSetup
95 public :: getNlocal
96 public :: setBox
97 public :: getDielect
98 public :: SimUsesPBC
99
100 public :: SimUsesDirectionalAtoms
101 public :: SimUsesLennardJones
102 public :: SimUsesElectrostatics
103 public :: SimUsesCharges
104 public :: SimUsesDipoles
105 public :: SimUsesSticky
106 public :: SimUsesStickyPower
107 public :: SimUsesGayBerne
108 public :: SimUsesEAM
109 public :: SimUsesShapes
110 public :: SimUsesFLARB
111 public :: SimUsesRF
112 public :: SimRequiresPrepairCalc
113 public :: SimRequiresPostpairCalc
114 public :: SimHasAtype
115
116 contains
117
118 subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
119 CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
120 CmolMembership, Cmfact, CnGroups, CglobalGroupMembership, &
121 status)
122
123 type (simtype) :: setThisSim
124 integer, intent(inout) :: CnGlobal, CnLocal
125 integer, dimension(CnLocal),intent(inout) :: c_idents
126
127 integer :: CnLocalExcludes
128 integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
129 integer :: CnGlobalExcludes
130 integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
131 integer, dimension(CnGlobal),intent(in) :: CmolMembership
132 !! Result status, success = 0, status = -1
133 integer, intent(out) :: status
134 integer :: i, j, me, thisStat, alloc_stat, myNode, id1, id2
135 integer :: ia
136
137 !! mass factors used for molecular cutoffs
138 real ( kind = dp ), dimension(CnLocal) :: Cmfact
139 integer, intent(in):: CnGroups
140 integer, dimension(CnGlobal), intent(in):: CglobalGroupMembership
141 integer :: maxSkipsForAtom, glPointer
142
143 #ifdef IS_MPI
144 integer, allocatable, dimension(:) :: c_idents_Row
145 integer, allocatable, dimension(:) :: c_idents_Col
146 integer :: nAtomsInRow, nGroupsInRow, aid
147 integer :: nAtomsInCol, nGroupsInCol, gid
148 #endif
149
150 simulation_setup_complete = .false.
151 status = 0
152
153 ! copy C struct into fortran type
154
155 nLocal = CnLocal
156 nGlobal = CnGlobal
157 nGroups = CnGroups
158
159 thisSim = setThisSim
160
161 nExcludes_Global = CnGlobalExcludes
162 nExcludes_Local = CnLocalExcludes
163
164 call InitializeForceGlobals(nLocal, thisStat)
165 if (thisStat /= 0) then
166 write(default_error,*) "SimSetup: InitializeForceGlobals error"
167 status = -1
168 return
169 endif
170
171 call InitializeSimGlobals(thisStat)
172 if (thisStat /= 0) then
173 write(default_error,*) "SimSetup: InitializeSimGlobals error"
174 status = -1
175 return
176 endif
177
178 #ifdef IS_MPI
179 ! We can only set up forces if mpiSimulation has been setup.
180 if (.not. isMPISimSet()) then
181 write(default_error,*) "MPI is not set"
182 status = -1
183 return
184 endif
185 nAtomsInRow = getNatomsInRow(plan_atom_row)
186 nAtomsInCol = getNatomsInCol(plan_atom_col)
187 nGroupsInRow = getNgroupsInRow(plan_group_row)
188 nGroupsInCol = getNgroupsInCol(plan_group_col)
189 mynode = getMyNode()
190
191 allocate(c_idents_Row(nAtomsInRow),stat=alloc_stat)
192 if (alloc_stat /= 0 ) then
193 status = -1
194 return
195 endif
196
197 allocate(c_idents_Col(nAtomsInCol),stat=alloc_stat)
198 if (alloc_stat /= 0 ) then
199 status = -1
200 return
201 endif
202
203 call gather(c_idents, c_idents_Row, plan_atom_row)
204 call gather(c_idents, c_idents_Col, plan_atom_col)
205
206 do i = 1, nAtomsInRow
207 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
208 atid_Row(i) = me
209 enddo
210
211 do i = 1, nAtomsInCol
212 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
213 atid_Col(i) = me
214 enddo
215
216 !! free temporary ident arrays
217 if (allocated(c_idents_Col)) then
218 deallocate(c_idents_Col)
219 end if
220 if (allocated(c_idents_Row)) then
221 deallocate(c_idents_Row)
222 endif
223
224 #endif
225
226 #ifdef IS_MPI
227 allocate(groupStartRow(nGroupsInRow+1),stat=alloc_stat)
228 if (alloc_stat /= 0 ) then
229 status = -1
230 return
231 endif
232 allocate(groupStartCol(nGroupsInCol+1),stat=alloc_stat)
233 if (alloc_stat /= 0 ) then
234 status = -1
235 return
236 endif
237 allocate(groupListRow(nAtomsInRow),stat=alloc_stat)
238 if (alloc_stat /= 0 ) then
239 status = -1
240 return
241 endif
242 allocate(groupListCol(nAtomsInCol),stat=alloc_stat)
243 if (alloc_stat /= 0 ) then
244 status = -1
245 return
246 endif
247 allocate(mfactRow(nAtomsInRow),stat=alloc_stat)
248 if (alloc_stat /= 0 ) then
249 status = -1
250 return
251 endif
252 allocate(mfactCol(nAtomsInCol),stat=alloc_stat)
253 if (alloc_stat /= 0 ) then
254 status = -1
255 return
256 endif
257 allocate(mfactLocal(nLocal),stat=alloc_stat)
258 if (alloc_stat /= 0 ) then
259 status = -1
260 return
261 endif
262
263 glPointer = 1
264
265 do i = 1, nGroupsInRow
266
267 gid = GroupRowToGlobal(i)
268 groupStartRow(i) = glPointer
269
270 do j = 1, nAtomsInRow
271 aid = AtomRowToGlobal(j)
272 if (CglobalGroupMembership(aid) .eq. gid) then
273 groupListRow(glPointer) = j
274 glPointer = glPointer + 1
275 endif
276 enddo
277 enddo
278 groupStartRow(nGroupsInRow+1) = nAtomsInRow + 1
279
280 glPointer = 1
281
282 do i = 1, nGroupsInCol
283
284 gid = GroupColToGlobal(i)
285 groupStartCol(i) = glPointer
286
287 do j = 1, nAtomsInCol
288 aid = AtomColToGlobal(j)
289 if (CglobalGroupMembership(aid) .eq. gid) then
290 groupListCol(glPointer) = j
291 glPointer = glPointer + 1
292 endif
293 enddo
294 enddo
295 groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1
296
297 mfactLocal = Cmfact
298
299 call gather(mfactLocal, mfactRow, plan_atom_row)
300 call gather(mfactLocal, mfactCol, plan_atom_col)
301
302 if (allocated(mfactLocal)) then
303 deallocate(mfactLocal)
304 end if
305 #else
306 allocate(groupStartRow(nGroups+1),stat=alloc_stat)
307 if (alloc_stat /= 0 ) then
308 status = -1
309 return
310 endif
311 allocate(groupStartCol(nGroups+1),stat=alloc_stat)
312 if (alloc_stat /= 0 ) then
313 status = -1
314 return
315 endif
316 allocate(groupListRow(nLocal),stat=alloc_stat)
317 if (alloc_stat /= 0 ) then
318 status = -1
319 return
320 endif
321 allocate(groupListCol(nLocal),stat=alloc_stat)
322 if (alloc_stat /= 0 ) then
323 status = -1
324 return
325 endif
326 allocate(mfactRow(nLocal),stat=alloc_stat)
327 if (alloc_stat /= 0 ) then
328 status = -1
329 return
330 endif
331 allocate(mfactCol(nLocal),stat=alloc_stat)
332 if (alloc_stat /= 0 ) then
333 status = -1
334 return
335 endif
336 allocate(mfactLocal(nLocal),stat=alloc_stat)
337 if (alloc_stat /= 0 ) then
338 status = -1
339 return
340 endif
341
342 glPointer = 1
343 do i = 1, nGroups
344 groupStartRow(i) = glPointer
345 groupStartCol(i) = glPointer
346 do j = 1, nLocal
347 if (CglobalGroupMembership(j) .eq. i) then
348 groupListRow(glPointer) = j
349 groupListCol(glPointer) = j
350 glPointer = glPointer + 1
351 endif
352 enddo
353 enddo
354 groupStartRow(nGroups+1) = nLocal + 1
355 groupStartCol(nGroups+1) = nLocal + 1
356
357 do i = 1, nLocal
358 mfactRow(i) = Cmfact(i)
359 mfactCol(i) = Cmfact(i)
360 end do
361
362 #endif
363
364
365 ! We build the local atid's for both mpi and nonmpi
366 do i = 1, nLocal
367
368 me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
369 atid(i) = me
370
371 enddo
372
373 do i = 1, nExcludes_Local
374 excludesLocal(1,i) = CexcludesLocal(1,i)
375 excludesLocal(2,i) = CexcludesLocal(2,i)
376 enddo
377
378 #ifdef IS_MPI
379 allocate(nSkipsForAtom(nAtomsInRow), stat=alloc_stat)
380 #else
381 allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
382 #endif
383 if (alloc_stat /= 0 ) then
384 thisStat = -1
385 write(*,*) 'Could not allocate nSkipsForAtom array'
386 return
387 endif
388
389 maxSkipsForAtom = 0
390 #ifdef IS_MPI
391 do j = 1, nAtomsInRow
392 #else
393 do j = 1, nLocal
394 #endif
395 nSkipsForAtom(j) = 0
396 #ifdef IS_MPI
397 id1 = AtomRowToGlobal(j)
398 #else
399 id1 = j
400 #endif
401 do i = 1, nExcludes_Local
402 if (excludesLocal(1,i) .eq. id1 ) then
403 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
404
405 if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
406 maxSkipsForAtom = nSkipsForAtom(j)
407 endif
408 endif
409 if (excludesLocal(2,i) .eq. id1 ) then
410 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
411
412 if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
413 maxSkipsForAtom = nSkipsForAtom(j)
414 endif
415 endif
416 end do
417 enddo
418
419 #ifdef IS_MPI
420 allocate(skipsForAtom(nAtomsInRow, maxSkipsForAtom), stat=alloc_stat)
421 #else
422 allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
423 #endif
424 if (alloc_stat /= 0 ) then
425 write(*,*) 'Could not allocate skipsForAtom array'
426 return
427 endif
428
429 #ifdef IS_MPI
430 do j = 1, nAtomsInRow
431 #else
432 do j = 1, nLocal
433 #endif
434 nSkipsForAtom(j) = 0
435 #ifdef IS_MPI
436 id1 = AtomRowToGlobal(j)
437 #else
438 id1 = j
439 #endif
440 do i = 1, nExcludes_Local
441 if (excludesLocal(1,i) .eq. id1 ) then
442 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
443 ! exclude lists have global ID's so this line is
444 ! the same in MPI and non-MPI
445 id2 = excludesLocal(2,i)
446 skipsForAtom(j, nSkipsForAtom(j)) = id2
447 endif
448 if (excludesLocal(2, i) .eq. id1 ) then
449 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
450 ! exclude lists have global ID's so this line is
451 ! the same in MPI and non-MPI
452 id2 = excludesLocal(1,i)
453 skipsForAtom(j, nSkipsForAtom(j)) = id2
454 endif
455 end do
456 enddo
457
458 do i = 1, nExcludes_Global
459 excludesGlobal(i) = CexcludesGlobal(i)
460 enddo
461
462 do i = 1, nGlobal
463 molMemberShipList(i) = CmolMembership(i)
464 enddo
465
466 call createSimHasAtype(alloc_stat)
467 if (alloc_stat /= 0) then
468 status = -1
469 end if
470
471 if (status == 0) simulation_setup_complete = .true.
472
473 end subroutine SimulationSetup
474
475 subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
476 real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
477 integer :: cBoxIsOrthorhombic
478 integer :: smallest, status, i
479
480 Hmat = cHmat
481 HmatInv = cHmatInv
482 if (cBoxIsOrthorhombic .eq. 0 ) then
483 boxIsOrthorhombic = .false.
484 else
485 boxIsOrthorhombic = .true.
486 endif
487
488 return
489 end subroutine setBox
490
491 function getDielect() result(dielect)
492 real( kind = dp ) :: dielect
493 dielect = thisSim%dielect
494 end function getDielect
495
496 function SimUsesPBC() result(doesit)
497 logical :: doesit
498 doesit = thisSim%SIM_uses_PBC
499 end function SimUsesPBC
500
501 function SimUsesDirectionalAtoms() result(doesit)
502 logical :: doesit
503 doesit = thisSim%SIM_uses_dipoles .or. thisSim%SIM_uses_Sticky .or. &
504 thisSim%SIM_uses_StickyPower .or. &
505 thisSim%SIM_uses_GayBerne .or. thisSim%SIM_uses_Shapes
506 end function SimUsesDirectionalAtoms
507
508 function SimUsesLennardJones() result(doesit)
509 logical :: doesit
510 doesit = thisSim%SIM_uses_LennardJones
511 end function SimUsesLennardJones
512
513 function SimUsesElectrostatics() result(doesit)
514 logical :: doesit
515 doesit = thisSim%SIM_uses_Electrostatics
516 end function SimUsesElectrostatics
517
518 function SimUsesCharges() result(doesit)
519 logical :: doesit
520 doesit = thisSim%SIM_uses_Charges
521 end function SimUsesCharges
522
523 function SimUsesDipoles() result(doesit)
524 logical :: doesit
525 doesit = thisSim%SIM_uses_Dipoles
526 end function SimUsesDipoles
527
528 function SimUsesSticky() result(doesit)
529 logical :: doesit
530 doesit = thisSim%SIM_uses_Sticky
531 end function SimUsesSticky
532
533 function SimUsesStickyPower() result(doesit)
534 logical :: doesit
535 doesit = thisSim%SIM_uses_StickyPower
536 end function SimUsesStickyPower
537
538 function SimUsesGayBerne() result(doesit)
539 logical :: doesit
540 doesit = thisSim%SIM_uses_GayBerne
541 end function SimUsesGayBerne
542
543 function SimUsesEAM() result(doesit)
544 logical :: doesit
545 doesit = thisSim%SIM_uses_EAM
546 end function SimUsesEAM
547
548 function SimUsesShapes() result(doesit)
549 logical :: doesit
550 doesit = thisSim%SIM_uses_Shapes
551 end function SimUsesShapes
552
553 function SimUsesFLARB() result(doesit)
554 logical :: doesit
555 doesit = thisSim%SIM_uses_FLARB
556 end function SimUsesFLARB
557
558 function SimUsesRF() result(doesit)
559 logical :: doesit
560 doesit = thisSim%SIM_uses_RF
561 end function SimUsesRF
562
563 function SimRequiresPrepairCalc() result(doesit)
564 logical :: doesit
565 doesit = thisSim%SIM_uses_EAM
566 end function SimRequiresPrepairCalc
567
568 function SimRequiresPostpairCalc() result(doesit)
569 logical :: doesit
570 doesit = thisSim%SIM_uses_RF
571 end function SimRequiresPostpairCalc
572
573 ! Function returns true if the simulation has this atype
574 function SimHasAtype(thisAtype) result(doesit)
575 logical :: doesit
576 integer :: thisAtype
577 doesit = .false.
578 if(.not.allocated(SimHasAtypeMap)) return
579
580 doesit = SimHasAtypeMap(thisAtype)
581
582 end function SimHasAtype
583
584 subroutine createSimHasAtype(status)
585 integer, intent(out) :: status
586 integer :: alloc_stat
587 integer :: me_i
588 integer :: mpiErrors
589 integer :: nAtypes
590 status = 0
591
592 nAtypes = getSize(atypes)
593 ! Setup logical map for atypes in simulation
594 if (.not.allocated(SimHasAtypeMap)) then
595 allocate(SimHasAtypeMap(nAtypes),stat=alloc_stat)
596 if (alloc_stat /= 0 ) then
597 status = -1
598 return
599 end if
600 SimHasAtypeMap = .false.
601 end if
602
603 ! Loop through the local atoms and grab the atypes present
604 do me_i = 1,nLocal
605 SimHasAtypeMap(atid(me_i)) = .true.
606 end do
607 ! For MPI, we need to know all possible atypes present in
608 ! simulation on all processors. Use LOR operation to set map.
609 #ifdef IS_MPI
610 if (.not.allocated(SimHasAtypeMapTemp)) then
611 allocate(SimHasAtypeMapTemp(nAtypes),stat=alloc_stat)
612 if (alloc_stat /= 0 ) then
613 status = -1
614 return
615 end if
616 end if
617 call mpi_allreduce(SimHasAtypeMap, SimHasAtypeMaptemp, nAtypes, &
618 mpi_logical, MPI_LOR, mpi_comm_world, mpiErrors)
619 simHasAtypeMap = simHasAtypeMapTemp
620 deallocate(simHasAtypeMapTemp)
621 #endif
622 end subroutine createSimHasAtype
623
624 subroutine InitializeSimGlobals(thisStat)
625 integer, intent(out) :: thisStat
626 integer :: alloc_stat
627
628 thisStat = 0
629
630 call FreeSimGlobals()
631
632 allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
633 if (alloc_stat /= 0 ) then
634 thisStat = -1
635 return
636 endif
637
638 allocate(excludesGlobal(nExcludes_Global), stat=alloc_stat)
639 if (alloc_stat /= 0 ) then
640 thisStat = -1
641 return
642 endif
643
644 allocate(molMembershipList(nGlobal), stat=alloc_stat)
645 if (alloc_stat /= 0 ) then
646 thisStat = -1
647 return
648 endif
649
650 end subroutine InitializeSimGlobals
651
652 subroutine FreeSimGlobals()
653
654 !We free in the opposite order in which we allocate in.
655
656 if (allocated(skipsForAtom)) deallocate(skipsForAtom)
657 if (allocated(nSkipsForAtom)) deallocate(nSkipsForAtom)
658 if (allocated(mfactLocal)) deallocate(mfactLocal)
659 if (allocated(mfactCol)) deallocate(mfactCol)
660 if (allocated(mfactRow)) deallocate(mfactRow)
661 if (allocated(groupListCol)) deallocate(groupListCol)
662 if (allocated(groupListRow)) deallocate(groupListRow)
663 if (allocated(groupStartCol)) deallocate(groupStartCol)
664 if (allocated(groupStartRow)) deallocate(groupStartRow)
665 if (allocated(molMembershipList)) deallocate(molMembershipList)
666 if (allocated(excludesGlobal)) deallocate(excludesGlobal)
667 if (allocated(excludesLocal)) deallocate(excludesLocal)
668
669 end subroutine FreeSimGlobals
670
671 pure function getNlocal() result(n)
672 integer :: n
673 n = nLocal
674 end function getNlocal
675
676
677
678
679
680 end module simulation