--- trunk/OOPSE/libmdtools/calc_dipole_dipole.F90	2004/01/05 22:49:14	898
+++ trunk/OOPSE/libmdtools/calc_dipole_dipole.F90	2004/01/06 19:49:18	901
@@ -5,6 +5,7 @@ module dipole_dipole
   use atype_module
   use vector_class
   use simulation
+  use status
 #ifdef IS_MPI
   use mpiSimulation
 #endif
@@ -14,12 +15,18 @@ module dipole_dipole
   real(kind=dp), save :: rrf = 0.0
   real(kind=dp), save :: rt  = 0.0
    real(kind=dp), save :: pre = 0.0
-  logical, save :: dipole_initialized = .false.
+  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 setCutoffsDipole(this_rrf, this_rt)
@@ -30,17 +37,54 @@ contains
       ! pre converts from mu in units of debye to kcal/mol
     pre = 14.38362_dp
 
-    dipole_initialized = .true.
+    haveCutoffs = .true.
     
     return
   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, 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
@@ -56,11 +100,20 @@ contains
     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)
@@ -83,10 +136,9 @@ 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.lt.rt) then
@@ -96,20 +148,20 @@ contains
           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)
        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