ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/simulation.F90
Revision: 2204
Committed: Fri Apr 15 22:04:00 2005 UTC (19 years, 3 months ago) by gezelter
File size: 18628 byte(s)
Log Message:
xemacs has been drafted to perform our indentation services

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