--- trunk/OOPSE/libmdtools/calc_dipole_dipole.F90	2003/07/15 17:10:50	611
+++ trunk/OOPSE/libmdtools/calc_dipole_dipole.F90	2004/01/15 01:14:55	945
@@ -5,39 +5,86 @@ module dipole_dipole
   use atype_module
   use vector_class
   use simulation
+  use status
 #ifdef IS_MPI
   use mpiSimulation
 #endif
   implicit none
 
-  real(kind=dp), save :: rrf = 0.0
+  PRIVATE
+  real(kind=dp), save :: ecr = 0.0
   real(kind=dp), save :: rt  = 0.0
-  real(kind=dp), save :: rrfsq = 0.0
-  real(kind=dp), save :: pre = 0.0
-  logical, save :: dipole_initialized = .false.
+   real(kind=dp), save :: pre = 0.0
+  logical, save :: haveCutoffs = .false.
+  logical, save :: haveMomentMap = .false.
 
+  public::setCutoffsDipole
+  public::do_dipole_pair
+
+  type :: MomentList
+     real(kind=DP) :: dipole_moment = 0.0_DP
+  end type MomentList
+
+  type(MomentList), dimension(:),allocatable :: MomentMap
+
 contains
     
-  subroutine initialize_dipole(this_rrf, this_rt)
-    real(kind=dp), intent(in) :: this_rrf, this_rt
-    rrf = this_rrf
+  subroutine setCutoffsDipole(this_ecr, this_rt)
+    real(kind=dp), intent(in) :: this_ecr, this_rt
+    ecr = this_ecr
     rt = this_rt    
-    rrfsq = rrf * rrf
-    ! pre converts from mu in units of debye to kcal/mol
+
+      ! pre converts from mu in units of debye to kcal/mol
     pre = 14.38362_dp
 
-    dipole_initialized = .true.
+    haveCutoffs = .true.
     
     return
-  end subroutine initialize_dipole
+  end subroutine setCutoffsDipole
 
+  subroutine createMomentMap(status)
+    integer :: nAtypes
+    integer :: status
+    integer :: i
+    real (kind=DP) :: thisDP
+    logical :: thisProperty
 
+    status = 0
+
+    nAtypes = getSize(atypes)
+    
+    if (nAtypes == 0) then
+       status = -1
+       return
+    end if
+    
+    if (.not. allocated(MomentMap)) then
+       allocate(MomentMap(nAtypes))
+    endif
+
+    do i = 1, nAtypes
+
+       call getElementProperty(atypes, i, "is_DP", thisProperty)
+
+       if (thisProperty) then
+          call getElementProperty(atypes, i, "dipole_moment", thisDP)
+          MomentMap(i)%dipole_moment = thisDP
+       endif
+       
+    end do
+    
+    haveMomentMap = .true.
+    
+  end subroutine createMomentMap
+  
+  
   subroutine do_dipole_pair(atom1, atom2, d, rij, r2, pot, u_l, f, t, &
        do_pot, do_stress)
     
     logical :: do_pot, do_stress
 
-    integer atom1, atom2, me1, me2
+    integer atom1, atom2, me1, me2, id1, id2
+    integer :: localError
     real(kind=dp) :: rij, mu1, mu2
     real(kind=dp) :: dfact1, dfact2, dip2, r2, r3, r5
     real(kind=dp) :: dudx, dudy, dudz, dudu1x, dudu1y, dudu1z
@@ -46,18 +93,27 @@ contains
 
     real( kind = dp ) :: pot
     real( kind = dp ), dimension(3) :: d
-    real( kind = dp ), dimension(3,getNlocal()) :: u_l
-    real( kind = dp ), dimension(3,getNlocal()) :: f
-    real( kind = dp ), dimension(3,getNlocal()) :: t
+    real( kind = dp ), dimension(3,nLocal) :: u_l
+    real( kind = dp ), dimension(3,nLocal) :: f
+    real( kind = dp ), dimension(3,nLocal) :: t
     
     real (kind = dp), dimension(3) :: ul1
     real (kind = dp), dimension(3) :: ul2
 
-    if (.not.dipole_initialized) then
-       write(default_error,*) 'Dipole-dipole not initialized!'
+    if (.not. haveCutoffs) then
+       write(default_error,*) 'Dipole-dipole does not have cutoffs set!'
        return
     endif
-    
+
+    if (.not.haveMomentMap) then
+       localError = 0
+       call createMomentMap(localError)
+       if ( localError .ne. 0 ) then
+          call handleError("dipole-dipole", "MomentMap creation failed!")
+          return
+       end if
+    endif
+
 #ifdef IS_MPI
     me1 = atid_Row(atom1)
     ul1(1) = u_l_Row(1,atom1)
@@ -80,33 +136,32 @@ contains
     ul2(3) = u_l(3,atom2)
 #endif
 
+    mu1 = MomentMap(me1)%dipole_moment
+    mu2 = MomentMap(me2)%dipole_moment
 
-    call getElementProperty(atypes, me1, "dipole_moment", mu1)
-    call getElementProperty(atypes, me2, "dipole_moment", mu2)
-
-    if (rij.le.rrf) then
+    if (rij.le.ecr) then
        
        if (rij.lt.rt) then
           taper = 1.0d0
           dtdr = 0.0d0
        else
-          taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3)
-          dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
+          taper = (ecr + 2.0d0*rij - 3.0d0*rt)*(ecr-rij)**2/ ((ecr-rt)**3)
+          dtdr = 6.0d0*(rij*rij - rij*rt - rij*ecr +ecr*rt)/((ecr-rt)**3)
        endif
-
+       
        r3 = r2*rij
        r5 = r3*r2
        
        rdotu1 = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
        rdotu2 = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
        u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
-
+       
        dip2 = pre * mu1 * mu2
        dfact1 = 3.0d0*dip2 / r2
        dfact2 = 3.0d0*dip2 / r5
-
+       
        vterm = dip2*((u1dotu2/r3) - 3.0d0*(rdotu1*rdotu2/r5))
-
+       
        if (do_pot) then
 #ifdef IS_MPI 
           pot_row(atom1) = pot_row(atom1) + 0.5d0*vterm*taper
@@ -173,10 +228,17 @@ contains
        t(2,atom2) = t(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
        t(3,atom2) = t(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
 #endif
-       
+
        if (do_stress) then          
 
-          if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then
+#ifdef IS_MPI
+          id1 = tagRow(atom1)
+          id2 = tagColumn(atom2)
+#else
+          id1 = atom1
+          id2 = atom2
+#endif                 
+          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
 
              ! because the d vector is the rj - ri vector, and 
              ! because dudx, dudy, dudz are the (positive) force on