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 216 by chuckv, Sun Dec 29 16:43:11 2002 UTC vs.
Revision 246 by chuckv, Fri Jan 24 21:36:52 2003 UTC

# Line 1 | Line 1
1 + !! Calculates Long Range forces Lennard-Jones interactions.
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 $
6 +
7 +
8 +
9   module lj_ff
10    use simulation
11 <  use definitions, ONLY : dp, ndim
11 >  use definitions
12 > #ifdef IS_MPI
13 >  use mpiSimulation
14 > #endif
15    implicit none
16 +  PRIVATE
17  
18 + !! Number of lj_atypes in lj_atype_list
19    integer, save :: n_lj_atypes = 0
20  
21 + !! Starting Size for ljMixed Array
22 +  integer, parameter :: ljMixed_blocksize = 10
23 +
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 <     character(len = 15) :: name = "NULL"
32 > !! Mass of Particle
33       real ( kind = dp )  :: mass = 0.0_dp
34 <     real ( kind = dp )  :: epslon = 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
45  
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()
53 + !! LJ mixing array  
54 +  type (lj_atype), dimension(:,:), pointer :: ljMixed => null()
55 + !! identity pointer list for force loop.
56 +  type (lj_atypePtr), dimension(:), pointer :: identPtrList => null()
57  
58 +
59 + !! Neighbor list and commom arrays
60 + !! This should eventally get moved to a neighbor list type
61 +  integer, allocatable, dimension(:) :: point
62 +  integer, allocatable, dimension(:) :: list
63 +  integer, parameter :: listMultiplier = 80
64 +  integer :: nListAllocs = 0
65 +  integer, parameter :: maxListAllocs = 5
66 +
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
74 +
75 + !! Row Force array.
76 +  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 +
86 + !! Row lj_atype pointer list
87 +  type (lj_atypePtr), dimension(:), pointer :: identPtrListRow => null()
88 + !! Column lj_atype pointer list
89 +  type (lj_atypePtr), dimension(:), pointer :: identPtrListColumn => null()
90 + #endif
91 +
92 +
93 + !! Logical has lj force field module been initialized?
94 +  logical, save :: isljFFinit = .false.
95 +
96 +
97 + !! Public methods and data
98    public :: new_lj_atype
99 +  public :: do_lj_ff
100 +  public :: getLjPot
101 +  public :: init_ljFF
102  
103 +  
104 +
105 +
106   contains
107  
108 <  subroutine new_lj_atype(name,mass,epslon,sigma,status)
26 <    character( len = 15 ) :: name
108 >  subroutine new_lj_atype(ident,mass,epsilon,sigma,status)
109      real( kind = dp ), intent(in) :: mass
110 <    real( kind = dp ), intent(in) :: epslon
111 <    real( kind = dp ), intent(in) :: sigam
110 >    real( kind = dp ), intent(in) :: epsilon
111 >    real( kind = dp ), intent(in) :: sigma
112 >    integer, intent(in) :: ident
113      integer, intent(out) :: status
114  
115      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
120      integer :: alloc_error
121      integer :: atype_counter = 0
122 +    integer :: alloc_size
123  
124      status = 0
125  
126 +
127 +
128 + ! allocate a new atype    
129      allocate(this_lj_atype,stat=alloc_error)
130      if (alloc_error /= 0 ) then
131         status = -1
# Line 44 | Line 133 | contains
133      end if
134  
135   ! assign our new lj_atype information
136 <    this_lj_atype%name       = name
137 <    this_lj_atype%mass       = mass
138 <    this_lj_atype%epslon     = epslon
139 <    this_lj_atype%sigma      = sigma
136 >    this_lj_atype%mass        = mass
137 >    this_lj_atype%epsilon     = epsilon
138 >    this_lj_atype%sigma       = sigma
139 >    this_lj_atype%sigma2      = sigma * sigma
140 >    this_lj_atype%sigma6      = this_lj_atype%sigma2 * this_lj_atype%sigma2 &
141 >         * this_lj_atype%sigma2
142 > ! assume that this atype will be successfully added
143 >    this_lj_atype%atype_ident = ident
144 >    this_lj_atype%atype_number = n_lj_atypes + 1
145  
146  
147 + ! First time through allocate a array of size ljMixed_blocksize
148 +    if(.not. associated(ljMixed)) then
149 +       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)
172 +    endif
173 +
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 <       lj_atype_ptr => this_lj_atype
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
57 –    else ! we need to find the bottom of the list
58 –       lj_atype_ptr => lj_atype_list%next
206         find_end: do
207 <          if (.not. associated(lj_atype_ptr%next)) then
208 <             exit find_end
209 <          end if
210 <          lj_atype_ptr => lj_atype_ptr%next
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 => this_lj_atype
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 +
249 +    n_lj_atypes = n_lj_atypes + 1
250 +
251 + ! Set up self mixing
252 +
253    end subroutine new_lj_atype
254  
71 –  subroutine add_lj_atype()
255  
256 +  subroutine init_ljFF(nComponents,ident, status)
257 + !! Number of components in ident array
258 +    integer, intent(inout) :: nComponents
259 + !! Array of identities nComponents long corresponding to
260 + !! ljatype ident.
261 +    integer, dimension(nComponents),intent(inout) :: ident
262 + !!  Result status, success = 0, error = -1
263 +    integer, intent(out) :: Status
264  
265 +    integer :: alloc_stat
266  
267 <  end subroutine add_lj_atype
267 >    integer :: thisStat
268 >    integer :: i
269 > #ifdef IS_MPI
270 >    integer, allocatable, dimension(:) :: identRow
271 >    integer, allocatable, dimension(:) :: identCol
272 >    integer :: nrow
273 >    integer :: ncol
274 >    integer :: alloc_stat
275 > #endif
276 >    status = 0
277 >
278 >    write(*,*) "Initializing ljFF"
279 > !! if were're not in MPI, we just update ljatypePtrList
280 > #ifndef IS_MPI
281 >    call new_ljatypePtrList(nComponents,ident,identPtrList,thisStat)
282 >    if ( thisStat /= 0 ) then
283 >       status = -1
284 >       return
285 >    endif
286 >    write(*,*) "component pointer list initialized"
287 > !! Allocate pointer lists
288 >    allocate(point(nComponents),stat=alloc_stat)
289 >    if (alloc_stat /=0) then
290 >       status = -1
291 >       return
292 >    endif
293 >
294 >    allocate(list(nComponents*listMultiplier),stat=alloc_stat)
295 >    if (alloc_stat /=0) then
296 >       status = -1
297 >       return
298 >    endif
299 >    
300 > ! if were're in MPI, we also have to worry about row and col lists    
301 > #else
302 > ! We can only set up forces if mpiSimulation has been setup.
303 >    if (.not. isMPISimSet()) then
304 >       status = -1
305 >       return
306 >    endif
307 >    nrow = getNrow()
308 >    ncol = getNcol()
309 > !! Allocate temperary arrays to hold gather information
310 >    allocate(identRow(nrow),stat=alloc_stat)
311 >    if (alloc_stat /= 0 ) then
312 >       status = -1
313 >       return
314 >    endif
315 >
316 >    allocate(identCol(ncol),stat=alloc_stat)
317 >    if (alloc_stat /= 0 ) then
318 >       status = -1
319 >       return
320 >    endif
321 > !! Gather idents into row and column idents
322 >    call gather(ident,identRow,plan_row)
323 >    call gather(ident,identCol,plan_col)
324 >
325 > !! Create row and col pointer lists
326 >    call new_ljatypePtrList(nrow,identRow,identPtrListRow,thisStat)
327 >    if (thisStat /= 0 ) then
328 >       status = -1
329 >       return
330 >    endif
331 >
332 >    call new_ljatypePtrList(ncol,identCol,identPtrListCol,thisStat)
333 >    if (thisStat /= 0 ) then
334 >       status = -1
335 >       return
336 >    endif
337 >
338 > !! free temporary ident arrays
339 >    deallocate(identCol)
340 >    deallocate(identRow)
341 >
342 > !! Allocate Simulation arrays
343 > !! NOTE: This bit of code should be fixed, it can cause large
344 > !! memory fragmentation if call repeatedly
345 >
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 >
452 > !! Allocate neighbor lists for mpi simulations.
453 >    if (.not. allocated(point)) then
454 >       allocate(point(nrow),stat=alloc_stat)
455 >       if (alloc_stat /=0) then
456 >          status = -1
457 >          return
458 >       endif
459 >
460 >       allocate(list(ncol*listMultiplier),stat=alloc_stat)
461 >       if (alloc_stat /=0) then
462 >          status = -1
463 >          return
464 >       endif
465 >    else
466 >       deallocate(list)
467 >       deallocate(point)
468 >       allocate(point(nrow),stat=alloc_stat)
469 >       if (alloc_stat /=0) then
470 >          status = -1
471 >          return
472 >       endif
473 >
474 >       allocate(list(ncol*listMultiplier),stat=alloc_stat)
475 >       if (alloc_stat /=0) then
476 >          status = -1
477 >          return
478 >       endif
479 >    endif
480 >
481 > #endif
482 >    
483 >    
484 >    do i = 1, nComponents
485 >       write(*,*) "Ident pointer list ident: ", identPtrList(i)%this%atype_ident
486 >    end do
487 >
488 >    isljFFinit = .true.
489 >
490 >  end subroutine init_ljFF
491 >
492 >
493 >
494 >
495 >
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
503 >
504 >    integer :: thisIdent
505 >    integer :: i
506 >    integer :: alloc_error
507 >    type (lj_atype), pointer :: tmpPtr
508 >
509 >    if (present(status)) status = 0
510 >
511 >    write(*,*) "Creating new ljatypePtrList...."
512 > ! First time through, allocate list
513 >    if (.not.associated(PtrList)) then
514 >       write(*,*) "allocating pointer list...."
515 >       allocate(PtrList(mysize))
516 >    else
517 > ! We want to creat a new ident list so free old list
518 >       deallocate(PtrList)
519 >       allocate(PtrList(mysize))
520 >    endif
521  
522 + ! Match pointer list
523 +    write(*,*) "Doing ident search"
524 +    do i = 1, mysize
525 +       thisIdent = ident(i)
526 +       write(*,*) "Calling getLJatype for ident ", thisIdent
527 +       call getLJatype(thisIdent,tmpPtr)
528 +
529 +      if (.not. associated(tmpPtr)) then
530 +         write(*,*) "tmpptr was not associated"
531 +          status = -1
532 +          return
533 +       endif
534 +      
535 +       PtrList(i)%this => tmpPtr
536 +    end do
537  
538 +  end subroutine new_ljatypePtrList
539  
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()
546  
547 +    write(*,*) "Finding ljAtype for ident ", ident
548 +    ljAtypePtr => null()
549 +
550 +    if(.not. associated(lj_atype_list)) return
551 +
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
565 +
566 +
567 +
568 +
569 +  end subroutine getLJatype
570 +
571 +
572 + !! FORCE routine Calculates Lennard Jones forces.
573 + !------------------------------------------------------------->
574 +  subroutine do_lj_ff(q,f,potE,do_pot)
575 + !! Position array provided by C, dimensioned by getNlocal
576 +    real ( kind = dp ), dimension(3,getNlocal()) :: q
577 + !! Force array provided by C, dimensioned by getNlocal
578 +    real ( kind = dp ), dimension(3,getNlocal()) :: f
579 +    real ( kind = dp ) :: potE
580 +    logical ( kind = 2) :: do_pot
581 +    
582 +    type(lj_atype), pointer :: ljAtype_i
583 +    type(lj_atype), pointer :: ljAtype_j
584 +
585 + #ifdef MPI
586 +  real( kind = DP ), dimension(3,getNcol()) :: efr
587 +  real( kind = DP ) :: pot_local
588 + #else
589 +  real( kind = DP ), dimension(3,getNlocal()) :: efr
590 + #endif
591 +  
592 +  real( kind = DP )   :: pe
593 +  logical             :: update_nlist
594 +
595 +
596 +  integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
597 +  integer :: nlist
598 +  integer :: j_start
599 +  integer :: tag_i,tag_j
600 +  real( kind = DP ) ::  r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
601 +  real( kind = DP ) ::  rxi, ryi, rzi, rxij, ryij, rzij, rijsq
602 +  real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
603 +
604 + ! a rig that need to be fixed.
605 + #ifdef IS_MPI
606 +  logical :: newtons_thrd = .true.
607 +  real( kind = dp ) :: pe_local
608 +  integer :: nlocal
609 + #endif
610 +  integer :: nrow
611 +  integer :: ncol
612 +  integer :: natoms
613 +
614 + ! Make sure we are properly initialized.
615 +  if (.not. isljFFInit) then
616 +     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
617 +     return
618 +  endif
619 + #ifdef IS_MPI
620 +    if (.not. isMPISimSet()) then
621 +     write(default_error,*) "ERROR: mpiSimulation has not been properly initialized"
622 +     return
623 +  endif
624 + #endif
625 +
626 + !! initialize local variables  
627 +  natoms = getNlocal()
628 +  call getRcut(rcut,rcut2=rcutsq)
629 +  call getRlist(rlist,rlistsq)
630 + #ifndef IS_MPI
631 +  nrow = natoms - 1
632 +  ncol = natoms
633 + #else
634 +  nrow = getNrow(plan_row)
635 +  ncol = getNcol(plan_col)
636 +  nlocal = natoms
637 +  j_start = 1
638 + #endif
639 +
640 +  
641 + !! See if we need to update neighbor lists
642 +  call check(q,update_nlist)
643 +
644 + !--------------WARNING...........................
645 + ! Zero variables, NOTE:::: Forces are zeroed in C
646 + ! Zeroing them here could delete previously computed
647 + ! Forces.
648 + !------------------------------------------------
649 + #ifndef IS_MPI
650 + !  nloops = nloops + 1
651 +  pe = 0.0E0_DP
652 +
653 + #else
654 +    f_row = 0.0E0_DP
655 +    f_col = 0.0E0_DP
656 +
657 +    pe_local = 0.0E0_DP
658 +
659 +    eRow = 0.0E0_DP
660 +    eCol = 0.0E0_DP
661 +    eTemp = 0.0E0_DP
662 + #endif
663 +    efr = 0.0E0_DP
664 +
665 + ! communicate MPI positions
666 + #ifdef MPI    
667 +    call gather(q,qRow,plan_row3)
668 +    call gather(q,qCol,plan_col3)
669 + #endif
670 +
671 +
672 +  if (update_nlist) then
673 +
674 +     ! save current configuration, contruct neighbor list,
675 +     ! and calculate forces
676 +     call save_nlist(q)
677 +    
678 +     nlist = 0
679 +    
680 +    
681 +
682 +     do i = 1, nrow
683 +        point(i) = nlist + 1
684 + #ifdef MPI
685 +        ljAtype_i => identPtrListRow(i)%this
686 +        tag_i = tagRow(i)
687 +        rxi = qRow(1,i)
688 +        ryi = qRow(2,i)
689 +        rzi = qRow(3,i)
690 + #else
691 +        ljAtype_i   => identPtrList(i)%this
692 +        j_start = i + 1
693 +        rxi = q(1,i)
694 +        ryi = q(2,i)
695 +        rzi = q(3,i)
696 + #endif
697 +
698 +        inner: do j = j_start, ncol
699 + #ifdef MPI
700 + ! Assign identity pointers and tags
701 +           ljAtype_j => identPtrListColumn(j)%this
702 +           tag_j = tagCol(j)
703 +           if (newtons_thrd) then
704 +              if (tag_i <= tag_j) then
705 +                 if (mod(tag_i + tag_j,2) == 0) cycle inner
706 +              else                
707 +                 if (mod(tag_i + tag_j,2) == 1) cycle inner
708 +              endif
709 +           endif
710 +
711 +           rxij = wrap(rxi - qCol(1,j), 1)
712 +           ryij = wrap(ryi - qCol(2,j), 2)
713 +           rzij = wrap(rzi - qCol(3,j), 3)
714 + #else          
715 +           ljAtype_j   => identPtrList(j)%this
716 +           rxij = wrap(rxi - q(1,j), 1)
717 +           ryij = wrap(ryi - q(2,j), 2)
718 +           rzij = wrap(rzi - q(3,j), 3)
719 +      
720 + #endif          
721 +           rijsq = rxij*rxij + ryij*ryij + rzij*rzij
722 +
723 + #ifdef MPI
724 +             if (rijsq <=  rlistsq .AND. &
725 +                  tag_j /= tag_i) then
726 + #else
727 +             if (rijsq <  rlistsq) then
728 + #endif
729 +            
730 +              nlist = nlist + 1
731 +              if (nlist > size(list)) then
732 + !!  "Change how nlist size is done"
733 +                 write(DEFAULT_ERROR,*) "ERROR: nlist > list size"
734 +              endif
735 +              list(nlist) = j
736 +
737 +              
738 +              if (rijsq <  rcutsq) then
739 +                
740 +                 r = dsqrt(rijsq)
741 +      
742 +                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
743 +      
744 + #ifdef MPI
745 +                eRow(i) = eRow(i) + pot*0.5
746 +                eCol(i) = eCol(i) + pot*0.5
747 + #else
748 +                pe = pe + pot
749 + #endif                
750 +            
751 +                 efr(1,j) = -rxij
752 +                 efr(2,j) = -ryij
753 +                 efr(3,j) = -rzij
754 +
755 +                 do dim = 1, 3  
756 +
757 +            
758 +                    drdx1 = efr(dim,j) / r
759 +                    ftmp = dudr * drdx1
760 +
761 +
762 + #ifdef MPI
763 +                    fCol(dim,j) = fCol(dim,j) - ftmp
764 +                    fRow(dim,i) = fRow(dim,i) + ftmp
765 + #else                    
766 +            
767 +                    f(dim,j) = f(dim,j) - ftmp
768 +                    f(dim,i) = f(dim,i) + ftmp
769 +
770 + #endif                    
771 +                 enddo
772 +              endif
773 +           endif
774 +        enddo inner
775 +     enddo
776 +
777 + #ifdef MPI
778 +     point(nrow + 1) = nlist + 1
779 + #else
780 +     point(natoms) = nlist + 1
781 + #endif
782 +
783 +  else
784 +
785 +     ! use the list to find the neighbors
786 +     do i = 1, nrow
787 +        JBEG = POINT(i)
788 +        JEND = POINT(i+1) - 1
789 +        ! check thiat molecule i has neighbors
790 +        if (jbeg .le. jend) then
791 + #ifdef MPI
792 +           ljAtype_i => identPtrListRow(i)%this
793 +           rxi = qRow(1,i)
794 +           ryi = qRow(2,i)
795 +           rzi = qRow(3,i)
796 + #else
797 +           ljAtype_i   => identPtrList(i)%this
798 +           rxi = q(1,i)
799 +           ryi = q(2,i)
800 +           rzi = q(3,i)
801 + #endif
802 +           do jnab = jbeg, jend
803 +              j = list(jnab)
804 + #ifdef MPI
805 +              ljAtype_j = identPtrListColumn(j)%this
806 +              rxij = wrap(rxi - q_col(1,j), 1)
807 +              ryij = wrap(ryi - q_col(2,j), 2)
808 +              rzij = wrap(rzi - q_col(3,j), 3)
809 + #else
810 +              ljAtype_j = identPtrList(j)%this
811 +              rxij = wrap(rxi - q(1,j), 1)
812 +              ryij = wrap(ryi - q(2,j), 2)
813 +              rzij = wrap(rzi - q(3,j), 3)
814 + #endif
815 +              rijsq = rxij*rxij + ryij*ryij + rzij*rzij
816 +              
817 +              if (rijsq <  rcutsq) then
818 +
819 +                 r = dsqrt(rijsq)
820 +                
821 +                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
822 + #ifdef MPI
823 +                eRow(i) = eRow(i) + pot*0.5
824 +                eCol(i) = eCol(i) + pot*0.5
825 + #else
826 +                pe = pe + pot
827 + #endif                
828 +
829 +                
830 +                 efr(1,j) = -rxij
831 +                 efr(2,j) = -ryij
832 +                 efr(3,j) = -rzij
833 +
834 +                 do dim = 1, 3                        
835 +                    
836 +                    drdx1 = efr(dim,j) / r
837 +                    ftmp = dudr * drdx1
838 + #ifdef MPI
839 +                    fCol(dim,j) = fCol(dim,j) - ftmp
840 +                    fRow(dim,i) = fRow(dim,i) + ftmp
841 + #else                    
842 +                    f(dim,j) = f(dim,j) - ftmp
843 +                    f(dim,i) = f(dim,i) + ftmp
844 + #endif                    
845 +                 enddo
846 +              endif
847 +           enddo
848 +        endif
849 +     enddo
850 +  endif
851 +
852 +
853 +
854 + #ifdef MPI
855 +    !!distribute forces
856 +    call scatter(fRow,f,plan_row3)
857 +
858 +    call scatter(fCol,f_tmp,plan_col3)
859 +
860 +    do i = 1,nlocal
861 +       f(1:3,i) = f(1:3,i) + f_tmp(1:3,i)
862 +    end do
863 +
864 +
865 +    
866 +    if (do_pot) then
867 + ! scatter/gather pot_row into the members of my column
868 +       call scatter(e_row,e_tmp,plan_row)
869 +      
870 +       ! scatter/gather pot_local into all other procs
871 +       ! add resultant to get total pot
872 +       do i = 1, nlocal
873 +          pe_local = pe_local + e_tmp(i)
874 +       enddo
875 +       if (newtons_thrd) then
876 +          e_tmp = 0.0E0_DP
877 +          call scatter(e_col,e_tmp,plan_col)
878 +          do i = 1, nlocal
879 +             pe_local = pe_local + e_tmp(i)
880 +          enddo
881 +       endif
882 +    endif
883 +    potE = pe_local
884 + #endif
885 +
886 +    potE = pe
887 +
888 +
889 +  end subroutine do_lj_ff
890 +
891 + !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
892 + !! derivatives.
893 +  subroutine getLjPot(r,pot,dudr,atype1,atype2,d2,status)
894 + ! arguments
895 + !! Length of vector between particles
896 +    real( kind = dp ), intent(in)  :: r
897 + !! Potential Energy
898 +    real( kind = dp ), intent(out) :: pot
899 + !! Derivatve wrt postion
900 +    real( kind = dp ), intent(out) :: dudr
901 + !! Second Derivative, optional, used mainly for normal mode calculations.
902 +    real( kind = dp ), intent(out), optional :: d2
903 +    
904 +    type (lj_atype), pointer :: atype1
905 +    type (lj_atype), pointer :: atype2
906 +
907 +    integer, intent(out), optional :: status
908 +
909 + ! local Variables
910 +    real( kind = dp ) :: sigma
911 +    real( kind = dp ) :: sigma2
912 +    real( kind = dp ) :: sigma6
913 +    real( kind = dp ) :: epsilon
914 +
915 +    real( kind = dp ) :: rcut
916 +    real( kind = dp ) :: rcut2
917 +    real( kind = dp ) :: rcut6
918 +    real( kind = dp ) :: r2
919 +    real( kind = dp ) :: r6
920 +
921 +    real( kind = dp ) :: t6
922 +    real( kind = dp ) :: t12
923 +    real( kind = dp ) :: tp6
924 +    real( kind = dp ) :: tp12
925 +    real( kind = dp ) :: delta
926 +
927 +    logical :: doSec = .false.
928 +
929 +    integer :: errorStat
930 +
931 + !! Optional Argument Checking
932 + ! Check to see if we need to do second derivatives
933 +    
934 +    if (present(d2))     doSec = .true.
935 +    if (present(status)) status = 0
936 +
937 + ! Look up the correct parameters in the mixing matrix
938 +    sigma    = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma
939 +    sigma2   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2
940 +    sigma6   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6
941 +    epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
942 +
943 +
944 +
945 +    call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
946 +    
947 +    r2 = r * r
948 +    r6 = r2 * r2 * r2
949 +
950 +    t6  = sigma6/ r6
951 +    t12 = t6 * t6
952 +
953 +
954 +                                                                              
955 +    tp6 = sigma6 / rcut6
956 +    tp12 = tp6*tp6
957 +                                                                              
958 +    delta = -4.0E0_DP*epsilon * (tp12 - tp6)
959 +                                                                              
960 +    if (r.le.rcut) then
961 +       pot = 4.0E0_DP * epsilon * (t12 - t6) + delta
962 +       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / r
963 +       if(doSec)  d2 = 24.0E0_DP * epsilon * (26.0E0_DP*t12 - 7.0E0_DP*t6)/r/r
964 +    else
965 +       pot = 0.0E0_DP
966 +       dudr = 0.0E0_DP
967 +       if(doSec) d2 = 0.0E0_DP
968 +    endif
969 +                                                                              
970 +  return
971 +
972 +
973 +
974 +  end subroutine getLjPot
975 +
976 +
977 + !! Calculates the mixing for sigma or epslon based on character string for initialzition of mixing array.
978 +  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
979 +    character(len=*) :: thisParam
980 +    real(kind = dp)  :: param1
981 +    real(kind = dp)  :: param2
982 +    real(kind = dp ) :: myMixParam
983 +    integer, optional :: status
984 +
985 +
986 +    myMixParam = 0.0_dp
987 +
988 +    if (present(status)) status = 0
989 +
990 +    select case (thisParam)
991 +
992 +    case ("sigma")
993 +       myMixParam = 0.5_dp * (param1 + param2)
994 +    case ("epsilon")
995 +       myMixParam = sqrt(param1 * param2)
996 +    case default
997 +       status = -1
998 +    end select
999 +
1000 +  end function calcLJMix
1001 +
1002 +
1003 +
1004   end module lj_ff

Diff Legend

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