ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/forceGlobals.F90
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/forceGlobals.F90 (file contents):
Revision 309 by gezelter, Mon Mar 10 23:19:23 2003 UTC vs.
Revision 318 by chuckv, Wed Mar 12 13:51:04 2003 UTC

# Line 1 | Line 1
1   module forceGlobals
2 <
2 >  
3    use definitions
4    use simulation
5    use atype_typedefs
6 <  use generic_atypes
6 >  use vector_class
7   #ifdef IS_MPI
8    use mpiSimulation
9   #endif
10 +    
11 +  logical, save :: atype_vector_initialized = .false.
12 +  logical, save :: ff_initialized = .false.
13  
14 < !! Number of lj_atypes in lj_atype_list
12 <  integer, save :: n_atypes = 0
13 <
14 < !! Global list of lj atypes in simulation
15 <  type (atype), pointer :: ListHead => null()
16 <  type (atype), pointer :: ListTail => null()
17 <
18 <
19 <
20 <
21 <  logical, save :: firstTime = .True.
22 <
23 < !! Atype identity pointer lists
14 >  
15   #ifdef IS_MPI
16 < !! Row lj_atype pointer list
17 <  type (identPtrList), dimension(:), pointer :: identPtrListRow => null()
18 < !! Column lj_atype pointer list
19 <  type (identPtrList), dimension(:), pointer :: identPtrListColumn => null()
16 >  real( kind = dp ), allocatable, dimension(:,:) :: q_Row
17 >  real( kind = dp ), allocatable, dimension(:,:) :: q_Col
18 >  real( kind = dp ), allocatable, dimension(:,:) :: u_l_Row
19 >  real( kind = dp ), allocatable, dimension(:,:) :: u_l_Col
20 >  real( kind = dp ), allocatable, dimension(:,:) :: A_Row
21 >  real( kind = dp ), allocatable, dimension(:,:) :: A_Col
22 >  
23 >  real( kind = dp ), allocatable, dimension(:) :: pot_Row
24 >  real( kind = dp ), allocatable, dimension(:) :: pot_Col
25 >  real( kind = dp ), allocatable, dimension(:) :: pot_Temp
26 >  real( kind = dp ), allocatable, dimension(:,:) :: f_Row
27 >  real( kind = dp ), allocatable, dimension(:,:) :: f_Col
28 >  real( kind = dp ), allocatable, dimension(:,:) :: f_Temp
29 >  real( kind = dp ), allocatable, dimension(:,:) :: t_Row
30 >  real( kind = dp ), allocatable, dimension(:,:) :: t_Col
31 >  real( kind = dp ), allocatable, dimension(:,:) :: t_Temp
32 >  integer, allocatable, dimension(:) :: atid_Row
33 >  integer, allocatable, dimension(:) :: atid_Col
34   #else
35 <  type(identPtrList ), dimension(:), pointer :: identPtrListGlobal => null()
35 >  integer, allocatable, dimension(:) :: atid
36   #endif
37 <
33 <
34 < !! Logical has lj force field module been initialized?
35 <  logical, save :: isFFinit = .false.
36 <
37 < !! Use periodic boundry conditions
38 <  logical :: wrap = .false.
39 <
40 < !! Potential energy global module variables
41 < #ifdef IS_MPI
42 <    real( kind = dp ), allocatable, dimension(:,:) :: q_Row
43 <    real( kind = dp ), allocatable, dimension(:,:) :: q_Col
44 <    real( kind = dp ), allocatable, dimension(:,:) :: u_l_Row
45 <    real( kind = dp ), allocatable, dimension(:,:) :: u_l_Col
46 <    real( kind = dp ), allocatable, dimension(:,:) :: A_Row
47 <    real( kind = dp ), allocatable, dimension(:,:) :: A_Col
48 <
49 <    real( kind = dp ), allocatable, dimension(:) :: pot_Row
50 <    real( kind = dp ), allocatable, dimension(:) :: pot_Col
51 <    real( kind = dp ), allocatable, dimension(:) :: pot_Temp
52 <    real( kind = dp ), allocatable, dimension(:,:) :: f_Row
53 <    real( kind = dp ), allocatable, dimension(:,:) :: f_Col
54 <    real( kind = dp ), allocatable, dimension(:,:) :: f_Temp
55 <    real( kind = dp ), allocatable, dimension(:,:) :: t_Row
56 <    real( kind = dp ), allocatable, dimension(:,:) :: t_Col
57 <    real( kind = dp ), allocatable, dimension(:,:) :: t_Temp
58 < #endif
59 <
37 >  
38    real(kind = dp) :: pot = 0.0_dp
39    real(kind = dp), dimension(9) :: tau_Temp = 0.0_dp
40    real(kind = dp) :: virial_Temp = 0.0_dp
41 <
64 <  logical :: do_preForce  = .false.
65 <  logical :: do_postForce = .false.
66 <
67 <
68 <
69 < !! Public methods and data
41 >  
42    public :: new_atype
43    public :: init_FF
44  
45 + contains
46    
47 +  subroutine new_atype(c_ident, is_LJ, is_Sticky, is_DP, is_GB, &
48 +       lj_epsilon, lj_sigma, &
49 +       dipole_moment, &
50 +       sticky_w0, sticky_v0, &
51 +       GB_sigma, GB_l2b_ratio, GB_eps, GB_eps_ratio, GB_mu, GB_nu, &
52 +       status)
53 +    
54 +    real( kind = dp ), intent(in) :: lj_epsilon
55 +    real( kind = dp ), intent(in) :: lj_sigma
56 +    real( kind = dp ), intent(in) :: sticky_w0
57 +    real( kind = dp ), intent(in) :: sticky_v0
58 +    real( kind = dp ), intent(in) :: dipole_moment
59 +    real( kind = dp ), intent(in) :: GB_sigma, GB_l2b_ratio, GB_eps
60 +    real( kind = dp ), intent(in) :: GB_eps_ratio, GB_mu, GB_nu
61  
62 <
76 < contains
77 <
78 < !! Adds a new lj_atype to the list.
79 <  subroutine new_atype(ident, epsilon, sigma, &
80 <       is_LJ,is_Sticky,is_DP,is_GB,w0,v0,dipoleMoment,status)
81 <
82 <    real( kind = dp ), intent(in) :: epsilon
83 <    real( kind = dp ), intent(in) :: sigma
84 <    real( kind = dp ), intent(in) :: w0
85 <    real( kind = dp ), intent(in) :: v0
86 <    real( kind = dp ), intent(in) :: dipoleMoment
87 <
88 <    integer, intent(in) :: ident
62 >    integer, intent(in)  :: c_ident
63      integer, intent(out) :: status
64 <    integer, intent(in) :: is_Sticky
65 <    integer, intent(in) :: is_DP
66 <    integer, intent(in) :: is_GB
67 <    integer, intent(in) :: is_LJ
64 >    integer, intent(in)  :: is_Sticky
65 >    integer, intent(in)  :: is_DP
66 >    integer, intent(in)  :: is_GB
67 >    integer, intent(in)  :: is_LJ
68 >    integer :: me
69 >    logical :: l_is_LJ, l_is_DP, l_is_Sticky, l_is_GB
70  
95 –
96 –    type (atype), pointer :: the_new_atype
97 –    integer :: alloc_error
98 –    integer :: atype_counter = 0
99 –    integer :: alloc_size
100 –    integer :: err_stat
71      status = 0
72  
73 +    if (.not.atype_vector_initialized) then
74 +       !! There are 16 properties to worry about for now.  
75 +       !! Fix this if needed for more atomic properties
76 +       atypes = initialize(16)
77 +       if (.not.associated(atypes)) then
78 +          status = -1
79 +          return
80 +       endif
81 +    endif
82  
83 +    me = addElement(atypes)
84  
85 < ! allocate a new atype    
106 <    allocate(the_new_atype,stat=alloc_error)
107 <    if (alloc_error /= 0 ) then
108 <       status = -1
109 <       return
110 <    end if
85 >    call setElementProperty(atypes, me, "c_ident", c_ident)
86  
87 < ! assign our new atype information
87 >    l_is_LJ = (is_LJ .ne. 0)
88 >    l_is_DP = (is_DP .ne. 0)
89 >    l_is_Sticky = (is_Sticky .ne. 0)
90 >    l_is_GB = (is_GB .ne. 0)
91  
92 <    the_new_atype%lj_epsilon     = epsilon
93 <    the_new_atype%lj_sigma       = sigma
94 <    the_new_atype%sticky_w0       = w0
95 <    the_new_atype%sticky_v0       = v0
118 <    the_new_atype%dipole_moment       = dipoleMoment
119 <
120 <    
121 < ! assume that this atype will be successfully added
122 <    the_new_atype%atype_ident = ident
123 <    the_new_atype%atype_fortran_index = atype_count + 1
124 <    
125 <    if ( is_Sticky /= 0 )    the_new_atype%is_Sticky   = .true.
126 <    if ( is_GB /= 0 )        the_new_atype%is_GB       = .true.
127 <    if ( is_LJ /= 0 )        the_new_atype%is_LJ       = .true.
128 <    if ( is_DP /= 0 )        the_new_atype%is_DP       = .true.
92 >    call setElementProperty(atypes, me, "is_LJ", l_is_LJ)
93 >    call setElementProperty(atypes, me, "is_DP", l_is_DP)
94 >    call setElementProperty(atypes, me, "is_Sticky", l_is_Sticky)
95 >    call setElementProperty(atypes, me, "is_GB", l_is_GB)
96  
97 <    call add_atype(the_new_atype,ListHead,ListTail,err_stat)
98 <    if (err_stat /= 0 ) then
99 <       status = -1
133 <       return
97 >    if (l_is_LJ) then
98 >       call setElementProperty(atypes, me, "lj_sigma", lj_sigma)
99 >       call setElementProperty(atypes, me, "lj_epsilon", lj_epsilon)
100      endif
101 <    
102 <    atype_count = atype_count + 1
103 <    
104 <    
101 >    if (l_is_DP) then
102 >       call setElementProperty(atypes, me, "dipole_moment", dipole_moment)
103 >    endif
104 >    if (l_is_Sticky) then
105 >       call setElementProperty(atypes, me, "sticky_w0", sticky_w0)
106 >       call setElementProperty(atypes, me, "sticky_v0", sticky_v0)
107 >    endif
108 >    if (l_is_GB) then
109 >       call setElementProperty(atypes, me, "GB_sigma", GB_sigma)
110 >       call setElementProperty(atypes, me, "GB_l2b_ratio", GB_l2b_ratio)
111 >       call setElementProperty(atypes, me, "GB_eps", GB_eps)
112 >       call setElementProperty(atypes, me, "GB_eps_ratio", GB_eps_ratio)
113 >       call setElementProperty(atypes, me, "GB_mu", GB_mu)
114 >       call setElementProperty(atypes, me, "GB_nu", GB_nu)
115 >    endif
116 >
117    end subroutine new_atype
118  
119 <
120 <  subroutine init_FF(nComponents,ident, nExcludes,excludesLocal,status)
121 < !! Number of components in ident array
119 >  
120 >  subroutine init_FF(nComponents, c_idents, nExcludes, excludesLocal, status)
121 >    !! Number of components in c_idents array
122      integer, intent(inout) :: nComponents
123 < !! Array of identities nComponents long corresponding to
124 < !! ljatype ident.
125 <    integer, dimension(nComponents),intent(inout) :: ident
123 >    !! Array of identities nComponents long corresponding to
124 >    !! ljatype ident.
125 >    integer, dimension(nComponents),intent(inout) :: c_idents
126      integer :: nExcludes
127      integer, dimension(nExcludes),intent(inout) :: excludesLocal
128 < !!  Result status, success = 0, error = -1
129 <    integer, intent(out) :: Status
130 <
153 <    integer :: alloc_stat
154 <
155 <    integer :: thisStat
156 <    integer :: i
157 <
128 >    !!  Result status, success = 0, error = -1
129 >    integer, intent(out) :: status    
130 >    integer :: i, me, thisStat
131      integer :: myNode
132 +
133   #ifdef IS_MPI
134 <    integer, allocatable, dimension(:) :: identRow
135 <    integer, allocatable, dimension(:) :: identCol
134 >    integer, allocatable, dimension(:) :: c_idents_Row
135 >    integer, allocatable, dimension(:) :: c_idents_Col
136      integer :: nrow
137      integer :: ncol
138   #endif
139 +
140      status = 0
166 –  
141  
168 –    
142  
143 < !! if were're not in MPI, we just update ljatypePtrList
144 < #ifndef IS_MPI
172 <    call create_IdentPtrList(ident,ListHead,identPtrListGlobal,thisStat)
173 <    if ( thisStat /= 0 ) then
174 <       status = -1
175 <       return
176 <    endif
177 <
178 <    
179 < ! if were're in MPI, we also have to worry about row and col lists    
180 < #else
181 <  
182 < ! We can only set up forces if mpiSimulation has been setup.
143 > #ifdef IS_MPI
144 >    ! We can only set up forces if mpiSimulation has been setup.
145      if (.not. isMPISimSet()) then
146         write(default_error,*) "MPI is not set"
147         status = -1
# Line 188 | Line 150 | contains
150      nrow = getNrow(plan_row)
151      ncol = getNcol(plan_col)
152      mynode = getMyNode()
153 < !! Allocate temperary arrays to hold gather information
154 <    allocate(identRow(nrow),stat=alloc_stat)
153 >    
154 >    allocate(c_idents_Row(nrow),stat=alloc_stat)
155      if (alloc_stat /= 0 ) then
156         status = -1
157         return
158      endif
159  
160 <    allocate(identCol(ncol),stat=alloc_stat)
160 >    allocate(c_idents_Col(ncol),stat=alloc_stat)
161      if (alloc_stat /= 0 ) then
162         status = -1
163         return
164      endif
165  
166 < !! Gather idents into row and column idents
167 <
168 <    call gather(ident,identRow,plan_row)
169 <    call gather(ident,identCol,plan_col)
170 <    
209 <  
210 < !! Create row and col pointer lists
211 <  
212 <    call create_IdentPtrlst(identRow,ListHead,identPtrListRow,thisStat)
213 <    if (thisStat /= 0 ) then
166 >    call gather(c_idents, c_idents_Row, plan_row)
167 >    call gather(c_idents, c_idents_Col, plan_col)
168 >
169 >    allocate(atid_Row(nrow),stat=alloc_stat)
170 >    if (alloc_stat /= 0 ) then
171         status = -1
172         return
173      endif
174 <  
175 <    call create_IdentPtrlst(identCol,ListHead,identPtrListColumn,thisStat)
176 <    if (thisStat /= 0 ) then
174 >
175 >    allocate(atid_Col(ncol),stat=alloc_stat)
176 >    if (alloc_stat /= 0 ) then
177         status = -1
178         return
179      endif
180  
181 < !! free temporary ident arrays
181 >    do i = 1, nrow
182 >       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Row(i))
183 >       atid_Row(i) = me
184 >    enddo
185 >
186 >    do i = 1, ncol
187 >       me = getFirstMatchingElement(atypes, "c_ident", c_idents_Col(i))
188 >       atid_Col(i) = me
189 >    enddo
190 >
191 >    !! free temporary ident arrays
192      if (allocated(identCol)) then
193         deallocate(identCol)
194      end if
195      if (allocated(identCol)) then
196         deallocate(identRow)
197      endif
198 +    
199 + #else
200 +    do i = 1, nComponents
201  
202 +       me = getFirstMatchingElement(atypes, "c_ident", c_idents(i))
203 +       atid(i) = me
204 +
205 +    enddo
206   #endif
207 <    
207 >
208      call initForce_Modules(thisStat)
209      if (thisStat /= 0) then
210         status = -1
211         return
212      endif
213 <
214 < !! Create neighbor lists
213 >    
214 >    !! Create neighbor lists
215      call expandList(thisStat)
216      if (thisStat /= 0) then
217         status = -1
218         return
219      endif
246 –
247 –    isFFinit = .true.
220      
221 +    call setupGlobals(thisStat)
222 +    if (thisStat /= 0) then
223 +       status = -1
224 +       return
225 +    endif
226      
227 +    ff_initialized = .true.
228 +    
229 +    
230    end subroutine init_FF
231 +  
232 +  subroutine setupGlobals(thisStat)
233 +    integer, intent(out) :: thisStat
234 +    integer :: nrow
235 +    integer :: ncol
236 +    integer :: nlocal
237 +    integer :: ndim = 3
238 +    integer :: alloc_stat
239 +
240 +    thisStat = 0
241 +
242 + #ifdef IS_MPI
243 +    nrow = getNrow(plan_row)
244 +    ncol = getNcol(plan_col)
245 + #endif
246 +    nlocal = getNlocal()
247 +
248 +    call freeGlobals()
249 +
250 +    allocate(q_Row(ndim,nrow),stat=alloc_stat)
251 +    if (alloc_stat /= 0 ) then
252 +       thisStat = 0
253 +       return
254 +    endif
255 +    
256 +    
257 +    allocate(q_Col(ndim,ncol),stat=alloc_stat)
258 +    if (alloc_stat /= 0 ) then
259 +       thisStat = 0
260 +       return
261 +    endif
262 +    
263 +    
264 +    allocate(u_l_Row(ndim,nrow),stat=alloc_stat)
265 +    if (alloc_stat /= 0 ) then
266 +       thisStat = 0
267 +       return
268 +    endif
269 +    
270  
271 +    allocate(u_l_Col(ndim,ncol),stat=alloc_stat)
272 +    if (alloc_stat /= 0 ) then
273 +       thisStat = 0
274 +       return
275 +    endif
276 +    
277  
278 +    allocate(A_row(9,nrow),stat=alloc_stat)
279 +    if (alloc_stat /= 0 ) then
280 +       thisStat = 0
281 +       return
282 +    endif
283 +    
284 +    
285 +    allocate(A_Col(9,ncol),stat=alloc_stat)
286 +    if (alloc_stat /= 0 ) then
287 +       thisStat = 0
288 +       return
289 +    endif
290 +    
291  
292 +    allocate(pot_row(nrow),stat=alloc_stat)
293 +    if (alloc_stat /= 0 ) then
294 +       thisStat = 0
295 +       return
296 +    endif
297 +    
298 +    allocate(pot_Col(ncol),stat=alloc_stat)
299 +    if (alloc_stat /= 0 ) then
300 +       thisStat = 0
301 +       return
302 +    endif
303 +
304 +    allocate(pot_Temp(nlocal),stat=alloc_stat)
305 +    if (alloc_stat /= 0 ) then
306 +       thisStat = 0
307 +       return
308 +    endif
309 +    
310 +    allocate(f_Row(ndim,nrow),stat=alloc_stat)
311 +    if (alloc_stat /= 0 ) then
312 +       thisStat = 0
313 +       return
314 +    endif
315 +    
316 +    allocate(f_Col(ndim,ncol),stat=alloc_stat)
317 +    if (alloc_stat /= 0 ) then
318 +       thisStat = 0
319 +       return
320 +    endif
321 +    
322 +    allocate(f_Temp(ndim,nlocal),stat=alloc_stat)
323 +    if (alloc_stat /= 0 ) then
324 +       thisStat = 0
325 +       return
326 +    endif
327 +    
328 +    allocate(t_Row(ndim,nrow),stat=alloc_stat)
329 +    if (alloc_stat /= 0 ) then
330 +       thisStat = 0
331 +       return
332 +    endif
333 +    
334 +    allocate(t_Col(ndim,ncol),stat=alloc_stat)
335 +    if (alloc_stat /= 0 ) then
336 +       thisStat = 0
337 +       return
338 +    endif
339 +
340 +    allocate(t_temp(ndim,nlocal),stat=alloc_stat)
341 +    if (alloc_stat /= 0 ) then
342 +       thisStat = 0
343 +       return
344 +    endif
345 +
346 +    allocate(atid_Row(nrow),stat=alloc_stat)
347 +    if (alloc_stat /= 0 ) then
348 +       thisStat = 0
349 +       return
350 +    endif
351 +
352 +    allocate(atid_Col(ncol),stat=alloc_stat)
353 +    if (alloc_stat /= 0 ) then
354 +       thisStat = 0
355 +       return
356 +    endif
357 +
358 + #else
359 +
360 +    allocate(atid(nlocal),stat=alloc_stat)
361 +    if (alloc_stat /= 0 ) then
362 +       thisStat = 0
363 +       return
364 +    end if
365 + #endif
366 +
367 +  end subroutine setupGlobals
368 +  
369 +  subroutine freeGlobals()
370 +
371 + !We free in the opposite order in which we allocate in.
372 + #ifdef IS_MPI
373 +
374 +    if (allocated(atid_Col))   deallocate(atid_Col)
375 +    if (allocated(atid_Row))   deallocate(atid_Row)
376 +    if (allocated(t_Temp))     deallocate(t_Temp)
377 +    if (allocated(t_Col))      deallocate(t_Col)
378 +    if (allocated(t_Row))      deallocate(t_Row)
379 +    if (allocated(f_Temp))     deallocate(f_Temp)
380 +    if (allocated(f_Col))      deallocate(f_Col)
381 +    if (allocated(f_Row))      deallocate(f_Row)
382 +    if (allocated(pot_Temp))   deallocate(pot_Temp)
383 +    if (allocated(pot_Col))    deallocate(pot_Col)
384 +    if (allocated(pot_Row))    deallocate(pot_Row)
385 +    if (allocated(A_Col))      deallocate(A_Col)
386 +    if (allocated(A_Row))      deallocate(A_Row)
387 +    if (allocated(u_l_Col))    deallocate(u_l_Col)
388 +    if (allocated(u_l_Row))    deallocate(u_l_Row)
389 +    if (allocated(q_Col))      deallocate(q_Col)
390 +    if (allocated(q_Row))      deallocate(q_Row)
391 +
392 + #else
393 +
394 +    if (allocated(atid))       deallocate(atid)
395 +
396 + #endif
397 +
398 +
399 +  end subroutine freeGlobals
400 +
401 +
402 +
403 +
404 +  
405    subroutine initForce_Modules(thisStat)
406      integer, intent(out) :: thisStat
407      integer :: my_status
408      
409      thisStat = 0
410 <    call init_lj_FF(ListHead,my_status)
410 >    call init_lj_FF(my_status)
411      if (my_status /= 0) then
412         thisStat = -1
413         return
414      end if
415 <
415 >    
416    end subroutine initForce_Modules
417 <
267 <
417 >  
418   end module forceGlobals

Diff Legend

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