--- trunk/mdtools/md_code/lj_FF.F90	2002/12/19 21:59:51	215
+++ trunk/mdtools/md_code/lj_FF.F90	2003/01/22 21:45:20	240
@@ -1,41 +1,893 @@
+!! 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.9 2003-01-22 21:45:20 chuckv Exp $, $Date: 2003-01-22 21:45:20 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $
+
+
+
 module lj_ff
   use simulation
   use definitions, ONLY : dp, ndim
+#ifdef IS_MPI
+  use mpiSimulation
+#endif
   implicit none
+  PRIVATE
 
+!! Number of lj_atypes in lj_atype_list
   integer, save :: n_lj_atypes = 0
 
+!! Starting Size for ljMixed Array
+  integer, parameter :: ljMixed_blocksize = 10
+
+!! Basic atom type for a Lennard-Jones Atom.
   type, public :: lj_atype
      private
      sequence
+!! Unique number for place in linked list
+     integer :: atype_number = 0
+!! Unique indentifier number (ie atomic no, etc)
      integer :: atype_ident = 0
-     character(len = 15) :: name = "NULL"
+!! Mass of Particle
      real ( kind = dp )  :: mass = 0.0_dp
+!! Lennard-Jones epslon
      real ( kind = dp )  :: epslon = 0.0_dp
+!! Lennard-Jones Sigma
      real ( kind = dp )  :: sigma = 0.0_dp
+!! Lennard-Jones Sigma Squared
+     real ( kind = dp )  :: sigma2 = 0.0_dp
+!! Lennard-Jones Sigma to sixth
+     real ( kind = dp )  :: sigma6 = 0.0_dp
+!! Pointer for linked list creation
      type (lj_atype), pointer :: next => null()
   end type lj_atype
 
+!! Pointer type for atype ident array
+  type, public :: lj_atypePtr
+     type (lj_atype), pointer :: this => null()
+  end type lj_atypePtr
+
+!! Global list of lj atypes in simulation
   type (lj_atype), pointer :: lj_atype_list => null()
+!! LJ mixing array  
+  type (lj_atype), dimension(:,:), allocatable, pointer :: ljMixed =>null()
+!! identity pointer list for force loop.
+  type (lj_atypePtr), dimension(:), allocatable :: identPtrList
 
+
+!! 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
+
+#ifdef IS_MPI
+! Universal routines: All types of force calculations will need these arrays
+! Arrays specific to a type of force calculation should be declared in that module.
+  real( kind = dp ), allocatable, dimension(:,:) :: qRow
+  real( kind = dp ), allocatable, dimension(:,:) :: qColumn
+
+  real( kind = dp ), allocatable, dimension(:,:) :: fRow
+  real( kind = dp ), allocatable, dimension(:,:) :: fColumn
+
+  type (lj_atypePtr), dimension(:), allocatable :: identPtrListRow
+  type (lj_atypePtr), dimension(:), allocatable :: identPtrListColumn
+#endif
+
+
+
+  logical :: isljFFinit = .false.
+
+
+!! Public methods and data
   public :: new_lj_atype 
+  public :: do_lj_ff
+  public :: getLjPot
+  
 
+  
+
+
 contains
 
-  subroutine new_lj_atype(name,mass,epslon,sigma,status)
-    character( len = 15 ) :: name
+  subroutine new_lj_atype(ident,mass,epslon,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) :: 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), allocatable, dimension(:,:), pointer :: thisMix
+    type (lj_atype), allocatable, dimension(:,:), pointer :: oldMix
+    integer :: alloc_error
+    integer :: atype_counter = 0
+    integer :: alloc_size
+
     status = 0
 
+
+
+! allocate a new atype    
+    allocate(this_lj_atype,stat=alloc_error)
+    if (alloc_error /= 0 ) then
+       status = -1
+       return
+    end if
+
+! assign our new lj_atype information
+    this_lj_atype%mass       = mass
+    this_lj_atype%epslon     = epslon
+    this_lj_atype%sigma      = sigma
+    this_lj_atype%sigma2     = sigma * sigma
+    this_lj_atype%sigma6     = this_lj_atype%sigma2 * this_lj_atype%sigma2 &
+         * this_lj_atype%sigma2
+! assume that this atype will be successfully added
+    this_lj_atype%atype_ident = ident
+    this_lj_atype%number = n_lj_atypes + 1
+
+
+! First time through allocate a array of size ljMixed_blocksize
+    if(.not. associated(ljMixed)) then
+       allocate(thisMix(ljMixed_blocksize,ljMixed_blocksize))
+       if (alloc_error /= 0 ) then
+          status = -1
+          return
+       end if
+       ljMixed => thisMix
+! If we have outgrown ljMixed_blocksize, allocate a new matrix twice the size and
+! point ljMix at the new matrix.
+    else if( (n_lj_atypes + 1) > size(ljMixed)) then
+       alloc_size = 2*size(ljMix)
+       allocate(thisMix(alloc_size,alloc_size))
+       if (alloc_error /= 0 ) then
+          status = -1
+          return
+       end if
+! point oldMix at old ljMixed array
+       oldMix => ljMixed
+! Copy oldMix into new Mixed array      
+       thisMix = oldMix
+! Point ljMixed at new array
+       ljMixed => thisMix
+! Free old array so we don't have a memory leak
+       deallocate(oldMix)
+    endif
+
+
+
+
+
+! Find bottom of atype master list
+! 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 to insert new atype
+       lj_atype_ptr => lj_atype_list%next
+       atype_counter = 1
+       find_end: do
+          if (.not. associated(lj_atype_ptr%next)) then
+             exit find_end
+          end if
+! Set up mixing for new atype and current atype in list
+       ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma  =  &
+            calcLJMix("sigma",this_lj_atype%sigma, &
+            lj_atype_prt%sigma)
+
+       ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2  = &
+            ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma & 
+            * ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma
+
+       ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma6 = &
+            ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 &
+            * ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2 &
+            * ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%sigma2
+
+       ljMix(this_lj_atype%atype_number,lj_atype_ptr%atype_number)%epslon = &
+            calcLJMix("epslon",this_lj_atype%epslon, &
+            lj_atype_prt%epslon)
+
+! Advance to next pointer
+       lj_atype_ptr => lj_atype_ptr%next
+       atype_counter = atype_counter + 1
+          
+       end do find_end
+    end if
+
+
+
+    
+! Insert new atype at end of list
+    lj_atype_ptr => this_lj_atype
+! Increment number of atypes
+
+    n_lj_atypes = n_lj_atypes + 1
+
+! Set up self mixing
+
+    ljMix(n_lj_atypes,n_lj_atypes)%sigma  = this_lj_atype%sigma
+
+    ljMix(n_lj_atypes,n_lj_atypes)%sigma2  = ljMix(n_lj_atypes,n_lj_atypes)%sigma & 
+            * ljMix(n_lj_atypes,n_lj_atypes)%sigma
+
+    ljMix(n_lj_atypes,n_lj_atypes)%sigma6 = ljMix(n_lj_atypes,n_lj_atypes)%sigma2 &
+            * ljMix(n_lj_atypes,n_lj_atypes)%sigma2 &
+            * ljMix(n_lj_atypes,n_lj_atypes)%sigma2
+
+    ljMix(n_lj_atypes,n_lj_atypes)%epslon = this_lj_atype%epslon
+
+
   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 :: thisStat
+#ifdef IS_MPI
+    integer, allocatable, dimension(:) :: identRow
+    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 new_ljatypePtrList(nComponents,ident,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
+       status = -1
+       return
+    endif
+    nrow = getNrow()
+    ncol = getNcol()
+!! 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 new_ljatypePtrList(nrow,identRow,identPtrListRow,thisStat)
+    if (thisStat /= 0 ) then
+       status = -1
+       return
+    endif
+
+    call new_ljatypePtrList(ncol,identCol,identPtrListCol,thisStat)
+    if (thisStat /= 0 ) then
+       status = -1
+       return
+    endif
+
+!! free temporary ident arrays
+    deallocate(identCol)
+    deallocate(identRow)
+
+!! Allocate Simulation arrays
+!! NOTE: This bit of code should be fixed, it can cause large
+!! memory fragmentation if call repeatedly
+
+    if (.not.allocated(qRow)) then
+       allocate(qRow(3,nrow),stat=alloc_stat)
+       if (alloc_stat /= 0 ) then
+          status = -1
+          return
+       endif
+    else
+       deallocate(qrow)
+       allocate(qRow(3,nrow),stat=alloc_stat)
+       if (alloc_stat /= 0 ) then
+          status = -1
+          return
+       endif
+    endif
+
+    if (.not.allocated(3,qCol)) then
+       allocate(qCol(ncol),stat=alloc_stat)
+       if (alloc_stat /= 0 ) then
+          status = -1
+          return
+       endif
+    else
+       deallocate(qCol)
+       allocate(qCol(3,ncol),stat=alloc_stat)
+       if (alloc_stat /= 0 ) then
+          status = -1
+          return
+       endif
+    endif
+
+    if (.not.allocated(fRow)) then
+       allocate(fRow(3,nrow),stat=alloc_stat)
+       if (alloc_stat /= 0 ) then
+          status = -1
+          return
+       endif
+    else
+       deallocate(fRow)
+       allocate(fRow(3,nrow),stat=alloc_stat)
+       if (alloc_stat /= 0 ) then
+          status = -1
+          return
+       endif
+    endif
+
+    if (.not.allocated(fCol)) then
+       allocate(fCol(3,ncol),stat=alloc_stat)
+       if (alloc_stat /= 0 ) then
+          status = -1
+          return
+       endif
+    else
+       deallocate(fCol)
+       allocate(fCol(3,ncol),stat=alloc_stat)
+       if (alloc_stat /= 0 ) then
+          status = -1
+          return
+       endif
+    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
+    
+
+    isljFFinit = .true.
+
+
+  end subroutine init_ljFF
+
+
+
+
+
+!! Takes an ident array and creates an atype pointer list
+!! based on those identities
+  subroutine new_ljatypePtrList(mysize,ident,PtrList,status)
+    integer, intent(in) :: mysize
+    integer, intent(in) :: ident
+    integer, optional :: status
+    type(lj_atypePtr), dimension(:) :: PtrList
+
+    integer :: thisIdent
+    integer :: i
+    integer :: alloc_error
+    type (lj_atype), pointer :: tmpPtr
+
+    if (present(status)) status = 0
+
+! First time through, allocate list
+    if (.not.(allocated)) then
+       allocate(PtrList(mysize))
+    else
+! We want to creat a new ident list so free old list
+       deallocate(PrtList)
+       allocate(PtrList(mysize))
+    endif
+
+! Match pointer list
+    do i = 1, mysize
+       thisIdent = ident(i)
+       call getLJatype(thisIdent,tmpPtr)
+ 
+      if (.not. associated(tmpPtr)) then
+          status = -1
+          return
+       endif
+       
+       PtrList(i)%this => tmpPtr
+    end do
+
+  end subroutine new_ljatypePtrList
+
+!! Finds a lj_atype based upon numerical ident 
+!! returns a null pointer if error
+  subroutine getLJatype(ident,ljAtypePtr)
+    integer, intent(in) :: ident
+    type (lj_atype), intent(out),pointer :: ljAtypePtr => null()
+    
+    type (lj_atype), pointer :: tmplj_atype_ptr => null()
 
+    if(.not. associated(lj_atype_list)) return
 
+! Point at head of list.
+    tmplj_atype_ptr => lj_atype_list
+    find_ident: do 
+       if (.not.associated(tmplj_atype_ptr)) then
+          exit find_ident
+       else if( lj_atype_ptr%atype_ident == ident)
+          ljAtypePtr => tmplj_atype_ptr
+          exit find_ident
+       endif
+       tmplj_atype_ptr => tmplj_atype_ptr%next
+    end do find_ident
 
+  end subroutine getLJatype
 
 
+!! FORCE routine Calculates Lennard Jones forces.
+!------------------------------------------------------------->
+  subroutine do_lj_ff(q,f,potE,do_pot)
+    real ( kind = dp ), dimension(ndim,) :: q
+    real ( kind = dp ), dimension(ndim,nLRparticles) :: f
+    real ( kind = dp ) :: potE
+    logical ( kind = 2) :: do_pot
+    
+    type(lj_atype), pointer :: ljAtype_i
+    type(lj_atype), pointer :: ljAtype_j
 
+#ifdef MPI
+  real( kind = DP ), dimension(3,ncol) :: efr
+  real( kind = DP ) :: pot_local
+#else
+!  real( kind = DP ), dimension(3,natoms) :: efr
+#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
+
+  integer :: nrow
+  integer :: ncol
+
+  if (.not. isljFFInit) then
+     write(default_error,*) "ERROR: lj_FF has not been properly initialized"
+     return
+  endif
+
+#ifndef IS_MPI
+  nrow = natoms - 1
+  ncol = natoms
+#else
+  nrow = getNrow(plan_row)
+  ncol = getNcol(plan_col)
+  j_start = 1
+#endif
+
+  
+
+  call check(update_nlist)
+
+!--------------WARNING...........................
+! Zero variables, NOTE:::: Forces are zeroed in C
+! Zeroing them here could delete previously computed
+! Forces.
+!------------------------------------------------
+#ifndef IS_MPI
+  nloops = nloops + 1
+  pot = 0.0E0_DP
+  e = 0.0E0_DP
+#else
+    f_row = 0.0E0_DP
+    f_col = 0.0E0_DP
+
+    pot_local = 0.0E0_DP
+
+    e_row = 0.0E0_DP
+    e_col = 0.0E0_DP
+    e_tmp = 0.0E0_DP
+#endif
+    efr = 0.0E0_DP
+
+! communicate MPI positions
+#ifdef MPI    
+    call gather(q,qRow,plan_row3)
+    call gather(q,qCol,plan_col3)
+#endif
+
+#ifndef MPI
+
+#endif
+
+  if (update_nlist) then
+
+     ! save current configuration, contruct neighbor list, 
+     ! and calculate forces
+     call save_nlist()
+     
+     nlist = 0
+     
+     
+
+     do i = 1, nrow
+        point(i) = nlist + 1
+#ifdef 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 MPI
+! Assign identity pointers and tags
+           ljAtype_j => identPtrListColumn(j)%this
+           tag_j = tagCol(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 MPI
+             if (rijsq <=  rlstsq .AND. &
+                  tag_j /= tag_i) then
+#else
+             if (rijsq <  rlstsq) then
+#endif
+             
+              nlist = nlist + 1
+              if (nlist > size(list)) then
+#warning "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 MPI
+                e_row(i) = e_row(i) + pot*0.5
+                e_col(i) = e_col(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 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 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 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 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)
+#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 MPI
+                e_row(i) = e_row(i) + pot*0.5
+                e_col(i) = e_col(i) + pot*0.5
+#else
+               if (do_pot)  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 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 MPI
+    !!distribute forces
+    call scatter(fRow,f,plan_row3)
+
+    call scatter(fCol,f_tmp,plan_col3)
+    do i = 1,nlocal
+       do dim = 1,3
+          f(dim,i) = f(dim,i) + f_tmp(dim,i)
+       end do
+    end do
+
+
+    
+    if (do_pot) then
+! scatter/gather pot_row into the members of my column
+       call scatter(e_row,e_tmp,plan_row)
+       
+       ! scatter/gather pot_local into all other procs
+       ! add resultant to get total pot
+       do i = 1, nlocal
+          pot_local = pot_local + e_tmp(i)
+       enddo
+       if (newtons_thrd) then
+          e_tmp = 0.0E0_DP
+          call scatter(e_col,e_tmp,plan_col)
+          do i = 1, nlocal
+             pot_local = pot_local + e_tmp(i)
+          enddo
+       endif
+    endif
+#endif
+
+
+
+
+  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), intent(in), pointer :: atype1
+    type (lj_atype), intent(in), pointer :: atype2
+
+    integer, intent(out), optional :: status
+
+! local Variables
+    real( kind = dp ) :: sigma
+    real( kind = dp ) :: sigma2
+    real( kind = dp ) :: sigma6
+    real( kind = dp ) :: epslon
+
+    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
+    epslon  = ljMixed(atype1%atype_ident,atype2_atype_ident)%epslon
+ 
+
+
+    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
+       u = 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
+       u = 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 ("epslon")
+       myMixParam = sqrt(param1 * param2)
+    case default
+       status = -1
+    end select
+
+  end function calcLJMix
+
+
+
 end module lj_ff