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 235 by mmeineke, Fri Jan 10 21:56:10 2003 UTC vs.
Revision 239 by chuckv, Mon Jan 20 22:36:12 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.8 2003-01-20 22:36:12 chuckv Exp $, $Date: 2003-01-20 22:36:12 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
6 +
7 +
8 +
9   module lj_ff
10    use simulation
11    use definitions, ONLY : dp, ndim
# Line 13 | Line 21 | module lj_ff
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
# Line 48 | Line 57 | module lj_ff
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
# Line 59 | Line 72 | module lj_ff
72  
73    real( kind = dp ), allocatable, dimension(:,:) :: fRow
74    real( kind = dp ), allocatable, dimension(:,:) :: fColumn
75 +
76 +  type (lj_atypePtr), dimension(:), allocatable :: identPtrListRow
77 +  type (lj_atypePtr), dimension(:), allocatable :: identPtrListColumn
78   #endif
79  
80  
81  
82 +  logical :: isljFFinit = .false.
83  
84  
68 –
85   !! Public methods and data
86    public :: new_lj_atype
87    public :: do_lj_ff
# Line 211 | Line 227 | contains
227    end subroutine new_lj_atype
228  
229  
230 <  subroutine init_ljFF()
230 >  subroutine init_ljFF(nComponents,ident, status)
231 > !! Number of components in ident array
232 >    integer, intent(inout) :: nComponents
233 > !! Array of identities nComponents long corresponding to
234 > !! ljatype ident.
235 >    integer, dimension(nComponents),intent(inout) :: ident
236 > !!  Result status, success = 0, error = -1
237 >    integer, intent(out) :: Status
238  
239 <
239 >    integer :: thisStat
240   #ifdef IS_MPI
241 +    integer, allocatable, dimension(:) :: identRow
242 +    integer, allocatable, dimension(:) :: identCol
243 +    integer :: nrow
244 +    integer :: ncol
245 +    integer :: alloc_stat
246 + #endif
247 +    status = 0
248  
249 + !! if were're not in MPI, we just update ljatypePtrList
250 + #ifndef IS_MPI
251 +    call new_ljatypePtrList(nComponents,ident,identPtrList,thisStat)
252 +    if ( thisStat /= 0 ) then
253 +       status = -1
254 +       return
255 +    endif
256 + !! Allocate pointer lists
257 +    allocate(point(nComponents),stat=alloc_stat)
258 +    if (alloc_stat /=0) then
259 +       status = -1
260 +       return
261 +    endif
262 +
263 +    allocate(list(nComponents*listMultiplier),stat=alloc_stat)
264 +    if (alloc_stat /=0) then
265 +       status = -1
266 +       return
267 +    endif
268      
269 + ! if were're in MPI, we also have to worry about row and col lists    
270 + #else
271 + ! We can only set up forces if mpiSimulation has been setup.
272 +    if (.not. isMPISimSet()) then
273 +       status = -1
274 +       return
275 +    endif
276 +    nrow = getNrow()
277 +    ncol = getNcol()
278 + !! Allocate temperary arrays to hold gather information
279 +    allocate(identRow(nrow),stat=alloc_stat)
280 +    if (alloc_stat /= 0 ) then
281 +       status = -1
282 +       return
283 +    endif
284 +
285 +    allocate(identCol(ncol),stat=alloc_stat)
286 +    if (alloc_stat /= 0 ) then
287 +       status = -1
288 +       return
289 +    endif
290 + !! Gather idents into row and column idents
291 +    call gather(ident,identRow,plan_row)
292 +    call gather(ident,identCol,plan_col)
293 +
294 + !! Create row and col pointer lists
295 +    call new_ljatypePtrList(nrow,identRow,identPtrListRow,thisStat)
296 +    if (thisStat /= 0 ) then
297 +       status = -1
298 +       return
299 +    endif
300 +
301 +    call new_ljatypePtrList(ncol,identCol,identPtrListCol,thisStat)
302 +    if (thisStat /= 0 ) then
303 +       status = -1
304 +       return
305 +    endif
306 +
307 + !! free temporary ident arrays
308 +    deallocate(identCol)
309 +    deallocate(identRow)
310 +
311 + !! Allocate Simulation arrays
312 + !! NOTE: This bit of code should be fixed, it can cause large
313 + !! memory fragmentation if call repeatedly
314 +
315 +    if (.not.allocated(qRow)) then
316 +       allocate(qRow(3,nrow),stat=alloc_stat)
317 +       if (alloc_stat /= 0 ) then
318 +          status = -1
319 +          return
320 +       endif
321 +    else
322 +       deallocate(qrow)
323 +       allocate(qRow(3,nrow),stat=alloc_stat)
324 +       if (alloc_stat /= 0 ) then
325 +          status = -1
326 +          return
327 +       endif
328 +    endif
329 +
330 +    if (.not.allocated(3,qCol)) then
331 +       allocate(qCol(ncol),stat=alloc_stat)
332 +       if (alloc_stat /= 0 ) then
333 +          status = -1
334 +          return
335 +       endif
336 +    else
337 +       deallocate(qCol)
338 +       allocate(qCol(3,ncol),stat=alloc_stat)
339 +       if (alloc_stat /= 0 ) then
340 +          status = -1
341 +          return
342 +       endif
343 +    endif
344 +
345 +    if (.not.allocated(fRow)) then
346 +       allocate(fRow(3,nrow),stat=alloc_stat)
347 +       if (alloc_stat /= 0 ) then
348 +          status = -1
349 +          return
350 +       endif
351 +    else
352 +       deallocate(fRow)
353 +       allocate(fRow(3,nrow),stat=alloc_stat)
354 +       if (alloc_stat /= 0 ) then
355 +          status = -1
356 +          return
357 +       endif
358 +    endif
359 +
360 +    if (.not.allocated(fCol)) then
361 +       allocate(fCol(3,ncol),stat=alloc_stat)
362 +       if (alloc_stat /= 0 ) then
363 +          status = -1
364 +          return
365 +       endif
366 +    else
367 +       deallocate(fCol)
368 +       allocate(fCol(3,ncol),stat=alloc_stat)
369 +       if (alloc_stat /= 0 ) then
370 +          status = -1
371 +          return
372 +       endif
373 +    endif
374 + !! Allocate neighbor lists for mpi simulations.
375 +    if (.not. allocated(point)) then
376 +       allocate(point(nrow),stat=alloc_stat)
377 +       if (alloc_stat /=0) then
378 +          status = -1
379 +          return
380 +       endif
381  
382 +       allocate(list(ncol*listMultiplier),stat=alloc_stat)
383 +       if (alloc_stat /=0) then
384 +          status = -1
385 +          return
386 +       endif
387 +    else
388 +       deallocate(list)
389 +       deallocate(point)
390 +       allocate(point(nrow),stat=alloc_stat)
391 +       if (alloc_stat /=0) then
392 +          status = -1
393 +          return
394 +       endif
395 +
396 +       allocate(list(ncol*listMultiplier),stat=alloc_stat)
397 +       if (alloc_stat /=0) then
398 +          status = -1
399 +          return
400 +       endif
401 +    endif
402 +
403   #endif
404 +    
405  
406 +    isljFFinit = .true.
407 +
408 +
409    end subroutine init_ljFF
410  
411 +
412 +
413 +
414 +
415   !! Takes an ident array and creates an atype pointer list
416   !! based on those identities
417    subroutine new_ljatypePtrList(mysize,ident,PtrList,status)
# Line 286 | Line 476 | contains
476    end subroutine getLJatype
477  
478  
479 < ! FORCE routine
479 > !! FORCE routine Calculates Lennard Jones forces.
480   !------------------------------------------------------------->
481    subroutine do_lj_ff(q,f,potE,do_pot)
482      real ( kind = dp ), dimension(ndim,) :: q
483      real ( kind = dp ), dimension(ndim,nLRparticles) :: f
484      real ( kind = dp ) :: potE
485      logical ( kind = 2) :: do_pot
486 +    
487 +    type(lj_atype), pointer :: ljAtype_i
488 +    type(lj_atype), pointer :: ljAtype_j
489 +
490   #ifdef MPI
491    real( kind = DP ), dimension(3,ncol) :: efr
492    real( kind = DP ) :: pot_local
# Line 314 | Line 508 | contains
508    integer :: nrow
509    integer :: ncol
510  
511 +  if (.not. isljFFInit) then
512 +     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
513 +     return
514 +  endif
515  
318 –
516   #ifndef IS_MPI
517    nrow = natoms - 1
518    ncol = natoms
519   #else
520 +  nrow = getNrow(plan_row)
521 +  ncol = getNcol(plan_col)
522    j_start = 1
523   #endif
524  
525 +  
526 +
527    call check(update_nlist)
528  
529 + !--------------WARNING...........................
530 + ! Zero variables, NOTE:::: Forces are zeroed in C
531 + ! Zeroing them here could delete previously computed
532 + ! Forces.
533 + !------------------------------------------------
534   #ifndef IS_MPI
535    nloops = nloops + 1
536    pot = 0.0E0_DP
331 –  f = 0.0E0_DP
537    e = 0.0E0_DP
538   #else
539      f_row = 0.0E0_DP
# Line 344 | Line 549 | contains
549  
550   ! communicate MPI positions
551   #ifdef MPI    
552 <    call gather(q,q_row,plan_row3)
553 <    call gather(q,q_col,plan_col3)
552 >    call gather(q,qRow,plan_row3)
553 >    call gather(q,qCol,plan_col3)
554   #endif
555  
556   #ifndef MPI
# Line 365 | Line 570 | contains
570       do i = 1, nrow
571          point(i) = nlist + 1
572   #ifdef MPI
573 <        tag_i = tag_row(i)
574 <        rxi = q_row(1,i)
575 <        ryi = q_row(2,i)
576 <        rzi = q_row(3,i)
573 >        ljAtype_i => identPtrListRow(i)%this
574 >        tag_i = tagRow(i)
575 >        rxi = qRow(1,i)
576 >        ryi = qRow(2,i)
577 >        rzi = qRow(3,i)
578   #else
579 +        ljAtype_i   => identPtrList(i)%this
580          j_start = i + 1
581          rxi = q(1,i)
582          ryi = q(2,i)
# Line 378 | Line 585 | contains
585  
586          inner: do j = j_start, ncol
587   #ifdef MPI
588 <           tag_j = tag_col(j)
588 > ! Assign identity pointers and tags
589 >           ljAtype_j => identPtrListColumn(j)%this
590 >           tag_j = tagCol(j)
591             if (newtons_thrd) then
592                if (tag_i <= tag_j) then
593                   if (mod(tag_i + tag_j,2) == 0) cycle inner
# Line 387 | Line 596 | contains
596                endif
597             endif
598  
599 <           rxij = wrap(rxi - q_col(1,j), 1)
600 <           ryij = wrap(ryi - q_col(2,j), 2)
601 <           rzij = wrap(rzi - q_col(3,j), 3)
599 >           rxij = wrap(rxi - qCol(1,j), 1)
600 >           ryij = wrap(ryi - qCol(2,j), 2)
601 >           rzij = wrap(rzi - qCol(3,j), 3)
602   #else          
603 <        
603 >           ljAtype_j   => identPtrList(j)%this
604             rxij = wrap(rxi - q(1,j), 1)
605             ryij = wrap(ryi - q(2,j), 2)
606             rzij = wrap(rzi - q(3,j), 3)
# Line 408 | Line 617 | contains
617              
618                nlist = nlist + 1
619                if (nlist > size(list)) then
620 <                 call info("FORCE_LJ","error size list smaller then nlist")
621 <                 write(msg,*) "nlist size(list)", nlist, size(list)
413 <                 call info("FORCE_LJ",msg)
414 < #ifdef MPI
415 <                 call mpi_abort(MPI_COMM_WORLD,mpi_err)
416 < #endif
417 <                 stop
620 > #warning "Change how nlist size is done"
621 >                 write(DEFAULT_ERROR,*) "ERROR: nlist > list size"
622                endif
623                list(nlist) = j
624  
# Line 423 | Line 627 | contains
627                  
628                   r = dsqrt(rijsq)
629        
630 <                 call LJ_mix(r,pot,dudr,d2,i,j)
630 >                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
631        
632   #ifdef MPI
633                  e_row(i) = e_row(i) + pot*0.5
# Line 444 | Line 648 | contains
648  
649  
650   #ifdef MPI
651 <                    f_col(dim,j) = f_col(dim,j) - ftmp
652 <                    f_row(dim,i) = f_row(dim,i) + ftmp
651 >                    fCol(dim,j) = fCol(dim,j) - ftmp
652 >                    fRow(dim,i) = fRow(dim,i) + ftmp
653   #else                    
654              
655                      f(dim,j) = f(dim,j) - ftmp
# Line 473 | Line 677 | contains
677          ! check thiat molecule i has neighbors
678          if (jbeg .le. jend) then
679   #ifdef MPI
680 <           rxi = q_row(1,i)
681 <           ryi = q_row(2,i)
682 <           rzi = q_row(3,i)
680 >           ljAtype_i => identPtrListRow(i)%this
681 >           rxi = qRow(1,i)
682 >           ryi = qRow(2,i)
683 >           rzi = qRow(3,i)
684   #else
685 +           ljAtype_i   => identPtrList(i)%this
686             rxi = q(1,i)
687             ryi = q(2,i)
688             rzi = q(3,i)
# Line 484 | Line 690 | contains
690             do jnab = jbeg, jend
691                j = list(jnab)
692   #ifdef MPI
693 <                rxij = wrap(rxi - q_col(1,j), 1)
694 <                ryij = wrap(ryi - q_col(2,j), 2)
695 <                rzij = wrap(rzi - q_col(3,j), 3)
693 >              ljAtype_j = identPtrListColumn(j)%this
694 >              rxij = wrap(rxi - q_col(1,j), 1)
695 >              ryij = wrap(ryi - q_col(2,j), 2)
696 >              rzij = wrap(rzi - q_col(3,j), 3)
697   #else
698 <                rxij = wrap(rxi - q(1,j), 1)
699 <                ryij = wrap(ryi - q(2,j), 2)
700 <                rzij = wrap(rzi - q(3,j), 3)
698 >              ljAtype_j = identPtrList(j)%this
699 >              rxij = wrap(rxi - q(1,j), 1)
700 >              ryij = wrap(ryi - q(2,j), 2)
701 >              rzij = wrap(rzi - q(3,j), 3)
702   #endif
703                rijsq = rxij*rxij + ryij*ryij + rzij*rzij
704                
# Line 498 | Line 706 | contains
706  
707                   r = dsqrt(rijsq)
708                  
709 <                 call LJ_mix(r,pot,dudr,d2,i,j)
709 >                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
710   #ifdef MPI
711                  e_row(i) = e_row(i) + pot*0.5
712                  e_col(i) = e_col(i) + pot*0.5
# Line 516 | Line 724 | contains
724                      drdx1 = efr(dim,j) / r
725                      ftmp = dudr * drdx1
726   #ifdef MPI
727 <                    f_col(dim,j) = f_col(dim,j) - ftmp
728 <                    f_row(dim,i) = f_row(dim,i) + ftmp
727 >                    fCol(dim,j) = fCol(dim,j) - ftmp
728 >                    fRow(dim,i) = fRow(dim,i) + ftmp
729   #else                    
730                      f(dim,j) = f(dim,j) - ftmp
731                      f(dim,i) = f(dim,i) + ftmp
# Line 533 | Line 741 | contains
741  
742   #ifdef MPI
743      !!distribute forces
744 <    call scatter(f_row,f,plan_row3)
744 >    call scatter(fRow,f,plan_row3)
745  
746 <    call scatter(f_col,f_tmp,plan_col3)
746 >    call scatter(fCol,f_tmp,plan_col3)
747      do i = 1,nlocal
748         do dim = 1,3
749            f(dim,i) = f(dim,i) + f_tmp(dim,i)
# Line 568 | Line 776 | contains
776  
777    end subroutine do_lj_ff
778  
779 < !! Calculates the potential between two lj particles, optionally returns second
779 > !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
780   !! derivatives.
781    subroutine getLjPot(r,pot,dudr,atype1,atype2,d2,status)
782   ! arguments
# Line 581 | Line 789 | contains
789   !! Second Derivative, optional, used mainly for normal mode calculations.
790      real( kind = dp ), intent(out), optional :: d2
791      
792 <    type (lj_atype), intent(in) :: atype1
793 <    type (lj_atype), intent(in) :: atype2
792 >    type (lj_atype), intent(in), pointer :: atype1
793 >    type (lj_atype), intent(in), pointer :: atype2
794  
795      integer, intent(out), optional :: status
796  
# Line 654 | Line 862 | contains
862    end subroutine getLjPot
863  
864  
865 < !! Calculates the mixing for sigma or epslon based on character string
865 > !! Calculates the mixing for sigma or epslon based on character string for initialzition of mixing array.
866    function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
867      character(len=*) :: thisParam
868      real(kind = dp)  :: param1

Diff Legend

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