ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/lj_FF.F90
(Generate patch)

Comparing trunk/mdtools/md_code/lj_FF.F90 (file contents):
Revision 246 by chuckv, Fri Jan 24 21:36:52 2003 UTC vs.
Revision 247 by chuckv, Mon Jan 27 18:28:11 2003 UTC

# Line 2 | Line 2
2   !! Corresponds to the force field defined in lj_FF.cpp
3   !! @author Charles F. Vardeman II
4   !! @author Matthew Meineke
5 < !! @version $Id: lj_FF.F90,v 1.10 2003-01-24 21:36:52 chuckv Exp $, $Date: 2003-01-24 21:36:52 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $
5 > !! @version $Id: lj_FF.F90,v 1.11 2003-01-27 18:28:11 chuckv Exp $, $Date: 2003-01-27 18:28:11 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $
6  
7  
8  
9   module lj_ff
10    use simulation
11    use definitions
12 +  use generic_atypes
13   #ifdef IS_MPI
14    use mpiSimulation
15   #endif
# Line 18 | Line 19 | module lj_ff
19   !! Number of lj_atypes in lj_atype_list
20    integer, save :: n_lj_atypes = 0
21  
22 < !! Starting Size for ljMixed Array
23 <  integer, parameter :: ljMixed_blocksize = 10
22 > !! Global list of lj atypes in simulation
23 >  type (lj_atype), pointer :: ljListHead => null()
24 >  type (lj_atype), pointer :: ljListTail => null()
25  
24 !! Basic atom type for a Lennard-Jones Atom.
25  type, public :: lj_atype
26     private
27     sequence
28 !! Unique number for place in linked list
29     integer :: atype_number = 0
30 !! Unique indentifier number (ie atomic no, etc)
31     integer :: atype_ident = 0
32 !! Mass of Particle
33     real ( kind = dp )  :: mass = 0.0_dp
34 !! Lennard-Jones epslon
35     real ( kind = dp )  :: epsilon = 0.0_dp
36 !! Lennard-Jones Sigma
37     real ( kind = dp )  :: sigma = 0.0_dp
38 !! Lennard-Jones Sigma Squared
39     real ( kind = dp )  :: sigma2 = 0.0_dp
40 !! Lennard-Jones Sigma to sixth
41     real ( kind = dp )  :: sigma6 = 0.0_dp
42 !! Pointer for linked list creation
43     type (lj_atype), pointer :: next => null()
44  end type lj_atype
26  
46 !! Pointer type for atype ident array
47  type, public :: lj_atypePtr
48     type (lj_atype), pointer :: this => null()
49  end type lj_atypePtr
50
51 !! Global list of lj atypes in simulation
52  type (lj_atype), pointer :: lj_atype_list => null()
27   !! LJ mixing array  
28    type (lj_atype), dimension(:,:), pointer :: ljMixed => null()
55 !! identity pointer list for force loop.
56  type (lj_atypePtr), dimension(:), pointer :: identPtrList => null()
29  
30  
31   !! Neighbor list and commom arrays
# Line 64 | Line 36 | module lj_ff
36    integer :: nListAllocs = 0
37    integer, parameter :: maxListAllocs = 5
38  
67 #ifdef IS_MPI
68 ! Universal routines: All types of force calculations will need these arrays
69 ! Arrays specific to a type of force calculation should be declared in that module.
70 !! Row position array.
71  real( kind = dp ), allocatable, dimension(:,:) :: qRow
72 !! Column position array.
73  real( kind = dp ), allocatable, dimension(:,:) :: qColumn
39  
40 < !! Row Force array.
41 <  real( kind = dp ), allocatable, dimension(:,:) :: fRow
77 < !! Column Force Array.
78 <  real( kind = dp ), allocatable, dimension(:,:) :: fColumn
79 <
80 < !! Row potential energy array
81 <  real( kind = dp ), allocatable, dimension(:) :: eRow
82 < !! Column potential energy array
83 <  real( kind = dp ), allocatable, dimension(:) :: eColumn
84 <
85 <
40 > !! Atype identity pointer lists
41 > #ifdef IS_MPI
42   !! Row lj_atype pointer list
43    type (lj_atypePtr), dimension(:), pointer :: identPtrListRow => null()
44   !! Column lj_atype pointer list
45    type (lj_atypePtr), dimension(:), pointer :: identPtrListColumn => null()
46 + #else
47 +  type( lj_identPtrList ), dimension(:), pointer :: identPtrList => null()
48   #endif
49  
50  
# Line 105 | Line 63 | contains
63  
64   contains
65  
66 + !! Adds a new lj_atype to the list.
67    subroutine new_lj_atype(ident,mass,epsilon,sigma,status)
68      real( kind = dp ), intent(in) :: mass
69      real( kind = dp ), intent(in) :: epsilon
# Line 112 | Line 71 | contains
71      integer, intent(in) :: ident
72      integer, intent(out) :: status
73  
74 <    type (lj_atype), pointer :: this_lj_atype
116 <    type (lj_atype), pointer :: lj_atype_ptr
117 <
118 <    type (lj_atype), dimension(:,:), pointer :: thisMix
119 <    type (lj_atype), dimension(:,:), pointer :: oldMix
74 >    type (lj_atype), pointer :: newLJ_atype
75      integer :: alloc_error
76      integer :: atype_counter = 0
77      integer :: alloc_size
78 <
78 >    integer :: err_stat
79      status = 0
80  
81  
82  
83   ! allocate a new atype    
84 <    allocate(this_lj_atype,stat=alloc_error)
84 >    allocate(newLJ_atype,stat=alloc_error)
85      if (alloc_error /= 0 ) then
86         status = -1
87         return
88      end if
89  
90   ! assign our new lj_atype information
91 <    this_lj_atype%mass        = mass
92 <    this_lj_atype%epsilon     = epsilon
93 <    this_lj_atype%sigma       = sigma
94 <    this_lj_atype%sigma2      = sigma * sigma
95 <    this_lj_atype%sigma6      = this_lj_atype%sigma2 * this_lj_atype%sigma2 &
96 <         * this_lj_atype%sigma2
91 >    newLJ_atype%mass        = mass
92 >    newLJ_atype%epsilon     = epsilon
93 >    newLJ_atype%sigma       = sigma
94 >    newLJ_atype%sigma2      = sigma * sigma
95 >    newLJ_atype%sigma6      = newLJ_atype%sigma2 * newLJ_atype%sigma2 &
96 >         * newLJ_atype%sigma2
97   ! assume that this atype will be successfully added
98 <    this_lj_atype%atype_ident = ident
99 <    this_lj_atype%atype_number = n_lj_atypes + 1
98 >    newLJ_atype%atype_ident = ident
99 >    newLJ_atype%atype_number = n_lj_atypes + 1
100  
101 <
102 < ! First time through allocate a array of size ljMixed_blocksize
103 <    if(.not. associated(ljMixed)) then
104 <       allocate(thisMix(ljMixed_blocksize,ljMixed_blocksize))
150 <       if (alloc_error /= 0 ) then
151 <          status = -1
152 <          return
153 <       end if
154 <       ljMixed => thisMix
155 < ! If we have outgrown ljMixed_blocksize, allocate a new matrix twice the size and
156 < ! point ljMix at the new matrix.
157 <    else if( (n_lj_atypes + 1) > size(ljMixed)) then
158 <       alloc_size = 2*size(ljMixed)
159 <       allocate(thisMix(alloc_size,alloc_size))
160 <       if (alloc_error /= 0 ) then
161 <          status = -1
162 <          return
163 <       end if
164 < ! point oldMix at old ljMixed array
165 <       oldMix => ljMixed
166 < ! Copy oldMix into new Mixed array      
167 <       thisMix = oldMix
168 < ! Point ljMixed at new array
169 <       ljMixed => thisMix
170 < ! Free old array so we don't have a memory leak
171 <       deallocate(oldMix)
101 >    call add_atype(newLJ_atype,ljListHead,ljListTail,err_stat)
102 >    if (err_stat /= 0 ) then
103 >       status = -1
104 >       return
105      endif
106  
174
175
176
177
178 ! Find bottom of atype master list
179 ! if lj_atype_list is null then we are at the top of the list.
180    if (.not. associated(lj_atype_list)) then
181
182       write(*,*) "Associating lists for first time"
183      
184
185       this_lj_atype%atype_number = 1
186      
187       ljMixed(n_lj_atypes,n_lj_atypes)%sigma  = this_lj_atype%sigma
188      
189       ljMixed(n_lj_atypes,n_lj_atypes)%sigma2  = ljMixed(n_lj_atypes,n_lj_atypes)%sigma &
190            * ljMixed(n_lj_atypes,n_lj_atypes)%sigma
191      
192       ljMixed(n_lj_atypes,n_lj_atypes)%sigma6 = ljMixed(n_lj_atypes,n_lj_atypes)%sigma2 &
193            * ljMixed(n_lj_atypes,n_lj_atypes)%sigma2 &
194            * ljMixed(n_lj_atypes,n_lj_atypes)%sigma2
195      
196       ljMixed(n_lj_atypes,n_lj_atypes)%epsilon = this_lj_atype%epsilon
197      
198       lj_atype_list => this_lj_atype
199      
200       write(*,*) "lj_atype_list first time through -- ident", lj_atype_list%atype_ident
201    else ! we need to find the bottom of the list to insert new atype
202       write(*,*) "Second time through"
203       lj_atype_ptr => lj_atype_list
204       write(*,*) "lj_atype_ptr to next"
205       atype_counter = 1
206       find_end: do
207          if (.not. associated(lj_atype_ptr)) exit find_end
208
209          write(*,*) "Inside find_end-this ident", lj_atype_ptr%atype_ident
210 ! Set up mixing for new atype and current atype in list
211       ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma  =  &
212            calcLJMix("sigma",this_lj_atype%sigma, &
213            lj_atype_ptr%sigma)
214
215       ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2  = &
216            ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma &
217            * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma
218
219       ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma6 = &
220            ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 &
221            * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 &
222            * ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2
223
224       ljMixed(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%epsilon = &
225            calcLJMix("epsilon",this_lj_atype%epsilon, &
226            lj_atype_ptr%epsilon)
227
228 ! Advance to next pointer
229       lj_atype_ptr => lj_atype_ptr%next          
230       atype_counter = atype_counter + 1          
231       end do find_end
232
233       lj_atype_ptr => this_lj_atype
234       write(*,*) "just added: ", lj_atype_ptr%atype_ident
235    end if
236
237    write(*,*) "lj_atype_list now contains.."
238    
239    lj_atype_ptr => lj_atype_list
240    do while (associated(lj_atype_ptr))
241       write(*,*) "lj_atype_list contains ident", lj_atype_ptr%atype_ident
242       lj_atype_ptr => lj_atype_ptr%next
243    end do
244    
245 ! Insert new atype at end of list
246
247 ! Increment number of atypes
248
107      n_lj_atypes = n_lj_atypes + 1
108  
251 ! Set up self mixing
109  
110    end subroutine new_lj_atype
111  
# Line 274 | Line 131 | contains
131      integer :: alloc_stat
132   #endif
133      status = 0
277
278    write(*,*) "Initializing ljFF"
134   !! if were're not in MPI, we just update ljatypePtrList
135   #ifndef IS_MPI
136 <    call new_ljatypePtrList(nComponents,ident,identPtrList,thisStat)
136 >    call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat)
137      if ( thisStat /= 0 ) then
138         status = -1
139         return
140      endif
141 <    write(*,*) "component pointer list initialized"
141 >
142   !! Allocate pointer lists
143      allocate(point(nComponents),stat=alloc_stat)
144      if (alloc_stat /=0) then
# Line 323 | Line 178 | contains
178      call gather(ident,identCol,plan_col)
179  
180   !! Create row and col pointer lists
181 <    call new_ljatypePtrList(nrow,identRow,identPtrListRow,thisStat)
181 >    call create_IdentPtrlst(identRow,ljListTail,identPtrListRow,thisStat)
182      if (thisStat /= 0 ) then
183         status = -1
184         return
185      endif
186  
187 <    call new_ljatypePtrList(ncol,identCol,identPtrListCol,thisStat)
187 >    call create_IdentPtrlst(identCol,ljListTail,identPtrListCol,thisStat)
188      if (thisStat /= 0 ) then
189         status = -1
190         return
# Line 339 | Line 194 | contains
194      deallocate(identCol)
195      deallocate(identRow)
196  
342 !! Allocate Simulation arrays
343 !! NOTE: This bit of code should be fixed, it can cause large
344 !! memory fragmentation if call repeatedly
197  
346    if (.not.allocated(qRow)) then
347       allocate(qRow(3,nrow),stat=alloc_stat)
348       if (alloc_stat /= 0 ) then
349          status = -1
350          return
351       endif
352    else
353       deallocate(qrow)
354       allocate(qRow(3,nrow),stat=alloc_stat)
355       if (alloc_stat /= 0 ) then
356          status = -1
357          return
358       endif
359    endif
360
361    if (.not.allocated(3,qCol)) then
362       allocate(qCol(ncol),stat=alloc_stat)
363       if (alloc_stat /= 0 ) then
364          status = -1
365          return
366       endif
367    else
368       deallocate(qCol)
369       allocate(qCol(3,ncol),stat=alloc_stat)
370       if (alloc_stat /= 0 ) then
371          status = -1
372          return
373       endif
374    endif
375
376    if (.not.allocated(fRow)) then
377       allocate(fRow(3,nrow),stat=alloc_stat)
378       if (alloc_stat /= 0 ) then
379          status = -1
380          return
381       endif
382    else
383       deallocate(fRow)
384       allocate(fRow(3,nrow),stat=alloc_stat)
385       if (alloc_stat /= 0 ) then
386          status = -1
387          return
388       endif
389    endif
390
391    if (.not.allocated(fCol)) then
392       allocate(fCol(3,ncol),stat=alloc_stat)
393       if (alloc_stat /= 0 ) then
394          status = -1
395          return
396       endif
397    else
398       deallocate(fCol)
399       allocate(fCol(3,ncol),stat=alloc_stat)
400       if (alloc_stat /= 0 ) then
401          status = -1
402          return
403       endif
404    endif
405
406    if (.not.allocated(eRow)) then
407       allocate(eRow(nrow),stat=alloc_stat)
408       if (alloc_stat /= 0 ) then
409          status = -1
410          return
411       endif
412    else
413       deallocate(eRow)
414       allocate(eRow(nrow),stat=alloc_stat)
415       if (alloc_stat /= 0 ) then
416          status = -1
417          return
418       endif
419    endif
420
421    if (.not.allocated(eCol)) then
422       allocate(eCol(ncol),stat=alloc_stat)
423       if (alloc_stat /= 0 ) then
424          status = -1
425          return
426       endif
427    else
428       deallocate(eCol)
429       allocate(eCol(ncol),stat=alloc_stat)
430       if (alloc_stat /= 0 ) then
431          status = -1
432          return
433       endif
434    endif
435
436    if (.not.allocated(eTemp)) then
437       allocate(eTemp(nComponents),stat=alloc_stat)
438       if (alloc_stat /= 0 ) then
439          status = -1
440          return
441       endif
442    else
443       deallocate(eTemp)
444       allocate(eTemp(nComponets),stat=alloc_stat)
445       if (alloc_stat /= 0 ) then
446          status = -1
447          return
448       endif
449    endif
450
451
198   !! Allocate neighbor lists for mpi simulations.
199      if (.not. allocated(point)) then
200         allocate(point(nrow),stat=alloc_stat)
# Line 480 | Line 226 | contains
226  
227   #endif
228      
229 <    
230 <    do i = 1, nComponents
231 <       write(*,*) "Ident pointer list ident: ", identPtrList(i)%this%atype_ident
232 <    end do
233 <
229 >    call createMixingList(thisStat)
230 >    if (thisStat /= 0) then
231 >       status = -1
232 >       return
233 >    endif
234      isljFFinit = .true.
235  
236    end subroutine init_ljFF
# Line 493 | Line 239 | contains
239  
240  
241  
496 !! Takes an ident array and creates an atype pointer list
497 !! based on those identities
498  subroutine new_ljatypePtrList(mysize,ident,PtrList,status)
499    integer, intent(in) :: mysize
500    integer,dimension(:), intent(in) :: ident
501    integer, optional :: status
502    type(lj_atypePtr), dimension(:), pointer :: PtrList
242  
243 <    integer :: thisIdent
243 >  subroutine createMixingList(status)
244 >    integer :: listSize
245 >    integer :: status
246      integer :: i
247 <    integer :: alloc_error
507 <    type (lj_atype), pointer :: tmpPtr
247 >    integer :: j
248  
249 <    if (present(status)) status = 0
249 >    integer :: outerCounter = 0
250 >    integer :: innerCounter = 0
251 >    type (lj_atype), pointer :: tmpPtrOuter => null()
252 >    type (lj_atype), pointer :: tmpPtrInner => null()
253 >    status = 0
254  
255 <    write(*,*) "Creating new ljatypePtrList...."
256 < ! First time through, allocate list
257 <    if (.not.associated(PtrList)) then
258 <       write(*,*) "allocating pointer list...."
259 <       allocate(PtrList(mysize))
260 <    else
517 < ! We want to creat a new ident list so free old list
518 <       deallocate(PtrList)
519 <       allocate(PtrList(mysize))
520 <    endif
255 >    listSize = getListLen(ljListHead)
256 >    if (listSize == 0) then
257 >       status = -1
258 >       return
259 >    end if
260 >  
261  
262 < ! Match pointer list
263 <    write(*,*) "Doing ident search"
264 <    do i = 1, mysize
265 <       thisIdent = ident(i)
266 <       write(*,*) "Calling getLJatype for ident ", thisIdent
267 <       call getLJatype(thisIdent,tmpPtr)
268 <
269 <      if (.not. associated(tmpPtr)) then
270 <         write(*,*) "tmpptr was not associated"
271 <          status = -1
272 <          return
262 >    if (.not. associated(ljMixed)) then
263 >       allocate(ljMixed(listSize,listSize))
264 >    else
265 >       status = -1
266 >       return
267 >    end if
268 >
269 >    
270 >    write(*,*) "Creating mixing list"
271 >    tmpPtrOuter => ljListHead
272 >    tmpPtrInner => tmpPtrOuter%next
273 >    do while (associated(tmpPtrOuter))
274 >       outerCounter = outerCounter + 1
275 >       write(*,*) "Outer index is: ", outerCounter
276 >       write(*,*) "For atype: ", tmpPtrOuter%atype_ident
277 > ! do self mixing rule
278 >       ljMixed(outerCounter,outerCounter)%sigma  = tmpPtrOuter%sigma
279 >                                                                                                  
280 >       ljMixed(outerCounter,outerCounter)%sigma2  = ljMixed(outerCounter,outerCounter)%sigma &
281 >            * ljMixed(outerCounter,outerCounter)%sigma
282 >                                                                                                  
283 >       ljMixed(outerCounter,outerCounter)%sigma6 = ljMixed(outerCounter,outerCounter)%sigma2 &
284 >            * ljMixed(outerCounter,outerCounter)%sigma2 &
285 >            * ljMixed(outerCounter,outerCounter)%sigma2
286 >                                                                                                  
287 >       ljMixed(outerCounter,outerCounter)%epsilon = tmpPtrOuter%epsilon
288 >
289 >       innerCounter = outerCounter + 1
290 >       do while (associated(tmpPtrInner))
291 >          
292 >          write(*,*) "Entered inner loop for: ", innerCounter
293 >          ljMixed(outerCounter,innerCounter)%sigma  =  &
294 >               calcLJMix("sigma",tmpPtrOuter%sigma, &
295 >               tmpPtrInner%sigma)
296 >          
297 >          ljMixed(outerCounter,innerCounter)%sigma2  = &
298 >               ljMixed(outerCounter,innerCounter)%sigma &
299 >               * ljMixed(outerCounter,innerCounter)%sigma
300 >          
301 >          ljMixed(outerCounter,innerCounter)%sigma6 = &
302 >               ljMixed(outerCounter,innerCounter)%sigma2 &
303 >               * ljMixed(outerCounter,innerCounter)%sigma2 &
304 >               * ljMixed(outerCounter,innerCounter)%sigma2
305 >          
306 >          ljMixed(outerCounter,innerCounter)%epsilon = &
307 >               calcLJMix("epsilon",tmpPtrOuter%epsilon, &
308 >               tmpPtrInner%epsilon)
309 >          ljMixed(innerCounter,outerCounter)%sigma = ljMixed(outerCounter,innerCounter)%sigma
310 >          ljMixed(innerCounter,outerCounter)%sigma2 = ljMixed(outerCounter,innerCounter)%sigma2
311 >          ljMixed(innerCounter,outerCounter)%sigma6 = ljMixed(outerCounter,innerCounter)%sigma6
312 >          ljMixed(innerCounter,outerCounter)%epsilon = ljMixed(outerCounter,innerCounter)%epsilon
313 >
314 >
315 >          tmpPtrInner => tmpPtrInner%next
316 >          innerCounter = innerCounter + 1
317 >       end do
318 > ! advance pointers
319 >       tmpPtrOuter => tmpPtrOuter%next
320 >       if (associated(tmpPtrOuter)) then
321 >          tmpPtrInner => tmpPtrOuter%next
322         endif
323        
535       PtrList(i)%this => tmpPtr
324      end do
325  
326 <  end subroutine new_ljatypePtrList
326 >  end subroutine createMixingList
327  
540 !! Finds a lj_atype based upon numerical ident
541 !! returns a null pointer if error
542  subroutine getLJatype(ident,ljAtypePtr)
543    integer, intent(in) :: ident
544    type (lj_atype), pointer :: ljAtypePtr    
545    type (lj_atype), pointer :: tmplj_atype_ptr => null()
328  
547    write(*,*) "Finding ljAtype for ident ", ident
548    ljAtypePtr => null()
329  
550    if(.not. associated(lj_atype_list)) return
330  
552 ! Point at head of list.
553    write(*,*) "Searching at head of list"
554    tmplj_atype_ptr => lj_atype_list
555    find_ident: do
556       if (.not.associated(tmplj_atype_ptr)) then
557          write(*,*) "Find_ident, pointer not associated"
558          exit find_ident
559       else if( tmplj_atype_ptr%atype_ident == ident) then
560          ljAtypePtr => tmplj_atype_ptr
561          exit find_ident
562       endif
563       tmplj_atype_ptr => tmplj_atype_ptr%next
564    end do find_ident
331  
332  
333  
334  
569  end subroutine getLJatype
570
571
335   !! FORCE routine Calculates Lennard Jones forces.
336   !------------------------------------------------------------->
337    subroutine do_lj_ff(q,f,potE,do_pot)
# Line 589 | Line 352 | contains
352    real( kind = DP ), dimension(3,getNlocal()) :: efr
353   #endif
354    
355 + !! Local arrays needed for MPI
356 + #ifdef IS_MPI
357 +  real(kind = dp), dimension(3:getNrow()) :: qRow
358 +  real(kind = dp), dimension(3:getNcol()) :: qCol
359 +
360 +  real(kind = dp), dimension(3:getNrow()) :: fRow
361 +  real(kind = dp), dimension(3:getNcol()) :: fCol
362 +
363 +  real(kind = dp), dimension(getNrow()) :: eRow
364 +  real(kind = dp), dimension(getNcol()) :: eCol
365 +
366 +  real(kind = dp), dimension(getNlocal()) :: eTemp
367 + #endif
368 +
369 +
370 +
371    real( kind = DP )   :: pe
372    logical             :: update_nlist
373  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines