OOPSE/libmdtools/calc_reaction_field.F90

module reaction_field
  use force_globals
  use definitions
  use atype_module
  use vector_class
  use simulation
#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 :: rf_initialized = .false., haveCuts = .false.
  logical, save :: haveDie = .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 = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) 
    
    haveDie = .true.
    if (haveCuts) rf_initialized = .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)
    
    haveCuts = .true.
    if (haveDie) rf_initialized = .true.

  end subroutine setCutoffsRF

  
  subroutine accumulate_rf(atom1, atom2, rij, u_l)

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

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

    if (.not.rf_initialized) then
       write(default_error,*) 'Reaction field not initialized!'
       return
    endif
    
    if (rij.le.rrf) then
          
       if (rij.lt.rt) then
          taper = 1.0d0
       else
          write(*,*) 'rf in taper region'
          taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((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
       
       call getElementProperty(atypes, me1, "dipole_moment", mu1)
       call getElementProperty(atypes, me2, "dipole_moment", mu2)
       
       
#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
       
    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,getNlocal()) :: 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,getNlocal()) :: u_l    
    real (kind = dp), dimension(3,getNlocal()) :: t

    if (.not.rf_initialized) 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*(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, f, do_stress)
    
    integer, intent(in) :: atom1, atom2
    real(kind=dp), dimension(3), intent(in) :: d
    real(kind=dp), intent(in) :: rij
    real( kind = dp ), dimension(3,getNlocal()) :: u_l
    real( kind = dp ), dimension(3,getNlocal()) :: f
    logical, intent(in) :: do_stress
    
    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
    
    if (.not.rf_initialized) 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) = 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
       
       call getElementProperty(atypes, me1, "dipole_moment", mu1)
       call getElementProperty(atypes, me2, "dipole_moment", mu2)
       
       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
       
       if (do_stress) 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, and dudz are the
             ! (positive) force on atom i (negative on atom j) we need
             ! a negative sign here:

             tau_Temp(1) = tau_Temp(1) - d(1) * dudx
             tau_Temp(2) = tau_Temp(2) - d(1) * dudy
             tau_Temp(3) = tau_Temp(3) - d(1) * dudz
             tau_Temp(4) = tau_Temp(4) - d(2) * dudx
             tau_Temp(5) = tau_Temp(5) - d(2) * dudy
             tau_Temp(6) = tau_Temp(6) - d(2) * dudz
             tau_Temp(7) = tau_Temp(7) - d(3) * dudx
             tau_Temp(8) = tau_Temp(8) - d(3) * dudy
             tau_Temp(9) = tau_Temp(9) - d(3) * dudz
             virial_Temp = virial_Temp + &
                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
          endif
       endif       
    endif
    
    return
  end subroutine rf_correct_forces
end module reaction_field
Revision:	730
Committed:	Wed Aug 27 16:25:11 2003 UTC (20 years, 10 months ago) by gezelter
File size:	8627 byte(s)
Log Message:	More fixes for stress tensor parallel bug.
#	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	mmeineke	377	#ifdef IS_MPI
8			use mpiSimulation
9			#endif
10			implicit none
11
12	mmeineke	626	PRIVATE
13
14			real(kind=dp), save :: rrf = 1.0_dp
15	mmeineke	377	real(kind=dp), save :: rt
16	mmeineke	626	real(kind=dp), save :: dielect = 1.0_dp
17			real(kind=dp), save :: rrfsq = 1.0_dp
18	mmeineke	377	real(kind=dp), save :: pre
19	mmeineke	626	logical, save :: rf_initialized = .false., haveCuts = .false.
20			logical, save :: haveDie = .false.
21	mmeineke	377
22	mmeineke	626	PUBLIC::initialize_rf
23			PUBLIC::setCutoffsRF
24			PUBLIC::accumulate_rf
25			PUBLIC::accumulate_self_rf
26			PUBLIC::reaction_field_final
27			PUBLIC::rf_correct_forces
28
29	mmeineke	377	contains
30
31	mmeineke	626	subroutine initialize_rf(this_dielect)
32			real(kind=dp), intent(in) :: this_dielect
33
34			dielect = this_dielect
35	mmeineke	377
36	mmeineke	626	pre = 14.38362d02.0d0(dielect-1.0d0)/((2.0d0dielect+1.0d0)rrfsq*rrf)
37
38			haveDie = .true.
39			if (haveCuts) rf_initialized = .true.
40
41			return
42			end subroutine initialize_rf
43
44			subroutine setCutoffsRF( this_rrf, this_rt )
45
46			real(kind=dp), intent(in) :: this_rrf, this_rt
47
48	mmeineke	377	rrf = this_rrf
49	mmeineke	626	rt = this_rt
50
51	mmeineke	377	rrfsq = rrf * rrf
52			pre = 14.38362d02.0d0(dielect-1.0d0)/((2.0d0dielect+1.0d0)rrfsq*rrf)
53
54	mmeineke	626	haveCuts = .true.
55			if (haveDie) rf_initialized = .true.
56	gezelter	394
57	mmeineke	626	end subroutine setCutoffsRF
58
59	mmeineke	377
60			subroutine accumulate_rf(atom1, atom2, rij, u_l)
61
62			integer, intent(in) :: atom1, atom2
63			real (kind = dp), intent(in) :: rij
64	chuckv	460	real (kind = dp), dimension(3,getNlocal()) :: u_l
65	mmeineke	377
66			integer :: me1, me2
67			real (kind = dp) :: taper, mu1, mu2
68			real (kind = dp), dimension(3) :: ul1
69			real (kind = dp), dimension(3) :: ul2
70
71			if (.not.rf_initialized) then
72			write(default_error,*) 'Reaction field not initialized!'
73			return
74			endif
75
76			if (rij.le.rrf) then
77
78			if (rij.lt.rt) then
79			taper = 1.0d0
80			else
81	gezelter	621	write(,) 'rf in taper region'
82	mmeineke	377	taper = (rrf + 2.0d0rij - 3.0d0rt)(rrf-rij)2/ ((rrf-rt)*3)
83			endif
84
85			#ifdef IS_MPI
86			me1 = atid_Row(atom1)
87			ul1(1) = u_l_Row(1,atom1)
88			ul1(2) = u_l_Row(2,atom1)
89			ul1(3) = u_l_Row(3,atom1)
90
91			me2 = atid_Col(atom2)
92			ul2(1) = u_l_Col(1,atom2)
93			ul2(2) = u_l_Col(2,atom2)
94			ul2(3) = u_l_Col(3,atom2)
95			#else
96			me1 = atid(atom1)
97			ul1(1) = u_l(1,atom1)
98			ul1(2) = u_l(2,atom1)
99			ul1(3) = u_l(3,atom1)
100
101			me2 = atid(atom2)
102			ul2(1) = u_l(1,atom2)
103			ul2(2) = u_l(2,atom2)
104			ul2(3) = u_l(3,atom2)
105			#endif
106
107			call getElementProperty(atypes, me1, "dipole_moment", mu1)
108			call getElementProperty(atypes, me2, "dipole_moment", mu2)
109
110
111			#ifdef IS_MPI
112			rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)mu2taper
113			rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)mu2taper
114			rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)mu2taper
115
116	gezelter	394	rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)mu1taper
117			rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)mu1taper
118			rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)mu1taper
119	mmeineke	377	#else
120			rf(1,atom1) = rf(1,atom1) + ul2(1)mu2taper
121			rf(2,atom1) = rf(2,atom1) + ul2(2)mu2taper
122			rf(3,atom1) = rf(3,atom1) + ul2(3)mu2taper
123
124	gezelter	394	rf(1,atom2) = rf(1,atom2) + ul1(1)mu1taper
125			rf(2,atom2) = rf(2,atom2) + ul1(2)mu1taper
126			rf(3,atom2) = rf(3,atom2) + ul1(3)mu1taper
127	mmeineke	377	#endif
128
129			endif
130			return
131			end subroutine accumulate_rf
132
133			subroutine accumulate_self_rf(atom1, mu1, u_l)
134
135			integer, intent(in) :: atom1
136			real(kind=dp), intent(in) :: mu1
137	chuckv	460	real(kind=dp), dimension(3,getNlocal()) :: u_l
138	mmeineke	377
139			!! should work for both MPI and non-MPI version since this is not pairwise.
140			rf(1,atom1) = rf(1,atom1) + u_l(1,atom1)*mu1
141			rf(2,atom1) = rf(2,atom1) + u_l(2,atom1)*mu1
142			rf(3,atom1) = rf(3,atom1) + u_l(3,atom1)*mu1
143
144			return
145			end subroutine accumulate_self_rf
146
147			subroutine reaction_field_final(a1, mu1, u_l, rfpot, t, do_pot)
148
149			integer, intent(in) :: a1
150			real (kind=dp), intent(in) :: mu1
151			real (kind=dp), intent(inout) :: rfpot
152			logical, intent(in) :: do_pot
153	chuckv	460	real (kind = dp), dimension(3,getNlocal()) :: u_l
154			real (kind = dp), dimension(3,getNlocal()) :: t
155	mmeineke	377
156			if (.not.rf_initialized) then
157			write(default_error,*) 'Reaction field not initialized!'
158			return
159			endif
160
161			! compute torques on dipoles:
162			! pre converts from mu in units of debye to kcal/mol
163
164	gezelter	394	! The torque contribution is dipole cross reaction_field
165
166	mmeineke	377	t(1,a1) = t(1,a1) + premu1(u_l(2,a1)rf(3,a1) - u_l(3,a1)rf(2,a1))
167			t(2,a1) = t(2,a1) + premu1(u_l(3,a1)rf(1,a1) - u_l(1,a1)rf(3,a1))
168			t(3,a1) = t(3,a1) + premu1(u_l(1,a1)rf(2,a1) - u_l(2,a1)rf(1,a1))
169
170			! the potential contribution is -1/2 dipole dot reaction_field
171
172			if (do_pot) then
173			rfpot = rfpot - 0.5d0 * pre * mu1 * &
174			(rf(1,a1)u_l(1,a1) + rf(2,a1)u_l(2,a1) + rf(3,a1)*u_l(3,a1))
175			endif
176
177			return
178			end subroutine reaction_field_final
179
180			subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, f, do_stress)
181
182			integer, intent(in) :: atom1, atom2
183			real(kind=dp), dimension(3), intent(in) :: d
184			real(kind=dp), intent(in) :: rij
185	chuckv	460	real( kind = dp ), dimension(3,getNlocal()) :: u_l
186			real( kind = dp ), dimension(3,getNlocal()) :: f
187	mmeineke	377	logical, intent(in) :: do_stress
188
189			real (kind = dp), dimension(3) :: ul1
190			real (kind = dp), dimension(3) :: ul2
191			real (kind = dp) :: dtdr
192			real (kind = dp) :: dudx, dudy, dudz, u1dotu2
193	gezelter	730	integer :: me1, me2, id1, id2
194	mmeineke	377	real (kind = dp) :: mu1, mu2
195
196			if (.not.rf_initialized) then
197			write(default_error,*) 'Reaction field not initialized!'
198			return
199			endif
200
201			if (rij.le.rrf) then
202
203			if (rij.lt.rt) then
204			dtdr = 0.0d0
205			else
206	gezelter	621	write(,) 'rf correct in taper region'
207	mmeineke	377	dtdr = 6.0d0(rijrij - rijrt - rijrrf +rrfrt)/((rrf-rt)*3)
208			endif
209
210			#ifdef IS_MPI
211			me1 = atid_Row(atom1)
212			ul1(1) = u_l_Row(1,atom1)
213			ul1(2) = u_l_Row(2,atom1)
214			ul1(3) = u_l_Row(3,atom1)
215
216			me2 = atid_Col(atom2)
217			ul2(1) = u_l_Col(1,atom2)
218			ul2(2) = u_l_Col(2,atom2)
219			ul2(3) = u_l_Col(3,atom2)
220			#else
221			me1 = atid(atom1)
222			ul1(1) = u_l(1,atom1)
223			ul1(2) = u_l(2,atom1)
224			ul1(3) = u_l(3,atom1)
225
226			me2 = atid(atom2)
227			ul2(1) = u_l(1,atom2)
228			ul2(2) = u_l(2,atom2)
229			ul2(3) = u_l(3,atom2)
230			#endif
231
232			call getElementProperty(atypes, me1, "dipole_moment", mu1)
233			call getElementProperty(atypes, me2, "dipole_moment", mu2)
234
235			u1dotu2 = ul1(1)ul2(1) + ul1(2)ul2(2) + ul1(3)*ul2(3)
236
237			dudx = - premu1mu2u1dotu2dtdr*d(1)/rij
238			dudy = - premu1mu2u1dotu2dtdr*d(2)/rij
239			dudz = - premu1mu2u1dotu2dtdr*d(3)/rij
240
241			#ifdef IS_MPI
242			f_Row(1,atom1) = f_Row(1,atom1) + dudx
243			f_Row(2,atom1) = f_Row(2,atom1) + dudy
244			f_Row(3,atom1) = f_Row(3,atom1) + dudz
245
246			f_Col(1,atom2) = f_Col(1,atom2) - dudx
247			f_Col(2,atom2) = f_Col(2,atom2) - dudy
248			f_Col(3,atom2) = f_Col(3,atom2) - dudz
249			#else
250			f(1,atom1) = f(1,atom1) + dudx
251			f(2,atom1) = f(2,atom1) + dudy
252			f(3,atom1) = f(3,atom1) + dudz
253
254			f(1,atom2) = f(1,atom2) - dudx
255			f(2,atom2) = f(2,atom2) - dudy
256			f(3,atom2) = f(3,atom2) - dudz
257			#endif
258
259			if (do_stress) then
260	gezelter	611
261	gezelter	730	#ifdef IS_MPI
262			id1 = tagRow(atom1)
263			id2 = tagColumn(atom2)
264			#else
265			id1 = atom1
266			id2 = atom2
267			#endif
268
269			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
270
271	gezelter	611	! because the d vector is the rj - ri vector, and
272			! because dudx, dudy, and dudz are the
273			! (positive) force on atom i (negative on atom j) we need
274			! a negative sign here:
275
276			tau_Temp(1) = tau_Temp(1) - d(1) * dudx
277			tau_Temp(2) = tau_Temp(2) - d(1) * dudy
278			tau_Temp(3) = tau_Temp(3) - d(1) * dudz
279			tau_Temp(4) = tau_Temp(4) - d(2) * dudx
280			tau_Temp(5) = tau_Temp(5) - d(2) * dudy
281			tau_Temp(6) = tau_Temp(6) - d(2) * dudz
282			tau_Temp(7) = tau_Temp(7) - d(3) * dudx
283			tau_Temp(8) = tau_Temp(8) - d(3) * dudy
284			tau_Temp(9) = tau_Temp(9) - d(3) * dudz
285	gezelter	483	virial_Temp = virial_Temp + &
286			(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
287			endif
288	mmeineke	377	endif
289			endif
290
291			return
292			end subroutine rf_correct_forces
293			end module reaction_field