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 394 by gezelter, Mon Mar 24 21:55:34 2003 UTC vs.
Revision 898 by chuckv, Mon Jan 5 22:49:14 2004 UTC

# Line 3 | Line 3 | module reaction_field
3    use definitions
4    use atype_module
5    use vector_class
6 +  use simulation
7   #ifdef IS_MPI
8    use mpiSimulation
9   #endif
10    implicit none
11  
12 <  real(kind=dp), save :: rrf
12 >  PRIVATE
13 >  
14 >  real(kind=dp), save :: rrf = 1.0_dp
15    real(kind=dp), save :: rt
16 <  real(kind=dp), save :: dielect
17 <  real(kind=dp), save :: rrfsq
16 >  real(kind=dp), save :: dielect = 1.0_dp
17 >  real(kind=dp), save :: rrfsq = 1.0_dp
18    real(kind=dp), save :: pre
19 <  logical, save :: rf_initialized = .false.
19 >  logical, save :: rf_initialized = .false., haveCuts = .false.
20 >  logical, save :: haveDie = .false.
21  
22 +  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   contains
30    
31 <  subroutine initialize_rf(this_rrf, this_rt, this_dielect)
32 <    real(kind=dp), intent(in) :: this_rrf, this_rt, this_dielect
33 <
23 <    rrf = this_rrf
24 <    rt = this_rt
31 >  subroutine initialize_rf(this_dielect)
32 >    real(kind=dp), intent(in) :: this_dielect
33 >    
34      dielect = this_dielect
35 +
36 +    pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+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 +    rrf = this_rrf
49 +    rt  = this_rt
50 +
51      rrfsq = rrf * rrf
52      pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf)
53      
54 +    haveCuts = .true.
55 +    if (haveDie) rf_initialized = .true.
56  
57 <    write(*,*) 'rrf = ', rrf
58 <    write(*,*) 'rt = ', rt
33 <    write(*,*) 'dielect = ', dielect
34 <    write(*,*) 'pre = ', pre
35 <    rf_initialized = .true.
36 <    
37 <    return
38 <  end subroutine initialize_rf
57 >  end subroutine setCutoffsRF
58 >
59    
60    subroutine accumulate_rf(atom1, atom2, rij, u_l)
61  
62      integer, intent(in) :: atom1, atom2
63      real (kind = dp), intent(in) :: rij
64 <    real (kind = dp), dimension(:,:) :: u_l    
64 >    real (kind = dp), dimension(3,nLocal) :: u_l    
65  
66      integer :: me1, me2
67      real (kind = dp) :: taper, mu1, mu2
# Line 58 | Line 78 | contains
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
84        
# Line 113 | Line 134 | contains
134      
135      integer, intent(in) :: atom1
136      real(kind=dp), intent(in) :: mu1
137 <    real(kind=dp), dimension(:,:) :: u_l
137 >    real(kind=dp), dimension(3,nLocal) :: u_l
138      
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
# Line 129 | Line 150 | contains
150      real (kind=dp), intent(in) :: mu1
151      real (kind=dp), intent(inout) :: rfpot
152      logical, intent(in) :: do_pot
153 <    real (kind = dp), dimension(:,:) :: u_l    
154 <    real (kind = dp), dimension(:,:) :: t
153 >    real (kind = dp), dimension(3,nLocal) :: u_l    
154 >    real (kind = dp), dimension(3,nLocal) :: t
155  
156      if (.not.rf_initialized) then
157         write(default_error,*) 'Reaction field not initialized!'
# Line 161 | Line 182 | contains
182      integer, intent(in) :: atom1, atom2
183      real(kind=dp), dimension(3), intent(in) :: d
184      real(kind=dp), intent(in) :: rij
185 <    real( kind = dp ), dimension(:,:) :: u_l
186 <    real( kind = dp ), dimension(:,:) :: f
185 >    real( kind = dp ), dimension(3,nLocal) :: u_l
186 >    real( kind = dp ), dimension(3,nLocal) :: f
187      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 <    integer :: me1, me2
193 >    integer :: me1, me2, id1, id2
194      real (kind = dp) :: mu1, mu2
195      
196      if (.not.rf_initialized) then
# Line 182 | Line 203 | contains
203         if (rij.lt.rt) then
204            dtdr = 0.0d0
205         else
206 +          write(*,*) 'rf correct in taper region'
207            dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
208         endif
209        
# Line 235 | Line 257 | contains
257   #endif
258        
259         if (do_stress) then          
260 <          tau_Temp(1) = tau_Temp(1) + dudx * d(1)
261 <          tau_Temp(2) = tau_Temp(2) + dudx * d(2)
262 <          tau_Temp(3) = tau_Temp(3) + dudx * d(3)
263 <          tau_Temp(4) = tau_Temp(4) + dudy * d(1)
264 <          tau_Temp(5) = tau_Temp(5) + dudy * d(2)
265 <          tau_Temp(6) = tau_Temp(6) + dudy * d(3)
266 <          tau_Temp(7) = tau_Temp(7) + dudz * d(1)
267 <          tau_Temp(8) = tau_Temp(8) + dudz * d(2)
268 <          tau_Temp(9) = tau_Temp(9) + dudz * d(3)
269 <          virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
260 >
261 > #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 >             ! 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 >             virial_Temp = virial_Temp + &
286 >                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
287 >          endif
288         endif      
289      endif
290      

Diff Legend

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