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

Diff Legend

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