ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/simulation_module.F90
Revision: 329
Committed: Wed Mar 12 22:27:59 2003 UTC (21 years, 5 months ago) by gezelter
File size: 12968 byte(s)
Log Message:
Stick a fork in it.  It's rare.

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     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 mmeineke 270
26 gezelter 325 #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 gezelter 329 real( kind = dp ), allocatable, dimension(:,:), public :: rf_Row
44     real( kind = dp ), allocatable, dimension(:,:), public :: rf_Col
45    
46 gezelter 325 integer, allocatable, dimension(:), public :: atid_Row
47     integer, allocatable, dimension(:), public :: atid_Col
48     #else
49 gezelter 329 real( kind = dp ), allocatable, dimension(:,:), public :: rf
50 gezelter 325 integer, allocatable, dimension(:), public :: atid
51     #endif
52 gezelter 329
53     real(kind = dp), dimension(9), public :: tau_Temp = 0.0_dp
54     real(kind = dp), public :: virial_Temp = 0.0_dp
55 gezelter 325
56     public :: SimulationSetup
57 mmeineke 270 public :: getBox
58     public :: getRcut
59     public :: getRlist
60 gezelter 312 public :: getRrf
61     public :: getRt
62 gezelter 329 public :: getDielect
63 mmeineke 270 public :: getNlocal
64 gezelter 325 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 mmeineke 285
74 mmeineke 270 interface getBox
75     module procedure getBox_3d
76     module procedure getBox_dim
77     end interface
78 gezelter 325
79 mmeineke 270 contains
80 gezelter 325
81     subroutine SimulationSetup(setThisSim, nComponents, c_idents, &
82     nExcludes_local, excludesLocal, nExcludes_global, excludesGlobal, &
83     status)
84    
85 chuckv 290 type (simtype) :: setThisSim
86 gezelter 325 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 chuckv 290
96 gezelter 325 #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 gezelter 312 ! copy C struct into fortran type
107 chuckv 290 thisSim = setThisSim
108 gezelter 325 natoms = nComponents
109     rcut2 = thisSim%rcut * thisSim%rcut
110     rcut6 = rcut2 * rcut2 * rcut2
111     rlist2 = thisSim%rlist * thisSim%rlist
112 mmeineke 270
113 gezelter 325 #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 mmeineke 270
130 gezelter 325 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 mmeineke 270 subroutine change_box_size(new_box_size)
196     real(kind=dp), dimension(3) :: new_box_size
197     thisSim%box = new_box_size
198     end subroutine change_box_size
199    
200     function getBox_3d() result(thisBox)
201     real( kind = dp ), dimension(3) :: thisBox
202     thisBox = thisSim%box
203     end function getBox_3d
204    
205     function getBox_dim(dim) result(thisBox)
206     integer, intent(in) :: dim
207     real( kind = dp ) :: thisBox
208 gezelter 312
209 mmeineke 270 thisBox = thisSim%box(dim)
210     end function getBox_dim
211 gezelter 312
212 gezelter 325 subroutine getRcut(thisrcut,rc2,rc6,status)
213 mmeineke 270 real( kind = dp ), intent(out) :: thisrcut
214 gezelter 325 real( kind = dp ), intent(out), optional :: rc2
215     real( kind = dp ), intent(out), optional :: rc6
216 mmeineke 270 integer, optional :: status
217    
218     if (present(status)) status = 0
219 gezelter 325
220     if (.not.simulation_setup_complete ) then
221 mmeineke 270 if (present(status)) status = -1
222     return
223     end if
224    
225     thisrcut = thisSim%rcut
226 gezelter 325 if(present(rc2)) rc2 = rcut2
227     if(present(rc6)) rc6 = rcut6
228 mmeineke 270 end subroutine getRcut
229    
230 gezelter 325 subroutine getRlist(thisrlist,rl2,status)
231 mmeineke 270 real( kind = dp ), intent(out) :: thisrlist
232 gezelter 325 real( kind = dp ), intent(out), optional :: rl2
233 mmeineke 270
234     integer, optional :: status
235    
236     if (present(status)) status = 0
237    
238 gezelter 325 if (.not.simulation_setup_complete ) then
239 mmeineke 270 if (present(status)) status = -1
240     return
241     end if
242    
243     thisrlist = thisSim%rlist
244 gezelter 325 if(present(rl2)) rl2 = rlist2
245 gezelter 312 end subroutine getRlist
246 mmeineke 270
247 gezelter 312 function getRrf() result(rrf)
248     real( kind = dp ) :: rrf
249     rrf = thisSim%rrf
250     end function getRrf
251 mmeineke 270
252 gezelter 312 function getRt() result(rt)
253     real( kind = dp ) :: rt
254     rt = thisSim%rt
255     end function getRt
256 gezelter 329
257     function getDielect() result(dielect)
258     real( kind = dp ) :: dielect
259     dielect = thisSim%dielect
260     end function getDielect
261 gezelter 312
262     pure function getNlocal() result(nlocal)
263 mmeineke 270 integer :: nlocal
264 gezelter 325 nlocal = natoms
265 mmeineke 270 end function getNlocal
266 gezelter 325
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 gezelter 329 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 gezelter 325 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 gezelter 312
318 gezelter 325 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 gezelter 329 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 gezelter 325 #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 gezelter 329 allocate(rf(ndim,nlocal),stat=alloc_stat)
462     if (alloc_stat /= 0 ) then
463     thisStat = 0
464     return
465     endif
466    
467 gezelter 325 #endif
468    
469     end subroutine setupGlobals
470 gezelter 312
471 gezelter 325 subroutine freeGlobals()
472    
473     !We free in the opposite order in which we allocate in.
474     #ifdef IS_MPI
475    
476 gezelter 329 if (allocated(rf_Col)) deallocate(rf_Col)
477     if (allocated(rf_Row)) deallocate(rf_Row)
478 gezelter 325 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 gezelter 329 if (allocated(rf)) deallocate(rf)
499 gezelter 325 if (allocated(atid)) deallocate(atid)
500    
501     #endif
502    
503     end subroutine freeGlobals
504 gezelter 309
505 mmeineke 270 end module simulation