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 230 by chuckv, Thu Jan 9 19:40:38 2003 UTC vs.
Revision 240 by chuckv, Wed Jan 22 21:45:20 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.9 2003-01-22 21:45:20 chuckv Exp $, $Date: 2003-01-22 21:45:20 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $
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
28   !! Unique number for place in linked list
29       integer :: atype_number = 0
21 – !! Name of atype
22 –     character(len = 15) :: name = "NULL"
30   !! Unique indentifier number (ie atomic no, etc)
31       integer :: atype_ident = 0
32   !! Mass of Particle
# Line 37 | Line 44 | module lj_ff
44    end type lj_atype
45  
46   !! Pointer type for atype ident array
47 <  type :: lj_atypePtr
47 >  type, public :: lj_atypePtr
48       type (lj_atype), pointer :: this => null()
49    end type lj_atypePtr
50  
# Line 50 | 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 61 | 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  
70 –
85   !! Public methods and data
86    public :: new_lj_atype
87    public :: do_lj_ff
# Line 79 | Line 93 | contains
93  
94   contains
95  
96 <  subroutine new_lj_atype(name,ident,mass,epslon,sigma,status)
83 <    character( len = 15 ) :: name
96 >  subroutine new_lj_atype(ident,mass,epslon,sigma,status)
97      real( kind = dp ), intent(in) :: mass
98      real( kind = dp ), intent(in) :: epslon
99      real( kind = dp ), intent(in) :: sigma
# Line 108 | Line 121 | contains
121      end if
122  
123   ! assign our new lj_atype information
111 –    this_lj_atype%name       = name
124      this_lj_atype%mass       = mass
125      this_lj_atype%epslon     = epslon
126      this_lj_atype%sigma      = sigma
# Line 215 | 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 290 | 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
# Line 306 | Line 496 | contains
496    
497    real( kind = DP )   :: pe
498    logical,            :: update_nlist
309 –  logical             :: do_pot
499  
500 +
501    integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
502    integer :: nlist
503    integer :: j_start
504    integer :: tag_i,tag_j
505    real( kind = DP ) ::  r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
506    real( kind = DP ) ::  rxi, ryi, rzi, rxij, ryij, rzij, rijsq
507 +  real( kind = DP ) ::  rlistsq, rcutsq
508  
509    integer :: nrow
510    integer :: ncol
511  
512 +  if (.not. isljFFInit) then
513 +     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
514 +     return
515 +  endif
516  
322 –
517   #ifndef IS_MPI
518    nrow = natoms - 1
519    ncol = natoms
520   #else
521 +  nrow = getNrow(plan_row)
522 +  ncol = getNcol(plan_col)
523    j_start = 1
524   #endif
525  
526 +  
527 +
528    call check(update_nlist)
529  
530 <  do_pot = .false.
531 <  if (present(pe)) do_pot = .true.
532 <
530 > !--------------WARNING...........................
531 > ! Zero variables, NOTE:::: Forces are zeroed in C
532 > ! Zeroing them here could delete previously computed
533 > ! Forces.
534 > !------------------------------------------------
535   #ifndef IS_MPI
536    nloops = nloops + 1
537    pot = 0.0E0_DP
338 –  f = 0.0E0_DP
538    e = 0.0E0_DP
539   #else
540      f_row = 0.0E0_DP
# Line 351 | Line 550 | contains
550  
551   ! communicate MPI positions
552   #ifdef MPI    
553 <    call gather(q,q_row,plan_row3)
554 <    call gather(q,q_col,plan_col3)
553 >    call gather(q,qRow,plan_row3)
554 >    call gather(q,qCol,plan_col3)
555   #endif
556  
557   #ifndef MPI
# Line 372 | Line 571 | contains
571       do i = 1, nrow
572          point(i) = nlist + 1
573   #ifdef MPI
574 <        tag_i = tag_row(i)
575 <        rxi = q_row(1,i)
576 <        ryi = q_row(2,i)
577 <        rzi = q_row(3,i)
574 >        ljAtype_i => identPtrListRow(i)%this
575 >        tag_i = tagRow(i)
576 >        rxi = qRow(1,i)
577 >        ryi = qRow(2,i)
578 >        rzi = qRow(3,i)
579   #else
580 +        ljAtype_i   => identPtrList(i)%this
581          j_start = i + 1
582          rxi = q(1,i)
583          ryi = q(2,i)
# Line 385 | Line 586 | contains
586  
587          inner: do j = j_start, ncol
588   #ifdef MPI
589 <           tag_j = tag_col(j)
589 > ! Assign identity pointers and tags
590 >           ljAtype_j => identPtrListColumn(j)%this
591 >           tag_j = tagCol(j)
592             if (newtons_thrd) then
593                if (tag_i <= tag_j) then
594                   if (mod(tag_i + tag_j,2) == 0) cycle inner
# Line 394 | Line 597 | contains
597                endif
598             endif
599  
600 <           rxij = wrap(rxi - q_col(1,j), 1)
601 <           ryij = wrap(ryi - q_col(2,j), 2)
602 <           rzij = wrap(rzi - q_col(3,j), 3)
600 >           rxij = wrap(rxi - qCol(1,j), 1)
601 >           ryij = wrap(ryi - qCol(2,j), 2)
602 >           rzij = wrap(rzi - qCol(3,j), 3)
603   #else          
604 <        
604 >           ljAtype_j   => identPtrList(j)%this
605             rxij = wrap(rxi - q(1,j), 1)
606             ryij = wrap(ryi - q(2,j), 2)
607             rzij = wrap(rzi - q(3,j), 3)
# Line 415 | Line 618 | contains
618              
619                nlist = nlist + 1
620                if (nlist > size(list)) then
621 <                 call info("FORCE_LJ","error size list smaller then nlist")
622 <                 write(msg,*) "nlist size(list)", nlist, size(list)
420 <                 call info("FORCE_LJ",msg)
421 < #ifdef MPI
422 <                 call mpi_abort(MPI_COMM_WORLD,mpi_err)
423 < #endif
424 <                 stop
621 > #warning "Change how nlist size is done"
622 >                 write(DEFAULT_ERROR,*) "ERROR: nlist > list size"
623                endif
624                list(nlist) = j
625  
# Line 430 | Line 628 | contains
628                  
629                   r = dsqrt(rijsq)
630        
631 <                 call LJ_mix(r,pot,dudr,d2,i,j)
631 >                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
632        
633   #ifdef MPI
634                  e_row(i) = e_row(i) + pot*0.5
# Line 451 | Line 649 | contains
649  
650  
651   #ifdef MPI
652 <                    f_col(dim,j) = f_col(dim,j) - ftmp
653 <                    f_row(dim,i) = f_row(dim,i) + ftmp
652 >                    fCol(dim,j) = fCol(dim,j) - ftmp
653 >                    fRow(dim,i) = fRow(dim,i) + ftmp
654   #else                    
655              
656                      f(dim,j) = f(dim,j) - ftmp
# Line 480 | Line 678 | contains
678          ! check thiat molecule i has neighbors
679          if (jbeg .le. jend) then
680   #ifdef MPI
681 <           rxi = q_row(1,i)
682 <           ryi = q_row(2,i)
683 <           rzi = q_row(3,i)
681 >           ljAtype_i => identPtrListRow(i)%this
682 >           rxi = qRow(1,i)
683 >           ryi = qRow(2,i)
684 >           rzi = qRow(3,i)
685   #else
686 +           ljAtype_i   => identPtrList(i)%this
687             rxi = q(1,i)
688             ryi = q(2,i)
689             rzi = q(3,i)
# Line 491 | Line 691 | contains
691             do jnab = jbeg, jend
692                j = list(jnab)
693   #ifdef MPI
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)
694 >              ljAtype_j = identPtrListColumn(j)%this
695 >              rxij = wrap(rxi - q_col(1,j), 1)
696 >              ryij = wrap(ryi - q_col(2,j), 2)
697 >              rzij = wrap(rzi - q_col(3,j), 3)
698   #else
699 <                rxij = wrap(rxi - q(1,j), 1)
700 <                ryij = wrap(ryi - q(2,j), 2)
701 <                rzij = wrap(rzi - q(3,j), 3)
699 >              ljAtype_j = identPtrList(j)%this
700 >              rxij = wrap(rxi - q(1,j), 1)
701 >              ryij = wrap(ryi - q(2,j), 2)
702 >              rzij = wrap(rzi - q(3,j), 3)
703   #endif
704                rijsq = rxij*rxij + ryij*ryij + rzij*rzij
705                
# Line 505 | Line 707 | contains
707  
708                   r = dsqrt(rijsq)
709                  
710 <                 call LJ_mix(r,pot,dudr,d2,i,j)
710 >                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
711   #ifdef MPI
712                  e_row(i) = e_row(i) + pot*0.5
713                  e_col(i) = e_col(i) + pot*0.5
# Line 523 | Line 725 | contains
725                      drdx1 = efr(dim,j) / r
726                      ftmp = dudr * drdx1
727   #ifdef MPI
728 <                    f_col(dim,j) = f_col(dim,j) - ftmp
729 <                    f_row(dim,i) = f_row(dim,i) + ftmp
728 >                    fCol(dim,j) = fCol(dim,j) - ftmp
729 >                    fRow(dim,i) = fRow(dim,i) + ftmp
730   #else                    
731                      f(dim,j) = f(dim,j) - ftmp
732                      f(dim,i) = f(dim,i) + ftmp
# Line 540 | Line 742 | contains
742  
743   #ifdef MPI
744      !!distribute forces
745 <    call scatter(f_row,f,plan_row3)
745 >    call scatter(fRow,f,plan_row3)
746  
747 <    call scatter(f_col,f_tmp,plan_col3)
747 >    call scatter(fCol,f_tmp,plan_col3)
748      do i = 1,nlocal
749         do dim = 1,3
750            f(dim,i) = f(dim,i) + f_tmp(dim,i)
# Line 553 | Line 755 | contains
755      
756      if (do_pot) then
757   ! scatter/gather pot_row into the members of my column
758 <  call scatter(e_row,e_tmp,plan_row)
758 >       call scatter(e_row,e_tmp,plan_row)
759        
760         ! scatter/gather pot_local into all other procs
761         ! add resultant to get total pot
# Line 575 | Line 777 | contains
777  
778    end subroutine do_lj_ff
779  
780 < !! Calculates the potential between two lj particles, optionally returns second
780 > !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
781   !! derivatives.
782    subroutine getLjPot(r,pot,dudr,atype1,atype2,d2,status)
783   ! arguments
# Line 588 | Line 790 | contains
790   !! Second Derivative, optional, used mainly for normal mode calculations.
791      real( kind = dp ), intent(out), optional :: d2
792      
793 <    type (lj_atype), intent(in) :: atype1
794 <    type (lj_atype), intent(in) :: atype2
793 >    type (lj_atype), intent(in), pointer :: atype1
794 >    type (lj_atype), intent(in), pointer :: atype2
795  
796      integer, intent(out), optional :: status
797  
# Line 661 | Line 863 | contains
863    end subroutine getLjPot
864  
865  
866 < !! Calculates the mixing for sigma or epslon based on character string
866 > !! Calculates the mixing for sigma or epslon based on character string for initialzition of mixing array.
867    function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
868      character(len=*) :: thisParam
869      real(kind = dp)  :: param1

Diff Legend

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