--- trunk/OOPSE/libmdtools/do_Forces.F90	2004/01/06 18:54:57	900
+++ trunk/OOPSE/libmdtools/do_Forces.F90	2004/04/15 16:18:26	1113
@@ -4,7 +4,7 @@
 
 !! @author Charles F. Vardeman II
 !! @author Matthew Meineke
-!! @version $Id: do_Forces.F90,v 1.45 2004-01-06 18:54:57 chuckv Exp $, $Date: 2004-01-06 18:54:57 $, $Name: not supported by cvs2svn $, $Revision: 1.45 $
+!! @version $Id: do_Forces.F90,v 1.49 2004-04-15 16:18:26 tim Exp $, $Date: 2004-04-15 16:18:26 $, $Name: not supported by cvs2svn $, $Revision: 1.49 $
 
 module do_Forces
   use force_globals
@@ -15,6 +15,7 @@ module do_Forces
   use lj
   use sticky_pair
   use dipole_dipole
+  use charge_charge
   use reaction_field
   use gb_pair
   use vector_class
@@ -38,12 +39,14 @@ module do_Forces
   logical, save :: haveSaneForceField = .false.
   logical, save :: FF_uses_LJ
   logical, save :: FF_uses_sticky
+  logical, save :: FF_uses_charges
   logical, save :: FF_uses_dipoles
   logical, save :: FF_uses_RF
   logical, save :: FF_uses_GB
   logical, save :: FF_uses_EAM
   logical, save :: SIM_uses_LJ
   logical, save :: SIM_uses_sticky
+  logical, save :: SIM_uses_charges
   logical, save :: SIM_uses_dipoles
   logical, save :: SIM_uses_RF
   logical, save :: SIM_uses_GB
@@ -73,6 +76,8 @@ module do_Forces
      logical :: is_dp     = .false.
      logical :: is_gb     = .false.
      logical :: is_eam    = .false.
+     logical :: is_charge = .false.
+     real(kind=DP) :: charge = 0.0_DP
      real(kind=DP) :: dipole_moment = 0.0_DP
   end type Properties
 
@@ -114,6 +119,15 @@ contains
     do i = 1, nAtypes
        call getElementProperty(atypes, i, "is_LJ", thisProperty)
        PropertyMap(i)%is_LJ = thisProperty
+
+       call getElementProperty(atypes, i, "is_Charge", thisProperty)
+       PropertyMap(i)%is_Charge = thisProperty
+       
+       if (thisProperty) then
+          call getElementProperty(atypes, i, "charge", thisDPproperty)
+          PropertyMap(i)%charge = thisDPproperty
+       endif
+
        call getElementProperty(atypes, i, "is_DP", thisProperty)
        PropertyMap(i)%is_DP = thisProperty
 
@@ -137,6 +151,7 @@ contains
   subroutine setSimVariables()
     SIM_uses_LJ = SimUsesLJ()
     SIM_uses_sticky = SimUsesSticky()
+    SIM_uses_charges = SimUsesCharges()
     SIM_uses_dipoles = SimUsesDipoles()
     SIM_uses_RF = SimUsesRF()
     SIM_uses_GB = SimUsesGB()
@@ -236,13 +251,17 @@ contains
   
     FF_uses_LJ = .false.
     FF_uses_sticky = .false.
+    FF_uses_charges = .false.
     FF_uses_dipoles = .false.
     FF_uses_GB = .false.
     FF_uses_EAM = .false.
     
     call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
     if (nMatches .gt. 0) FF_uses_LJ = .true.
-    
+
+    call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
+    if (nMatches .gt. 0) FF_uses_charges = .true.   
+
     call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
     if (nMatches .gt. 0) FF_uses_dipoles = .true.
     
@@ -786,7 +805,7 @@ contains
     if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
        
        if (FF_uses_RF .and. SIM_uses_RF) then
-          
+ 
 #ifdef IS_MPI
           call scatter(rf_Row,rf,plan_row3d)
           call scatter(rf_Col,rf_Temp,plan_col3d)
@@ -867,38 +886,37 @@ contains
     real ( kind = dp ), intent(inout)    :: rijsq
     real ( kind = dp )                :: r
     real ( kind = dp ), intent(inout) :: d(3)
-    logical :: is_LJ_i, is_LJ_j
-    logical :: is_DP_i, is_DP_j
-    logical :: is_GB_i, is_GB_j
-    logical :: is_EAM_i,is_EAM_j
-    logical :: is_Sticky_i, is_Sticky_j
     integer :: me_i, me_j
-    integer :: propPack_i
-    integer :: propPack_j
+
     r = sqrt(rijsq)
 
 #ifdef IS_MPI
     if (tagRow(i) .eq. tagColumn(j)) then
        write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
     endif
-
     me_i = atid_row(i)
     me_j = atid_col(j)
-
 #else
-
     me_i = atid(i)
     me_j = atid(j)
-
 #endif
     
     if (FF_uses_LJ .and. SIM_uses_LJ) then
-
-       if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) &
-            call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
-
+       
+       if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
+          call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
+       endif
+       
     endif
-
+    
+    if (FF_uses_charges .and. SIM_uses_charges) then
+       
+       if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
+          call do_charge_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
+       endif
+       
+    endif
+    
     if (FF_uses_dipoles .and. SIM_uses_dipoles) then
        
        if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
@@ -907,9 +925,9 @@ contains
           if (FF_uses_RF .and. SIM_uses_RF) then
              call accumulate_rf(i, j, r, u_l)
              call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
-          endif
-          
+          endif          
        endif
+
     endif
 
     if (FF_uses_Sticky .and. SIM_uses_sticky) then
@@ -918,6 +936,7 @@ contains
           call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
                do_pot, do_stress)
        endif
+
     endif
 
 
@@ -929,16 +948,14 @@ contains
        endif
 
     endif
-     
-
-  
-   if (FF_uses_EAM .and. SIM_uses_EAM) then
-      
-      if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
-         call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
-      endif
-
-   endif
+       
+    if (FF_uses_EAM .and. SIM_uses_EAM) then
+       
+       if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
+          call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
+       endif
+       
+    endif
 
   end subroutine do_pair