ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/simulation.F90
Revision: 1948
Committed: Fri Jan 14 20:31:16 2005 UTC (19 years, 6 months ago) by gezelter
File size: 17783 byte(s)
Log Message:
separating modules and C/Fortran interface subroutines

File Contents

# User Rev Content
1 gezelter 1930 !!
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 gezelter 1608 !! 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 gezelter 1633
94     public :: SimUsesDirectionalAtoms
95     public :: SimUsesLennardJones
96     public :: SimUsesElectrostatics
97 gezelter 1608 public :: SimUsesCharges
98     public :: SimUsesDipoles
99     public :: SimUsesSticky
100 gezelter 1633 public :: SimUsesGayBerne
101     public :: SimUsesEAM
102     public :: SimUsesShapes
103     public :: SimUsesFLARB
104 gezelter 1608 public :: SimUsesRF
105     public :: SimRequiresPrepairCalc
106     public :: SimRequiresPostpairCalc
107 gezelter 1633
108 gezelter 1608
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 gezelter 1633 function SimUsesDirectionalAtoms() result(doesit)
490 gezelter 1608 logical :: doesit
491 gezelter 1633 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 gezelter 1608
495 gezelter 1633 function SimUsesLennardJones() result(doesit)
496 gezelter 1608 logical :: doesit
497 gezelter 1633 doesit = thisSim%SIM_uses_LennardJones
498     end function SimUsesLennardJones
499 gezelter 1608
500 gezelter 1633 function SimUsesElectrostatics() result(doesit)
501     logical :: doesit
502     doesit = thisSim%SIM_uses_Electrostatics
503     end function SimUsesElectrostatics
504    
505 gezelter 1608 function SimUsesCharges() result(doesit)
506     logical :: doesit
507 gezelter 1633 doesit = thisSim%SIM_uses_Charges
508 gezelter 1608 end function SimUsesCharges
509    
510     function SimUsesDipoles() result(doesit)
511     logical :: doesit
512 gezelter 1633 doesit = thisSim%SIM_uses_Dipoles
513 gezelter 1608 end function SimUsesDipoles
514    
515 gezelter 1633 function SimUsesSticky() result(doesit)
516 gezelter 1608 logical :: doesit
517 gezelter 1633 doesit = thisSim%SIM_uses_Sticky
518     end function SimUsesSticky
519 gezelter 1608
520 gezelter 1633 function SimUsesGayBerne() result(doesit)
521 gezelter 1608 logical :: doesit
522 gezelter 1633 doesit = thisSim%SIM_uses_GayBerne
523     end function SimUsesGayBerne
524    
525 gezelter 1608 function SimUsesEAM() result(doesit)
526     logical :: doesit
527     doesit = thisSim%SIM_uses_EAM
528     end function SimUsesEAM
529    
530 gezelter 1633 function SimUsesShapes() result(doesit)
531 gezelter 1608 logical :: doesit
532 gezelter 1633 doesit = thisSim%SIM_uses_Shapes
533     end function SimUsesShapes
534 gezelter 1608
535 gezelter 1633 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 gezelter 1608 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