ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/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

# 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 gezelter 2204
88 gezelter 1608 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 chrisfen 2220 public :: SimUsesStickyPower
101 gezelter 1633 public :: SimUsesGayBerne
102     public :: SimUsesEAM
103     public :: SimUsesShapes
104     public :: SimUsesFLARB
105 gezelter 1608 public :: SimUsesRF
106     public :: SimRequiresPrepairCalc
107     public :: SimRequiresPostpairCalc
108 gezelter 1633
109 gezelter 2204
110 gezelter 1608 contains
111 gezelter 2204
112 gezelter 1608 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 gezelter 2204
185 gezelter 1608 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 gezelter 2204
218 gezelter 1608 #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 gezelter 2204
257 gezelter 1608 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 gezelter 2204
296 gezelter 1608 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 gezelter 2204
356 gezelter 1608 #endif
357    
358    
359 gezelter 2204 ! We build the local atid's for both mpi and nonmpi
360 gezelter 1608 do i = 1, nLocal
361 gezelter 2204
362 gezelter 1608 me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
363     atid(i) = me
364 gezelter 2204
365 gezelter 1608 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 gezelter 2204 do j = 1, nLocal
388 gezelter 1608 #endif
389 gezelter 2204 nSkipsForAtom(j) = 0
390 gezelter 1608 #ifdef IS_MPI
391 gezelter 2204 id1 = AtomRowToGlobal(j)
392 gezelter 1608 #else
393 gezelter 2204 id1 = j
394 gezelter 1608 #endif
395 gezelter 2204 do i = 1, nExcludes_Local
396     if (excludesLocal(1,i) .eq. id1 ) then
397     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
398 gezelter 1608
399 gezelter 2204 if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
400     maxSkipsForAtom = nSkipsForAtom(j)
401     endif
402 gezelter 1608 endif
403 gezelter 2204 if (excludesLocal(2,i) .eq. id1 ) then
404     nSkipsForAtom(j) = nSkipsForAtom(j) + 1
405 gezelter 1608
406 gezelter 2204 if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
407     maxSkipsForAtom = nSkipsForAtom(j)
408     endif
409 gezelter 1608 endif
410 gezelter 2204 end do
411     enddo
412 gezelter 1608
413     #ifdef IS_MPI
414 gezelter 2204 allocate(skipsForAtom(nAtomsInRow, maxSkipsForAtom), stat=alloc_stat)
415 gezelter 1608 #else
416 gezelter 2204 allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
417 gezelter 1608 #endif
418 gezelter 2204 if (alloc_stat /= 0 ) then
419     write(*,*) 'Could not allocate skipsForAtom array'
420     return
421     endif
422 gezelter 1608
423     #ifdef IS_MPI
424 gezelter 2204 do j = 1, nAtomsInRow
425 gezelter 1608 #else
426 gezelter 2204 do j = 1, nLocal
427 gezelter 1608 #endif
428 gezelter 2204 nSkipsForAtom(j) = 0
429 gezelter 1608 #ifdef IS_MPI
430 gezelter 2204 id1 = AtomRowToGlobal(j)
431 gezelter 1608 #else
432 gezelter 2204 id1 = j
433 gezelter 1608 #endif
434 gezelter 2204 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 gezelter 1608 endif
476    
477 gezelter 2204 return
478     end subroutine setBox
479 gezelter 1608
480 gezelter 2204 function getDielect() result(dielect)
481     real( kind = dp ) :: dielect
482     dielect = thisSim%dielect
483     end function getDielect
484 gezelter 1608
485 gezelter 2204 function SimUsesPBC() result(doesit)
486     logical :: doesit
487     doesit = thisSim%SIM_uses_PBC
488     end function SimUsesPBC
489 gezelter 1608
490 gezelter 2204 function SimUsesDirectionalAtoms() result(doesit)
491     logical :: doesit
492 chrisfen 2220 doesit = thisSim%SIM_uses_dipoles .or. thisSim%SIM_uses_Sticky .or. &
493     thisSim%SIM_uses_StickyPower .or. &
494 gezelter 2204 thisSim%SIM_uses_GayBerne .or. thisSim%SIM_uses_Shapes
495     end function SimUsesDirectionalAtoms
496 gezelter 1608
497 gezelter 2204 function SimUsesLennardJones() result(doesit)
498     logical :: doesit
499     doesit = thisSim%SIM_uses_LennardJones
500     end function SimUsesLennardJones
501 gezelter 1633
502 gezelter 2204 function SimUsesElectrostatics() result(doesit)
503     logical :: doesit
504     doesit = thisSim%SIM_uses_Electrostatics
505     end function SimUsesElectrostatics
506 gezelter 1608
507 gezelter 2204 function SimUsesCharges() result(doesit)
508     logical :: doesit
509     doesit = thisSim%SIM_uses_Charges
510     end function SimUsesCharges
511 gezelter 1608
512 gezelter 2204 function SimUsesDipoles() result(doesit)
513     logical :: doesit
514     doesit = thisSim%SIM_uses_Dipoles
515     end function SimUsesDipoles
516 gezelter 1608
517 gezelter 2204 function SimUsesSticky() result(doesit)
518     logical :: doesit
519     doesit = thisSim%SIM_uses_Sticky
520     end function SimUsesSticky
521 gezelter 1608
522 chrisfen 2220 function SimUsesStickyPower() result(doesit)
523     logical :: doesit
524     doesit = thisSim%SIM_uses_StickyPower
525     end function SimUsesStickyPower
526    
527 gezelter 2204 function SimUsesGayBerne() result(doesit)
528     logical :: doesit
529     doesit = thisSim%SIM_uses_GayBerne
530     end function SimUsesGayBerne
531 gezelter 1608
532 gezelter 2204 function SimUsesEAM() result(doesit)
533     logical :: doesit
534     doesit = thisSim%SIM_uses_EAM
535     end function SimUsesEAM
536 gezelter 1633
537 gezelter 2204 function SimUsesShapes() result(doesit)
538     logical :: doesit
539     doesit = thisSim%SIM_uses_Shapes
540     end function SimUsesShapes
541 gezelter 1633
542 gezelter 2204 function SimUsesFLARB() result(doesit)
543     logical :: doesit
544     doesit = thisSim%SIM_uses_FLARB
545     end function SimUsesFLARB
546 gezelter 1608
547 gezelter 2204 function SimUsesRF() result(doesit)
548     logical :: doesit
549     doesit = thisSim%SIM_uses_RF
550     end function SimUsesRF
551 gezelter 1608
552 gezelter 2204 function SimRequiresPrepairCalc() result(doesit)
553     logical :: doesit
554     doesit = thisSim%SIM_uses_EAM
555     end function SimRequiresPrepairCalc
556 gezelter 1608
557 gezelter 2204 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