ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/simulation.F90
Revision: 2274
Committed: Wed Aug 17 15:26:42 2005 UTC (18 years, 11 months ago) by gezelter
File size: 20569 byte(s)
Log Message:
added fCutoffPolicy.h

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
62 type (simtype), public, save :: thisSim
63
64 logical, save :: simulation_setup_complete = .false.
65
66 integer, public, save :: nLocal, nGlobal
67 integer, public, save :: nGroups, nGroupGlobal
68 integer, public, save :: nExcludes_Global = 0
69 integer, public, save :: nExcludes_Local = 0
70 integer, allocatable, dimension(:,:), public :: excludesLocal
71 integer, allocatable, dimension(:), public :: excludesGlobal
72 integer, allocatable, dimension(:), public :: molMembershipList
73 integer, allocatable, dimension(:), public :: groupListRow
74 integer, allocatable, dimension(:), public :: groupStartRow
75 integer, allocatable, dimension(:), public :: groupListCol
76 integer, allocatable, dimension(:), public :: groupStartCol
77 integer, allocatable, dimension(:), public :: groupListLocal
78 integer, allocatable, dimension(:), public :: groupStartLocal
79 integer, allocatable, dimension(:), public :: nSkipsForAtom
80 integer, allocatable, dimension(:,:), public :: skipsForAtom
81 real(kind=dp), allocatable, dimension(:), public :: mfactRow
82 real(kind=dp), allocatable, dimension(:), public :: mfactCol
83 real(kind=dp), allocatable, dimension(:), public :: mfactLocal
84
85 logical, allocatable, dimension(:) :: simHasAtypeMap
86 real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
87 logical, public, save :: boxIsOrthorhombic
88
89 public :: SimulationSetup
90 public :: getNlocal
91 public :: setBox
92 public :: getDielect
93 public :: SimUsesPBC
94
95 public :: SimUsesDirectionalAtoms
96 public :: SimUsesLennardJones
97 public :: SimUsesElectrostatics
98 public :: SimUsesCharges
99 public :: SimUsesDipoles
100 public :: SimUsesSticky
101 public :: SimUsesStickyPower
102 public :: SimUsesGayBerne
103 public :: SimUsesEAM
104 public :: SimUsesShapes
105 public :: SimUsesFLARB
106 public :: SimUsesRF
107 public :: SimRequiresPrepairCalc
108 public :: SimRequiresPostpairCalc
109 public :: SimHasAtype
110
111 contains
112
113 subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
114 CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
115 CmolMembership, Cmfact, CnGroups, CglobalGroupMembership, &
116 status)
117
118 type (simtype) :: setThisSim
119 integer, intent(inout) :: CnGlobal, CnLocal
120 integer, dimension(CnLocal),intent(inout) :: c_idents
121
122 integer :: CnLocalExcludes
123 integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
124 integer :: CnGlobalExcludes
125 integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
126 integer, dimension(CnGlobal),intent(in) :: CmolMembership
127 !! Result status, success = 0, status = -1
128 integer, intent(out) :: status
129 integer :: i, j, me, thisStat, alloc_stat, myNode, id1, id2
130 integer :: ia
131
132 !! mass factors used for molecular cutoffs
133 real ( kind = dp ), dimension(CnLocal) :: Cmfact
134 integer, intent(in):: CnGroups
135 integer, dimension(CnGlobal), intent(in):: CglobalGroupMembership
136 integer :: maxSkipsForAtom, glPointer
137
138 #ifdef IS_MPI
139 integer, allocatable, dimension(:) :: c_idents_Row
140 integer, allocatable, dimension(:) :: c_idents_Col
141 integer :: nAtomsInRow, nGroupsInRow, aid
142 integer :: nAtomsInCol, nGroupsInCol, gid
143 #endif
144
145 simulation_setup_complete = .false.
146 status = 0
147
148 ! copy C struct into fortran type
149
150 nLocal = CnLocal
151 nGlobal = CnGlobal
152 nGroups = CnGroups
153
154 thisSim = setThisSim
155
156 nExcludes_Global = CnGlobalExcludes
157 nExcludes_Local = CnLocalExcludes
158
159 call InitializeForceGlobals(nLocal, thisStat)
160 if (thisStat /= 0) then
161 write(default_error,*) "SimSetup: InitializeForceGlobals error"
162 status = -1
163 return
164 endif
165
166 call InitializeSimGlobals(thisStat)
167 if (thisStat /= 0) then
168 write(default_error,*) "SimSetup: InitializeSimGlobals error"
169 status = -1
170 return
171 endif
172
173 #ifdef IS_MPI
174 ! We can only set up forces if mpiSimulation has been setup.
175 if (.not. isMPISimSet()) then
176 write(default_error,*) "MPI is not set"
177 status = -1
178 return
179 endif
180 nAtomsInRow = getNatomsInRow(plan_atom_row)
181 nAtomsInCol = getNatomsInCol(plan_atom_col)
182 nGroupsInRow = getNgroupsInRow(plan_group_row)
183 nGroupsInCol = getNgroupsInCol(plan_group_col)
184 mynode = getMyNode()
185
186 allocate(c_idents_Row(nAtomsInRow),stat=alloc_stat)
187 if (alloc_stat /= 0 ) then
188 status = -1
189 return
190 endif
191
192 allocate(c_idents_Col(nAtomsInCol),stat=alloc_stat)
193 if (alloc_stat /= 0 ) then
194 status = -1
195 return
196 endif
197
198 call gather(c_idents, c_idents_Row, plan_atom_row)
199 call gather(c_idents, c_idents_Col, plan_atom_col)
200
201 do i = 1, nAtomsInRow
202 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
203 atid_Row(i) = me
204 enddo
205
206 do i = 1, nAtomsInCol
207 me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
208 atid_Col(i) = me
209 enddo
210
211 !! free temporary ident arrays
212 if (allocated(c_idents_Col)) then
213 deallocate(c_idents_Col)
214 end if
215 if (allocated(c_idents_Row)) then
216 deallocate(c_idents_Row)
217 endif
218
219 #endif
220
221 #ifdef IS_MPI
222 allocate(groupStartRow(nGroupsInRow+1),stat=alloc_stat)
223 if (alloc_stat /= 0 ) then
224 status = -1
225 return
226 endif
227 allocate(groupStartCol(nGroupsInCol+1),stat=alloc_stat)
228 if (alloc_stat /= 0 ) then
229 status = -1
230 return
231 endif
232 allocate(groupListRow(nAtomsInRow),stat=alloc_stat)
233 if (alloc_stat /= 0 ) then
234 status = -1
235 return
236 endif
237 allocate(groupListCol(nAtomsInCol),stat=alloc_stat)
238 if (alloc_stat /= 0 ) then
239 status = -1
240 return
241 endif
242 allocate(mfactRow(nAtomsInRow),stat=alloc_stat)
243 if (alloc_stat /= 0 ) then
244 status = -1
245 return
246 endif
247 allocate(mfactCol(nAtomsInCol),stat=alloc_stat)
248 if (alloc_stat /= 0 ) then
249 status = -1
250 return
251 endif
252 allocate(mfactLocal(nLocal),stat=alloc_stat)
253 if (alloc_stat /= 0 ) then
254 status = -1
255 return
256 endif
257
258 glPointer = 1
259
260 do i = 1, nGroupsInRow
261
262 gid = GroupRowToGlobal(i)
263 groupStartRow(i) = glPointer
264
265 do j = 1, nAtomsInRow
266 aid = AtomRowToGlobal(j)
267 if (CglobalGroupMembership(aid) .eq. gid) then
268 groupListRow(glPointer) = j
269 glPointer = glPointer + 1
270 endif
271 enddo
272 enddo
273 groupStartRow(nGroupsInRow+1) = nAtomsInRow + 1
274
275 glPointer = 1
276
277 do i = 1, nGroupsInCol
278
279 gid = GroupColToGlobal(i)
280 groupStartCol(i) = glPointer
281
282 do j = 1, nAtomsInCol
283 aid = AtomColToGlobal(j)
284 if (CglobalGroupMembership(aid) .eq. gid) then
285 groupListCol(glPointer) = j
286 glPointer = glPointer + 1
287 endif
288 enddo
289 enddo
290 groupStartCol(nGroupsInCol+1) = nAtomsInCol + 1
291
292 mfactLocal = Cmfact
293
294 call gather(mfactLocal, mfactRow, plan_atom_row)
295 call gather(mfactLocal, mfactCol, plan_atom_col)
296
297 if (allocated(mfactLocal)) then
298 deallocate(mfactLocal)
299 end if
300 #else
301 allocate(groupStartRow(nGroups+1),stat=alloc_stat)
302 if (alloc_stat /= 0 ) then
303 status = -1
304 return
305 endif
306 allocate(groupStartCol(nGroups+1),stat=alloc_stat)
307 if (alloc_stat /= 0 ) then
308 status = -1
309 return
310 endif
311 allocate(groupListRow(nLocal),stat=alloc_stat)
312 if (alloc_stat /= 0 ) then
313 status = -1
314 return
315 endif
316 allocate(groupListCol(nLocal),stat=alloc_stat)
317 if (alloc_stat /= 0 ) then
318 status = -1
319 return
320 endif
321 allocate(mfactRow(nLocal),stat=alloc_stat)
322 if (alloc_stat /= 0 ) then
323 status = -1
324 return
325 endif
326 allocate(mfactCol(nLocal),stat=alloc_stat)
327 if (alloc_stat /= 0 ) then
328 status = -1
329 return
330 endif
331 allocate(mfactLocal(nLocal),stat=alloc_stat)
332 if (alloc_stat /= 0 ) then
333 status = -1
334 return
335 endif
336
337 glPointer = 1
338 do i = 1, nGroups
339 groupStartRow(i) = glPointer
340 groupStartCol(i) = glPointer
341 do j = 1, nLocal
342 if (CglobalGroupMembership(j) .eq. i) then
343 groupListRow(glPointer) = j
344 groupListCol(glPointer) = j
345 glPointer = glPointer + 1
346 endif
347 enddo
348 enddo
349 groupStartRow(nGroups+1) = nLocal + 1
350 groupStartCol(nGroups+1) = nLocal + 1
351
352 do i = 1, nLocal
353 mfactRow(i) = Cmfact(i)
354 mfactCol(i) = Cmfact(i)
355 end do
356
357 #endif
358
359
360 ! We build the local atid's for both mpi and nonmpi
361 do i = 1, nLocal
362
363 me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
364 atid(i) = me
365
366 enddo
367
368 do i = 1, nExcludes_Local
369 excludesLocal(1,i) = CexcludesLocal(1,i)
370 excludesLocal(2,i) = CexcludesLocal(2,i)
371 enddo
372
373 #ifdef IS_MPI
374 allocate(nSkipsForAtom(nAtomsInRow), stat=alloc_stat)
375 #else
376 allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
377 #endif
378 if (alloc_stat /= 0 ) then
379 thisStat = -1
380 write(*,*) 'Could not allocate nSkipsForAtom array'
381 return
382 endif
383
384 maxSkipsForAtom = 0
385 #ifdef IS_MPI
386 do j = 1, nAtomsInRow
387 #else
388 do j = 1, nLocal
389 #endif
390 nSkipsForAtom(j) = 0
391 #ifdef IS_MPI
392 id1 = AtomRowToGlobal(j)
393 #else
394 id1 = j
395 #endif
396 do i = 1, nExcludes_Local
397 if (excludesLocal(1,i) .eq. id1 ) then
398 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
399
400 if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
401 maxSkipsForAtom = nSkipsForAtom(j)
402 endif
403 endif
404 if (excludesLocal(2,i) .eq. id1 ) then
405 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
406
407 if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
408 maxSkipsForAtom = nSkipsForAtom(j)
409 endif
410 endif
411 end do
412 enddo
413
414 #ifdef IS_MPI
415 allocate(skipsForAtom(nAtomsInRow, maxSkipsForAtom), stat=alloc_stat)
416 #else
417 allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
418 #endif
419 if (alloc_stat /= 0 ) then
420 write(*,*) 'Could not allocate skipsForAtom array'
421 return
422 endif
423
424 #ifdef IS_MPI
425 do j = 1, nAtomsInRow
426 #else
427 do j = 1, nLocal
428 #endif
429 nSkipsForAtom(j) = 0
430 #ifdef IS_MPI
431 id1 = AtomRowToGlobal(j)
432 #else
433 id1 = j
434 #endif
435 do i = 1, nExcludes_Local
436 if (excludesLocal(1,i) .eq. id1 ) then
437 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
438 ! exclude lists have global ID's so this line is
439 ! the same in MPI and non-MPI
440 id2 = excludesLocal(2,i)
441 skipsForAtom(j, nSkipsForAtom(j)) = id2
442 endif
443 if (excludesLocal(2, i) .eq. id1 ) then
444 nSkipsForAtom(j) = nSkipsForAtom(j) + 1
445 ! exclude lists have global ID's so this line is
446 ! the same in MPI and non-MPI
447 id2 = excludesLocal(1,i)
448 skipsForAtom(j, nSkipsForAtom(j)) = id2
449 endif
450 end do
451 enddo
452
453 do i = 1, nExcludes_Global
454 excludesGlobal(i) = CexcludesGlobal(i)
455 enddo
456
457 do i = 1, nGlobal
458 molMemberShipList(i) = CmolMembership(i)
459 enddo
460
461 call createSimHasAtype(alloc_stat)
462 if (alloc_stat /= 0) then
463 status = -1
464 end if
465
466 if (status == 0) simulation_setup_complete = .true.
467
468 end subroutine SimulationSetup
469
470 subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
471 real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
472 integer :: cBoxIsOrthorhombic
473 integer :: smallest, status, i
474
475 Hmat = cHmat
476 HmatInv = cHmatInv
477 if (cBoxIsOrthorhombic .eq. 0 ) then
478 boxIsOrthorhombic = .false.
479 else
480 boxIsOrthorhombic = .true.
481 endif
482
483 return
484 end subroutine setBox
485
486 function getDielect() result(dielect)
487 real( kind = dp ) :: dielect
488 dielect = thisSim%dielect
489 end function getDielect
490
491 function SimUsesPBC() result(doesit)
492 logical :: doesit
493 doesit = thisSim%SIM_uses_PBC
494 end function SimUsesPBC
495
496 function SimUsesDirectionalAtoms() result(doesit)
497 logical :: doesit
498 doesit = thisSim%SIM_uses_dipoles .or. thisSim%SIM_uses_Sticky .or. &
499 thisSim%SIM_uses_StickyPower .or. &
500 thisSim%SIM_uses_GayBerne .or. thisSim%SIM_uses_Shapes
501 end function SimUsesDirectionalAtoms
502
503 function SimUsesLennardJones() result(doesit)
504 logical :: doesit
505 doesit = thisSim%SIM_uses_LennardJones
506 end function SimUsesLennardJones
507
508 function SimUsesElectrostatics() result(doesit)
509 logical :: doesit
510 doesit = thisSim%SIM_uses_Electrostatics
511 end function SimUsesElectrostatics
512
513 function SimUsesCharges() result(doesit)
514 logical :: doesit
515 doesit = thisSim%SIM_uses_Charges
516 end function SimUsesCharges
517
518 function SimUsesDipoles() result(doesit)
519 logical :: doesit
520 doesit = thisSim%SIM_uses_Dipoles
521 end function SimUsesDipoles
522
523 function SimUsesSticky() result(doesit)
524 logical :: doesit
525 doesit = thisSim%SIM_uses_Sticky
526 end function SimUsesSticky
527
528 function SimUsesStickyPower() result(doesit)
529 logical :: doesit
530 doesit = thisSim%SIM_uses_StickyPower
531 end function SimUsesStickyPower
532
533 function SimUsesGayBerne() result(doesit)
534 logical :: doesit
535 doesit = thisSim%SIM_uses_GayBerne
536 end function SimUsesGayBerne
537
538 function SimUsesEAM() result(doesit)
539 logical :: doesit
540 doesit = thisSim%SIM_uses_EAM
541 end function SimUsesEAM
542
543 function SimUsesShapes() result(doesit)
544 logical :: doesit
545 doesit = thisSim%SIM_uses_Shapes
546 end function SimUsesShapes
547
548 function SimUsesFLARB() result(doesit)
549 logical :: doesit
550 doesit = thisSim%SIM_uses_FLARB
551 end function SimUsesFLARB
552
553 function SimUsesRF() result(doesit)
554 logical :: doesit
555 doesit = thisSim%SIM_uses_RF
556 end function SimUsesRF
557
558 function SimRequiresPrepairCalc() result(doesit)
559 logical :: doesit
560 doesit = thisSim%SIM_uses_EAM
561 end function SimRequiresPrepairCalc
562
563 function SimRequiresPostpairCalc() result(doesit)
564 logical :: doesit
565 doesit = thisSim%SIM_uses_RF
566 end function SimRequiresPostpairCalc
567
568 ! Function returns true if the simulation has this atype
569 function SimHasAtype(thisAtype) result(doesit)
570 logical :: doesit
571 integer :: thisAtype
572 doesit = .false.
573 if(.not.allocated(SimHasAtypeMap)) return
574
575 doesit = SimHasAtypeMap(thisAtype)
576
577 end function SimHasAtype
578
579 subroutine createSimHasAtype(status)
580 integer, intent(out) :: status
581 integer :: alloc_stat
582 integer :: me_i
583 integer :: mpiErrors
584 integer :: nAtypes
585 status = 0
586
587 nAtypes = getSize(atypes)
588 ! Setup logical map for atypes in simulation
589 if (.not.allocated(SimHasAtypeMap)) then
590 allocate(SimHasAtypeMap(nAtypes),stat=alloc_stat)
591 if (alloc_stat /= 0 ) then
592 status = -1
593 return
594 end if
595 SimHasAtypeMap = .false.
596 end if
597 ! Loop through the local atoms and grab the atypes present
598 do me_i = 1,nLocal
599 SimHasAtypeMap(atid(me_i)) = .true.
600 end do
601 ! For MPI, we need to know all possible atypes present in
602 ! simulation on all processors. Use LOR operation to set map.
603 #ifdef IS_MPI
604 call mpi_allreduce(SimHasAtypeMap, SimHasAtypeMap, nAtypes, &
605 mpi_logical, MPI_LOR, mpi_comm_world, mpiErrors)
606 #endif
607 end subroutine createSimHasAtype
608
609 subroutine InitializeSimGlobals(thisStat)
610 integer, intent(out) :: thisStat
611 integer :: alloc_stat
612
613 thisStat = 0
614
615 call FreeSimGlobals()
616
617 allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
618 if (alloc_stat /= 0 ) then
619 thisStat = -1
620 return
621 endif
622
623 allocate(excludesGlobal(nExcludes_Global), stat=alloc_stat)
624 if (alloc_stat /= 0 ) then
625 thisStat = -1
626 return
627 endif
628
629 allocate(molMembershipList(nGlobal), stat=alloc_stat)
630 if (alloc_stat /= 0 ) then
631 thisStat = -1
632 return
633 endif
634
635 end subroutine InitializeSimGlobals
636
637 subroutine FreeSimGlobals()
638
639 !We free in the opposite order in which we allocate in.
640
641 if (allocated(skipsForAtom)) deallocate(skipsForAtom)
642 if (allocated(nSkipsForAtom)) deallocate(nSkipsForAtom)
643 if (allocated(mfactLocal)) deallocate(mfactLocal)
644 if (allocated(mfactCol)) deallocate(mfactCol)
645 if (allocated(mfactRow)) deallocate(mfactRow)
646 if (allocated(groupListCol)) deallocate(groupListCol)
647 if (allocated(groupListRow)) deallocate(groupListRow)
648 if (allocated(groupStartCol)) deallocate(groupStartCol)
649 if (allocated(groupStartRow)) deallocate(groupStartRow)
650 if (allocated(molMembershipList)) deallocate(molMembershipList)
651 if (allocated(excludesGlobal)) deallocate(excludesGlobal)
652 if (allocated(excludesLocal)) deallocate(excludesLocal)
653
654 end subroutine FreeSimGlobals
655
656 pure function getNlocal() result(n)
657 integer :: n
658 n = nLocal
659 end function getNlocal
660
661
662
663
664
665 end module simulation