ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_reaction_field.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/calc_reaction_field.F90 (file contents):
Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC vs.
Revision 1192 by gezelter, Mon May 24 21:03:30 2004 UTC

# Line 4 | Line 4 | module reaction_field
4    use atype_module
5    use vector_class
6    use simulation
7 +  use status
8   #ifdef IS_MPI
9    use mpiSimulation
10   #endif
# Line 16 | Line 17 | module reaction_field
17    real(kind=dp), save :: dielect = 1.0_dp
18    real(kind=dp), save :: rrfsq = 1.0_dp
19    real(kind=dp), save :: pre
20 <  logical, save :: rf_initialized = .false., haveCuts = .false.
21 <  logical, save :: haveDie = .false.
20 >  logical, save :: haveCutoffs = .false.
21 >  logical, save :: haveMomentMap = .false.
22 >  logical, save :: haveDielectric = .false.
23  
24 +  type :: MomentList
25 +     real(kind=DP) :: dipole_moment = 0.0_DP
26 +  end type MomentList
27 +
28 +  type(MomentList), dimension(:),allocatable :: MomentMap
29 +
30    PUBLIC::initialize_rf
31    PUBLIC::setCutoffsRF
32    PUBLIC::accumulate_rf
# Line 35 | Line 43 | contains
43  
44      pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
45      
46 <    haveDie = .true.
39 <    if (haveCuts) rf_initialized = .true.
46 >    haveDielectric = .true.
47  
48      return
49    end subroutine initialize_rf
# Line 51 | Line 58 | contains
58      rrfsq = rrf * rrf
59      pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
60      
61 <    haveCuts = .true.
55 <    if (haveDie) rf_initialized = .true.
61 >    haveCutoffs = .true.
62  
63    end subroutine setCutoffsRF
64  
65 <  
66 <  subroutine accumulate_rf(atom1, atom2, rij, u_l)
65 >  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 +  subroutine accumulate_rf(atom1, atom2, rij, u_l, taper)
101 +
102      integer, intent(in) :: atom1, atom2
103      real (kind = dp), intent(in) :: rij
104      real (kind = dp), dimension(3,nLocal) :: u_l    
105  
106      integer :: me1, me2
107 <    real (kind = dp) :: taper, mu1, mu2
107 >    real (kind = dp), intent(in) :: taper
108 >    real (kind = dp):: mu1, mu2
109      real (kind = dp), dimension(3) :: ul1
110      real (kind = dp), dimension(3) :: ul2  
111  
112 <    if (.not.rf_initialized) then
112 >    integer :: localError
113 >
114 >    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
115         write(default_error,*) 'Reaction field not initialized!'
116         return
117      endif
118 +
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      
76 –    if (rij.le.rrf) then
77 –          
78 –       if (rij.lt.rt) then
79 –          taper = 1.0d0
80 –       else
81 –          write(*,*) 'rf in taper region'
82 –          taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3)
83 –       endif
128        
129   #ifdef IS_MPI
130 <       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)
130 >    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   #else
140 <       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)
140 >    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   #endif
150 <      
151 <       call getElementProperty(atypes, me1, "dipole_moment", mu1)
152 <       call getElementProperty(atypes, me2, "dipole_moment", mu2)
153 <      
110 <      
150 >    
151 >    mu1 = MomentMap(me1)%dipole_moment
152 >    mu2 = MomentMap(me2)%dipole_moment
153 >    
154   #ifdef IS_MPI
155 <       rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
156 <       rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
157 <       rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
158 <      
159 <       rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
160 <       rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
161 <       rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
155 >    rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
156 >    rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
157 >    rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
158 >    
159 >    rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
160 >    rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
161 >    rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
162   #else
163 <       rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
164 <       rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
165 <       rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
166 <      
167 <       rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
168 <       rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
169 <       rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper    
163 >    rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
164 >    rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
165 >    rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
166 >    
167 >    rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
168 >    rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
169 >    rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper    
170   #endif
171 <      
172 <    endif
171 >    
172 >    
173      return  
174    end subroutine accumulate_rf
175  
# Line 153 | Line 196 | contains
196      real (kind = dp), dimension(3,nLocal) :: u_l    
197      real (kind = dp), dimension(3,nLocal) :: t
198  
199 <    if (.not.rf_initialized) then
199 >    integer :: localError
200 >
201 >    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
202         write(default_error,*) 'Reaction field not initialized!'
203         return
204      endif
205  
206 +    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      ! compute torques on dipoles:
216      ! pre converts from mu in units of debye to kcal/mol
217      
# Line 173 | Line 227 | contains
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 <    
230 >
231      return
232    end subroutine reaction_field_final
233    
234 <  subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, f, do_stress)
234 >  subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, taper, f, fpair)
235      
236      integer, intent(in) :: atom1, atom2
237      real(kind=dp), dimension(3), intent(in) :: d
238 <    real(kind=dp), intent(in) :: rij
238 >    real(kind=dp), intent(in) :: rij, taper
239      real( kind = dp ), dimension(3,nLocal) :: u_l
240      real( kind = dp ), dimension(3,nLocal) :: f
241 <    logical, intent(in) :: do_stress
241 >    real( kind = dp ), dimension(3), intent(inout) :: fpair
242      
243      real (kind = dp), dimension(3) :: ul1
244      real (kind = dp), dimension(3) :: ul2
# Line 193 | Line 247 | contains
247      integer :: me1, me2, id1, id2
248      real (kind = dp) :: mu1, mu2
249      
250 <    if (.not.rf_initialized) then
250 >    integer :: localError
251 >
252 >    if ((.not.haveDielectric).or.(.not.haveCutoffs)) then
253         write(default_error,*) 'Reaction field not initialized!'
254         return
255      endif
256  
257 +    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      if (rij.le.rrf) then
267        
268         if (rij.lt.rt) then
269            dtdr = 0.0d0
270         else
271 <          write(*,*) 'rf correct in taper region'
271 > !         write(*,*) 'rf correct in taper region'
272            dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
273         endif
274        
# Line 229 | Line 294 | contains
294         ul2(3) = u_l(3,atom2)
295   #endif
296        
297 <       call getElementProperty(atypes, me1, "dipole_moment", mu1)
298 <       call getElementProperty(atypes, me2, "dipole_moment", mu2)
297 >       mu1 = MomentMap(me1)%dipole_moment
298 >       mu2 = MomentMap(me2)%dipole_moment
299        
300         u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
301        
# Line 255 | Line 320 | contains
320         f(2,atom2) = f(2,atom2) - dudy
321         f(3,atom2) = f(3,atom2) - dudz
322   #endif
258 –      
259 –       if (do_stress) then          
323  
324   #ifdef IS_MPI
325 <          id1 = tagRow(atom1)
326 <          id2 = tagColumn(atom2)
325 >       id1 = tagRow(atom1)
326 >       id2 = tagColumn(atom2)
327   #else
328 <          id1 = atom1
329 <          id2 = atom2
328 >       id1 = atom1
329 >       id2 = atom2
330   #endif
331 <
332 <          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
333 <
334 <             ! because the d vector is the rj - ri vector, and
335 <             ! because dudx, dudy, and dudz are the
336 <             ! (positive) force on atom i (negative on atom j) we need
337 <             ! a negative sign here:
338 <
339 <             tau_Temp(1) = tau_Temp(1) - d(1) * dudx
340 <             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 <             virial_Temp = virial_Temp + &
286 <                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
287 <          endif
288 <       endif      
289 <    endif
290 <    
331 >      
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      return
342    end subroutine rf_correct_forces
343   end module reaction_field

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines