--- trunk/mdtools/md_code/lj_FF.F90	2003/01/10 21:56:10	235
+++ trunk/mdtools/md_code/lj_FF.F90	2003/01/20 22:36:12	239
@@ -1,3 +1,11 @@
+!! 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.8 2003-01-20 22:36:12 chuckv Exp $, $Date: 2003-01-20 22:36:12 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
+
+
+
 module lj_ff
   use simulation
   use definitions, ONLY : dp, ndim
@@ -13,6 +21,7 @@ module lj_ff
 !! Starting Size for ljMixed Array
   integer, parameter :: ljMixed_blocksize = 10
 
+!! Basic atom type for a Lennard-Jones Atom.
   type, public :: lj_atype
      private
      sequence
@@ -48,8 +57,12 @@ module lj_ff
 
 
 !! 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
@@ -59,13 +72,16 @@ module lj_ff
 
   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
@@ -211,17 +227,191 @@ contains
   end subroutine new_lj_atype
 
 
-  subroutine init_ljFF()
+  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)
@@ -286,13 +476,17 @@ contains
   end subroutine getLJatype
 
 
-! FORCE routine
+!! 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
@@ -314,21 +508,32 @@ contains
   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
-  f = 0.0E0_DP
   e = 0.0E0_DP
 #else
     f_row = 0.0E0_DP
@@ -344,8 +549,8 @@ contains
 
 ! communicate MPI positions
 #ifdef MPI    
-    call gather(q,q_row,plan_row3)
-    call gather(q,q_col,plan_col3)
+    call gather(q,qRow,plan_row3)
+    call gather(q,qCol,plan_col3)
 #endif
 
 #ifndef MPI
@@ -365,11 +570,13 @@ contains
      do i = 1, nrow
         point(i) = nlist + 1
 #ifdef MPI
-        tag_i = tag_row(i)
-        rxi = q_row(1,i)
-        ryi = q_row(2,i)
-        rzi = q_row(3,i)
+        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)
@@ -378,7 +585,9 @@ contains
 
         inner: do j = j_start, ncol
 #ifdef MPI
-           tag_j = tag_col(j)
+! 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
@@ -387,11 +596,11 @@ contains
               endif
            endif
 
-           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)
            ryij = wrap(ryi - q(2,j), 2)
            rzij = wrap(rzi - q(3,j), 3)
@@ -408,13 +617,8 @@ contains
              
               nlist = nlist + 1
               if (nlist > size(list)) then
-                 call info("FORCE_LJ","error size list smaller then nlist") 
-                 write(msg,*) "nlist size(list)", nlist, size(list)
-                 call info("FORCE_LJ",msg)
-#ifdef MPI
-                 call mpi_abort(MPI_COMM_WORLD,mpi_err)
-#endif
-                 stop
+#warning "Change how nlist size is done"
+                 write(DEFAULT_ERROR,*) "ERROR: nlist > list size"
               endif
               list(nlist) = j
 
@@ -423,7 +627,7 @@ contains
                 
                  r = dsqrt(rijsq)
        
-                 call LJ_mix(r,pot,dudr,d2,i,j)
+                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
        
 #ifdef MPI
                 e_row(i) = e_row(i) + pot*0.5
@@ -444,8 +648,8 @@ contains
 
 
 #ifdef MPI
-                    f_col(dim,j) = f_col(dim,j) - ftmp
-                    f_row(dim,i) = f_row(dim,i) + ftmp
+                    fCol(dim,j) = fCol(dim,j) - ftmp
+                    fRow(dim,i) = fRow(dim,i) + ftmp
 #else                    
              
                     f(dim,j) = f(dim,j) - ftmp
@@ -473,10 +677,12 @@ contains
         ! check thiat molecule i has neighbors
         if (jbeg .le. jend) then
 #ifdef MPI
-           rxi = q_row(1,i)
-           ryi = q_row(2,i)
-           rzi = q_row(3,i)
+           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)
@@ -484,13 +690,15 @@ contains
            do jnab = jbeg, jend
               j = list(jnab)
 #ifdef MPI
-                rxij = wrap(rxi - q_col(1,j), 1)
-                ryij = wrap(ryi - q_col(2,j), 2)
-                rzij = wrap(rzi - q_col(3,j), 3)
+              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
-                rxij = wrap(rxi - q(1,j), 1)
-                ryij = wrap(ryi - q(2,j), 2)
-                rzij = wrap(rzi - q(3,j), 3)
+              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
               
@@ -498,7 +706,7 @@ contains
 
                  r = dsqrt(rijsq)
                  
-                 call LJ_mix(r,pot,dudr,d2,i,j)
+                 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
@@ -516,8 +724,8 @@ contains
                     drdx1 = efr(dim,j) / r
                     ftmp = dudr * drdx1
 #ifdef MPI
-                    f_col(dim,j) = f_col(dim,j) - ftmp
-                    f_row(dim,i) = f_row(dim,i) + ftmp
+                    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
@@ -533,9 +741,9 @@ contains
 
 #ifdef MPI
     !!distribute forces
-    call scatter(f_row,f,plan_row3)
+    call scatter(fRow,f,plan_row3)
 
-    call scatter(f_col,f_tmp,plan_col3)
+    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)
@@ -568,7 +776,7 @@ contains
 
   end subroutine do_lj_ff
 
-!! Calculates the potential between two lj particles, optionally returns second
+!! 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
@@ -581,8 +789,8 @@ contains
 !! Second Derivative, optional, used mainly for normal mode calculations.
     real( kind = dp ), intent(out), optional :: d2
     
-    type (lj_atype), intent(in) :: atype1
-    type (lj_atype), intent(in) :: atype2
+    type (lj_atype), intent(in), pointer :: atype1
+    type (lj_atype), intent(in), pointer :: atype2
 
     integer, intent(out), optional :: status
 
@@ -654,7 +862,7 @@ contains
   end subroutine getLjPot
 
 
-!! Calculates the mixing for sigma or epslon based on character string
+!! 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