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 258 by chuckv, Thu Jan 30 22:29: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.15 2003-01-30 22:29:37 chuckv Exp $, $Date: 2003-01-30 22:29:37 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $
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 >
171   !! Allocate temperary arrays to hold gather information
172      allocate(identRow(nrow),stat=alloc_stat)
173      if (alloc_stat /= 0 ) then
# Line 173 | Line 180 | contains
180         status = -1
181         return
182      endif
183 +
184   !! Gather idents into row and column idents
185 +
186      call gather(ident,identRow,plan_row)
187      call gather(ident,identCol,plan_col)
188 <
188 >    
189 >
190   !! Create row and col pointer lists
191 <    call create_IdentPtrlst(identRow,ljListTail,identPtrListRow,thisStat)
191 >
192 >    call create_IdentPtrlst(identRow,ljListHead,identPtrListRow,thisStat)
193      if (thisStat /= 0 ) then
194         status = -1
195         return
196      endif
197  
198 <    call create_IdentPtrlst(identCol,ljListTail,identPtrListCol,thisStat)
198 >    call create_IdentPtrlst(identCol,ljListHead,identPtrListColumn,thisStat)
199      if (thisStat /= 0 ) then
200         status = -1
201         return
202      endif
203 <
203 >    write(*,*) "After createIdent row-col lists"
204   !! free temporary ident arrays
205 <    deallocate(identCol)
206 <    deallocate(identRow)
207 <
205 >    if (allocated(identCol)) then
206 >       deallocate(identCol)
207 >    end if
208 >    if (allocated(identCol)) then
209 >       deallocate(identRow)
210 >    endif
211  
212   !! Allocate neighbor lists for mpi simulations.
213      if (.not. allocated(point)) then
# Line 342 | Line 356 | contains
356      type(lj_atype), pointer :: ljAtype_i
357      type(lj_atype), pointer :: ljAtype_j
358  
359 < #ifdef MPI
360 <  real( kind = DP ), dimension(3,getNcol()) :: efr
359 > #ifdef IS_MPI
360 >  real( kind = DP ), dimension(3,getNcol(plan_col)) :: efr
361    real( kind = DP ) :: pot_local
362   #else
363    real( kind = DP ), dimension(3,getNlocal()) :: efr
# Line 351 | Line 365 | contains
365    
366   !! Local arrays needed for MPI
367   #ifdef IS_MPI
368 <  real(kind = dp), dimension(3:getNrow()) :: qRow
369 <  real(kind = dp), dimension(3:getNcol()) :: qCol
368 >  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow
369 >  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol
370  
371 <  real(kind = dp), dimension(3:getNrow()) :: fRow
372 <  real(kind = dp), dimension(3:getNcol()) :: fCol
371 >  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow
372 >  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol
373 >  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp
374  
375 <  real(kind = dp), dimension(getNrow()) :: eRow
376 <  real(kind = dp), dimension(getNcol()) :: eCol
375 >  real(kind = dp), dimension(getNrow(plan_row)) :: eRow
376 >  real(kind = dp), dimension(getNcol(plan_col)) :: eCol
377  
378    real(kind = dp), dimension(getNlocal()) :: eTemp
379   #endif
# Line 387 | Line 402 | contains
402    integer :: ncol
403    integer :: natoms
404  
405 +
406 +
407 +
408   ! Make sure we are properly initialized.
409    if (.not. isljFFInit) then
410       write(default_error,*) "ERROR: lj_FF has not been properly initialized"
# Line 403 | Line 421 | contains
421    natoms = getNlocal()
422    call getRcut(rcut,rcut2=rcutsq)
423    call getRlist(rlist,rlistsq)
424 +
425   #ifndef IS_MPI
426    nrow = natoms - 1
427    ncol = natoms
# Line 416 | Line 435 | contains
435    
436   !! See if we need to update neighbor lists
437    call check(q,update_nlist)
438 +  if (firstTime) then
439 +     update_nlist = .true.
440 +     firstTime = .false.
441 +  endif
442  
443   !--------------WARNING...........................
444   ! Zero variables, NOTE:::: Forces are zeroed in C
# Line 427 | Line 450 | contains
450    pe = 0.0E0_DP
451  
452   #else
453 <    f_row = 0.0E0_DP
454 <    f_col = 0.0E0_DP
453 >    fRow = 0.0E0_DP
454 >    fCol = 0.0E0_DP
455  
456      pe_local = 0.0E0_DP
457  
# Line 439 | Line 462 | contains
462      efr = 0.0E0_DP
463  
464   ! communicate MPI positions
465 < #ifdef MPI    
466 <    call gather(q,qRow,plan_row3)
467 <    call gather(q,qCol,plan_col3)
465 > #ifdef IS_MPI    
466 >    call gather(q,qRow,plan_row3d)
467 >    call gather(q,qCol,plan_col3d)
468   #endif
469  
470  
# Line 453 | Line 476 | contains
476      
477       nlist = 0
478      
479 <    
479 >    
480  
481       do i = 1, nrow
482          point(i) = nlist + 1
483 < #ifdef MPI
483 > #ifdef IS_MPI
484          ljAtype_i => identPtrListRow(i)%this
485          tag_i = tagRow(i)
486          rxi = qRow(1,i)
# Line 472 | Line 495 | contains
495   #endif
496  
497          inner: do j = j_start, ncol
498 < #ifdef MPI
498 > #ifdef IS_MPI
499   ! Assign identity pointers and tags
500             ljAtype_j => identPtrListColumn(j)%this
501 <           tag_j = tagCol(j)
501 >           tag_j = tagColumn(j)
502             if (newtons_thrd) then
503                if (tag_i <= tag_j) then
504                   if (mod(tag_i + tag_j,2) == 0) cycle inner
# Line 496 | Line 519 | contains
519   #endif          
520             rijsq = rxij*rxij + ryij*ryij + rzij*rzij
521  
522 < #ifdef MPI
522 > #ifdef IS_MPI
523               if (rijsq <=  rlistsq .AND. &
524                    tag_j /= tag_i) then
525   #else
526 +          
527               if (rijsq <  rlistsq) then
528   #endif
529              
# Line 510 | Line 534 | contains
534                endif
535                list(nlist) = j
536  
537 <              
537 >    
538                if (rijsq <  rcutsq) then
539                  
540                   r = dsqrt(rijsq)
541        
542                   call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
543        
544 < #ifdef MPI
544 > #ifdef IS_MPI
545                  eRow(i) = eRow(i) + pot*0.5
546                  eCol(i) = eCol(i) + pot*0.5
547   #else
548 <                pe = pe + pot
548 >                    pe = pe + pot
549   #endif                
550              
551                   efr(1,j) = -rxij
# Line 535 | Line 559 | contains
559                      ftmp = dudr * drdx1
560  
561  
562 < #ifdef MPI
562 > #ifdef IS_MPI
563                      fCol(dim,j) = fCol(dim,j) - ftmp
564                      fRow(dim,i) = fRow(dim,i) + ftmp
565   #else                    
# Line 550 | Line 574 | contains
574          enddo inner
575       enddo
576  
577 < #ifdef MPI
577 > #ifdef IS_MPI
578       point(nrow + 1) = nlist + 1
579   #else
580       point(natoms) = nlist + 1
# Line 564 | Line 588 | contains
588          JEND = POINT(i+1) - 1
589          ! check thiat molecule i has neighbors
590          if (jbeg .le. jend) then
591 < #ifdef MPI
591 > #ifdef IS_MPI
592             ljAtype_i => identPtrListRow(i)%this
593             rxi = qRow(1,i)
594             ryi = qRow(2,i)
# Line 577 | Line 601 | contains
601   #endif
602             do jnab = jbeg, jend
603                j = list(jnab)
604 < #ifdef MPI
604 > #ifdef IS_MPI
605                ljAtype_j = identPtrListColumn(j)%this
606 <              rxij = wrap(rxi - q_col(1,j), 1)
607 <              ryij = wrap(ryi - q_col(2,j), 2)
608 <              rzij = wrap(rzi - q_col(3,j), 3)
606 >              rxij = wrap(rxi - qCol(1,j), 1)
607 >              ryij = wrap(ryi - qCol(2,j), 2)
608 >              rzij = wrap(rzi - qCol(3,j), 3)
609   #else
610                ljAtype_j = identPtrList(j)%this
611                rxij = wrap(rxi - q(1,j), 1)
# Line 595 | Line 619 | contains
619                   r = dsqrt(rijsq)
620                  
621                   call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
622 < #ifdef MPI
622 > #ifdef IS_MPI
623                  eRow(i) = eRow(i) + pot*0.5
624                  eCol(i) = eCol(i) + pot*0.5
625   #else
# Line 611 | Line 635 | contains
635                      
636                      drdx1 = efr(dim,j) / r
637                      ftmp = dudr * drdx1
638 < #ifdef MPI
638 > #ifdef IS_MPI
639                      fCol(dim,j) = fCol(dim,j) - ftmp
640                      fRow(dim,i) = fRow(dim,i) + ftmp
641   #else                    
# Line 627 | Line 651 | contains
651  
652  
653  
654 < #ifdef MPI
654 > #ifdef IS_MPI
655      !!distribute forces
656 <    call scatter(fRow,f,plan_row3)
656 >    call scatter(fRow,f,plan_row3d)
657  
658 <    call scatter(fCol,f_tmp,plan_col3)
658 >    call scatter(fCol,fMPITemp,plan_col3d)
659  
660      do i = 1,nlocal
661 <       f(1:3,i) = f(1:3,i) + f_tmp(1:3,i)
661 >       f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
662      end do
663  
664  
665      
666      if (do_pot) then
667   ! scatter/gather pot_row into the members of my column
668 <       call scatter(e_row,e_tmp,plan_row)
668 >       call scatter(eRow,eTemp,plan_row)
669        
670         ! scatter/gather pot_local into all other procs
671         ! add resultant to get total pot
672         do i = 1, nlocal
673 <          pe_local = pe_local + e_tmp(i)
673 >          pe_local = pe_local + eTemp(i)
674         enddo
675         if (newtons_thrd) then
676 <          e_tmp = 0.0E0_DP
677 <          call scatter(e_col,e_tmp,plan_col)
676 >          eTemp = 0.0E0_DP
677 >          call scatter(eCol,eTemp,plan_col)
678            do i = 1, nlocal
679 <             pe_local = pe_local + e_tmp(i)
679 >             pe_local = pe_local + eTemp(i)
680            enddo
681         endif
682      endif
# Line 661 | Line 685 | contains
685  
686      potE = pe
687  
688 +  
689  
690 +
691    end subroutine do_lj_ff
692  
693   !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
# Line 717 | Line 743 | contains
743      epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
744  
745  
746 <
746 >    
747 >
748      call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
749      
750      r2 = r * r

Diff Legend

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