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 378 by mmeineke, Fri Mar 21 17:42:12 2003 UTC vs.
Revision 1150 by gezelter, Fri May 7 21:35:05 2004 UTC

# Line 6 | Line 6 | module simulation
6    use force_globals
7    use vector_class
8    use atype_module
9 <  use lj
9 >  use switcheroo
10   #ifdef IS_MPI
11    use mpiSimulation
12   #endif
# Line 16 | 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
31 +  integer, allocatable, dimension(:), public :: molMembershipList
32 +  integer, allocatable, dimension(:), public :: groupList
33 +  integer, allocatable, dimension(:), public :: groupStart
34 +  real(kind=dp), allocatable, dimension(:), public :: mfact
35  
36 <  real(kind=dp), save :: rcut2 = 0.0_DP
37 <  real(kind=dp), save :: rcut6 = 0.0_DP
32 <  real(kind=dp), save :: rlist2 = 0.0_DP
33 <  real(kind=dp), public, dimension(3), save :: box
34 <
36 >  real(kind=dp), public, dimension(3,3), save :: Hmat, HmatInv
37 >  logical, public, save :: boxIsOrthorhombic
38    
39    public :: SimulationSetup
40    public :: getNlocal
41    public :: setBox
39 –  public :: setBox_3d
40 –  public :: getBox
41 –  public :: setRcut
42 –  public :: getRcut
43 –  public :: getRlist
44 –  public :: getRrf
45 –  public :: getRt
42    public :: getDielect
43    public :: SimUsesPBC
44    public :: SimUsesLJ
45 +  public :: SimUsesCharges
46    public :: SimUsesDipoles
47    public :: SimUsesSticky
48    public :: SimUsesRF
# Line 54 | Line 51 | module simulation
51    public :: SimRequiresPrepairCalc
52    public :: SimRequiresPostpairCalc
53    public :: SimUsesDirectionalAtoms
57 –
58 –  interface getBox
59 –     module procedure getBox_3d
60 –     module procedure getBox_1d
61 –  end interface
54    
63 –  interface setBox
64 –     module procedure setBox_3d
65 –     module procedure setBox_1d
66 –  end interface
67 –  
55   contains
56    
57 <  subroutine SimulationSetup(setThisSim, nComponents, c_idents, &
57 >  subroutine SimulationSetup(setThisSim, CnGlobal, CnLocal, c_idents, &
58         CnLocalExcludes, CexcludesLocal, CnGlobalExcludes, CexcludesGlobal, &
59 +       CmolMembership, Cmfact, CnGroup, CgroupList, CgroupStart, &
60         status)    
61  
62      type (simtype) :: setThisSim
63 <    integer, intent(inout) :: nComponents
64 <    integer, dimension(nComponents),intent(inout) :: c_idents
63 >    integer, intent(inout) :: CnGlobal, CnLocal
64 >    integer, dimension(CnLocal),intent(inout) :: c_idents
65  
66      integer :: CnLocalExcludes
67      integer, dimension(2,CnLocalExcludes), intent(in) :: CexcludesLocal
68      integer :: CnGlobalExcludes
69      integer, dimension(CnGlobalExcludes), intent(in) :: CexcludesGlobal
70 +    integer, dimension(CnGlobal),intent(in) :: CmolMembership
71      !!  Result status, success = 0, status = -1
72      integer, intent(out) :: status
73      integer :: i, me, thisStat, alloc_stat, myNode
74 +    integer :: gStart, gEnd, groupSize, biggestGroupSize, ia
75 +
76 +    !! mass factors used for molecular cutoffs
77 +    real ( kind = dp ), dimension(CnLocal) :: Cmfact
78 +    integer, intent(in):: CnGroup
79 +    integer, dimension(CnLocal),intent(in) :: CgroupList
80 +    integer, dimension(CnGroup),intent(in) :: CgroupStart
81 +
82   #ifdef IS_MPI
83      integer, allocatable, dimension(:) :: c_idents_Row
84      integer, allocatable, dimension(:) :: c_idents_Col
# Line 93 | Line 90 | contains
90      status = 0
91  
92      ! copy C struct into fortran type
93 +
94 +    nLocal = CnLocal
95 +    nGlobal = CnGlobal
96 +    nGroup = CnGroup
97 +
98      thisSim = setThisSim
97 –    natoms = nComponents
98 –    rcut2 = thisSim%rcut * thisSim%rcut
99 –    rcut6 = rcut2 * rcut2 * rcut2
100 –    rlist2 = thisSim%rlist * thisSim%rlist
101 –    box = thisSim%box
99  
100      nExcludes_Global = CnGlobalExcludes
101      nExcludes_Local = CnLocalExcludes
102  
103 <    call InitializeForceGlobals(natoms, thisStat)
103 >    call InitializeForceGlobals(nLocal, thisStat)
104      if (thisStat /= 0) then
105 +       write(default_error,*) "SimSetup: InitializeForceGlobals error"
106         status = -1
107         return
108      endif
109  
110      call InitializeSimGlobals(thisStat)
111      if (thisStat /= 0) then
112 +       write(default_error,*) "SimSetup: InitializeSimGlobals error"
113         status = -1
114         return
115      endif
# Line 159 | Line 158 | contains
158         deallocate(c_idents_Row)
159      endif
160      
161 < #else
162 <    do i = 1, nComponents
161 > #endif
162 >
163 > ! We build the local atid's for both mpi and nonmpi
164 >    do i = 1, nLocal
165        
166         me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
167         atid(i) = me
168    
169      enddo
169 – #endif
170  
171 –    !! Create neighbor lists
172 –    call expandNeighborList(nComponents, thisStat)
173 –    if (thisStat /= 0) then
174 –       status = -1
175 –       return
176 –    endif
177 –
171      do i = 1, nExcludes_Local
172         excludesLocal(1,i) = CexcludesLocal(1,i)
173         excludesLocal(2,i) = CexcludesLocal(2,i)
# Line 183 | Line 176 | contains
176      do i = 1, nExcludes_Global
177         excludesGlobal(i) = CexcludesGlobal(i)
178      enddo
179 +
180 +    do i = 1, nGlobal
181 +       molMemberShipList(i) = CmolMembership(i)
182 +    enddo
183      
184 +    biggestGroupSize = 0
185 +    do i = 1, nGroup
186 +       groupStart(i) = CgroupStart(i)
187 +       groupSize = 0
188 +       gStart = groupStart(i)      
189 +       if (i .eq. nGroup) then
190 +          gEnd = nLocal
191 +       else
192 +          gEnd = CgroupStart(i+1) - 1
193 +       endif
194 +       do ia = gStart, gEnd
195 +          groupList(ia) = CgroupList(ia)
196 +          groupSize = groupSize + 1
197 +       enddo
198 +       if (groupSize .gt. biggestGroupSize) biggestGroupSize = groupSize      
199 +    enddo
200 +    groupStart(nGroup+1) = nLocal+1
201 +
202 +    do i = 1, nLocal
203 +       mfact(i) = Cmfact(i)
204 +    enddo    
205 +    
206      if (status == 0) simulation_setup_complete = .true.
207      
208    end subroutine SimulationSetup
209    
210 <  subroutine setBox_3d(new_box_size)
211 <    real(kind=dp), dimension(3) :: new_box_size
210 >  subroutine setBox(cHmat, cHmatInv, cBoxIsOrthorhombic)
211 >    real(kind=dp), dimension(3,3) :: cHmat, cHmatInv
212 >    integer :: cBoxIsOrthorhombic
213      integer :: smallest, status, i
214 <
215 <    thisSim%box = new_box_size
216 <    box = thisSim%box
217 <
218 <    smallest = 1
219 <    do i = 2, 3
220 <       if (new_box_size(i) .lt. new_box_size(smallest)) smallest = i
221 <    end do
222 <    if (thisSim%rcut .gt. 0.5_dp * new_box_size(smallest)) &
203 <         call setRcut(0.5_dp * new_box_size(smallest), status)
214 >    
215 >    Hmat = cHmat
216 >    HmatInv = cHmatInv
217 >    if (cBoxIsOrthorhombic .eq. 0 ) then
218 >       boxIsOrthorhombic = .false.
219 >    else
220 >       boxIsOrthorhombic = .true.
221 >    endif
222 >    
223      return    
224 <  end subroutine setBox_3d
206 <
207 <  subroutine setBox_1d(dim, new_box_size)
208 <    integer :: dim, status
209 <    real(kind=dp) :: new_box_size
210 <    thisSim%box(dim) = new_box_size
211 <    box(dim) = thisSim%box(dim)
212 <    if (thisSim%rcut .gt. 0.5_dp * new_box_size) &
213 <         call setRcut(0.5_dp * new_box_size, status)
214 <  end subroutine setBox_1d
224 >  end subroutine setBox
225  
216 –  subroutine setRcut(new_rcut, status)
217 –    real(kind = dp) :: new_rcut
218 –    integer :: myStatus, status
219 –    thisSim%rcut = new_rcut
220 –    rcut2 = thisSim%rcut * thisSim%rcut
221 –    rcut6 = rcut2 * rcut2 * rcut2
222 –    myStatus = 0
223 –    call LJ_new_rcut(new_rcut, myStatus)
224 –    if (myStatus .ne. 0) then
225 –       write(default_error, *) 'LJ module refused our rcut!'
226 –       status = -1
227 –       return
228 –    endif
229 –    status = 0
230 –    return
231 –  end subroutine setRcut
232 –
233 –  function getBox_3d() result(thisBox)
234 –    real( kind = dp ), dimension(3) :: thisBox
235 –    thisBox = thisSim%box
236 –  end function getBox_3d
237 –  
238 –  function getBox_1d(dim) result(thisBox)
239 –    integer, intent(in) :: dim
240 –    real( kind = dp ) :: thisBox
241 –    
242 –    thisBox = thisSim%box(dim)
243 –  end function getBox_1d
244 –    
245 –  subroutine getRcut(thisrcut,rc2,rc6,status)
246 –    real( kind = dp ), intent(out) :: thisrcut
247 –    real( kind = dp ), intent(out), optional :: rc2
248 –    real( kind = dp ), intent(out), optional :: rc6
249 –    integer, optional :: status
250 –
251 –    if (present(status)) status = 0
252 –    
253 –    if (.not.simulation_setup_complete ) then
254 –       if (present(status)) status = -1
255 –       return
256 –    end if
257 –    
258 –    thisrcut = thisSim%rcut
259 –    if(present(rc2)) rc2 = rcut2
260 –    if(present(rc6)) rc6 = rcut6
261 –  end subroutine getRcut
262 –  
263 –  subroutine getRlist(thisrlist,rl2,status)
264 –    real( kind = dp ), intent(out) :: thisrlist
265 –    real( kind = dp ), intent(out), optional :: rl2
266 –
267 –    integer, optional :: status
268 –
269 –    if (present(status)) status = 0
270 –
271 –    if (.not.simulation_setup_complete ) then
272 –       if (present(status)) status = -1
273 –       return
274 –    end if
275 –    
276 –    thisrlist = thisSim%rlist
277 –    if(present(rl2)) rl2 = rlist2
278 –  end subroutine getRlist
279 –
280 –  function getRrf() result(rrf)
281 –    real( kind = dp ) :: rrf
282 –    rrf = thisSim%rrf
283 –  end function getRrf
284 –  
285 –  function getRt() result(rt)
286 –    real( kind = dp ) :: rt
287 –    rt = thisSim%rt
288 –  end function getRt
289 –
226    function getDielect() result(dielect)
227      real( kind = dp ) :: dielect
228      dielect = thisSim%dielect
# Line 307 | Line 243 | contains
243      doesit = thisSim%SIM_uses_sticky
244    end function SimUsesSticky
245  
246 +  function SimUsesCharges() result(doesit)
247 +    logical :: doesit
248 +    doesit = thisSim%SIM_uses_charges
249 +  end function SimUsesCharges
250 +
251    function SimUsesDipoles() result(doesit)
252      logical :: doesit
253      doesit = thisSim%SIM_uses_dipoles
# Line 362 | Line 303 | contains
303         thisStat = -1
304         return
305      endif
306 +
307 +    allocate(molMembershipList(nGlobal), stat=alloc_stat)
308 +    if (alloc_stat /= 0 ) then
309 +       thisStat = -1
310 +       return
311 +    endif
312 +
313 +    allocate(groupStart(nGroup+1), stat=alloc_stat)
314 +    if (alloc_stat /= 0 ) then
315 +       thisStat = -1
316 +       return
317 +    endif
318 +
319 +    allocate(groupList(nLocal), stat=alloc_stat)
320 +    if (alloc_stat /= 0 ) then
321 +       thisStat = -1
322 +       return
323 +    endif
324 +
325 +    allocate(mfact(nLocal), stat=alloc_stat)
326 +    if (alloc_stat /= 0 ) then
327 +       thisStat = -1
328 +       return
329 +    endif
330      
331    end subroutine InitializeSimGlobals
332    
333    subroutine FreeSimGlobals()
334      
335      !We free in the opposite order in which we allocate in.
336 <    
336 >
337 >    if (allocated(mfact)) deallocate(mfact)
338 >    if (allocated(groupList)) deallocate(groupList)    
339 >    if (allocated(groupStart)) deallocate(groupStart)    
340 >    if (allocated(molMembershipList)) deallocate(molMembershipList)    
341      if (allocated(excludesGlobal)) deallocate(excludesGlobal)
342      if (allocated(excludesLocal)) deallocate(excludesLocal)
343 <
343 >    
344    end subroutine FreeSimGlobals
345 <
346 <  pure function getNlocal() result(nlocal)
347 <    integer :: nlocal
348 <    nlocal = natoms
345 >  
346 >  pure function getNlocal() result(n)
347 >    integer :: n
348 >    n = nLocal
349    end function getNlocal
381 –
350    
351 +  
352   end module simulation

Diff Legend

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