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

Comparing trunk/OOPSE/libmdtools/calc_charge_charge.F90 (file contents):
Revision 1133 by gezelter, Mon Apr 26 14:29:18 2004 UTC vs.
Revision 1138 by gezelter, Wed Apr 28 21:39:22 2004 UTC

# Line 16 | Line 16 | module charge_charge
16    real(kind=dp), save :: rt  = 0.0_DP
17    real(kind=dp), save :: dielect  = 0.0_DP
18    real(kind=dp), save :: pre = 0.0_DP
19 –  real(kind=dp), save :: pre_rf = 0.0_DP
20 –  real(kind=dp), save :: ecr3 = 0.0_DP
21 –  real(kind=dp), save :: b0 = 0.0_DP
22 –  real(kind=dp), save :: pre_correction = 0.0_DP
19    
24 –  
20    logical, save :: haveCutoffs = .false.
21    logical, save :: haveChargeMap = .false.
22    logical, save :: haveDielectric = .false.
# Line 64 | Line 59 | contains
59      !   write(default_error,*) 'Dielect constant in charge module is not set!'
60      !endif
61      
67 –    ecr3 = ecr * ecr * ecr
68 –    b0 = 2.0d0*(dielect-1.0d0)/(2.0d0*dielect+1.0d0)
69 –    pre_rf = pre * b0 /(2.0d0* ecr3)
70 –    pre_correction = pre* (1 + 0.5d0 * b0) /ecr
62      haveCutoffs = .true.
63      
64      return
# Line 109 | Line 100 | contains
100    end subroutine createChargeMap
101    
102    
103 <  subroutine do_charge_pair(atom1, atom2, d, rij, r2, pot, f, &
104 <       do_pot, do_stress)
103 >  subroutine do_charge_pair(atom1, atom2, d, rij, r2, dc, rcij, rc2, mfact, &
104 >       pot, f, do_pot, do_stress, molecular_cutoffs)
105      
106 <    logical :: do_pot, do_stress
106 >    logical :: do_pot, do_stress, molecular_cutoffs
107  
108      integer atom1, atom2, me1, me2, id1, id2
109      integer :: localError
110 <    real(kind=dp) :: rij, r2, q1, q2
110 >    real(kind=dp) :: rij, r2, q1, q2, rcij, rc2
111      real(kind=dp) :: drdx, drdy, drdz, dudr, fx, fy, fz
112      real(kind=dp) :: taper, dtdr, vterm
113  
114      real( kind = dp ) :: pot
115 <    real( kind = dp ), dimension(3) :: d
115 >    real( kind = dp ), dimension(3) :: d, dc
116      real( kind = dp ), dimension(3,nLocal) :: f
117 <    real( kind = dp) :: ecr3, b0
118 <    real( kind = dp) :: vcoulomb, vrf, vcorrection
117 >    real( kind = dp ), dimension(nLocal) :: mfact
118 >    real( kind = dp ) :: theR, term1, term2, fxij, fyij, fzij
119 >    real( kind = dp ) :: fxab, fyab, fzab, fxba, fyba, fzba
120      
121      if (.not. haveCutoffs) then
122         write(default_error,*) 'charge-charge does not have cutoffs set!'
# Line 148 | Line 140 | contains
140      me2 = atid(atom2)
141   #endif
142  
151 –
143      q1 = ChargeMap(me1)%charge
144      q2 = ChargeMap(me2)%charge
145 <    
146 <    if (rij.le.ecr) then
145 >
146 >    if (molecular_cutoffs) then
147 >       theR = rcij
148 >    else
149 >       theR = rij
150 >    endif
151 >
152 >    if (theR.le.ecr) then
153        
154 <       if (rij.lt.rt) then
154 >       if (theR.lt.rt) then
155            taper = 1.0d0
156            dtdr = 0.0d0
157         else
158 <          taper = (ecr + 2.0d0*rij - 3.0d0*rt)*(ecr-rij)**2/ ((ecr-rt)**3)
159 <          dtdr = 6.0d0*(rij*rij - rij*rt - rij*ecr +ecr*rt)/((ecr-rt)**3)
158 >          taper = (ecr + 2.0d0*theR - 3.0d0*rt)*(ecr-theR)**2/ ((ecr-rt)**3)
159 >          dtdr = 6.0d0*(theR*theR - theR*rt - theR*ecr +ecr*rt)/((ecr-rt)**3)
160         endif
161  
162 <       vcoulomb = pre * q1 * q2 /rij
166 <       vrf = pre_rf * q1 * q2 * r2
167 <       ! vcorrection = pre_correction * q1 * q2
168 <      
169 <       vterm = vcoulomb + vrf
162 >       vterm = pre * q1 * q2 / rij
163  
164 <       ! dudr = vterm * ( dtdr  - taper / rij )
172 <      
173 <       dudr =  vterm * dtdr  + taper * ( 2.0 * vrf - vcoulomb ) / rij
164 >       if (molecular_cutoffs) then
165  
166 <       drdx = d(1) / rij
167 <       drdy = d(2) / rij
168 <       drdz = d(3) / rij
166 >          term1 = vterm * (taper / rij)
167 >        
168 >          fxij = term1 * d(1) / rij
169 >          fyij = term1 * d(2) / rij
170 >          fzij = term1 * d(3) / rij
171  
172 <       fx = dudr * drdx
180 <       fy = dudr * drdy
181 <       fz = dudr * drdz          
172 >          term2 = vterm * dtdr * mfact(atom1)
173  
174 +          fxab = term2 * dc(1) / rcij
175 +          fyab = term2 * dc(2) / rcij
176 +          fzab = term2 * dc(3) / rcij
177 +
178 +          term2 = vterm * dtdr * mfact(atom2)
179 +
180 +          fxba = term2 * dc(1) / rcij
181 +          fyba = term2 * dc(2) / rcij
182 +          fzba = term2 * dc(3) / rcij
183 +
184 +       else
185 +
186 +          dudr =  vterm * (dtdr  - taper / rij)
187 +          drdx = d(1) / rij
188 +          drdy = d(2) / rij
189 +          drdz = d(3) / rij
190 +          
191 +          fx = dudr * drdx
192 +          fy = dudr * drdy
193 +          fz = dudr * drdz          
194 +
195 +       endif
196 +
197 +
198   #ifdef IS_MPI
199         if (do_pot) then
200            pot_Row(atom1) = pot_Row(atom1) + vterm*taper*0.5
# Line 194 | Line 209 | contains
209         f_Col(2,atom2) = f_Col(2,atom2) - fy
210         f_Col(3,atom2) = f_Col(3,atom2) - fz
211  
212 +       if (molecular_cutoffs) then
213 +          f_Row(1,atom1) = f_Row(1,atom1) + fxab
214 +          f_Row(2,atom1) = f_Row(2,atom1) + fyab
215 +          f_Row(3,atom1) = f_Row(3,atom1) + fzab
216 +          
217 +          f_Col(1,atom2) = f_Col(1,atom2) - fxba
218 +          f_Col(2,atom2) = f_Col(2,atom2) - fyba
219 +          f_Col(3,atom2) = f_Col(3,atom2) - fzba
220 +       endif
221 +
222   #else
223  
224         if (do_pot) pot = pot + vterm*taper
# Line 206 | Line 231 | contains
231         f(2,atom2) = f(2,atom2) - fy
232         f(3,atom2) = f(3,atom2) - fz
233  
234 +       if (molecular_cutoffs) then
235 +          f(1,atom1) = f(1,atom1) + fxab
236 +          f(2,atom1) = f(2,atom1) + fyab
237 +          f(3,atom1) = f(3,atom1) + fzab
238 +          
239 +          f(1,atom2) = f(1,atom2) - fxba
240 +          f(2,atom2) = f(2,atom2) - fyba
241 +          f(3,atom2) = f(3,atom2) - fzba
242 +       endif
243 +
244   #endif
245  
246         if (do_stress) then
# Line 233 | Line 268 | contains
268               tau_Temp(6) = tau_Temp(6) - d(2) * fz
269               tau_Temp(8) = tau_Temp(8) - d(3) * fy
270               tau_Temp(9) = tau_Temp(9) - d(3) * fz
271 +
272 +             if (molecular_cutoffs) then
273 +                
274 +                tau_Temp(1) = tau_Temp(1) - d(1) * fxab
275 +                tau_Temp(2) = tau_Temp(2) - d(1) * fyab
276 +                tau_Temp(3) = tau_Temp(3) - d(1) * fzab
277 +                tau_Temp(4) = tau_Temp(4) - d(2) * fxab
278 +                tau_Temp(5) = tau_Temp(5) - d(2) * fyab
279 +                tau_Temp(6) = tau_Temp(6) - d(2) * fzab
280 +                tau_Temp(8) = tau_Temp(8) - d(3) * fyab
281 +                tau_Temp(9) = tau_Temp(9) - d(3) * fzab
282 +             endif
283              
284               virial_Temp = virial_Temp + &
285                    (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))

Diff Legend

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