--- trunk/OOPSE/libmdtools/calc_reaction_field.F90	2004/01/05 22:49:14	898
+++ trunk/OOPSE/libmdtools/calc_reaction_field.F90	2004/01/06 19:49:18	901
@@ -4,6 +4,7 @@ module reaction_field
   use atype_module
   use vector_class
   use simulation
+  use status
 #ifdef IS_MPI
   use mpiSimulation
 #endif
@@ -16,9 +17,16 @@ module reaction_field
   real(kind=dp), save :: dielect = 1.0_dp
   real(kind=dp), save :: rrfsq = 1.0_dp
   real(kind=dp), save :: pre
-  logical, save :: rf_initialized = .false., haveCuts = .false.
-  logical, save :: haveDie = .false.
+  logical, save :: haveCutoffs = .false.
+  logical, save :: haveMomentMap = .false.
+  logical, save :: haveDielectric = .false.
 
+  type :: MomentList
+     real(kind=DP) :: dipole_moment = 0.0_DP
+  end type MomentList
+
+  type(MomentList), dimension(:),allocatable :: MomentMap
+
   PUBLIC::initialize_rf
   PUBLIC::setCutoffsRF
   PUBLIC::accumulate_rf
@@ -35,8 +43,7 @@ contains
 
     pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) 
     
-    haveDie = .true.
-    if (haveCuts) rf_initialized = .true.
+    haveDielectric = .true.
 
     return
   end subroutine initialize_rf
@@ -51,12 +58,45 @@ contains
     rrfsq = rrf * rrf
     pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
     
-    haveCuts = .true.
-    if (haveDie) rf_initialized = .true.
+    haveCutoffs = .true.
 
   end subroutine setCutoffsRF
 
-  
+  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 accumulate_rf(atom1, atom2, rij, u_l)
 
     integer, intent(in) :: atom1, atom2
@@ -68,10 +108,21 @@ contains
     real (kind = dp), dimension(3) :: ul1
     real (kind = dp), dimension(3) :: ul2   
 
-    if (.not.rf_initialized) then
+    integer :: localError
+
+    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
        write(default_error,*) 'Reaction field not initialized!'
        return
     endif
+
+    if (.not.haveMomentMap) then
+       localError = 0
+       call createMomentMap(localError)
+       if ( localError .ne. 0 ) then
+          call handleError("reaction-field", "MomentMap creation failed!")
+          return
+       end if
+    endif
     
     if (rij.le.rrf) then
           
@@ -104,10 +155,9 @@ contains
        ul2(3) = u_l(3,atom2)
 #endif
        
-       call getElementProperty(atypes, me1, "dipole_moment", mu1)
-       call getElementProperty(atypes, me2, "dipole_moment", mu2)
+       mu1 = MomentMap(me1)%dipole_moment
+       mu2 = MomentMap(me2)%dipole_moment
        
-       
 #ifdef IS_MPI
        rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
        rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
@@ -153,11 +203,22 @@ contains
     real (kind = dp), dimension(3,nLocal) :: u_l    
     real (kind = dp), dimension(3,nLocal) :: t
 
-    if (.not.rf_initialized) then
+    integer :: localError
+
+    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
        write(default_error,*) 'Reaction field not initialized!'
        return
     endif
 
+    if (.not.haveMomentMap) then
+       localError = 0
+       call createMomentMap(localError)
+       if ( localError .ne. 0 ) then
+          call handleError("reaction-field", "MomentMap creation failed!")
+          return
+       end if
+    endif
+
     ! compute torques on dipoles:
     ! pre converts from mu in units of debye to kcal/mol
     
@@ -193,11 +254,22 @@ contains
     integer :: me1, me2, id1, id2
     real (kind = dp) :: mu1, mu2
     
-    if (.not.rf_initialized) then
+    integer :: localError
+
+    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
        write(default_error,*) 'Reaction field not initialized!'
        return
     endif
 
+    if (.not.haveMomentMap) then
+       localError = 0
+       call createMomentMap(localError)
+       if ( localError .ne. 0 ) then
+          call handleError("reaction-field", "MomentMap creation failed!")
+          return
+       end if
+    endif
+
     if (rij.le.rrf) then
        
        if (rij.lt.rt) then
@@ -229,8 +301,8 @@ contains
        ul2(3) = u_l(3,atom2)
 #endif
        
-       call getElementProperty(atypes, me1, "dipole_moment", mu1)
-       call getElementProperty(atypes, me2, "dipole_moment", mu2)
+       mu1 = MomentMap(me1)%dipole_moment
+       mu2 = MomentMap(me2)%dipole_moment
        
        u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)