ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/simulation_module.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/simulation_module.F90 (file contents):
Revision 480 by chuckv, Tue Apr 8 17:16:22 2003 UTC vs.
Revision 1183 by gezelter, Fri May 21 15:58:48 2004 UTC

# Line 6 | Line 6 | module simulation
6    use force_globals
7    use vector_class
8    use atype_module
9 +  use switcheroo
10   #ifdef IS_MPI
11    use mpiSimulation
12   #endif
# Line 15 | Line 16 | module simulation
16  
17   #define __FORTRAN90
18   #include "fSimulation.h"
19 + #include "fSwitchingFunction.h"
20  
21 <  type (simtype), public :: thisSim
21 >  type (simtype), public, save :: thisSim
22  
23    logical, save :: simulation_setup_complete = .false.
24  
25 <  integer, public, save :: natoms
25 >  integer, public, save :: nLocal, nGlobal
26 >  integer, public, save :: nGroup, 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
30 >  integer, allocatable, dimension(:),   public :: excludesGlobal
31 >  integer, allocatable, dimension(:),   public :: molMembershipList
32 >  integer, allocatable, dimension(:),   public :: groupList
33 >  integer, allocatable, dimension(:),   public :: groupStart
34 >  integer, allocatable, dimension(:),   public :: nSkipsForAtom
35 >  integer, allocatable, dimension(:,:), public :: skipsForAtom
36 >  real(kind=dp), allocatable, dimension(:), public :: mfact
37  
38 <  real(kind=dp), save :: rcut2 = 0.0_DP
39 <  real(kind=dp), save :: rcut6 = 0.0_DP
31 <  real(kind=dp), save :: rlist2 = 0.0_DP
32 <  real(kind=dp), public, dimension(3), save :: box
33 <
38 >  real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
39 >  logical, public, save :: boxIsOrthorhombic
40    
41    public :: SimulationSetup
42    public :: getNlocal
43    public :: setBox
38 –  public :: setBox_3d
39 –  public :: getBox
40 –  public :: setRcut
41 –  public :: getRcut
42 –  public :: getRlist
43 –  public :: getRrf
44 –  public :: getRt
44    public :: getDielect
45    public :: SimUsesPBC
46    public :: SimUsesLJ
47 +  public :: SimUsesCharges
48    public :: SimUsesDipoles
49    public :: SimUsesSticky
50    public :: SimUsesRF
# Line 53 | Line 53 | module simulation
53    public :: SimRequiresPrepairCalc
54    public :: SimRequiresPostpairCalc
55    public :: SimUsesDirectionalAtoms
56 –
57 –  interface getBox
58 –     module procedure getBox_3d
59 –     module procedure getBox_1d
60 –  end interface
56    
62 –  interface setBox
63 –     module procedure setBox_3d
64 –     module procedure setBox_1d
65 –  end interface
66 –  
57   contains
58    
59 <  subroutine SimulationSetup(setThisSim, nComponents, c_idents, &
59 >  subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
60         CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
61 +       CmolMembership, Cmfact, CnGroup, CgroupList, CgroupStart, &
62         status)    
63  
64      type (simtype) :: setThisSim
65 <    integer, intent(inout) :: nComponents
66 <    integer, dimension(nComponents),intent(inout) :: c_idents
65 >    integer, intent(inout) :: CnGlobal, CnLocal
66 >    integer, dimension(CnLocal),intent(inout) :: c_idents
67  
68      integer :: CnLocalExcludes
69      integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
70      integer :: CnGlobalExcludes
71      integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
72 +    integer, dimension(CnGlobal),intent(in) :: CmolMembership
73      !!  Result status, success = 0, status = -1
74      integer, intent(out) :: status
75 <    integer :: i, me, thisStat, alloc_stat, myNode
75 >    integer :: i, j, me, thisStat, alloc_stat, myNode, id1, id2
76 >    integer :: gStart, gEnd, groupSize, biggestGroupSize, ia
77 >
78 >    !! mass factors used for molecular cutoffs
79 >    real ( kind = dp ), dimension(CnLocal) :: Cmfact
80 >    integer, intent(in):: CnGroup
81 >    integer, dimension(CnLocal),intent(in) :: CgroupList
82 >    integer, dimension(CnGroup),intent(in) :: CgroupStart
83 >    integer :: maxSkipsForAtom
84 >
85   #ifdef IS_MPI
86      integer, allocatable, dimension(:) :: c_idents_Row
87      integer, allocatable, dimension(:) :: c_idents_Col
# Line 92 | Line 93 | contains
93      status = 0
94  
95      ! copy C struct into fortran type
96 +
97 +    nLocal = CnLocal
98 +    nGlobal = CnGlobal
99 +    nGroup = CnGroup
100 +
101      thisSim = setThisSim
96 –    natoms = nComponents
97 –    rcut2 = thisSim%rcut * thisSim%rcut
98 –    rcut6 = rcut2 * rcut2 * rcut2
99 –    rlist2 = thisSim%rlist * thisSim%rlist
100 –    box = thisSim%box
102  
103      nExcludes_Global = CnGlobalExcludes
104      nExcludes_Local = CnLocalExcludes
105  
106 <    call InitializeForceGlobals(natoms, thisStat)
106 >    call InitializeForceGlobals(nLocal, thisStat)
107      if (thisStat /= 0) then
108         write(default_error,*) "SimSetup: InitializeForceGlobals error"
109         status = -1
# Line 160 | Line 161 | contains
161         deallocate(c_idents_Row)
162      endif
163      
164 < #else
165 <    do i = 1, nComponents
164 > #endif
165 >
166 > ! We build the local atid's for both mpi and nonmpi
167 >    do i = 1, nLocal
168        
169         me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
170         atid(i) = me
171    
172      enddo
170 – #endif
173  
172 –
173 –
174      do i = 1, nExcludes_Local
175         excludesLocal(1,i) = CexcludesLocal(1,i)
176         excludesLocal(2,i) = CexcludesLocal(2,i)
177      enddo
178 +
179 +    maxSkipsForAtom = 0
180 +    do j = 1, nLocal
181 +       nSkipsForAtom(j) = 0
182 + #ifdef IS_MPI
183 +       id1 = tagRow(j)
184 + #else
185 +       id1 = j
186 + #endif
187 +       do i = 1, nExcludes_Local
188 +          if (excludesLocal(1,i) .eq. id1 ) then
189 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
190 +
191 +             if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
192 +                maxSkipsForAtom = nSkipsForAtom(j)
193 +             endif
194 +          endif
195 +          if (excludesLocal(2,i) .eq. id1 ) then
196 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
197 +
198 +             if (nSkipsForAtom(j) .gt. maxSkipsForAtom) then
199 +                maxSkipsForAtom = nSkipsForAtom(j)
200 +             endif
201 +          endif
202 +       end do
203 +    enddo
204 +
205 +    allocate(skipsForAtom(nLocal, maxSkipsForAtom), stat=alloc_stat)
206 +    if (alloc_stat /= 0 ) then
207 +       write(*,*) 'Could not allocate skipsForAtom array'
208 +       return
209 +    endif
210 +
211 +    do j = 1, nLocal
212 +       nSkipsForAtom(j) = 0
213 + #ifdef IS_MPI
214 +       id1 = tagRow(j)
215 + #else
216 +       id1 = j
217 + #endif
218 +       do i = 1, nExcludes_Local
219 +          if (excludesLocal(1,i) .eq. id1 ) then
220 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
221 + #ifdef IS_MPI
222 +             id2 = tagColumn(excludesLocal(2,i))
223 + #else
224 +             id2 = excludesLocal(2,i)
225 + #endif
226 +             skipsForAtom(j, nSkipsForAtom(j)) = id2
227 +          endif
228 +          if (excludesLocal(2, i) .eq. id2 ) then
229 +             nSkipsForAtom(j) = nSkipsForAtom(j) + 1
230 + #ifdef IS_MPI
231 +             id2 = tagColumn(excludesLocal(1,i))
232 + #else
233 +             id2 = excludesLocal(1,i)
234 + #endif
235 +             skipsForAtom(j, nSkipsForAtom(j)) = id2
236 +          endif
237 +       end do
238 +    enddo
239      
240      do i = 1, nExcludes_Global
241         excludesGlobal(i) = CexcludesGlobal(i)
242 +    enddo
243 +
244 +    do i = 1, nGlobal
245 +       molMemberShipList(i) = CmolMembership(i)
246 +    enddo
247 +    
248 +    biggestGroupSize = 0
249 +    do i = 1, nGroup
250 +       groupStart(i) = CgroupStart(i)
251 +       groupSize = 0
252 +       gStart = groupStart(i)      
253 +       if (i .eq. nGroup) then
254 +          gEnd = nLocal
255 +       else
256 +          gEnd = CgroupStart(i+1) - 1
257 +       endif
258 +       do ia = gStart, gEnd
259 +          groupList(ia) = CgroupList(ia)
260 +          groupSize = groupSize + 1
261 +       enddo
262 +       if (groupSize .gt. biggestGroupSize) biggestGroupSize = groupSize      
263      enddo
264 +    groupStart(nGroup+1) = nLocal+1
265  
266 +    do i = 1, nLocal
267 +       mfact(i) = Cmfact(i)
268 +    enddo    
269 +    
270      if (status == 0) simulation_setup_complete = .true.
271      
272    end subroutine SimulationSetup
273    
274 <  subroutine setBox_3d(new_box_size)
275 <    real(kind=dp), dimension(3) :: new_box_size
274 >  subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
275 >    real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
276 >    integer :: cBoxIsOrthorhombic
277      integer :: smallest, status, i
190 –
191 –    thisSim%box = new_box_size
192 –    box = thisSim%box
193 –
194 –    return    
195 –  end subroutine setBox_3d
196 –
197 –  subroutine setBox_1d(dim, new_box_size)
198 –    integer :: dim, status
199 –    real(kind=dp) :: new_box_size
200 –    thisSim%box(dim) = new_box_size
201 –    box(dim) = thisSim%box(dim)
202 –  end subroutine setBox_1d
203 –
204 –  subroutine setRcut(new_rcut, status)
205 –    real(kind = dp) :: new_rcut
206 –    integer :: myStatus, status
207 –    thisSim%rcut = new_rcut
208 –    rcut2 = thisSim%rcut * thisSim%rcut
209 –    rcut6 = rcut2 * rcut2 * rcut2
210 –    status = 0
211 –    return
212 –  end subroutine setRcut
213 –
214 –  function getBox_3d() result(thisBox)
215 –    real( kind = dp ), dimension(3) :: thisBox
216 –    thisBox = thisSim%box
217 –  end function getBox_3d
218 –  
219 –  function getBox_1d(dim) result(thisBox)
220 –    integer, intent(in) :: dim
221 –    real( kind = dp ) :: thisBox
278      
279 <    thisBox = thisSim%box(dim)
280 <  end function getBox_1d
281 <    
282 <  subroutine getRcut(thisrcut,rc2,rc6,status)
283 <    real( kind = dp ), intent(out) :: thisrcut
284 <    real( kind = dp ), intent(out), optional :: rc2
285 <    real( kind = dp ), intent(out), optional :: rc6
230 <    integer, optional :: status
231 <
232 <    if (present(status)) status = 0
279 >    Hmat = cHmat
280 >    HmatInv = cHmatInv
281 >    if (cBoxIsOrthorhombic .eq. 0 ) then
282 >       boxIsOrthorhombic = .false.
283 >    else
284 >       boxIsOrthorhombic = .true.
285 >    endif
286      
287 <    if (.not.simulation_setup_complete ) then
288 <       if (present(status)) status = -1
236 <       return
237 <    end if
238 <    
239 <    thisrcut = thisSim%rcut
240 <    if(present(rc2)) rc2 = rcut2
241 <    if(present(rc6)) rc6 = rcut6
242 <  end subroutine getRcut
243 <  
244 <  subroutine getRlist(thisrlist,rl2,status)
245 <    real( kind = dp ), intent(out) :: thisrlist
246 <    real( kind = dp ), intent(out), optional :: rl2
287 >    return    
288 >  end subroutine setBox
289  
248 –    integer, optional :: status
249 –
250 –    if (present(status)) status = 0
251 –
252 –    if (.not.simulation_setup_complete ) then
253 –       if (present(status)) status = -1
254 –       return
255 –    end if
256 –    
257 –    thisrlist = thisSim%rlist
258 –    if(present(rl2)) rl2 = rlist2
259 –  end subroutine getRlist
260 –
261 –  function getRrf() result(rrf)
262 –    real( kind = dp ) :: rrf
263 –    rrf = thisSim%rrf
264 –    write(*,*) 'getRrf = ', rrf, thisSim%rrf
265 –  end function getRrf
266 –  
267 –  function getRt() result(rt)
268 –    real( kind = dp ) :: rt
269 –    rt = thisSim%rt
270 –  end function getRt
271 –
290    function getDielect() result(dielect)
291      real( kind = dp ) :: dielect
292      dielect = thisSim%dielect
# Line 289 | Line 307 | contains
307      doesit = thisSim%SIM_uses_sticky
308    end function SimUsesSticky
309  
310 +  function SimUsesCharges() result(doesit)
311 +    logical :: doesit
312 +    doesit = thisSim%SIM_uses_charges
313 +  end function SimUsesCharges
314 +
315    function SimUsesDipoles() result(doesit)
316      logical :: doesit
317      doesit = thisSim%SIM_uses_dipoles
# Line 332 | Line 355 | contains
355      thisStat = 0
356      
357      call FreeSimGlobals()    
358 +
359 +    allocate(nSkipsForAtom(nLocal), stat=alloc_stat)
360 +    if (alloc_stat /= 0 ) then
361 +       thisStat = -1
362 +       return
363 +    endif
364      
365      allocate(excludesLocal(2,nExcludes_Local), stat=alloc_stat)
366      if (alloc_stat /= 0 ) then
# Line 344 | Line 373 | contains
373         thisStat = -1
374         return
375      endif
376 +
377 +    allocate(molMembershipList(nGlobal), stat=alloc_stat)
378 +    if (alloc_stat /= 0 ) then
379 +       thisStat = -1
380 +       return
381 +    endif
382 +
383 +    allocate(groupStart(nGroup+1), stat=alloc_stat)
384 +    if (alloc_stat /= 0 ) then
385 +       thisStat = -1
386 +       return
387 +    endif
388 +
389 +    allocate(groupList(nLocal), stat=alloc_stat)
390 +    if (alloc_stat /= 0 ) then
391 +       thisStat = -1
392 +       return
393 +    endif
394 +
395 +    allocate(mfact(nLocal), stat=alloc_stat)
396 +    if (alloc_stat /= 0 ) then
397 +       thisStat = -1
398 +       return
399 +    endif
400      
401    end subroutine InitializeSimGlobals
402    
403    subroutine FreeSimGlobals()
404      
405      !We free in the opposite order in which we allocate in.
406 <    
406 >
407 >    if (allocated(skipsForAtom)) deallocate(skipsForAtom)
408 >    if (allocated(mfact)) deallocate(mfact)
409 >    if (allocated(groupList)) deallocate(groupList)    
410 >    if (allocated(groupStart)) deallocate(groupStart)    
411 >    if (allocated(molMembershipList)) deallocate(molMembershipList)    
412      if (allocated(excludesGlobal)) deallocate(excludesGlobal)
413      if (allocated(excludesLocal)) deallocate(excludesLocal)
414 <
414 >    
415    end subroutine FreeSimGlobals
416 <
417 <  pure function getNlocal() result(nlocal)
418 <    integer :: nlocal
419 <    nlocal = natoms
416 >  
417 >  pure function getNlocal() result(n)
418 >    integer :: n
419 >    n = nLocal
420    end function getNlocal
363 –
421    
422 +  
423   end module simulation

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines