ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/simulation_module.F90
Revision: 358
Committed: Mon Mar 17 21:07:50 2003 UTC (21 years, 4 months ago) by gezelter
File size: 14086 byte(s)
Log Message:
bug fixes for new fForceField.h

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