ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/DarkSide/simulation.F90
Revision: 2402
Committed: Tue Nov 1 19:09:30 2005 UTC (18 years, 10 months ago) by chrisfen
File size: 21304 byte(s)
Log Message:
still fixing up wolf...

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