UseTheForce/DarkSide/reactionField.F90

!!
!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
!!
!! The University of Notre Dame grants you ("Licensee") a
!! non-exclusive, royalty free, license to use, modify and
!! redistribute this software in source and binary code form, provided
!! that the following conditions are met:
!!
!! 1. Acknowledgement of the program authors must be made in any
!!    publication of scientific results based in part on use of the
!!    program.  An acceptable form of acknowledgement is citation of
!!    the article in which the program was described (Matthew
!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
!!    Parallel Simulation Engine for Molecular Dynamics,"
!!    J. Comput. Chem. 26, pp. 252-271 (2005))
!!
!! 2. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 3. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer in the
!!    documentation and/or other materials provided with the
!!    distribution.
!!
!! This software is provided "AS IS," without a warranty of any
!! kind. All express or implied conditions, representations and
!! warranties, including any implied warranty of merchantability,
!! fitness for a particular purpose or non-infringement, are hereby
!! excluded.  The University of Notre Dame and its licensors shall not
!! be liable for any damages suffered by licensee as a result of
!! using, modifying or distributing the software or its
!! derivatives. In no event will the University of Notre Dame or its
!! licensors be liable for any lost revenue, profit or data, or for
!! direct, indirect, special, consequential, incidental or punitive
!! damages, however caused and regardless of the theory of liability,
!! arising out of the use of or inability to use software, even if the
!! University of Notre Dame has been advised of the possibility of
!! such damages.
!!

module reaction_field
  use force_globals
  use definitions
  use atype_module
  use vector_class
  use simulation
  use status
  use electrostatic_module
#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 :: haveDielectric = .false.

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

    haveCutoffs = .true.

  end subroutine setCutoffsRF

  subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper)

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

    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

#ifdef IS_MPI
    me1 = atid_Row(atom1)
    ul1(1) = eFrame_Row(3,atom1)
    ul1(2) = eFrame_Row(6,atom1)
    ul1(3) = eFrame_Row(9,atom1)

    me2 = atid_Col(atom2)
    ul2(1) = eFrame_Col(3,atom2)
    ul2(2) = eFrame_Col(6,atom2)
    ul2(3) = eFrame_Col(9,atom2)
#else
    me1 = atid(atom1)
    ul1(1) = eFrame(3,atom1)
    ul1(2) = eFrame(6,atom1)
    ul1(3) = eFrame(9,atom1)

    me2 = atid(atom2)
    ul2(1) = eFrame(3,atom2)
    ul2(2) = eFrame(6,atom2)
    ul2(3) = eFrame(9,atom2)
#endif

    mu1 = getDipoleMoment(me1)
    mu2 = getDipoleMoment(me2)

#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, eFrame)

    integer, intent(in) :: atom1
    real(kind=dp), intent(in) :: mu1
    real(kind=dp), dimension(9,nLocal) :: eFrame

    !! should work for both MPI and non-MPI version since this is not pairwise.
    rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
    rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
    rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1

    return
  end subroutine accumulate_self_rf

  subroutine reaction_field_final(a1, mu1, eFrame, 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(9,nLocal) :: eFrame
    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

    ! 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*(eFrame(6,a1)*rf(3,a1) - eFrame(9,a1)*rf(2,a1))
    t(2,a1) = t(2,a1) + pre*mu1*(eFrame(9,a1)*rf(1,a1) - eFrame(3,a1)*rf(3,a1))
    t(3,a1) = t(3,a1) + pre*mu1*(eFrame(3,a1)*rf(2,a1) - eFrame(6,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)*eFrame(3,a1) + rf(2,a1)*eFrame(6,a1) + rf(3,a1)*eFrame(9,a1))
    endif

    return
  end subroutine reaction_field_final

  subroutine rf_correct_forces(atom1, atom2, d, rij, eFrame, 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(9,nLocal) :: eFrame
    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 (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) = eFrame_Row(3,atom1)
       ul1(2) = eFrame_Row(6,atom1)
       ul1(3) = eFrame_Row(9,atom1)

       me2 = atid_Col(atom2)
       ul2(1) = eFrame_Col(3,atom2)
       ul2(2) = eFrame_Col(6,atom2)
       ul2(3) = eFrame_Col(9,atom2)
#else
       me1 = atid(atom1)
       ul1(1) = eFrame(3,atom1)
       ul1(2) = eFrame(6,atom1)
       ul1(3) = eFrame(9,atom1)

       me2 = atid(atom2)
       ul2(1) = eFrame(3,atom2)
       ul2(2) = eFrame(6,atom2)
       ul2(3) = eFrame(9,atom2)
#endif

       mu1 = getDipoleMoment(me1)
       mu2 = getDipoleMoment(me2)

       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:	2204
Committed:	Fri Apr 15 22:04:00 2005 UTC (19 years, 4 months ago) by gezelter
File size:	9540 byte(s)
Log Message:	xemacs has been drafted to perform our indentation services
#	User	Rev	Content
1	gezelter	1930	!!
2			!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			!!
4			!! The University of Notre Dame grants you ("Licensee") a
5			!! non-exclusive, royalty free, license to use, modify and
6			!! redistribute this software in source and binary code form, provided
7			!! that the following conditions are met:
8			!!
9			!! 1. Acknowledgement of the program authors must be made in any
10			!! publication of scientific results based in part on use of the
11			!! program. An acceptable form of acknowledgement is citation of
12			!! the article in which the program was described (Matthew
13			!! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			!! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			!! Parallel Simulation Engine for Molecular Dynamics,"
16			!! J. Comput. Chem. 26, pp. 252-271 (2005))
17			!!
18			!! 2. Redistributions of source code must retain the above copyright
19			!! notice, this list of conditions and the following disclaimer.
20			!!
21			!! 3. Redistributions in binary form must reproduce the above copyright
22			!! notice, this list of conditions and the following disclaimer in the
23			!! documentation and/or other materials provided with the
24			!! distribution.
25			!!
26			!! This software is provided "AS IS," without a warranty of any
27			!! kind. All express or implied conditions, representations and
28			!! warranties, including any implied warranty of merchantability,
29			!! fitness for a particular purpose or non-infringement, are hereby
30			!! excluded. The University of Notre Dame and its licensors shall not
31			!! be liable for any damages suffered by licensee as a result of
32			!! using, modifying or distributing the software or its
33			!! derivatives. In no event will the University of Notre Dame or its
34			!! licensors be liable for any lost revenue, profit or data, or for
35			!! direct, indirect, special, consequential, incidental or punitive
36			!! damages, however caused and regardless of the theory of liability,
37			!! arising out of the use of or inability to use software, even if the
38			!! University of Notre Dame has been advised of the possibility of
39			!! such damages.
40			!!
41
42	gezelter	1608	module reaction_field
43			use force_globals
44			use definitions
45			use atype_module
46			use vector_class
47			use simulation
48			use status
49	chrisfen	2129	use electrostatic_module
50	gezelter	1608	#ifdef IS_MPI
51			use mpiSimulation
52			#endif
53			implicit none
54
55			PRIVATE
56	gezelter	2204
57	gezelter	1608	real(kind=dp), save :: rrf = 1.0_dp
58			real(kind=dp), save :: rt
59			real(kind=dp), save :: dielect = 1.0_dp
60			real(kind=dp), save :: rrfsq = 1.0_dp
61			real(kind=dp), save :: pre
62	chrisfen	2129
63	gezelter	1608	logical, save :: haveCutoffs = .false.
64			logical, save :: haveDielectric = .false.
65
66			PUBLIC::initialize_rf
67			PUBLIC::setCutoffsRF
68			PUBLIC::accumulate_rf
69			PUBLIC::accumulate_self_rf
70			PUBLIC::reaction_field_final
71			PUBLIC::rf_correct_forces
72
73			contains
74	gezelter	2204
75	gezelter	1608	subroutine initialize_rf(this_dielect)
76			real(kind=dp), intent(in) :: this_dielect
77	gezelter	2204
78	gezelter	1608	dielect = this_dielect
79
80	chrisfen	2129	pre = pre22 * 2.0d0(dielect-1.0d0)/((2.0d0dielect+1.0d0)rrfsqrrf)
81	gezelter	2204
82	gezelter	1608	haveDielectric = .true.
83
84			return
85			end subroutine initialize_rf
86
87			subroutine setCutoffsRF( this_rrf, this_rt )
88
89			real(kind=dp), intent(in) :: this_rrf, this_rt
90
91			rrf = this_rrf
92			rt = this_rt
93
94			rrfsq = rrf * rrf
95	chrisfen	2129	pre = pre22 * 2.0d0(dielect-1.0d0)/((2.0d0dielect+1.0d0)rrfsqrrf)
96	gezelter	2204
97	gezelter	1608	haveCutoffs = .true.
98
99			end subroutine setCutoffsRF
100
101	gezelter	1930	subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper)
102	gezelter	1608
103			integer, intent(in) :: atom1, atom2
104			real (kind = dp), intent(in) :: rij
105	gezelter	1930	real (kind = dp), dimension(9,nLocal) :: eFrame
106	gezelter	1608
107			integer :: me1, me2
108			real (kind = dp), intent(in) :: taper
109			real (kind = dp):: mu1, mu2
110			real (kind = dp), dimension(3) :: ul1
111			real (kind = dp), dimension(3) :: ul2
112
113			integer :: localError
114
115			if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
116			write(default_error,*) 'Reaction field not initialized!'
117			return
118	gezelter	2204	endif
119
120	gezelter	1608	#ifdef IS_MPI
121			me1 = atid_Row(atom1)
122	gezelter	1930	ul1(1) = eFrame_Row(3,atom1)
123			ul1(2) = eFrame_Row(6,atom1)
124			ul1(3) = eFrame_Row(9,atom1)
125	gezelter	2204
126	gezelter	1608	me2 = atid_Col(atom2)
127	gezelter	1930	ul2(1) = eFrame_Col(3,atom2)
128			ul2(2) = eFrame_Col(6,atom2)
129			ul2(3) = eFrame_Col(9,atom2)
130	gezelter	1608	#else
131			me1 = atid(atom1)
132	gezelter	1930	ul1(1) = eFrame(3,atom1)
133			ul1(2) = eFrame(6,atom1)
134			ul1(3) = eFrame(9,atom1)
135	gezelter	2204
136	gezelter	1608	me2 = atid(atom2)
137	gezelter	1930	ul2(1) = eFrame(3,atom2)
138			ul2(2) = eFrame(6,atom2)
139			ul2(3) = eFrame(9,atom2)
140	gezelter	1608	#endif
141	gezelter	2204
142	chrisfen	2129	mu1 = getDipoleMoment(me1)
143			mu2 = getDipoleMoment(me2)
144	gezelter	2204
145	gezelter	1608	#ifdef IS_MPI
146			rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)mu2taper
147			rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)mu2taper
148			rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)mu2taper
149	gezelter	2204
150	gezelter	1608	rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)mu1taper
151			rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)mu1taper
152			rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)mu1taper
153			#else
154			rf(1,atom1) = rf(1,atom1) + ul2(1)mu2taper
155			rf(2,atom1) = rf(2,atom1) + ul2(2)mu2taper
156			rf(3,atom1) = rf(3,atom1) + ul2(3)mu2taper
157	gezelter	2204
158	gezelter	1608	rf(1,atom2) = rf(1,atom2) + ul1(1)mu1taper
159			rf(2,atom2) = rf(2,atom2) + ul1(2)mu1taper
160			rf(3,atom2) = rf(3,atom2) + ul1(3)mu1taper
161			#endif
162	gezelter	2204
163
164	gezelter	1608	return
165			end subroutine accumulate_rf
166
167	gezelter	1930	subroutine accumulate_self_rf(atom1, mu1, eFrame)
168	gezelter	2204
169	gezelter	1608	integer, intent(in) :: atom1
170			real(kind=dp), intent(in) :: mu1
171	gezelter	1930	real(kind=dp), dimension(9,nLocal) :: eFrame
172	gezelter	2204
173	gezelter	1608	!! should work for both MPI and non-MPI version since this is not pairwise.
174	gezelter	1930	rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
175			rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
176			rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1
177	gezelter	2204
178	gezelter	1608	return
179			end subroutine accumulate_self_rf
180	gezelter	2204
181	gezelter	1930	subroutine reaction_field_final(a1, mu1, eFrame, rfpot, t, do_pot)
182	gezelter	2204
183	gezelter	1608	integer, intent(in) :: a1
184			real (kind=dp), intent(in) :: mu1
185			real (kind=dp), intent(inout) :: rfpot
186			logical, intent(in) :: do_pot
187	gezelter	1930	real (kind = dp), dimension(9,nLocal) :: eFrame
188	gezelter	1608	real (kind = dp), dimension(3,nLocal) :: t
189
190			integer :: localError
191
192			if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
193			write(default_error,*) 'Reaction field not initialized!'
194			return
195			endif
196
197			! compute torques on dipoles:
198			! pre converts from mu in units of debye to kcal/mol
199	gezelter	2204
200	gezelter	1608	! The torque contribution is dipole cross reaction_field
201
202	gezelter	1930	t(1,a1) = t(1,a1) + premu1(eFrame(6,a1)rf(3,a1) - eFrame(9,a1)rf(2,a1))
203			t(2,a1) = t(2,a1) + premu1(eFrame(9,a1)rf(1,a1) - eFrame(3,a1)rf(3,a1))
204			t(3,a1) = t(3,a1) + premu1(eFrame(3,a1)rf(2,a1) - eFrame(6,a1)rf(1,a1))
205	gezelter	2204
206	gezelter	1608	! the potential contribution is -1/2 dipole dot reaction_field
207	gezelter	2204
208	gezelter	1608	if (do_pot) then
209			rfpot = rfpot - 0.5d0 * pre * mu1 * &
210	gezelter	1930	(rf(1,a1)eFrame(3,a1) + rf(2,a1)eFrame(6,a1) + rf(3,a1)*eFrame(9,a1))
211	gezelter	1608	endif
212	gezelter	2204
213	gezelter	1608	return
214			end subroutine reaction_field_final
215	gezelter	2204
216	gezelter	1930	subroutine rf_correct_forces(atom1, atom2, d, rij, eFrame, taper, f, fpair)
217	gezelter	2204
218	gezelter	1608	integer, intent(in) :: atom1, atom2
219			real(kind=dp), dimension(3), intent(in) :: d
220			real(kind=dp), intent(in) :: rij, taper
221	gezelter	1930	real( kind = dp ), dimension(9,nLocal) :: eFrame
222	gezelter	1608	real( kind = dp ), dimension(3,nLocal) :: f
223			real( kind = dp ), dimension(3), intent(inout) :: fpair
224	gezelter	2204
225	gezelter	1608	real (kind = dp), dimension(3) :: ul1
226			real (kind = dp), dimension(3) :: ul2
227			real (kind = dp) :: dtdr
228			real (kind = dp) :: dudx, dudy, dudz, u1dotu2
229			integer :: me1, me2, id1, id2
230			real (kind = dp) :: mu1, mu2
231	gezelter	2204
232	gezelter	1608	integer :: localError
233
234			if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
235			write(default_error,*) 'Reaction field not initialized!'
236			return
237			endif
238
239			if (rij.le.rrf) then
240	gezelter	2204
241	gezelter	1608	if (rij.lt.rt) then
242			dtdr = 0.0d0
243			else
244	gezelter	2204	! write(,) 'rf correct in taper region'
245	gezelter	1608	dtdr = 6.0d0(rijrij - rijrt - rijrrf +rrfrt)/((rrf-rt)*3)
246			endif
247	gezelter	2204
248	gezelter	1608	#ifdef IS_MPI
249			me1 = atid_Row(atom1)
250	gezelter	1930	ul1(1) = eFrame_Row(3,atom1)
251			ul1(2) = eFrame_Row(6,atom1)
252			ul1(3) = eFrame_Row(9,atom1)
253	gezelter	2204
254	gezelter	1608	me2 = atid_Col(atom2)
255	gezelter	1930	ul2(1) = eFrame_Col(3,atom2)
256			ul2(2) = eFrame_Col(6,atom2)
257			ul2(3) = eFrame_Col(9,atom2)
258	gezelter	1608	#else
259			me1 = atid(atom1)
260	gezelter	1930	ul1(1) = eFrame(3,atom1)
261			ul1(2) = eFrame(6,atom1)
262			ul1(3) = eFrame(9,atom1)
263	gezelter	2204
264	gezelter	1608	me2 = atid(atom2)
265	gezelter	1930	ul2(1) = eFrame(3,atom2)
266			ul2(2) = eFrame(6,atom2)
267			ul2(3) = eFrame(9,atom2)
268	gezelter	1608	#endif
269	gezelter	2204
270	chrisfen	2129	mu1 = getDipoleMoment(me1)
271			mu2 = getDipoleMoment(me2)
272	gezelter	2204
273	gezelter	1608	u1dotu2 = ul1(1)ul2(1) + ul1(2)ul2(2) + ul1(3)*ul2(3)
274	gezelter	2204
275	gezelter	1608	dudx = - premu1mu2u1dotu2dtdr*d(1)/rij
276			dudy = - premu1mu2u1dotu2dtdr*d(2)/rij
277			dudz = - premu1mu2u1dotu2dtdr*d(3)/rij
278	gezelter	2204
279	gezelter	1608	#ifdef IS_MPI
280			f_Row(1,atom1) = f_Row(1,atom1) + dudx
281			f_Row(2,atom1) = f_Row(2,atom1) + dudy
282			f_Row(3,atom1) = f_Row(3,atom1) + dudz
283	gezelter	2204
284	gezelter	1608	f_Col(1,atom2) = f_Col(1,atom2) - dudx
285			f_Col(2,atom2) = f_Col(2,atom2) - dudy
286			f_Col(3,atom2) = f_Col(3,atom2) - dudz
287			#else
288			f(1,atom1) = f(1,atom1) + dudx
289			f(2,atom1) = f(2,atom1) + dudy
290			f(3,atom1) = f(3,atom1) + dudz
291	gezelter	2204
292	gezelter	1608	f(1,atom2) = f(1,atom2) - dudx
293			f(2,atom2) = f(2,atom2) - dudy
294			f(3,atom2) = f(3,atom2) - dudz
295			#endif
296
297			#ifdef IS_MPI
298			id1 = AtomRowToGlobal(atom1)
299			id2 = AtomColToGlobal(atom2)
300			#else
301			id1 = atom1
302			id2 = atom2
303			#endif
304	gezelter	2204
305	gezelter	1608	if (molMembershipList(id1) .ne. molMembershipList(id2)) then
306	gezelter	2204
307	gezelter	1608	fpair(1) = fpair(1) + dudx
308			fpair(2) = fpair(2) + dudy
309			fpair(3) = fpair(3) + dudz
310	gezelter	2204
311	gezelter	1608	endif
312	gezelter	2204
313	gezelter	1608	end if
314			return
315			end subroutine rf_correct_forces
316			end module reaction_field