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 285 by mmeineke, Wed Feb 26 18:45:57 2003 UTC vs.
Revision 328 by gezelter, Wed Mar 12 20:00:58 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 8 | Line 12 | module simulation
12    PRIVATE
13  
14   #define __FORTRAN90
15 < #include "../headers/fsimulation.h"
15 > #include "fSimulation.h"
16  
17    type (simtype), public :: thisSim
14 !! Tag for MPI calculations  
15  integer, allocatable, dimension(:) :: tag
18  
19 < #ifdef IS_MPI
18 <  integer, allocatable, dimension(:) :: tag_row
19 <  integer, allocatable, dimension(:) :: tag_column
20 < #endif
19 >  logical, save :: simulation_setup_complete = .false.
20  
21 < !! WARNING: use_pbc hardcoded, fixme
22 <   logical :: setSim = .false.
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 < !! array for saving previous positions for neighbor lists.  
27 <  real( kind = dp ), allocatable,dimension(:,:),save :: q0
28 <
29 <
30 <  public :: wrap
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
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  
40 !  public :: setRcut
41
42  interface wrap
43     module procedure wrap_1d
44     module procedure wrap_3d
45  end interface
46
66    interface getBox
67       module procedure getBox_3d
68       module procedure getBox_dim
69    end interface
70 +  
71 + contains
72 +  
73 +  subroutine SimulationSetup(setThisSim, nComponents, c_idents, &
74 +       nExcludes_local, excludesLocal, nExcludes_global, excludesGlobal, &
75 +       status)
76 +    
77 +    type (simtype) :: setThisSim
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 + #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 +    natoms = nComponents
101 +    rcut2 = thisSim%rcut * thisSim%rcut
102 +    rcut6 = rcut2 * rcut2 * rcut2
103 +    rlist2 = thisSim%rlist * thisSim%rlist
104  
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 < contains
148 >    do i = 1, ncol
149 >       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
150 >       atid_Col(i) = me
151 >    enddo
152  
153 <  subroutine setSimulation(nLRParticles,box,rlist,rcut,ensemble,mixingRule,use_pbc)
154 <    integer, intent(in) :: nLRParticles
155 <    real(kind = dp ), intent(in), dimension(3) :: box
156 <    real(kind = dp ), intent(in) :: rlist
157 <    real(kind = dp ), intent(in) :: rcut
158 <    character( len = stringLen), intent(in)  :: ensemble
159 <    character( len = stringLen), intent(in)  :: mixingRule
70 <    logical, intent(in) :: use_pbc
71 <    integer :: alloc_stat
72 <    if( setsim ) return  ! simulation is already initialized
73 <    setSim = .true.
74 <
75 <    thisSim%nLRParticles = nLRParticles
76 <    thisSim%box          = box
77 <    thisSim%rlist        = rlist
78 <    thisSIm%rlistsq      = rlist * rlist
79 <    thisSim%rcut         = rcut
80 <    thisSim%rcutsq       = rcut * rcut
81 <    thisSim%rcut6        = thisSim%rcutsq * thisSim%rcutsq * thisSim%rcutsq
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 <    thisSim%ensemble = ensemble
162 <    thisSim%mixingRule = mixingRule
163 <    thisSim%use_pbc = use_pbc
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 <    if (.not. allocated(q0)) then
171 <       allocate(q0(3,nLRParticles),stat=alloc_stat)
170 >    call setupGlobals(thisStat)
171 >    if (thisStat /= 0) then
172 >       status = -1
173 >       return
174      endif
90  end subroutine setSimulation
175  
176 <  function getNparticles() result(nparticles)
177 <    integer :: nparticles
178 <    nparticles = thisSim%nLRparticles
179 <  end function getNparticles
180 <
181 <
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
100
189      thisSim%box = new_box_size
102
190    end subroutine change_box_size
191  
105
192    function getBox_3d() result(thisBox)
193      real( kind = dp ), dimension(3) :: thisBox
194      thisBox = thisSim%box
# Line 111 | Line 197 | contains
197    function getBox_dim(dim) result(thisBox)
198      integer, intent(in) :: dim
199      real( kind = dp ) :: thisBox
200 <
200 >    
201      thisBox = thisSim%box(dim)
202    end function getBox_dim
203 <  
204 <
119 <  function wrap_1d(r,dim) result(this_wrap)
120 <    
121 <    
122 <    real( kind = DP ) :: r
123 <    real( kind = DP ) :: this_wrap
124 <    integer           :: dim
125 <    
126 <    if (use_pbc) then
127 <       !     this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP)
128 <       this_wrap = r - thisSim%box(dim)*nint(r/thisSim%box(dim))
129 <    else
130 <       this_wrap = r
131 <    endif
132 <    
133 <    return
134 <  end function wrap_1d
135 <
136 <  function wrap_3d(r) result(this_wrap)
137 <    real( kind = dp ), dimension(3), intent(in) :: r
138 <    real( kind = dp ), dimension(3) :: this_wrap
139 <
140 <    
141 <    if (this_sim%use_pbc) then
142 <       !     this_wrap = r - box(dim)*dsign(1.0E0_DP,r)*int(abs(r/box(dim)) + 0.5E0_DP)
143 <       this_wrap = r - thisSim%box*nint(r/thisSim%box)
144 <    else
145 <       this_wrap = r
146 <    endif
147 <  end function wrap_3d
148 <
149 <  
150 <
151 <  subroutine getRcut(thisrcut,rcut2,rcut6,status)
203 >    
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
167 <
218 >    if(present(rc2)) rc2 = rcut2
219 >    if(present(rc6)) rc6 = rcut6
220    end subroutine getRcut
221    
222 <  
171 <  
172 <
173 <  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
188 <
189 <
236 >    if(present(rl2)) rl2 = rlist2
237    end subroutine getRlist
191  
192  
238  
239 < pure function getNlocal() result(nlocal)
239 >  function getRrf() result(rrf)
240 >    real( kind = dp ) :: rrf
241 >    rrf = thisSim%rrf
242 >  end function getRrf
243 >  
244 >  function getRt() result(rt)
245 >    real( kind = dp ) :: rt
246 >    rt = thisSim%rt
247 >  end function getRt
248 >  
249 >  pure function getNlocal() result(nlocal)
250      integer :: nlocal
251 <    nlocal = thisSim%nLRparticles
251 >    nlocal = natoms
252    end function getNlocal
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 isEnsemble(this_ensemble) result(is_this_ensemble)
265 <    character(len = *) :: this_ensemble
266 <    logical :: is_this_enemble
267 <    is_this_ensemble = .false.
204 <    if (this_ensemble == thisSim%ensemble) is_this_ensemble = .true.
205 <  end function isEnsemble
264 >  function SimUsesSticky() result(doesit)
265 >    logical :: doesit
266 >    doesit = thisSim%SIM_uses_sticky
267 >  end function SimUsesSticky
268  
269 <  function returnEnsemble() result(thisEnsemble)
270 <    character (len = len(thisSim%ensemble)) :: thisEnsemble
271 <    thisEnsemble = thisSim%ensemble
272 <  end function returnEnsemble
269 >  function SimUsesDipoles() result(doesit)
270 >    logical :: doesit
271 >    doesit = thisSim%SIM_uses_dipoles
272 >  end function SimUsesDipoles
273  
274 <  function returnMixingRules() result(thisMixingRule)
275 <    character (len = len(thisSim%ensemble)) :: thisMixingRule
276 <    thisMixingRule = thisSim%MixingRule
277 <  end function returnMixingRules
274 >  function SimUsesRF() result(doesit)
275 >    logical :: doesit
276 >    doesit = thisSim%SIM_uses_RF
277 >  end function SimUsesRF
278  
279 <  function isPBC() result(PBCset)
280 <    logical :: PBCset
281 <    PBCset = .false.
282 <    if (thisSim%use_pbc) PBCset = .true.
221 <  end function isPBC
279 >  function SimUsesGB() result(doesit)
280 >    logical :: doesit
281 >    doesit = thisSim%SIM_uses_GB
282 >  end function SimUsesGB
283  
284 <  pure function getStringLen() result (thislen)
285 <    integer :: thislen    
286 <    thislen = string_len
287 <  end function setStringLen
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 +  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 +  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