--- trunk/mdtools/md_code/lj_FF.F90 2003/01/27 18:28:11 247 +++ trunk/mdtools/md_code/lj_FF.F90 2003/01/30 20:03:37 254 @@ -2,7 +2,7 @@ !! Corresponds to the force field defined in lj_FF.cpp !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: lj_FF.F90,v 1.11 2003-01-27 18:28:11 chuckv Exp $, $Date: 2003-01-27 18:28:11 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $ +!! @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 $ @@ -36,13 +36,14 @@ module lj_ff integer :: nListAllocs = 0 integer, parameter :: maxListAllocs = 5 + logical, save :: firstTime = .True. !! Atype identity pointer lists #ifdef IS_MPI !! Row lj_atype pointer list - type (lj_atypePtr), dimension(:), pointer :: identPtrListRow => null() + type (lj_identPtrList), dimension(:), pointer :: identPtrListRow => null() !! Column lj_atype pointer list - type (lj_atypePtr), dimension(:), pointer :: identPtrListColumn => null() + type (lj_identPtrList), dimension(:), pointer :: identPtrListColumn => null() #else type( lj_identPtrList ), dimension(:), pointer :: identPtrList => null() #endif @@ -128,9 +129,12 @@ contains integer, allocatable, dimension(:) :: identCol integer :: nrow integer :: ncol - integer :: alloc_stat #endif status = 0 + + + + !! if were're not in MPI, we just update ljatypePtrList #ifndef IS_MPI call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat) @@ -154,13 +158,15 @@ contains ! if were're in MPI, we also have to worry about row and col lists #else + ! We can only set up forces if mpiSimulation has been setup. if (.not. isMPISimSet()) then + write(default_error,*) "MPI is not set" status = -1 return endif - nrow = getNrow() - ncol = getNcol() + nrow = getNrow(plan_row) + ncol = getNcol(plan_col) !! Allocate temperary arrays to hold gather information allocate(identRow(nrow),stat=alloc_stat) if (alloc_stat /= 0 ) then @@ -173,18 +179,22 @@ contains status = -1 return endif + !! Gather idents into row and column idents + call gather(ident,identRow,plan_row) call gather(ident,identCol,plan_col) + !! Create row and col pointer lists - call create_IdentPtrlst(identRow,ljListTail,identPtrListRow,thisStat) + + call create_IdentPtrlst(identRow,ljListHead,identPtrListRow,thisStat) if (thisStat /= 0 ) then status = -1 return endif - call create_IdentPtrlst(identCol,ljListTail,identPtrListCol,thisStat) + call create_IdentPtrlst(identCol,ljListHead,identPtrListColumn,thisStat) if (thisStat /= 0 ) then status = -1 return @@ -267,13 +277,11 @@ contains end if - write(*,*) "Creating mixing list" + tmpPtrOuter => ljListHead tmpPtrInner => tmpPtrOuter%next do while (associated(tmpPtrOuter)) outerCounter = outerCounter + 1 - write(*,*) "Outer index is: ", outerCounter - write(*,*) "For atype: ", tmpPtrOuter%atype_ident ! do self mixing rule ljMixed(outerCounter,outerCounter)%sigma = tmpPtrOuter%sigma @@ -289,7 +297,6 @@ contains innerCounter = outerCounter + 1 do while (associated(tmpPtrInner)) - write(*,*) "Entered inner loop for: ", innerCounter ljMixed(outerCounter,innerCounter)%sigma = & calcLJMix("sigma",tmpPtrOuter%sigma, & tmpPtrInner%sigma) @@ -345,8 +352,8 @@ contains type(lj_atype), pointer :: ljAtype_i type(lj_atype), pointer :: ljAtype_j -#ifdef MPI - real( kind = DP ), dimension(3,getNcol()) :: efr +#ifdef IS_MPI + real( kind = DP ), dimension(3,getNcol(plan_col)) :: efr real( kind = DP ) :: pot_local #else real( kind = DP ), dimension(3,getNlocal()) :: efr @@ -354,14 +361,15 @@ contains !! Local arrays needed for MPI #ifdef IS_MPI - real(kind = dp), dimension(3:getNrow()) :: qRow - real(kind = dp), dimension(3:getNcol()) :: qCol + real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow + real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol - real(kind = dp), dimension(3:getNrow()) :: fRow - real(kind = dp), dimension(3:getNcol()) :: fCol + real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow + real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol + real(kind = dp), dimension(3,getNlocal()) :: fMPITemp - real(kind = dp), dimension(getNrow()) :: eRow - real(kind = dp), dimension(getNcol()) :: eCol + real(kind = dp), dimension(getNrow(plan_row)) :: eRow + real(kind = dp), dimension(getNcol(plan_col)) :: eCol real(kind = dp), dimension(getNlocal()) :: eTemp #endif @@ -390,6 +398,9 @@ contains integer :: ncol integer :: natoms + + + ! Make sure we are properly initialized. if (.not. isljFFInit) then write(default_error,*) "ERROR: lj_FF has not been properly initialized" @@ -406,6 +417,7 @@ contains natoms = getNlocal() call getRcut(rcut,rcut2=rcutsq) call getRlist(rlist,rlistsq) + #ifndef IS_MPI nrow = natoms - 1 ncol = natoms @@ -419,6 +431,10 @@ contains !! See if we need to update neighbor lists call check(q,update_nlist) + if (firstTime) then + update_nlist = .true. + firstTime = .false. + endif !--------------WARNING........................... ! Zero variables, NOTE:::: Forces are zeroed in C @@ -430,8 +446,8 @@ contains pe = 0.0E0_DP #else - f_row = 0.0E0_DP - f_col = 0.0E0_DP + fRow = 0.0E0_DP + fCol = 0.0E0_DP pe_local = 0.0E0_DP @@ -442,9 +458,9 @@ contains efr = 0.0E0_DP ! communicate MPI positions -#ifdef MPI - call gather(q,qRow,plan_row3) - call gather(q,qCol,plan_col3) +#ifdef IS_MPI + call gather(q,qRow,plan_row3d) + call gather(q,qCol,plan_col3d) #endif @@ -456,11 +472,11 @@ contains nlist = 0 - + do i = 1, nrow point(i) = nlist + 1 -#ifdef MPI +#ifdef IS_MPI ljAtype_i => identPtrListRow(i)%this tag_i = tagRow(i) rxi = qRow(1,i) @@ -475,10 +491,10 @@ contains #endif inner: do j = j_start, ncol -#ifdef MPI +#ifdef IS_MPI ! Assign identity pointers and tags ljAtype_j => identPtrListColumn(j)%this - tag_j = tagCol(j) + tag_j = tagColumn(j) if (newtons_thrd) then if (tag_i <= tag_j) then if (mod(tag_i + tag_j,2) == 0) cycle inner @@ -499,10 +515,11 @@ contains #endif rijsq = rxij*rxij + ryij*ryij + rzij*rzij -#ifdef MPI +#ifdef IS_MPI if (rijsq <= rlistsq .AND. & tag_j /= tag_i) then #else + if (rijsq < rlistsq) then #endif @@ -513,18 +530,18 @@ contains endif list(nlist) = j - + if (rijsq < rcutsq) then r = dsqrt(rijsq) call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j) -#ifdef MPI +#ifdef IS_MPI eRow(i) = eRow(i) + pot*0.5 eCol(i) = eCol(i) + pot*0.5 #else - pe = pe + pot + pe = pe + pot #endif efr(1,j) = -rxij @@ -538,7 +555,7 @@ contains ftmp = dudr * drdx1 -#ifdef MPI +#ifdef IS_MPI fCol(dim,j) = fCol(dim,j) - ftmp fRow(dim,i) = fRow(dim,i) + ftmp #else @@ -553,7 +570,7 @@ contains enddo inner enddo -#ifdef MPI +#ifdef IS_MPI point(nrow + 1) = nlist + 1 #else point(natoms) = nlist + 1 @@ -567,7 +584,7 @@ contains JEND = POINT(i+1) - 1 ! check thiat molecule i has neighbors if (jbeg .le. jend) then -#ifdef MPI +#ifdef IS_MPI ljAtype_i => identPtrListRow(i)%this rxi = qRow(1,i) ryi = qRow(2,i) @@ -580,11 +597,11 @@ contains #endif do jnab = jbeg, jend j = list(jnab) -#ifdef MPI +#ifdef IS_MPI ljAtype_j = identPtrListColumn(j)%this - rxij = wrap(rxi - q_col(1,j), 1) - ryij = wrap(ryi - q_col(2,j), 2) - rzij = wrap(rzi - q_col(3,j), 3) + rxij = wrap(rxi - qCol(1,j), 1) + ryij = wrap(ryi - qCol(2,j), 2) + rzij = wrap(rzi - qCol(3,j), 3) #else ljAtype_j = identPtrList(j)%this rxij = wrap(rxi - q(1,j), 1) @@ -598,7 +615,7 @@ contains r = dsqrt(rijsq) call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j) -#ifdef MPI +#ifdef IS_MPI eRow(i) = eRow(i) + pot*0.5 eCol(i) = eCol(i) + pot*0.5 #else @@ -614,7 +631,7 @@ contains drdx1 = efr(dim,j) / r ftmp = dudr * drdx1 -#ifdef MPI +#ifdef IS_MPI fCol(dim,j) = fCol(dim,j) - ftmp fRow(dim,i) = fRow(dim,i) + ftmp #else @@ -630,32 +647,32 @@ contains -#ifdef MPI +#ifdef IS_MPI !!distribute forces - call scatter(fRow,f,plan_row3) + call scatter(fRow,f,plan_row3d) - call scatter(fCol,f_tmp,plan_col3) + call scatter(fCol,fMPITemp,plan_col3d) do i = 1,nlocal - f(1:3,i) = f(1:3,i) + f_tmp(1:3,i) + f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i) end do if (do_pot) then ! scatter/gather pot_row into the members of my column - call scatter(e_row,e_tmp,plan_row) + call scatter(eRow,eTemp,plan_row) ! scatter/gather pot_local into all other procs ! add resultant to get total pot do i = 1, nlocal - pe_local = pe_local + e_tmp(i) + pe_local = pe_local + eTemp(i) enddo if (newtons_thrd) then - e_tmp = 0.0E0_DP - call scatter(e_col,e_tmp,plan_col) + eTemp = 0.0E0_DP + call scatter(eCol,eTemp,plan_col) do i = 1, nlocal - pe_local = pe_local + e_tmp(i) + pe_local = pe_local + eTemp(i) enddo endif endif @@ -664,7 +681,9 @@ contains potE = pe + + end subroutine do_lj_ff !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second @@ -720,7 +739,8 @@ contains epsilon = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon - + + call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat) r2 = r * r