--- trunk/OOPSE/libmdtools/calc_LJ_FF.F90	2003/04/11 15:16:59	490
+++ trunk/OOPSE/libmdtools/calc_LJ_FF.F90	2004/05/07 21:35:05	1150
@@ -2,11 +2,12 @@
 !! Corresponds to the force field defined in lj_FF.cpp 
 !! @author Charles F. Vardeman II
 !! @author Matthew Meineke
-!! @version $Id: calc_LJ_FF.F90,v 1.5 2003-04-11 15:16:59 gezelter Exp $, $Date: 2003-04-11 15:16:59 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $
+!! @version $Id: calc_LJ_FF.F90,v 1.18 2004-05-07 21:35:04 gezelter Exp $, $Date: 2004-05-07 21:35:04 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
 
 module lj
   use definitions
   use atype_module
+  use switcheroo
   use vector_class
   use simulation
 #ifdef IS_MPI
@@ -21,7 +22,9 @@ module lj
 #include "fForceField.h"
 
   integer, save :: LJ_Mixing_Policy
-  integer, save :: LJ_rcut
+  real(kind=DP), save :: LJ_rcut
+  logical, save :: havePolicy = .false., haveCut = .false.
+
   
   !! Logical has lj force field module been initialized?
   
@@ -29,7 +32,7 @@ module lj
   
   !! Public methods and data
   public :: init_LJ_FF
-  public :: LJ_new_rcut
+  public :: setCutoffLJ
   public :: do_lj_pair
   
   type :: lj_mixed_params
@@ -45,17 +48,17 @@ module lj
 
      real ( kind = dp )  :: delta  = 0.0_dp
 
+
   end type lj_mixed_params
   
   type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
   
 contains
 
-  subroutine init_LJ_FF(mix_Policy, rcut, status)
+  subroutine init_LJ_FF(mix_Policy, status)
     integer, intent(in) :: mix_Policy
     integer, intent(out) :: status
     integer :: myStatus
-    real(kind=dp) :: rcut
     
     if (mix_Policy == LB_MIXING_RULE) then
        LJ_Mixing_Policy = LB_MIXING_RULE
@@ -68,34 +71,51 @@ contains
           return
        endif
     endif
-    
-    LJ_rcut = rcut
 
-    status = 0
-    call createMixingList(myStatus)
-    if (myStatus /= 0) then
-       status = -1
-       return
+    havePolicy = .true.
+
+    if (haveCut) then
+       status = 0
+       call createMixingList(myStatus)
+       if (myStatus /= 0) then
+          status = -1
+          return
+       end if
+       
+       LJ_FF_initialized = .true.
     end if
-    
-    LJ_FF_initialized = .true.
-    
+  
   end subroutine init_LJ_FF
   
-  subroutine LJ_new_rcut(rcut, status)
+  subroutine setCutoffLJ(rcut, status)
     integer :: status, myStatus
     real(kind=dp) :: rcut
 
-    LJ_rcut = rcut
+#define __FORTRAN90
+#include "fSwitchingFunction.h"
+
     status = 0
-    call createMixingList(myStatus)
-    if (myStatus /= 0) then
-       status = -1
-       return
-    end if
 
+    LJ_rcut = rcut
+    call set_switch(LJ_SWITCH, rcut, rcut)
+
+!!$    ! ATTENTION! This is a hardwiring of rcut!
+!!$    LJ_rcut = 9.0d0
+    haveCut = .true.
+
+    if (havePolicy) then
+       status = 0
+       call createMixingList(myStatus)
+       if (myStatus /= 0) then
+          status = -1
+          return
+       end if
+       
+       LJ_FF_initialized = .true.
+    end if    
+    
     return
-  end subroutine LJ_new_rcut
+  end subroutine setCutoffLJ
   
   subroutine createMixingList(status)
     integer :: nAtypes
@@ -115,13 +135,11 @@ contains
         
     if (.not. associated(ljMixed)) then
        allocate(ljMixed(nAtypes, nAtypes))
-    else
-       status = -1
-       return
-    end if
+    endif
 
     rcut6 = LJ_rcut**6
 
+! This loops through all atypes, even those that don't support LJ forces.
     do i = 1, nAtypes
 
        call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
@@ -129,11 +147,13 @@ contains
        ! do self mixing rule
        ljMixed(i,i)%sigma   = mySigma_i
        
-       ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6                                                                                                  
-       ljMixed(i,i)%tp6     = ljMixed(i,i)%sigma6/rcut6
+       ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
 
+       ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
+
        ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
 
+
        ljMixed(i,i)%epsilon = myEpsilon_i
 
        ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
@@ -176,12 +196,13 @@ contains
     
   end subroutine createMixingList
   
-  subroutine do_lj_pair(atom1, atom2, d, rij, r2, pot, f, do_pot, do_stress)
+  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, f, &
+       do_pot, do_stress)
 
     integer, intent(in) ::  atom1, atom2
     real( kind = dp ), intent(in) :: rij, r2
-    real( kind = dp ) :: pot
-    real( kind = dp ), dimension(3,getNlocal()) :: f    
+    real( kind = dp ) :: pot, sw, vpair
+    real( kind = dp ), dimension(3,nLocal) :: f    
     real( kind = dp ), intent(in), dimension(3) :: d
     logical, intent(in) :: do_pot, do_stress
 
@@ -197,18 +218,17 @@ contains
     real( kind = dp ) :: delta
     integer :: id1, id2
 
-
     if (rij.lt.LJ_rcut)  then
 
        ! Look up the correct parameters in the mixing matrix
 #ifdef IS_MPI
        sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
        epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
-       delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
+       !delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
 #else
        sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
        epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
-       delta    = ljMixed(atid(atom1),atid(atom2))%delta
+       !delta    = ljMixed(atid(atom1),atid(atom2))%delta
 #endif
        
        r6 = r2 * r2 * r2
@@ -216,9 +236,11 @@ contains
        t6  = sigma6/ r6
        t12 = t6 * t6     
              
-       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta       
+       !pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta       
+       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) * sw
+       vpair = vpair + pot_temp
        
-       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
+       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
        
        drdx = d(1) / rij
        drdy = d(2) / rij
@@ -228,6 +250,7 @@ contains
        fy = dudr * drdy
        fz = dudr * drdz
        
+       
 #ifdef IS_MPI
        if (do_pot) then
           pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
@@ -265,21 +288,25 @@ contains
 #endif
 
           if (molMembershipList(id1) .ne. molMembershipList(id2)) then
-             tau_Temp(1) = tau_Temp(1) + fx * d(1)
-             tau_Temp(2) = tau_Temp(2) + fx * d(2)
-             tau_Temp(3) = tau_Temp(3) + fx * d(3)
-             tau_Temp(4) = tau_Temp(4) + fy * d(1)
-             tau_Temp(5) = tau_Temp(5) + fy * d(2)
-             tau_Temp(6) = tau_Temp(6) + fy * d(3)
-             tau_Temp(7) = tau_Temp(7) + fz * d(1)
-             tau_Temp(8) = tau_Temp(8) + fz * d(2)
-             tau_Temp(9) = tau_Temp(9) + fz * d(3)
+ 
+             ! because the d vector is the rj - ri vector, and 
+             ! because fx, fy, fz are the force on atom i, we need a
+             ! negative sign here:
+
+             tau_Temp(1) = tau_Temp(1) - d(1) * fx
+             tau_Temp(2) = tau_Temp(2) - d(1) * fy
+             tau_Temp(3) = tau_Temp(3) - d(1) * fz
+             tau_Temp(4) = tau_Temp(4) - d(2) * fx
+             tau_Temp(5) = tau_Temp(5) - d(2) * fy
+             tau_Temp(6) = tau_Temp(6) - d(2) * fz
+             tau_Temp(7) = tau_Temp(7) - d(3) * fx
+             tau_Temp(8) = tau_Temp(8) - d(3) * fy
+             tau_Temp(9) = tau_Temp(9) - d(3) * fz
+             
              virial_Temp = virial_Temp + &
                   (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
-          else
-             write(0,*) 'skipping ', id1, id2, molMembershipList(id1), molMembershipList(id2)
-          endif
 
+          endif
        endif
        
     endif
@@ -302,7 +329,7 @@ contains
     
     if (present(status)) status = 0
     select case (LJ_Mixing_Policy)
-    case (LB_MIXING_RULE)
+    case (1)
        select case (thisParam)
        case ("sigma")
           myMixParam = 0.5_dp * (param1 + param2)