[ViewVC] Diff of: group/trunk/mdtools/md_code/lj

Comparing trunk/mdtools/md_code/lj_FF.F90 (file contents):
Revision 215 by chuckv, Thu Dec 19 21:59:51 2002 UTC vs.
Revision 258 by chuckv, Thu Jan 30 22:29:37 2003 UTC

-+
+!! 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 :: newLJ_atype
-+
+    integer :: alloc_error
-+
+    integer :: atype_counter = 0
-+
+    integer :: alloc_size
-+
+    integer :: err_stat
+    status = 0
-+
-+
-+
+! 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
-+
+    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
-+
-+
+    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
-+
-+
+    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
-+
-+
-+
-+
+    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
-+
-+
-+
+          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

Diff Legend

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

Comparing trunk/mdtools/md_code/lj_FF.F90 (file contents): Revision 215 by chuckv, Thu Dec 19 21:59:51 2002 UTC vs. Revision 258 by chuckv, Thu Jan 30 22:29:37 2003 UTC

Diff Legend

Comparing trunk/mdtools/md_code/lj_FF.F90 (file contents):
Revision 215 by chuckv, Thu Dec 19 21:59:51 2002 UTC vs.
Revision 258 by chuckv, Thu Jan 30 22:29:37 2003 UTC