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 261 by chuckv, Mon Feb 3 21:15:59 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.17 2003-02-03 21:15:59 chuckv Exp $, $Date: 2003-02-03 21:15:59 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $
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 123 | Line 124 | contains
124  
125      integer :: thisStat
126      integer :: i
127 +
128 +    integer :: myNode
129   #ifdef IS_MPI
130      integer, allocatable, dimension(:) :: identRow
131      integer, allocatable, dimension(:) :: identCol
132      integer :: nrow
133      integer :: ncol
131 –    integer :: alloc_stat
134   #endif
135      status = 0
136 +  
137 +
138 +    
139 +
140   !! if were're not in MPI, we just update ljatypePtrList
141   #ifndef IS_MPI
142      call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat)
# Line 154 | Line 160 | contains
160      
161   ! if were're in MPI, we also have to worry about row and col lists    
162   #else
163 +  
164   ! We can only set up forces if mpiSimulation has been setup.
165      if (.not. isMPISimSet()) then
166 +       write(default_error,*) "MPI is not set"
167         status = -1
168         return
169      endif
170 <    nrow = getNrow()
171 <    ncol = getNcol()
170 >    nrow = getNrow(plan_row)
171 >    ncol = getNcol(plan_col)
172 >    mynode = getMyNode()
173   !! Allocate temperary arrays to hold gather information
174      allocate(identRow(nrow),stat=alloc_stat)
175      if (alloc_stat /= 0 ) then
# Line 173 | Line 182 | contains
182         status = -1
183         return
184      endif
185 +
186   !! Gather idents into row and column idents
187 +
188      call gather(ident,identRow,plan_row)
189      call gather(ident,identCol,plan_col)
190 <
190 >    
191 >  
192   !! Create row and col pointer lists
193 <    call create_IdentPtrlst(identRow,ljListTail,identPtrListRow,thisStat)
193 >  
194 >    call create_IdentPtrlst(identRow,ljListHead,identPtrListRow,thisStat)
195      if (thisStat /= 0 ) then
196         status = -1
197         return
198      endif
199 <
200 <    call create_IdentPtrlst(identCol,ljListTail,identPtrListCol,thisStat)
199 >  
200 >    call create_IdentPtrlst(identCol,ljListHead,identPtrListColumn,thisStat)
201      if (thisStat /= 0 ) then
202         status = -1
203         return
204      endif
205  
206   !! free temporary ident arrays
207 <    deallocate(identCol)
208 <    deallocate(identRow)
209 <
207 >    if (allocated(identCol)) then
208 >       deallocate(identCol)
209 >    end if
210 >    if (allocated(identCol)) then
211 >       deallocate(identRow)
212 >    endif
213  
214   !! Allocate neighbor lists for mpi simulations.
215      if (.not. allocated(point)) then
# Line 233 | Line 249 | contains
249      endif
250      isljFFinit = .true.
251  
252 +
253    end subroutine init_ljFF
254  
255  
# Line 342 | Line 359 | contains
359      type(lj_atype), pointer :: ljAtype_i
360      type(lj_atype), pointer :: ljAtype_j
361  
362 < #ifdef MPI
363 <  real( kind = DP ), dimension(3,getNcol()) :: efr
362 > #ifdef IS_MPI
363 >  real( kind = DP ), dimension(3,getNcol(plan_col)) :: efr
364    real( kind = DP ) :: pot_local
365   #else
366    real( kind = DP ), dimension(3,getNlocal()) :: efr
# Line 351 | Line 368 | contains
368    
369   !! Local arrays needed for MPI
370   #ifdef IS_MPI
371 <  real(kind = dp), dimension(3:getNrow()) :: qRow
372 <  real(kind = dp), dimension(3:getNcol()) :: qCol
371 >  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow
372 >  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol
373  
374 <  real(kind = dp), dimension(3:getNrow()) :: fRow
375 <  real(kind = dp), dimension(3:getNcol()) :: fCol
374 >  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow
375 >  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol
376 >  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp
377  
378 <  real(kind = dp), dimension(getNrow()) :: eRow
379 <  real(kind = dp), dimension(getNcol()) :: eCol
378 >  real(kind = dp), dimension(getNrow(plan_row)) :: eRow
379 >  real(kind = dp), dimension(getNcol(plan_col)) :: eCol
380  
381    real(kind = dp), dimension(getNlocal()) :: eTemp
382   #endif
# Line 387 | Line 405 | contains
405    integer :: ncol
406    integer :: natoms
407  
408 +
409 +
410 +
411   ! Make sure we are properly initialized.
412    if (.not. isljFFInit) then
413       write(default_error,*) "ERROR: lj_FF has not been properly initialized"
# Line 403 | Line 424 | contains
424    natoms = getNlocal()
425    call getRcut(rcut,rcut2=rcutsq)
426    call getRlist(rlist,rlistsq)
427 +
428   #ifndef IS_MPI
429    nrow = natoms - 1
430    ncol = natoms
# Line 416 | Line 438 | contains
438    
439   !! See if we need to update neighbor lists
440    call check(q,update_nlist)
441 +  if (firstTime) then
442 +     update_nlist = .true.
443 +     firstTime = .false.
444 +  endif
445  
446   !--------------WARNING...........................
447   ! Zero variables, NOTE:::: Forces are zeroed in C
# Line 427 | Line 453 | contains
453    pe = 0.0E0_DP
454  
455   #else
456 <    f_row = 0.0E0_DP
457 <    f_col = 0.0E0_DP
456 >    fRow = 0.0E0_DP
457 >    fCol = 0.0E0_DP
458  
459      pe_local = 0.0E0_DP
460  
# Line 439 | Line 465 | contains
465      efr = 0.0E0_DP
466  
467   ! communicate MPI positions
468 < #ifdef MPI    
469 <    call gather(q,qRow,plan_row3)
470 <    call gather(q,qCol,plan_col3)
468 > #ifdef IS_MPI    
469 >    call gather(q,qRow,plan_row3d)
470 >    call gather(q,qCol,plan_col3d)
471   #endif
472  
473  
# Line 453 | Line 479 | contains
479      
480       nlist = 0
481      
482 <    
482 >    
483  
484       do i = 1, nrow
485          point(i) = nlist + 1
486 < #ifdef MPI
486 > #ifdef IS_MPI
487          ljAtype_i => identPtrListRow(i)%this
488          tag_i = tagRow(i)
489          rxi = qRow(1,i)
# Line 472 | Line 498 | contains
498   #endif
499  
500          inner: do j = j_start, ncol
501 < #ifdef MPI
501 > #ifdef IS_MPI
502   ! Assign identity pointers and tags
503             ljAtype_j => identPtrListColumn(j)%this
504 <           tag_j = tagCol(j)
504 >           tag_j = tagColumn(j)
505             if (newtons_thrd) then
506                if (tag_i <= tag_j) then
507                   if (mod(tag_i + tag_j,2) == 0) cycle inner
# Line 496 | Line 522 | contains
522   #endif          
523             rijsq = rxij*rxij + ryij*ryij + rzij*rzij
524  
525 < #ifdef MPI
525 > #ifdef IS_MPI
526               if (rijsq <=  rlistsq .AND. &
527                    tag_j /= tag_i) then
528   #else
529 +          
530               if (rijsq <  rlistsq) then
531   #endif
532              
# Line 510 | Line 537 | contains
537                endif
538                list(nlist) = j
539  
540 <              
540 >    
541                if (rijsq <  rcutsq) then
542                  
543                   r = dsqrt(rijsq)
544        
545                   call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
546        
547 < #ifdef MPI
547 > #ifdef IS_MPI
548                  eRow(i) = eRow(i) + pot*0.5
549                  eCol(i) = eCol(i) + pot*0.5
550   #else
551 <                pe = pe + pot
551 >                    pe = pe + pot
552   #endif                
553              
554                   efr(1,j) = -rxij
# Line 535 | Line 562 | contains
562                      ftmp = dudr * drdx1
563  
564  
565 < #ifdef MPI
565 > #ifdef IS_MPI
566                      fCol(dim,j) = fCol(dim,j) - ftmp
567                      fRow(dim,i) = fRow(dim,i) + ftmp
568   #else                    
# Line 550 | Line 577 | contains
577          enddo inner
578       enddo
579  
580 < #ifdef MPI
580 > #ifdef IS_MPI
581       point(nrow + 1) = nlist + 1
582   #else
583       point(natoms) = nlist + 1
# Line 564 | Line 591 | contains
591          JEND = POINT(i+1) - 1
592          ! check thiat molecule i has neighbors
593          if (jbeg .le. jend) then
594 < #ifdef MPI
594 > #ifdef IS_MPI
595             ljAtype_i => identPtrListRow(i)%this
596             rxi = qRow(1,i)
597             ryi = qRow(2,i)
# Line 577 | Line 604 | contains
604   #endif
605             do jnab = jbeg, jend
606                j = list(jnab)
607 < #ifdef MPI
607 > #ifdef IS_MPI
608                ljAtype_j = identPtrListColumn(j)%this
609 <              rxij = wrap(rxi - q_col(1,j), 1)
610 <              ryij = wrap(ryi - q_col(2,j), 2)
611 <              rzij = wrap(rzi - q_col(3,j), 3)
609 >              rxij = wrap(rxi - qCol(1,j), 1)
610 >              ryij = wrap(ryi - qCol(2,j), 2)
611 >              rzij = wrap(rzi - qCol(3,j), 3)
612   #else
613                ljAtype_j = identPtrList(j)%this
614                rxij = wrap(rxi - q(1,j), 1)
# Line 595 | Line 622 | contains
622                   r = dsqrt(rijsq)
623                  
624                   call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
625 < #ifdef MPI
625 > #ifdef IS_MPI
626                  eRow(i) = eRow(i) + pot*0.5
627                  eCol(i) = eCol(i) + pot*0.5
628   #else
# Line 611 | Line 638 | contains
638                      
639                      drdx1 = efr(dim,j) / r
640                      ftmp = dudr * drdx1
641 < #ifdef MPI
641 > #ifdef IS_MPI
642                      fCol(dim,j) = fCol(dim,j) - ftmp
643                      fRow(dim,i) = fRow(dim,i) + ftmp
644   #else                    
# Line 627 | Line 654 | contains
654  
655  
656  
657 < #ifdef MPI
657 > #ifdef IS_MPI
658      !!distribute forces
632 –    call scatter(fRow,f,plan_row3)
659  
660 <    call scatter(fCol,f_tmp,plan_col3)
660 >    call scatter(fRow,f,plan_row3d)
661  
662 +    call scatter(fCol,fMPITemp,plan_col3d)
663 +
664      do i = 1,nlocal
665 <       f(1:3,i) = f(1:3,i) + f_tmp(1:3,i)
665 >       f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
666      end do
667  
668  
669      
670      if (do_pot) then
671   ! scatter/gather pot_row into the members of my column
672 <       call scatter(e_row,e_tmp,plan_row)
672 >       call scatter(eRow,eTemp,plan_row)
673        
674         ! scatter/gather pot_local into all other procs
675         ! add resultant to get total pot
676         do i = 1, nlocal
677 <          pe_local = pe_local + e_tmp(i)
677 >          pe_local = pe_local + eTemp(i)
678         enddo
679         if (newtons_thrd) then
680 <          e_tmp = 0.0E0_DP
681 <          call scatter(e_col,e_tmp,plan_col)
680 >          eTemp = 0.0E0_DP
681 >          call scatter(eCol,eTemp,plan_col)
682            do i = 1, nlocal
683 <             pe_local = pe_local + e_tmp(i)
683 >             pe_local = pe_local + eTemp(i)
684            enddo
685         endif
686      endif
# Line 661 | Line 689 | contains
689  
690      potE = pe
691  
664 –
692    end subroutine do_lj_ff
693  
694   !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
# Line 717 | Line 744 | contains
744      epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
745  
746  
747 <
747 >    
748 >
749      call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
750      
751      r2 = r * r

Diff Legend

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