ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/simulation_module.F90
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/simulation_module.F90 (file contents):
Revision 312 by gezelter, Tue Mar 11 17:46:18 2003 UTC vs.
Revision 325 by gezelter, Wed Mar 12 19:10:54 2003 UTC

# Line 1 | Line 1
1   !! Fortran interface to C entry plug.
2  
3   module simulation
4 <  use definitions, ONLY :dp
4 >  use definitions
5 >  use neighborLists
6 >  use atype_typedefs
7   #ifdef IS_MPI
8    use mpiSimulation
9   #endif
# Line 14 | Line 16 | module simulation
16  
17    type (simtype), public :: thisSim
18  
19 <  logical :: setSim = .false.
19 >  logical, save :: simulation_setup_complete = .false.
20  
21 <  integer,public :: natoms
21 >  integer :: natoms
22 >  real(kind=dp), save :: rcut2 = 0.0_DP
23 >  real(kind=dp), save :: rcut6 = 0.0_DP
24 >  real(kind=dp), save :: rlist2 = 0.0_DP
25  
26 + #ifdef IS_MPI
27 +  real( kind = dp ), allocatable, dimension(:,:), public :: q_Row
28 +  real( kind = dp ), allocatable, dimension(:,:), public :: q_Col
29 +  real( kind = dp ), allocatable, dimension(:,:), public :: u_l_Row
30 +  real( kind = dp ), allocatable, dimension(:,:), public :: u_l_Col
31 +  real( kind = dp ), allocatable, dimension(:,:), public :: A_Row
32 +  real( kind = dp ), allocatable, dimension(:,:), public :: A_Col
33 +  
34 +  real( kind = dp ), allocatable, dimension(:), public :: pot_Row
35 +  real( kind = dp ), allocatable, dimension(:), public :: pot_Col
36 +  real( kind = dp ), allocatable, dimension(:), public :: pot_Temp
37 +  real( kind = dp ), allocatable, dimension(:,:), public :: f_Row
38 +  real( kind = dp ), allocatable, dimension(:,:), public :: f_Col
39 +  real( kind = dp ), allocatable, dimension(:,:), public :: f_Temp
40 +  real( kind = dp ), allocatable, dimension(:,:), public :: t_Row
41 +  real( kind = dp ), allocatable, dimension(:,:), public :: t_Col
42 +  real( kind = dp ), allocatable, dimension(:,:), public :: t_Temp
43 +  integer, allocatable, dimension(:), public :: atid_Row
44 +  integer, allocatable, dimension(:), public :: atid_Col
45 + #else
46 +  integer, allocatable, dimension(:), public :: atid
47 + #endif
48 +  
49 +  public :: SimulationSetup
50    public :: getBox
51    public :: getRcut
52    public :: getRlist
53    public :: getRrf
54    public :: getRt
55    public :: getNlocal
56 <  public :: setSimulation
57 <  public :: isEnsemble
58 <  public :: isPBC
59 <  public :: getStringLen
60 <  public :: returnMixingRules
61 <  public :: doStress
56 >  public :: SimUsesPBC
57 >  public :: SimUsesLJ
58 >  public :: SimUsesDipoles
59 >  public :: SimUsesSticky
60 >  public :: SimUsesRF
61 >  public :: SimUsesGB
62 >  public :: SimUsesEAM
63 >  public :: SimRequiresPrepairCalc
64 >  public :: SimRequiresPostpairCalc
65  
34 – !  public :: setRcut
66    interface getBox
67       module procedure getBox_3d
68       module procedure getBox_dim
69    end interface
70 <
70 >  
71   contains
72 <
73 <  subroutine setSimulation(setThisSim,error)
72 >  
73 >  subroutine SimulationSetup(setThisSim, nComponents, c_idents, &
74 >       nExcludes_local, excludesLocal, nExcludes_global, excludesGlobal, &
75 >       status)
76 >    
77      type (simtype) :: setThisSim
78 <    integer :: error
79 <    integer :: alloc_stat
78 >    integer, intent(inout) :: nComponents
79 >    integer, dimension(nComponents),intent(inout) :: c_idents
80 >    integer :: nExcludes_local
81 >    integer, dimension(nExcludes_local),intent(inout) :: excludesLocal
82 >    integer :: nExcludes_global
83 >    integer, dimension(nExcludes_global),intent(inout) :: excludesGlobal
84 >    !!  Result status, success = 0, status = -1
85 >    integer, intent(out) :: status
86 >    integer :: i, me, thisStat, alloc_stat, myNode
87  
88 <    error = 0
89 <    setSim = .true.
88 > #ifdef IS_MPI
89 >    integer, allocatable, dimension(:) :: c_idents_Row
90 >    integer, allocatable, dimension(:) :: c_idents_Col
91 >    integer :: nrow
92 >    integer :: ncol
93 > #endif
94 >
95 >    simulation_setup_complete = .false.
96 >    status = 0
97 >
98      ! copy C struct into fortran type
99      thisSim = setThisSim
100 <  end subroutine setSimulation
100 >    natoms = nComponents
101 >    rcut2 = thisSim%rcut * thisSim%rcut
102 >    rcut6 = rcut2 * rcut2 * rcut2
103 >    rlist2 = thisSim%rlist * thisSim%rlist
104  
105 <  function getNparticles() result(nparticles)
106 <    integer :: nparticles
107 <    nparticles = thisSim%nLRparticles
108 <  end function getNparticles
105 > #ifdef IS_MPI
106 >    ! We can only set up forces if mpiSimulation has been setup.
107 >    if (.not. isMPISimSet()) then
108 >       write(default_error,*) "MPI is not set"
109 >       status = -1
110 >       return
111 >    endif
112 >    nrow = getNrow(plan_row)
113 >    ncol = getNcol(plan_col)
114 >    mynode = getMyNode()
115 >    
116 >    allocate(c_idents_Row(nrow),stat=alloc_stat)
117 >    if (alloc_stat /= 0 ) then
118 >       status = -1
119 >       return
120 >    endif
121  
122 +    allocate(c_idents_Col(ncol),stat=alloc_stat)
123 +    if (alloc_stat /= 0 ) then
124 +       status = -1
125 +       return
126 +    endif
127 +
128 +    call gather(c_idents, c_idents_Row, plan_row)
129 +    call gather(c_idents, c_idents_Col, plan_col)
130 +
131 +    allocate(atid_Row(nrow),stat=alloc_stat)
132 +    if (alloc_stat /= 0 ) then
133 +       status = -1
134 +       return
135 +    endif
136 +
137 +    allocate(atid_Col(ncol),stat=alloc_stat)
138 +    if (alloc_stat /= 0 ) then
139 +       status = -1
140 +       return
141 +    endif
142 +
143 +    do i = 1, nrow
144 +       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
145 +       atid_Row(i) = me
146 +    enddo
147 +
148 +    do i = 1, ncol
149 +       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
150 +       atid_Col(i) = me
151 +    enddo
152 +
153 +    !! free temporary ident arrays
154 +    if (allocated(c_idents_Col)) then
155 +       deallocate(c_idents_Col)
156 +    end if
157 +    if (allocated(c_idents_Row)) then
158 +       deallocate(c_idents_Row)
159 +    endif
160 +    
161 + #else
162 +    do i = 1, nComponents
163 +      
164 +       me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
165 +       atid(i) = me
166 +      
167 +    enddo
168 + #endif
169 +
170 +    call setupGlobals(thisStat)
171 +    if (thisStat /= 0) then
172 +       status = -1
173 +       return
174 +    endif
175 +
176 +    !! Create neighbor lists
177 +    call expandNeighborList(nComponents, thisStat)
178 +    if (thisStat /= 0) then
179 +       status = -1
180 +       return
181 +    endif
182 +    
183 +    if (status == 0) simulation_setup_complete = .true.
184 +      
185 +  end subroutine SimulationSetup
186 +  
187    subroutine change_box_size(new_box_size)
188      real(kind=dp), dimension(3) :: new_box_size
189      thisSim%box = new_box_size
# Line 72 | Line 201 | contains
201      thisBox = thisSim%box(dim)
202    end function getBox_dim
203      
204 <  subroutine getRcut(thisrcut,rcut2,rcut6,status)
204 >  subroutine getRcut(thisrcut,rc2,rc6,status)
205      real( kind = dp ), intent(out) :: thisrcut
206 <    real( kind = dp ), intent(out), optional :: rcut2
207 <    real( kind = dp ), intent(out), optional :: rcut6
206 >    real( kind = dp ), intent(out), optional :: rc2
207 >    real( kind = dp ), intent(out), optional :: rc6
208      integer, optional :: status
209  
210      if (present(status)) status = 0
211 <
212 <    if (.not.setSim ) then
211 >    
212 >    if (.not.simulation_setup_complete ) then
213         if (present(status)) status = -1
214         return
215      end if
216      
217      thisrcut = thisSim%rcut
218 <    if(present(rcut2)) rcut2 = thisSim%rcutsq
219 <    if(present(rcut6)) rcut6 = thisSim%rcut6
218 >    if(present(rc2)) rc2 = rcut2
219 >    if(present(rc6)) rc6 = rcut6
220    end subroutine getRcut
221    
222 <  subroutine getRlist(thisrlist,rlist2,status)
222 >  subroutine getRlist(thisrlist,rl2,status)
223      real( kind = dp ), intent(out) :: thisrlist
224 <    real( kind = dp ), intent(out), optional :: rlist2
224 >    real( kind = dp ), intent(out), optional :: rl2
225  
226      integer, optional :: status
227  
228      if (present(status)) status = 0
229  
230 <    if (.not.setSim ) then
230 >    if (.not.simulation_setup_complete ) then
231         if (present(status)) status = -1
232         return
233      end if
234      
235      thisrlist = thisSim%rlist
236 <    if(present(rlist2)) rlist2 = thisSim%rlistsq
236 >    if(present(rl2)) rl2 = rlist2
237    end subroutine getRlist
238  
239    function getRrf() result(rrf)
# Line 119 | Line 248 | contains
248    
249    pure function getNlocal() result(nlocal)
250      integer :: nlocal
251 <    nlocal = thisSim%nLRparticles
251 >    nlocal = natoms
252    end function getNlocal
253 <  
254 <  function doStress() result(do_stress)
255 <    logical :: do_stress
256 <    do_stress = thisSim%do_stress
257 <  end function doStress
258 <  
259 <  function isEnsemble(this_ensemble) result(is_this_ensemble)
260 <    character(len = *) :: this_ensemble
261 <    logical :: is_this_ensemble
262 <    is_this_ensemble = .false.
263 <    if (this_ensemble == thisSim%ensemble) is_this_ensemble = .true.
264 <  end function isEnsemble
265 <  
266 <  function returnEnsemble() result(thisEnsemble)
267 <    character (len = len(thisSim%ensemble)) :: thisEnsemble
268 <    thisEnsemble = thisSim%ensemble
269 <  end function returnEnsemble
270 <  
271 <  function returnMixingRules() result(thisMixingRule)
272 <    character (len = len(thisSim%ensemble)) :: thisMixingRule
273 <    thisMixingRule = thisSim%MixingRule
274 <  end function returnMixingRules
253 >    
254 >  function SimUsesPBC() result(doesit)
255 >    logical :: doesit
256 >    doesit = thisSim%SIM_uses_PBC
257 >  end function SimUsesPBC
258 >
259 >  function SimUsesLJ() result(doesit)
260 >    logical :: doesit
261 >    doesit = thisSim%SIM_uses_LJ
262 >  end function SimUsesLJ
263 >
264 >  function SimUsesSticky() result(doesit)
265 >    logical :: doesit
266 >    doesit = thisSim%SIM_uses_sticky
267 >  end function SimUsesSticky
268 >
269 >  function SimUsesDipoles() result(doesit)
270 >    logical :: doesit
271 >    doesit = thisSim%SIM_uses_dipoles
272 >  end function SimUsesDipoles
273 >
274 >  function SimUsesRF() result(doesit)
275 >    logical :: doesit
276 >    doesit = thisSim%SIM_uses_RF
277 >  end function SimUsesRF
278 >
279 >  function SimUsesGB() result(doesit)
280 >    logical :: doesit
281 >    doesit = thisSim%SIM_uses_GB
282 >  end function SimUsesGB
283 >
284 >  function SimUsesEAM() result(doesit)
285 >    logical :: doesit
286 >    doesit = thisSim%SIM_uses_EAM
287 >  end function SimUsesEAM
288 >
289 >  function SimRequiresPrepairCalc() result(doesit)
290 >    logical :: doesit
291 >    doesit = thisSim%SIM_uses_EAM
292 >  end function SimRequiresPrepairCalc
293 >
294 >  function SimRequiresPostpairCalc() result(doesit)
295 >    logical :: doesit
296 >    doesit = thisSim%SIM_uses_RF
297 >  end function SimRequiresPostpairCalc
298    
299 <  function isPBC() result(PBCset)
300 <    logical :: PBCset
301 <    PBCset = .false.
302 <    if (thisSim%use_pbc) PBCset = .true.
303 <  end function isPBC
299 >  subroutine setupGlobals(thisStat)
300 >    integer, intent(out) :: thisStat
301 >    integer :: nrow
302 >    integer :: ncol
303 >    integer :: nlocal
304 >    integer :: ndim = 3
305 >    integer :: alloc_stat
306 >    
307 >    thisStat = 0
308 >    
309 > #ifdef IS_MPI
310 >    nrow = getNrow(plan_row)
311 >    ncol = getNcol(plan_col)
312 > #endif
313 >    nlocal = getNlocal()
314 >    
315 >    call freeGlobals()
316 >    
317 > #ifdef IS_MPI
318 >
319 >    allocate(q_Row(ndim,nrow),stat=alloc_stat)
320 >    if (alloc_stat /= 0 ) then
321 >       thisStat = 0
322 >       return
323 >    endif
324 >        
325 >    allocate(q_Col(ndim,ncol),stat=alloc_stat)
326 >    if (alloc_stat /= 0 ) then
327 >       thisStat = 0
328 >       return
329 >    endif
330 >      
331 >    allocate(u_l_Row(ndim,nrow),stat=alloc_stat)
332 >    if (alloc_stat /= 0 ) then
333 >       thisStat = 0
334 >       return
335 >    endif
336 >    
337 >    allocate(u_l_Col(ndim,ncol),stat=alloc_stat)
338 >    if (alloc_stat /= 0 ) then
339 >       thisStat = 0
340 >       return
341 >    endif
342 >    
343 >    allocate(A_row(9,nrow),stat=alloc_stat)
344 >    if (alloc_stat /= 0 ) then
345 >       thisStat = 0
346 >       return
347 >    endif
348 >        
349 >    allocate(A_Col(9,ncol),stat=alloc_stat)
350 >    if (alloc_stat /= 0 ) then
351 >       thisStat = 0
352 >       return
353 >    endif
354 >    
355 >    allocate(pot_row(nrow),stat=alloc_stat)
356 >    if (alloc_stat /= 0 ) then
357 >       thisStat = 0
358 >       return
359 >    endif
360 >    
361 >    allocate(pot_Col(ncol),stat=alloc_stat)
362 >    if (alloc_stat /= 0 ) then
363 >       thisStat = 0
364 >       return
365 >    endif
366 >
367 >    allocate(pot_Temp(nlocal),stat=alloc_stat)
368 >    if (alloc_stat /= 0 ) then
369 >       thisStat = 0
370 >       return
371 >    endif
372 >    
373 >    allocate(f_Row(ndim,nrow),stat=alloc_stat)
374 >    if (alloc_stat /= 0 ) then
375 >       thisStat = 0
376 >       return
377 >    endif
378 >    
379 >    allocate(f_Col(ndim,ncol),stat=alloc_stat)
380 >    if (alloc_stat /= 0 ) then
381 >       thisStat = 0
382 >       return
383 >    endif
384 >    
385 >    allocate(f_Temp(ndim,nlocal),stat=alloc_stat)
386 >    if (alloc_stat /= 0 ) then
387 >       thisStat = 0
388 >       return
389 >    endif
390 >    
391 >    allocate(t_Row(ndim,nrow),stat=alloc_stat)
392 >    if (alloc_stat /= 0 ) then
393 >       thisStat = 0
394 >       return
395 >    endif
396 >    
397 >    allocate(t_Col(ndim,ncol),stat=alloc_stat)
398 >    if (alloc_stat /= 0 ) then
399 >       thisStat = 0
400 >       return
401 >    endif
402 >
403 >    allocate(t_temp(ndim,nlocal),stat=alloc_stat)
404 >    if (alloc_stat /= 0 ) then
405 >       thisStat = 0
406 >       return
407 >    endif
408 >
409 >    allocate(atid_Row(nrow),stat=alloc_stat)
410 >    if (alloc_stat /= 0 ) then
411 >       thisStat = 0
412 >       return
413 >    endif
414 >
415 >    allocate(atid_Col(ncol),stat=alloc_stat)
416 >    if (alloc_stat /= 0 ) then
417 >       thisStat = 0
418 >       return
419 >    endif
420 >
421 > #else
422 >
423 >    allocate(atid(nlocal),stat=alloc_stat)
424 >    if (alloc_stat /= 0 ) then
425 >       thisStat = 0
426 >       return
427 >    end if
428 >
429 > #endif
430 >
431 >  end subroutine setupGlobals
432    
433 <  pure function getStringLen() result (thislen)
434 <    integer :: thislen    
435 <    thislen = string_len
436 <  end function getStringLen
433 >  subroutine freeGlobals()
434 >    
435 >    !We free in the opposite order in which we allocate in.
436 > #ifdef IS_MPI
437 >
438 >    if (allocated(atid_Col))   deallocate(atid_Col)
439 >    if (allocated(atid_Row))   deallocate(atid_Row)
440 >    if (allocated(t_Temp))     deallocate(t_Temp)
441 >    if (allocated(t_Col))      deallocate(t_Col)
442 >    if (allocated(t_Row))      deallocate(t_Row)
443 >    if (allocated(f_Temp))     deallocate(f_Temp)
444 >    if (allocated(f_Col))      deallocate(f_Col)
445 >    if (allocated(f_Row))      deallocate(f_Row)
446 >    if (allocated(pot_Temp))   deallocate(pot_Temp)
447 >    if (allocated(pot_Col))    deallocate(pot_Col)
448 >    if (allocated(pot_Row))    deallocate(pot_Row)
449 >    if (allocated(A_Col))      deallocate(A_Col)
450 >    if (allocated(A_Row))      deallocate(A_Row)
451 >    if (allocated(u_l_Col))    deallocate(u_l_Col)
452 >    if (allocated(u_l_Row))    deallocate(u_l_Row)
453 >    if (allocated(q_Col))      deallocate(q_Col)
454 >    if (allocated(q_Row))      deallocate(q_Row)
455 >    
456 > #else
457 >    
458 >    if (allocated(atid))       deallocate(atid)
459 >    
460 > #endif
461 >        
462 >  end subroutine freeGlobals
463    
464   end module simulation

Diff Legend

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