OOPSE/libmdtools/calc_reaction_field.F90

module reaction_field
  use force_globals
  use definitions
  use atype_module
  use vector_class
  use simulation
  use status
#ifdef IS_MPI
  use mpiSimulation
#endif
  implicit none

  PRIVATE
  
  real(kind=dp), save :: rrf = 1.0_dp
  real(kind=dp), save :: rt
  real(kind=dp), save :: dielect = 1.0_dp
  real(kind=dp), save :: rrfsq = 1.0_dp
  real(kind=dp), save :: pre
  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
  PUBLIC::accumulate_self_rf
  PUBLIC::reaction_field_final
  PUBLIC::rf_correct_forces

contains
  
  subroutine initialize_rf(this_dielect)
    real(kind=dp), intent(in) :: this_dielect
    
    dielect = this_dielect

    pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) 
    
    haveDielectric = .true.

    return
  end subroutine initialize_rf

  subroutine setCutoffsRF( this_rrf, this_rt )

    real(kind=dp), intent(in) :: this_rrf, this_rt

    rrf = this_rrf
    rt  = this_rt

    rrfsq = rrf * rrf
    pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
    
    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, taper)

    integer, intent(in) :: atom1, atom2
    real (kind = dp), intent(in) :: rij
    real (kind = dp), dimension(3,nLocal) :: u_l    

    integer :: me1, me2
    real (kind = dp), intent(in) :: taper
    real (kind = dp):: mu1, mu2
    real (kind = dp), dimension(3) :: ul1
    real (kind = dp), dimension(3) :: ul2   

    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
    
       
#ifdef IS_MPI
    me1 = atid_Row(atom1)
    ul1(1) = u_l_Row(1,atom1)
    ul1(2) = u_l_Row(2,atom1)
    ul1(3) = u_l_Row(3,atom1)
    
    me2 = atid_Col(atom2)
    ul2(1) = u_l_Col(1,atom2)
    ul2(2) = u_l_Col(2,atom2)
    ul2(3) = u_l_Col(3,atom2)
#else
    me1 = atid(atom1)
    ul1(1) = u_l(1,atom1)
    ul1(2) = u_l(2,atom1)
    ul1(3) = u_l(3,atom1)
    
    me2 = atid(atom2)
    ul2(1) = u_l(1,atom2)
    ul2(2) = u_l(2,atom2)
    ul2(3) = u_l(3,atom2)
#endif
    
    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
    rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
    
    rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
    rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
    rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
#else
    rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
    rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
    rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
    
    rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
    rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
    rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper     
#endif
    
    
    return  
  end subroutine accumulate_rf

  subroutine accumulate_self_rf(atom1, mu1, u_l)
    
    integer, intent(in) :: atom1
    real(kind=dp), intent(in) :: mu1
    real(kind=dp), dimension(3,nLocal) :: u_l
    
    !! should work for both MPI and non-MPI version since this is not pairwise.
    rf(1,atom1) = rf(1,atom1) + u_l(1,atom1)*mu1
    rf(2,atom1) = rf(2,atom1) + u_l(2,atom1)*mu1
    rf(3,atom1) = rf(3,atom1) + u_l(3,atom1)*mu1
        
    return
  end subroutine accumulate_self_rf
  
  subroutine reaction_field_final(a1, mu1, u_l, rfpot, t, do_pot)
            
    integer, intent(in) :: a1
    real (kind=dp), intent(in) :: mu1
    real (kind=dp), intent(inout) :: rfpot
    logical, intent(in) :: do_pot
    real (kind = dp), dimension(3,nLocal) :: u_l    
    real (kind = dp), dimension(3,nLocal) :: t

    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
    
    ! The torque contribution is dipole cross reaction_field   

    t(1,a1) = t(1,a1) + pre*mu1*(u_l(2,a1)*rf(3,a1) - u_l(3,a1)*rf(2,a1))
    t(2,a1) = t(2,a1) + pre*mu1*(u_l(3,a1)*rf(1,a1) - u_l(1,a1)*rf(3,a1))
    t(3,a1) = t(3,a1) + pre*mu1*(u_l(1,a1)*rf(2,a1) - u_l(2,a1)*rf(1,a1))
    
    ! the potential contribution is -1/2 dipole dot reaction_field
    
    if (do_pot) then
       rfpot = rfpot - 0.5d0 * pre * mu1 * &
            (rf(1,a1)*u_l(1,a1) + rf(2,a1)*u_l(2,a1) + rf(3,a1)*u_l(3,a1))
    endif
 
    return
  end subroutine reaction_field_final
  
  subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, taper, f, fpair)
    
    integer, intent(in) :: atom1, atom2
    real(kind=dp), dimension(3), intent(in) :: d
    real(kind=dp), intent(in) :: rij, taper
    real( kind = dp ), dimension(3,nLocal) :: u_l
    real( kind = dp ), dimension(3,nLocal) :: f
    real( kind = dp ), dimension(3), intent(inout) :: fpair
    
    real (kind = dp), dimension(3) :: ul1
    real (kind = dp), dimension(3) :: ul2
    real (kind = dp) :: dtdr
    real (kind = dp) :: dudx, dudy, dudz, u1dotu2
    integer :: me1, me2, id1, id2
    real (kind = dp) :: mu1, mu2
    
    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
          dtdr = 0.0d0
       else
 !         write(*,*) 'rf correct in taper region'
          dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
       endif
       
#ifdef IS_MPI
       me1 = atid_Row(atom1)
       ul1(1) = u_l_Row(1,atom1)
       ul1(2) = u_l_Row(2,atom1)
       ul1(3) = u_l_Row(3,atom1)
       
       me2 = atid_Col(atom2)
       ul2(1) = u_l_Col(1,atom2)
       ul2(2) = u_l_Col(2,atom2)
       ul2(3) = u_l_Col(3,atom2)
#else
       me1 = atid(atom1)
       ul1(1) = u_l(1,atom1)
       ul1(2) = u_l(2,atom1)
       ul1(3) = u_l(3,atom1)
       
       me2 = atid(atom2)
       ul2(1) = u_l(1,atom2)
       ul2(2) = u_l(2,atom2)
       ul2(3) = u_l(3,atom2)
#endif
       
       mu1 = MomentMap(me1)%dipole_moment
       mu2 = MomentMap(me2)%dipole_moment
       
       u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
       
       dudx = - pre*mu1*mu2*u1dotu2*dtdr*d(1)/rij
       dudy = - pre*mu1*mu2*u1dotu2*dtdr*d(2)/rij
       dudz = - pre*mu1*mu2*u1dotu2*dtdr*d(3)/rij
       
#ifdef IS_MPI
       f_Row(1,atom1) = f_Row(1,atom1) + dudx
       f_Row(2,atom1) = f_Row(2,atom1) + dudy
       f_Row(3,atom1) = f_Row(3,atom1) + dudz
       
       f_Col(1,atom2) = f_Col(1,atom2) - dudx
       f_Col(2,atom2) = f_Col(2,atom2) - dudy
       f_Col(3,atom2) = f_Col(3,atom2) - dudz
#else
       f(1,atom1) = f(1,atom1) + dudx
       f(2,atom1) = f(2,atom1) + dudy
       f(3,atom1) = f(3,atom1) + dudz
       
       f(1,atom2) = f(1,atom2) - dudx
       f(2,atom2) = f(2,atom2) - dudy
       f(3,atom2) = f(3,atom2) - dudz
#endif

#ifdef IS_MPI
       id1 = AtomRowToGlobal(atom1)
       id2 = AtomColToGlobal(atom2)
#else
       id1 = atom1
       id2 = atom2
#endif
       
       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
          
          fpair(1) = fpair(1) + dudx
          fpair(2) = fpair(2) + dudy
          fpair(3) = fpair(3) + dudz
          
       endif
       
    end if
    return
  end subroutine rf_correct_forces
end module reaction_field
Revision:	1198
Committed:	Thu May 27 00:48:12 2004 UTC (20 years, 1 month ago) by tim
File size:	9214 byte(s)
Log Message:	in the progress of fixing MPI version of cutoff group
#	User	Rev	Content
1	mmeineke	377	module reaction_field
2			use force_globals
3			use definitions
4			use atype_module
5			use vector_class
6	chuckv	460	use simulation
7	chuckv	901	use status
8	mmeineke	377	#ifdef IS_MPI
9			use mpiSimulation
10			#endif
11			implicit none
12
13	mmeineke	626	PRIVATE
14
15			real(kind=dp), save :: rrf = 1.0_dp
16	mmeineke	377	real(kind=dp), save :: rt
17	mmeineke	626	real(kind=dp), save :: dielect = 1.0_dp
18			real(kind=dp), save :: rrfsq = 1.0_dp
19	mmeineke	377	real(kind=dp), save :: pre
20	chuckv	901	logical, save :: haveCutoffs = .false.
21			logical, save :: haveMomentMap = .false.
22			logical, save :: haveDielectric = .false.
23	mmeineke	377
24	chuckv	901	type :: MomentList
25			real(kind=DP) :: dipole_moment = 0.0_DP
26			end type MomentList
27
28			type(MomentList), dimension(:),allocatable :: MomentMap
29
30	mmeineke	626	PUBLIC::initialize_rf
31			PUBLIC::setCutoffsRF
32			PUBLIC::accumulate_rf
33			PUBLIC::accumulate_self_rf
34			PUBLIC::reaction_field_final
35			PUBLIC::rf_correct_forces
36
37	mmeineke	377	contains
38
39	mmeineke	626	subroutine initialize_rf(this_dielect)
40			real(kind=dp), intent(in) :: this_dielect
41
42			dielect = this_dielect
43	mmeineke	377
44	mmeineke	626	pre = 14.38362d02.0d0(dielect-1.0d0)/((2.0d0dielect+1.0d0)rrfsq*rrf)
45
46	chuckv	901	haveDielectric = .true.
47	mmeineke	626
48			return
49			end subroutine initialize_rf
50
51			subroutine setCutoffsRF( this_rrf, this_rt )
52
53			real(kind=dp), intent(in) :: this_rrf, this_rt
54
55	mmeineke	377	rrf = this_rrf
56	mmeineke	626	rt = this_rt
57
58	mmeineke	377	rrfsq = rrf * rrf
59			pre = 14.38362d02.0d0(dielect-1.0d0)/((2.0d0dielect+1.0d0)rrfsq*rrf)
60
61	chuckv	901	haveCutoffs = .true.
62	gezelter	394
63	mmeineke	626	end subroutine setCutoffsRF
64
65	chuckv	901	subroutine createMomentMap(status)
66			integer :: nAtypes
67			integer :: status
68			integer :: i
69			real (kind=DP) :: thisDP
70			logical :: thisProperty
71
72			status = 0
73
74			nAtypes = getSize(atypes)
75
76			if (nAtypes == 0) then
77			status = -1
78			return
79			end if
80
81			if (.not. allocated(MomentMap)) then
82			allocate(MomentMap(nAtypes))
83			endif
84
85			do i = 1, nAtypes
86
87			call getElementProperty(atypes, i, "is_DP", thisProperty)
88
89			if (thisProperty) then
90			call getElementProperty(atypes, i, "dipole_moment", thisDP)
91			MomentMap(i)%dipole_moment = thisDP
92			endif
93
94			end do
95
96			haveMomentMap = .true.
97
98			end subroutine createMomentMap
99
100	gezelter	1160	subroutine accumulate_rf(atom1, atom2, rij, u_l, taper)
101	mmeineke	377
102			integer, intent(in) :: atom1, atom2
103			real (kind = dp), intent(in) :: rij
104	chuckv	898	real (kind = dp), dimension(3,nLocal) :: u_l
105	mmeineke	377
106			integer :: me1, me2
107	gezelter	1160	real (kind = dp), intent(in) :: taper
108			real (kind = dp):: mu1, mu2
109	mmeineke	377	real (kind = dp), dimension(3) :: ul1
110			real (kind = dp), dimension(3) :: ul2
111
112	chuckv	901	integer :: localError
113
114			if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
115	mmeineke	377	write(default_error,*) 'Reaction field not initialized!'
116			return
117			endif
118	chuckv	901
119			if (.not.haveMomentMap) then
120			localError = 0
121			call createMomentMap(localError)
122			if ( localError .ne. 0 ) then
123			call handleError("reaction-field", "MomentMap creation failed!")
124			return
125			end if
126			endif
127	mmeineke	377
128
129			#ifdef IS_MPI
130	gezelter	1160	me1 = atid_Row(atom1)
131			ul1(1) = u_l_Row(1,atom1)
132			ul1(2) = u_l_Row(2,atom1)
133			ul1(3) = u_l_Row(3,atom1)
134
135			me2 = atid_Col(atom2)
136			ul2(1) = u_l_Col(1,atom2)
137			ul2(2) = u_l_Col(2,atom2)
138			ul2(3) = u_l_Col(3,atom2)
139	mmeineke	377	#else
140	gezelter	1160	me1 = atid(atom1)
141			ul1(1) = u_l(1,atom1)
142			ul1(2) = u_l(2,atom1)
143			ul1(3) = u_l(3,atom1)
144
145			me2 = atid(atom2)
146			ul2(1) = u_l(1,atom2)
147			ul2(2) = u_l(2,atom2)
148			ul2(3) = u_l(3,atom2)
149	mmeineke	377	#endif
150	gezelter	1160
151			mu1 = MomentMap(me1)%dipole_moment
152			mu2 = MomentMap(me2)%dipole_moment
153
154	mmeineke	377	#ifdef IS_MPI
155	gezelter	1160	rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)mu2taper
156			rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)mu2taper
157			rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)mu2taper
158
159			rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)mu1taper
160			rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)mu1taper
161			rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)mu1taper
162	mmeineke	377	#else
163	gezelter	1160	rf(1,atom1) = rf(1,atom1) + ul2(1)mu2taper
164			rf(2,atom1) = rf(2,atom1) + ul2(2)mu2taper
165			rf(3,atom1) = rf(3,atom1) + ul2(3)mu2taper
166
167			rf(1,atom2) = rf(1,atom2) + ul1(1)mu1taper
168			rf(2,atom2) = rf(2,atom2) + ul1(2)mu1taper
169			rf(3,atom2) = rf(3,atom2) + ul1(3)mu1taper
170	mmeineke	377	#endif
171	gezelter	1160
172
173	mmeineke	377	return
174			end subroutine accumulate_rf
175
176			subroutine accumulate_self_rf(atom1, mu1, u_l)
177
178			integer, intent(in) :: atom1
179			real(kind=dp), intent(in) :: mu1
180	chuckv	898	real(kind=dp), dimension(3,nLocal) :: u_l
181	mmeineke	377
182			!! should work for both MPI and non-MPI version since this is not pairwise.
183			rf(1,atom1) = rf(1,atom1) + u_l(1,atom1)*mu1
184			rf(2,atom1) = rf(2,atom1) + u_l(2,atom1)*mu1
185			rf(3,atom1) = rf(3,atom1) + u_l(3,atom1)*mu1
186
187			return
188			end subroutine accumulate_self_rf
189
190			subroutine reaction_field_final(a1, mu1, u_l, rfpot, t, do_pot)
191
192			integer, intent(in) :: a1
193			real (kind=dp), intent(in) :: mu1
194			real (kind=dp), intent(inout) :: rfpot
195			logical, intent(in) :: do_pot
196	chuckv	898	real (kind = dp), dimension(3,nLocal) :: u_l
197			real (kind = dp), dimension(3,nLocal) :: t
198	mmeineke	377
199	chuckv	901	integer :: localError
200
201			if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
202	mmeineke	377	write(default_error,*) 'Reaction field not initialized!'
203			return
204			endif
205
206	chuckv	901	if (.not.haveMomentMap) then
207			localError = 0
208			call createMomentMap(localError)
209			if ( localError .ne. 0 ) then
210			call handleError("reaction-field", "MomentMap creation failed!")
211			return
212			end if
213			endif
214
215	mmeineke	377	! compute torques on dipoles:
216			! pre converts from mu in units of debye to kcal/mol
217
218	gezelter	394	! The torque contribution is dipole cross reaction_field
219
220	mmeineke	377	t(1,a1) = t(1,a1) + premu1(u_l(2,a1)rf(3,a1) - u_l(3,a1)rf(2,a1))
221			t(2,a1) = t(2,a1) + premu1(u_l(3,a1)rf(1,a1) - u_l(1,a1)rf(3,a1))
222			t(3,a1) = t(3,a1) + premu1(u_l(1,a1)rf(2,a1) - u_l(2,a1)rf(1,a1))
223
224			! the potential contribution is -1/2 dipole dot reaction_field
225
226			if (do_pot) then
227			rfpot = rfpot - 0.5d0 * pre * mu1 * &
228			(rf(1,a1)u_l(1,a1) + rf(2,a1)u_l(2,a1) + rf(3,a1)*u_l(3,a1))
229			endif
230	chrisfen	999
231	mmeineke	377	return
232			end subroutine reaction_field_final
233
234	gezelter	1192	subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, taper, f, fpair)
235	mmeineke	377
236			integer, intent(in) :: atom1, atom2
237			real(kind=dp), dimension(3), intent(in) :: d
238	gezelter	1160	real(kind=dp), intent(in) :: rij, taper
239	chuckv	898	real( kind = dp ), dimension(3,nLocal) :: u_l
240			real( kind = dp ), dimension(3,nLocal) :: f
241	gezelter	1192	real( kind = dp ), dimension(3), intent(inout) :: fpair
242	mmeineke	377
243			real (kind = dp), dimension(3) :: ul1
244			real (kind = dp), dimension(3) :: ul2
245			real (kind = dp) :: dtdr
246			real (kind = dp) :: dudx, dudy, dudz, u1dotu2
247	gezelter	730	integer :: me1, me2, id1, id2
248	mmeineke	377	real (kind = dp) :: mu1, mu2
249
250	chuckv	901	integer :: localError
251
252			if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
253	mmeineke	377	write(default_error,*) 'Reaction field not initialized!'
254			return
255			endif
256
257	chuckv	901	if (.not.haveMomentMap) then
258			localError = 0
259			call createMomentMap(localError)
260			if ( localError .ne. 0 ) then
261			call handleError("reaction-field", "MomentMap creation failed!")
262			return
263			end if
264			endif
265
266	mmeineke	377	if (rij.le.rrf) then
267
268			if (rij.lt.rt) then
269			dtdr = 0.0d0
270			else
271	chrisfen	999	! write(,) 'rf correct in taper region'
272	mmeineke	377	dtdr = 6.0d0(rijrij - rijrt - rijrrf +rrfrt)/((rrf-rt)*3)
273			endif
274
275			#ifdef IS_MPI
276			me1 = atid_Row(atom1)
277			ul1(1) = u_l_Row(1,atom1)
278			ul1(2) = u_l_Row(2,atom1)
279			ul1(3) = u_l_Row(3,atom1)
280
281			me2 = atid_Col(atom2)
282			ul2(1) = u_l_Col(1,atom2)
283			ul2(2) = u_l_Col(2,atom2)
284			ul2(3) = u_l_Col(3,atom2)
285			#else
286			me1 = atid(atom1)
287			ul1(1) = u_l(1,atom1)
288			ul1(2) = u_l(2,atom1)
289			ul1(3) = u_l(3,atom1)
290
291			me2 = atid(atom2)
292			ul2(1) = u_l(1,atom2)
293			ul2(2) = u_l(2,atom2)
294			ul2(3) = u_l(3,atom2)
295			#endif
296
297	chuckv	901	mu1 = MomentMap(me1)%dipole_moment
298			mu2 = MomentMap(me2)%dipole_moment
299	mmeineke	377
300			u1dotu2 = ul1(1)ul2(1) + ul1(2)ul2(2) + ul1(3)*ul2(3)
301
302			dudx = - premu1mu2u1dotu2dtdr*d(1)/rij
303			dudy = - premu1mu2u1dotu2dtdr*d(2)/rij
304			dudz = - premu1mu2u1dotu2dtdr*d(3)/rij
305
306			#ifdef IS_MPI
307			f_Row(1,atom1) = f_Row(1,atom1) + dudx
308			f_Row(2,atom1) = f_Row(2,atom1) + dudy
309			f_Row(3,atom1) = f_Row(3,atom1) + dudz
310
311			f_Col(1,atom2) = f_Col(1,atom2) - dudx
312			f_Col(2,atom2) = f_Col(2,atom2) - dudy
313			f_Col(3,atom2) = f_Col(3,atom2) - dudz
314			#else
315			f(1,atom1) = f(1,atom1) + dudx
316			f(2,atom1) = f(2,atom1) + dudy
317			f(3,atom1) = f(3,atom1) + dudz
318
319			f(1,atom2) = f(1,atom2) - dudx
320			f(2,atom2) = f(2,atom2) - dudy
321			f(3,atom2) = f(3,atom2) - dudz
322			#endif
323	gezelter	611
324	gezelter	730	#ifdef IS_MPI
325	tim	1198	id1 = AtomRowToGlobal(atom1)
326			id2 = AtomColToGlobal(atom2)
327	gezelter	730	#else
328	gezelter	1192	id1 = atom1
329			id2 = atom2
330	gezelter	730	#endif
331	gezelter	1192
332			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
333
334			fpair(1) = fpair(1) + dudx
335			fpair(2) = fpair(2) + dudy
336			fpair(3) = fpair(3) + dudz
337
338			endif
339
340			end if
341	mmeineke	377	return
342			end subroutine rf_correct_forces
343			end module reaction_field