ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/simulation_module.F90
Revision: 332
Committed: Thu Mar 13 15:28:43 2003 UTC (21 years, 5 months ago) by gezelter
File size: 14067 byte(s)
Log Message:
initialization routines, bug fixes, exclude lists

File Contents

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