ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/UseTheForce/DarkSide/simulation.F90
Revision: 1609
Committed: Wed Oct 20 04:12:01 2004 UTC (19 years, 8 months ago) by gezelter
File size: 16717 byte(s)
Log Message:
fixing some broken fortran stuff

File Contents

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