ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-1.0/libmdtools/simulation_module.F90
Revision: 1334
Committed: Fri Jul 16 18:58:03 2004 UTC (19 years, 11 months ago) by gezelter
File size: 15116 byte(s)
Log Message:
Initial import of OOPSE-1.0 source tree

File Contents

# User Rev Content
1 gezelter 1334 !! 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 "fSimulation.h"
19     #include "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