ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/DarkSide/simulation.F90
Revision: 2220
Committed: Thu May 5 14:47:35 2005 UTC (19 years, 4 months ago) by chrisfen
File size: 18880 byte(s)
Log Message:
OOPSE setup for TAP water.  It's not parametrized, but OOPSE will now let me run it...

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