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 252 by chuckv, Tue Jan 28 22:16:55 2003 UTC vs.
Revision 254 by chuckv, Thu Jan 30 20:03:37 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.12 2003-01-28 22:16:55 chuckv Exp $, $Date: 2003-01-28 22:16:55 $, $Name: not supported by cvs2svn $, $Revision: 1.12 $
5 > !! @version $Id: lj_FF.F90,v 1.14 2003-01-30 20:03:36 chuckv Exp $, $Date: 2003-01-30 20:03:36 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $
6  
7  
8  
# Line 36 | Line 36 | module lj_ff
36    integer :: nListAllocs = 0
37    integer, parameter :: maxListAllocs = 5
38  
39 +  logical, save :: firstTime = .True.
40  
41   !! Atype identity pointer lists
42   #ifdef IS_MPI
43   !! Row lj_atype pointer list
44 <  type (lj_atypePtr), dimension(:), pointer :: identPtrListRow => null()
44 >  type (lj_identPtrList), dimension(:), pointer :: identPtrListRow => null()
45   !! Column lj_atype pointer list
46 <  type (lj_atypePtr), dimension(:), pointer :: identPtrListColumn => null()
46 >  type (lj_identPtrList), dimension(:), pointer :: identPtrListColumn => null()
47   #else
48    type( lj_identPtrList ), dimension(:), pointer :: identPtrList => null()
49   #endif
# Line 128 | Line 129 | contains
129      integer, allocatable, dimension(:) :: identCol
130      integer :: nrow
131      integer :: ncol
131 –    integer :: alloc_stat
132   #endif
133      status = 0
134 +  
135 +
136 +
137 +
138   !! if were're not in MPI, we just update ljatypePtrList
139   #ifndef IS_MPI
140      call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat)
# Line 154 | Line 158 | contains
158      
159   ! if were're in MPI, we also have to worry about row and col lists    
160   #else
161 +  
162   ! We can only set up forces if mpiSimulation has been setup.
163      if (.not. isMPISimSet()) then
164 +       write(default_error,*) "MPI is not set"
165         status = -1
166         return
167      endif
168 <    nrow = getNrow()
169 <    ncol = getNcol()
168 >    nrow = getNrow(plan_row)
169 >    ncol = getNcol(plan_col)
170   !! Allocate temperary arrays to hold gather information
171      allocate(identRow(nrow),stat=alloc_stat)
172      if (alloc_stat /= 0 ) then
# Line 173 | Line 179 | contains
179         status = -1
180         return
181      endif
182 +
183   !! Gather idents into row and column idents
184 +
185      call gather(ident,identRow,plan_row)
186      call gather(ident,identCol,plan_col)
187 +
188  
189   !! Create row and col pointer lists
190 <    call create_IdentPtrlst(identRow,ljListTail,identPtrListRow,thisStat)
190 >
191 >    call create_IdentPtrlst(identRow,ljListHead,identPtrListRow,thisStat)
192      if (thisStat /= 0 ) then
193         status = -1
194         return
195      endif
196  
197 <    call create_IdentPtrlst(identCol,ljListTail,identPtrListCol,thisStat)
197 >    call create_IdentPtrlst(identCol,ljListHead,identPtrListColumn,thisStat)
198      if (thisStat /= 0 ) then
199         status = -1
200         return
# Line 342 | Line 352 | contains
352      type(lj_atype), pointer :: ljAtype_i
353      type(lj_atype), pointer :: ljAtype_j
354  
355 < #ifdef MPI
356 <  real( kind = DP ), dimension(3,getNcol()) :: efr
355 > #ifdef IS_MPI
356 >  real( kind = DP ), dimension(3,getNcol(plan_col)) :: efr
357    real( kind = DP ) :: pot_local
358   #else
359    real( kind = DP ), dimension(3,getNlocal()) :: efr
# Line 351 | Line 361 | contains
361    
362   !! Local arrays needed for MPI
363   #ifdef IS_MPI
364 <  real(kind = dp), dimension(3:getNrow()) :: qRow
365 <  real(kind = dp), dimension(3:getNcol()) :: qCol
364 >  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow
365 >  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol
366  
367 <  real(kind = dp), dimension(3:getNrow()) :: fRow
368 <  real(kind = dp), dimension(3:getNcol()) :: fCol
367 >  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow
368 >  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol
369 >  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp
370  
371 <  real(kind = dp), dimension(getNrow()) :: eRow
372 <  real(kind = dp), dimension(getNcol()) :: eCol
371 >  real(kind = dp), dimension(getNrow(plan_row)) :: eRow
372 >  real(kind = dp), dimension(getNcol(plan_col)) :: eCol
373  
374    real(kind = dp), dimension(getNlocal()) :: eTemp
375   #endif
# Line 387 | Line 398 | contains
398    integer :: ncol
399    integer :: natoms
400  
401 +
402 +
403 +
404   ! Make sure we are properly initialized.
405    if (.not. isljFFInit) then
406       write(default_error,*) "ERROR: lj_FF has not been properly initialized"
# Line 403 | Line 417 | contains
417    natoms = getNlocal()
418    call getRcut(rcut,rcut2=rcutsq)
419    call getRlist(rlist,rlistsq)
420 +
421   #ifndef IS_MPI
422    nrow = natoms - 1
423    ncol = natoms
# Line 416 | Line 431 | contains
431    
432   !! See if we need to update neighbor lists
433    call check(q,update_nlist)
434 +  if (firstTime) then
435 +     update_nlist = .true.
436 +     firstTime = .false.
437 +  endif
438  
439   !--------------WARNING...........................
440   ! Zero variables, NOTE:::: Forces are zeroed in C
# Line 427 | Line 446 | contains
446    pe = 0.0E0_DP
447  
448   #else
449 <    f_row = 0.0E0_DP
450 <    f_col = 0.0E0_DP
449 >    fRow = 0.0E0_DP
450 >    fCol = 0.0E0_DP
451  
452      pe_local = 0.0E0_DP
453  
# Line 439 | Line 458 | contains
458      efr = 0.0E0_DP
459  
460   ! communicate MPI positions
461 < #ifdef MPI    
462 <    call gather(q,qRow,plan_row3)
463 <    call gather(q,qCol,plan_col3)
461 > #ifdef IS_MPI    
462 >    call gather(q,qRow,plan_row3d)
463 >    call gather(q,qCol,plan_col3d)
464   #endif
465  
466  
# Line 453 | Line 472 | contains
472      
473       nlist = 0
474      
475 <    
475 >    
476  
477       do i = 1, nrow
478          point(i) = nlist + 1
479 < #ifdef MPI
479 > #ifdef IS_MPI
480          ljAtype_i => identPtrListRow(i)%this
481          tag_i = tagRow(i)
482          rxi = qRow(1,i)
# Line 472 | Line 491 | contains
491   #endif
492  
493          inner: do j = j_start, ncol
494 < #ifdef MPI
494 > #ifdef IS_MPI
495   ! Assign identity pointers and tags
496             ljAtype_j => identPtrListColumn(j)%this
497 <           tag_j = tagCol(j)
497 >           tag_j = tagColumn(j)
498             if (newtons_thrd) then
499                if (tag_i <= tag_j) then
500                   if (mod(tag_i + tag_j,2) == 0) cycle inner
# Line 496 | Line 515 | contains
515   #endif          
516             rijsq = rxij*rxij + ryij*ryij + rzij*rzij
517  
518 < #ifdef MPI
518 > #ifdef IS_MPI
519               if (rijsq <=  rlistsq .AND. &
520                    tag_j /= tag_i) then
521   #else
522 +          
523               if (rijsq <  rlistsq) then
524   #endif
525              
# Line 510 | Line 530 | contains
530                endif
531                list(nlist) = j
532  
533 <              
533 >    
534                if (rijsq <  rcutsq) then
535                  
536                   r = dsqrt(rijsq)
537        
538                   call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
539        
540 < #ifdef MPI
540 > #ifdef IS_MPI
541                  eRow(i) = eRow(i) + pot*0.5
542                  eCol(i) = eCol(i) + pot*0.5
543   #else
544 <                pe = pe + pot
544 >                    pe = pe + pot
545   #endif                
546              
547                   efr(1,j) = -rxij
# Line 535 | Line 555 | contains
555                      ftmp = dudr * drdx1
556  
557  
558 < #ifdef MPI
558 > #ifdef IS_MPI
559                      fCol(dim,j) = fCol(dim,j) - ftmp
560                      fRow(dim,i) = fRow(dim,i) + ftmp
561   #else                    
# Line 550 | Line 570 | contains
570          enddo inner
571       enddo
572  
573 < #ifdef MPI
573 > #ifdef IS_MPI
574       point(nrow + 1) = nlist + 1
575   #else
576       point(natoms) = nlist + 1
# Line 564 | Line 584 | contains
584          JEND = POINT(i+1) - 1
585          ! check thiat molecule i has neighbors
586          if (jbeg .le. jend) then
587 < #ifdef MPI
587 > #ifdef IS_MPI
588             ljAtype_i => identPtrListRow(i)%this
589             rxi = qRow(1,i)
590             ryi = qRow(2,i)
# Line 577 | Line 597 | contains
597   #endif
598             do jnab = jbeg, jend
599                j = list(jnab)
600 < #ifdef MPI
600 > #ifdef IS_MPI
601                ljAtype_j = identPtrListColumn(j)%this
602 <              rxij = wrap(rxi - q_col(1,j), 1)
603 <              ryij = wrap(ryi - q_col(2,j), 2)
604 <              rzij = wrap(rzi - q_col(3,j), 3)
602 >              rxij = wrap(rxi - qCol(1,j), 1)
603 >              ryij = wrap(ryi - qCol(2,j), 2)
604 >              rzij = wrap(rzi - qCol(3,j), 3)
605   #else
606                ljAtype_j = identPtrList(j)%this
607                rxij = wrap(rxi - q(1,j), 1)
# Line 595 | Line 615 | contains
615                   r = dsqrt(rijsq)
616                  
617                   call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
618 < #ifdef MPI
618 > #ifdef IS_MPI
619                  eRow(i) = eRow(i) + pot*0.5
620                  eCol(i) = eCol(i) + pot*0.5
621   #else
# Line 611 | Line 631 | contains
631                      
632                      drdx1 = efr(dim,j) / r
633                      ftmp = dudr * drdx1
634 < #ifdef MPI
634 > #ifdef IS_MPI
635                      fCol(dim,j) = fCol(dim,j) - ftmp
636                      fRow(dim,i) = fRow(dim,i) + ftmp
637   #else                    
# Line 627 | Line 647 | contains
647  
648  
649  
650 < #ifdef MPI
650 > #ifdef IS_MPI
651      !!distribute forces
652 <    call scatter(fRow,f,plan_row3)
652 >    call scatter(fRow,f,plan_row3d)
653  
654 <    call scatter(fCol,f_tmp,plan_col3)
654 >    call scatter(fCol,fMPITemp,plan_col3d)
655  
656      do i = 1,nlocal
657 <       f(1:3,i) = f(1:3,i) + f_tmp(1:3,i)
657 >       f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
658      end do
659  
660  
661      
662      if (do_pot) then
663   ! scatter/gather pot_row into the members of my column
664 <       call scatter(e_row,e_tmp,plan_row)
664 >       call scatter(eRow,eTemp,plan_row)
665        
666         ! scatter/gather pot_local into all other procs
667         ! add resultant to get total pot
668         do i = 1, nlocal
669 <          pe_local = pe_local + e_tmp(i)
669 >          pe_local = pe_local + eTemp(i)
670         enddo
671         if (newtons_thrd) then
672 <          e_tmp = 0.0E0_DP
673 <          call scatter(e_col,e_tmp,plan_col)
672 >          eTemp = 0.0E0_DP
673 >          call scatter(eCol,eTemp,plan_col)
674            do i = 1, nlocal
675 <             pe_local = pe_local + e_tmp(i)
675 >             pe_local = pe_local + eTemp(i)
676            enddo
677         endif
678      endif
# Line 661 | Line 681 | contains
681  
682      potE = pe
683  
684 +  
685  
686 +
687    end subroutine do_lj_ff
688  
689   !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
# Line 717 | Line 739 | contains
739      epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
740  
741  
742 <
742 >    
743 >
744      call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
745      
746      r2 = r * r

Diff Legend

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