--- trunk/mdtools/md_code/lj_FF.F90 2002/12/29 16:43:11 216 +++ trunk/mdtools/md_code/lj_FF.F90 2003/01/30 22:29:37 258 @@ -1,80 +1,807 @@ +!! Calculates Long Range forces Lennard-Jones interactions. +!! 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.15 2003-01-30 22:29:37 chuckv Exp $, $Date: 2003-01-30 22:29:37 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $ + + + module lj_ff use simulation - use definitions, ONLY : dp, ndim + use definitions + use generic_atypes +#ifdef IS_MPI + use mpiSimulation +#endif implicit none + PRIVATE +!! Number of lj_atypes in lj_atype_list integer, save :: n_lj_atypes = 0 - type, public :: lj_atype - private - sequence - integer :: atype_ident = 0 - character(len = 15) :: name = "NULL" - real ( kind = dp ) :: mass = 0.0_dp - real ( kind = dp ) :: epslon = 0.0_dp - real ( kind = dp ) :: sigma = 0.0_dp - type (lj_atype), pointer :: next => null() - end type lj_atype +!! Global list of lj atypes in simulation + type (lj_atype), pointer :: ljListHead => null() + type (lj_atype), pointer :: ljListTail => null() - type (lj_atype), pointer :: lj_atype_list => null() +!! LJ mixing array + type (lj_atype), dimension(:,:), pointer :: ljMixed => null() + + +!! Neighbor list and commom arrays +!! This should eventally get moved to a neighbor list type + integer, allocatable, dimension(:) :: point + integer, allocatable, dimension(:) :: list + integer, parameter :: listMultiplier = 80 + 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_identPtrList), dimension(:), pointer :: identPtrListRow => null() +!! Column lj_atype pointer list + type (lj_identPtrList), dimension(:), pointer :: identPtrListColumn => null() +#else + type( lj_identPtrList ), dimension(:), pointer :: identPtrList => null() +#endif + + +!! Logical has lj force field module been initialized? + logical, save :: isljFFinit = .false. + + +!! Public methods and data public :: new_lj_atype + public :: do_lj_ff + public :: getLjPot + public :: init_ljFF + + + contains - subroutine new_lj_atype(name,mass,epslon,sigma,status) - character( len = 15 ) :: name +!! Adds a new lj_atype to the list. + subroutine new_lj_atype(ident,mass,epsilon,sigma,status) real( kind = dp ), intent(in) :: mass - real( kind = dp ), intent(in) :: epslon - real( kind = dp ), intent(in) :: sigam + real( kind = dp ), intent(in) :: epsilon + real( kind = dp ), intent(in) :: sigma + integer, intent(in) :: ident integer, intent(out) :: status - type (lj_atype), pointer :: this_lj_atype - type (lj_atype), pointer :: lj_atype_ptr - + type (lj_atype), pointer :: newLJ_atype integer :: alloc_error integer :: atype_counter = 0 - + integer :: alloc_size + integer :: err_stat status = 0 - allocate(this_lj_atype,stat=alloc_error) + + +! allocate a new atype + allocate(newLJ_atype,stat=alloc_error) if (alloc_error /= 0 ) then status = -1 return end if ! assign our new lj_atype information - this_lj_atype%name = name - this_lj_atype%mass = mass - this_lj_atype%epslon = epslon - this_lj_atype%sigma = sigma + newLJ_atype%mass = mass + newLJ_atype%epsilon = epsilon + newLJ_atype%sigma = sigma + newLJ_atype%sigma2 = sigma * sigma + newLJ_atype%sigma6 = newLJ_atype%sigma2 * newLJ_atype%sigma2 & + * newLJ_atype%sigma2 +! assume that this atype will be successfully added + newLJ_atype%atype_ident = ident + newLJ_atype%atype_number = n_lj_atypes + 1 + call add_atype(newLJ_atype,ljListHead,ljListTail,err_stat) + if (err_stat /= 0 ) then + status = -1 + return + endif -! if lj_atype_list is null then we are at the top of the list. - if (.not. associated(lj_atype_list)) then - lj_atype_ptr => this_lj_atype - atype_counter = 1 - else ! we need to find the bottom of the list - lj_atype_ptr => lj_atype_list%next - find_end: do - if (.not. associated(lj_atype_ptr%next)) then - exit find_end - end if - lj_atype_ptr => lj_atype_ptr%next - end do find_end + n_lj_atypes = n_lj_atypes + 1 + + + end subroutine new_lj_atype + + + subroutine init_ljFF(nComponents,ident, status) + !! Number of components in ident array + integer, intent(inout) :: nComponents +!! Array of identities nComponents long corresponding to +!! ljatype ident. + integer, dimension(nComponents),intent(inout) :: ident +!! Result status, success = 0, error = -1 + integer, intent(out) :: Status + + integer :: alloc_stat + + integer :: thisStat + integer :: i +#ifdef IS_MPI + integer, allocatable, dimension(:) :: identRow + integer, allocatable, dimension(:) :: identCol + integer :: nrow + integer :: ncol +#endif + status = 0 + + + + +!! if were're not in MPI, we just update ljatypePtrList +#ifndef IS_MPI + call create_IdentPtrlst(ident,ljListHead,identPtrList,thisStat) + if ( thisStat /= 0 ) then + status = -1 + return + endif + +!! Allocate pointer lists + allocate(point(nComponents),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + + allocate(list(nComponents*listMultiplier),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + +! 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(plan_row) + ncol = getNcol(plan_col) + +!! Allocate temperary arrays to hold gather information + allocate(identRow(nrow),stat=alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif + + allocate(identCol(ncol),stat=alloc_stat) + if (alloc_stat /= 0 ) then + 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,ljListHead,identPtrListRow,thisStat) + if (thisStat /= 0 ) then + status = -1 + return + endif + + call create_IdentPtrlst(identCol,ljListHead,identPtrListColumn,thisStat) + if (thisStat /= 0 ) then + status = -1 + return + endif + write(*,*) "After createIdent row-col lists" +!! free temporary ident arrays + if (allocated(identCol)) then + deallocate(identCol) end if + if (allocated(identCol)) then + deallocate(identRow) + endif + +!! Allocate neighbor lists for mpi simulations. + if (.not. allocated(point)) then + allocate(point(nrow),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + + allocate(list(ncol*listMultiplier),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + else + deallocate(list) + deallocate(point) + allocate(point(nrow),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + + allocate(list(ncol*listMultiplier),stat=alloc_stat) + if (alloc_stat /=0) then + status = -1 + return + endif + endif + +#endif - lj_atype_ptr => this_lj_atype + call createMixingList(thisStat) + if (thisStat /= 0) then + status = -1 + return + endif + isljFFinit = .true. + + end subroutine init_ljFF + + + + + + + subroutine createMixingList(status) + integer :: listSize + integer :: status + integer :: i + integer :: j + + integer :: outerCounter = 0 + integer :: innerCounter = 0 + type (lj_atype), pointer :: tmpPtrOuter => null() + type (lj_atype), pointer :: tmpPtrInner => null() + status = 0 + + listSize = getListLen(ljListHead) + if (listSize == 0) then + status = -1 + return + end if + + + if (.not. associated(ljMixed)) then + allocate(ljMixed(listSize,listSize)) + else + status = -1 + return + end if + - end subroutine new_lj_atype - subroutine add_lj_atype() + tmpPtrOuter => ljListHead + tmpPtrInner => tmpPtrOuter%next + do while (associated(tmpPtrOuter)) + outerCounter = outerCounter + 1 +! do self mixing rule + ljMixed(outerCounter,outerCounter)%sigma = tmpPtrOuter%sigma + + ljMixed(outerCounter,outerCounter)%sigma2 = ljMixed(outerCounter,outerCounter)%sigma & + * ljMixed(outerCounter,outerCounter)%sigma + + ljMixed(outerCounter,outerCounter)%sigma6 = ljMixed(outerCounter,outerCounter)%sigma2 & + * ljMixed(outerCounter,outerCounter)%sigma2 & + * ljMixed(outerCounter,outerCounter)%sigma2 + + ljMixed(outerCounter,outerCounter)%epsilon = tmpPtrOuter%epsilon + innerCounter = outerCounter + 1 + do while (associated(tmpPtrInner)) + + ljMixed(outerCounter,innerCounter)%sigma = & + calcLJMix("sigma",tmpPtrOuter%sigma, & + tmpPtrInner%sigma) + + ljMixed(outerCounter,innerCounter)%sigma2 = & + ljMixed(outerCounter,innerCounter)%sigma & + * ljMixed(outerCounter,innerCounter)%sigma + + ljMixed(outerCounter,innerCounter)%sigma6 = & + ljMixed(outerCounter,innerCounter)%sigma2 & + * ljMixed(outerCounter,innerCounter)%sigma2 & + * ljMixed(outerCounter,innerCounter)%sigma2 + + ljMixed(outerCounter,innerCounter)%epsilon = & + calcLJMix("epsilon",tmpPtrOuter%epsilon, & + tmpPtrInner%epsilon) + ljMixed(innerCounter,outerCounter)%sigma = ljMixed(outerCounter,innerCounter)%sigma + ljMixed(innerCounter,outerCounter)%sigma2 = ljMixed(outerCounter,innerCounter)%sigma2 + ljMixed(innerCounter,outerCounter)%sigma6 = ljMixed(outerCounter,innerCounter)%sigma6 + ljMixed(innerCounter,outerCounter)%epsilon = ljMixed(outerCounter,innerCounter)%epsilon - end subroutine add_lj_atype + tmpPtrInner => tmpPtrInner%next + innerCounter = innerCounter + 1 + end do +! advance pointers + tmpPtrOuter => tmpPtrOuter%next + if (associated(tmpPtrOuter)) then + tmpPtrInner => tmpPtrOuter%next + endif + + end do + end subroutine createMixingList + + + + + +!! FORCE routine Calculates Lennard Jones forces. +!-------------------------------------------------------------> + subroutine do_lj_ff(q,f,potE,do_pot) +!! Position array provided by C, dimensioned by getNlocal + real ( kind = dp ), dimension(3,getNlocal()) :: q +!! Force array provided by C, dimensioned by getNlocal + real ( kind = dp ), dimension(3,getNlocal()) :: f + real ( kind = dp ) :: potE + logical ( kind = 2) :: do_pot + + type(lj_atype), pointer :: ljAtype_i + type(lj_atype), pointer :: ljAtype_j + +#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 +#endif + +!! Local arrays needed for MPI +#ifdef IS_MPI + real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow + real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol + + 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(plan_row)) :: eRow + real(kind = dp), dimension(getNcol(plan_col)) :: eCol + + real(kind = dp), dimension(getNlocal()) :: eTemp +#endif + + + + real( kind = DP ) :: pe + logical :: update_nlist + + + integer :: i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2 + integer :: nlist + integer :: j_start + integer :: tag_i,tag_j + real( kind = DP ) :: r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp + real( kind = DP ) :: rxi, ryi, rzi, rxij, ryij, rzij, rijsq + real( kind = DP ) :: rlistsq, rcutsq,rlist,rcut + +! a rig that need to be fixed. +#ifdef IS_MPI + logical :: newtons_thrd = .true. + real( kind = dp ) :: pe_local + integer :: nlocal +#endif + integer :: nrow + 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" + return + endif +#ifdef IS_MPI + if (.not. isMPISimSet()) then + write(default_error,*) "ERROR: mpiSimulation has not been properly initialized" + return + endif +#endif + +!! initialize local variables + natoms = getNlocal() + call getRcut(rcut,rcut2=rcutsq) + call getRlist(rlist,rlistsq) + +#ifndef IS_MPI + nrow = natoms - 1 + ncol = natoms +#else + nrow = getNrow(plan_row) + ncol = getNcol(plan_col) + nlocal = natoms + j_start = 1 +#endif + + +!! 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 +! Zeroing them here could delete previously computed +! Forces. +!------------------------------------------------ +#ifndef IS_MPI +! nloops = nloops + 1 + pe = 0.0E0_DP + +#else + fRow = 0.0E0_DP + fCol = 0.0E0_DP + + pe_local = 0.0E0_DP + + eRow = 0.0E0_DP + eCol = 0.0E0_DP + eTemp = 0.0E0_DP +#endif + efr = 0.0E0_DP + +! communicate MPI positions +#ifdef IS_MPI + call gather(q,qRow,plan_row3d) + call gather(q,qCol,plan_col3d) +#endif + + + if (update_nlist) then + + ! save current configuration, contruct neighbor list, + ! and calculate forces + call save_nlist(q) + + nlist = 0 + + + + do i = 1, nrow + point(i) = nlist + 1 +#ifdef IS_MPI + ljAtype_i => identPtrListRow(i)%this + tag_i = tagRow(i) + rxi = qRow(1,i) + ryi = qRow(2,i) + rzi = qRow(3,i) +#else + ljAtype_i => identPtrList(i)%this + j_start = i + 1 + rxi = q(1,i) + ryi = q(2,i) + rzi = q(3,i) +#endif + + inner: do j = j_start, ncol +#ifdef IS_MPI +! Assign identity pointers and tags + ljAtype_j => identPtrListColumn(j)%this + tag_j = tagColumn(j) + if (newtons_thrd) then + if (tag_i <= tag_j) then + if (mod(tag_i + tag_j,2) == 0) cycle inner + else + if (mod(tag_i + tag_j,2) == 1) cycle inner + endif + endif + + 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) + ryij = wrap(ryi - q(2,j), 2) + rzij = wrap(rzi - q(3,j), 3) + +#endif + rijsq = rxij*rxij + ryij*ryij + rzij*rzij + +#ifdef IS_MPI + if (rijsq <= rlistsq .AND. & + tag_j /= tag_i) then +#else + + if (rijsq < rlistsq) then +#endif + + nlist = nlist + 1 + if (nlist > size(list)) then +!! "Change how nlist size is done" + write(DEFAULT_ERROR,*) "ERROR: nlist > list size" + endif + list(nlist) = j + + + if (rijsq < rcutsq) then + + r = dsqrt(rijsq) + + call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j) + +#ifdef IS_MPI + eRow(i) = eRow(i) + pot*0.5 + eCol(i) = eCol(i) + pot*0.5 +#else + pe = pe + pot +#endif + + efr(1,j) = -rxij + efr(2,j) = -ryij + efr(3,j) = -rzij + + do dim = 1, 3 + + + drdx1 = efr(dim,j) / r + ftmp = dudr * drdx1 + + +#ifdef IS_MPI + fCol(dim,j) = fCol(dim,j) - ftmp + fRow(dim,i) = fRow(dim,i) + ftmp +#else + + f(dim,j) = f(dim,j) - ftmp + f(dim,i) = f(dim,i) + ftmp + +#endif + enddo + endif + endif + enddo inner + enddo + +#ifdef IS_MPI + point(nrow + 1) = nlist + 1 +#else + point(natoms) = nlist + 1 +#endif + + else + + ! use the list to find the neighbors + do i = 1, nrow + JBEG = POINT(i) + JEND = POINT(i+1) - 1 + ! check thiat molecule i has neighbors + if (jbeg .le. jend) then +#ifdef IS_MPI + ljAtype_i => identPtrListRow(i)%this + rxi = qRow(1,i) + ryi = qRow(2,i) + rzi = qRow(3,i) +#else + ljAtype_i => identPtrList(i)%this + rxi = q(1,i) + ryi = q(2,i) + rzi = q(3,i) +#endif + do jnab = jbeg, jend + j = list(jnab) +#ifdef IS_MPI + ljAtype_j = identPtrListColumn(j)%this + 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) + ryij = wrap(ryi - q(2,j), 2) + rzij = wrap(rzi - q(3,j), 3) +#endif + rijsq = rxij*rxij + ryij*ryij + rzij*rzij + + if (rijsq < rcutsq) then + + r = dsqrt(rijsq) + + call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j) +#ifdef IS_MPI + eRow(i) = eRow(i) + pot*0.5 + eCol(i) = eCol(i) + pot*0.5 +#else + pe = pe + pot +#endif + + + efr(1,j) = -rxij + efr(2,j) = -ryij + efr(3,j) = -rzij + + do dim = 1, 3 + + drdx1 = efr(dim,j) / r + ftmp = dudr * drdx1 +#ifdef IS_MPI + fCol(dim,j) = fCol(dim,j) - ftmp + fRow(dim,i) = fRow(dim,i) + ftmp +#else + f(dim,j) = f(dim,j) - ftmp + f(dim,i) = f(dim,i) + ftmp +#endif + enddo + endif + enddo + endif + enddo + endif + + + +#ifdef IS_MPI + !!distribute forces + call scatter(fRow,f,plan_row3d) + + call scatter(fCol,fMPITemp,plan_col3d) + + do i = 1,nlocal + 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(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 + eTemp(i) + enddo + if (newtons_thrd) then + eTemp = 0.0E0_DP + call scatter(eCol,eTemp,plan_col) + do i = 1, nlocal + pe_local = pe_local + eTemp(i) + enddo + endif + endif + potE = pe_local +#endif + + potE = pe + + + + + end subroutine do_lj_ff + +!! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second +!! derivatives. + subroutine getLjPot(r,pot,dudr,atype1,atype2,d2,status) +! arguments +!! Length of vector between particles + real( kind = dp ), intent(in) :: r +!! Potential Energy + real( kind = dp ), intent(out) :: pot +!! Derivatve wrt postion + real( kind = dp ), intent(out) :: dudr +!! Second Derivative, optional, used mainly for normal mode calculations. + real( kind = dp ), intent(out), optional :: d2 + + type (lj_atype), pointer :: atype1 + type (lj_atype), pointer :: atype2 + + integer, intent(out), optional :: status + +! local Variables + real( kind = dp ) :: sigma + real( kind = dp ) :: sigma2 + real( kind = dp ) :: sigma6 + real( kind = dp ) :: epsilon + + real( kind = dp ) :: rcut + real( kind = dp ) :: rcut2 + real( kind = dp ) :: rcut6 + real( kind = dp ) :: r2 + real( kind = dp ) :: r6 + + real( kind = dp ) :: t6 + real( kind = dp ) :: t12 + real( kind = dp ) :: tp6 + real( kind = dp ) :: tp12 + real( kind = dp ) :: delta + + logical :: doSec = .false. + + integer :: errorStat + +!! Optional Argument Checking +! Check to see if we need to do second derivatives + + if (present(d2)) doSec = .true. + if (present(status)) status = 0 + +! Look up the correct parameters in the mixing matrix + sigma = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma + sigma2 = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2 + sigma6 = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6 + epsilon = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon + + + + + call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat) + + r2 = r * r + r6 = r2 * r2 * r2 + + t6 = sigma6/ r6 + t12 = t6 * t6 + + + + tp6 = sigma6 / rcut6 + tp12 = tp6*tp6 + + delta = -4.0E0_DP*epsilon * (tp12 - tp6) + + if (r.le.rcut) then + pot = 4.0E0_DP * epsilon * (t12 - t6) + delta + dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / r + if(doSec) d2 = 24.0E0_DP * epsilon * (26.0E0_DP*t12 - 7.0E0_DP*t6)/r/r + else + pot = 0.0E0_DP + dudr = 0.0E0_DP + if(doSec) d2 = 0.0E0_DP + endif + + return + + + + end subroutine getLjPot + + +!! Calculates the mixing for sigma or epslon based on character string for initialzition of mixing array. + function calcLJMix(thisParam,param1,param2,status) result(myMixParam) + character(len=*) :: thisParam + real(kind = dp) :: param1 + real(kind = dp) :: param2 + real(kind = dp ) :: myMixParam + integer, optional :: status + + + myMixParam = 0.0_dp + + if (present(status)) status = 0 + + select case (thisParam) + + case ("sigma") + myMixParam = 0.5_dp * (param1 + param2) + case ("epsilon") + myMixParam = sqrt(param1 * param2) + case default + status = -1 + end select + + end function calcLJMix + + + end module lj_ff