ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/DarkSide/simulation.F90
Revision: 2310
Committed: Mon Sep 19 23:21:46 2005 UTC (19 years ago) by chrisfen
File size: 20633 byte(s)
Log Message:
Fixed bugs in reaction field, now it appears as though it really is working...

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