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 1138 by gezelter, Wed Apr 28 21:39:22 2004 UTC vs.
Revision 1160 by gezelter, Tue May 11 21:31:15 2004 UTC

# Line 12 | Line 12 | module charge_charge
12    implicit none
13  
14    PRIVATE
15 <  real(kind=dp), save :: ecr = 0.0_DP
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 <  
20 <  logical, save :: haveCutoffs = .false.
15 >
16 >  real(kind=dp), parameter :: pre = 332.06508_DP  
17    logical, save :: haveChargeMap = .false.
22 –  logical, save :: haveDielectric = .false.
18  
24 –  public::setCutoffsCharge
19    public::do_charge_pair
26 –  public::initialize_charge
20  
21    type :: ChargeList
22       real(kind=DP) :: charge = 0.0_DP
# Line 32 | Line 25 | contains
25    type(ChargeList), dimension(:), allocatable :: ChargeMap
26  
27   contains
28 <
36 <  subroutine initialize_charge(this_dielect)
37 <    real(kind=dp), intent(in) :: this_dielect
38 <    
39 <    dielect = this_dielect
40 <    haveDielectric = .true.
41 <    
42 <    ! because setCutoffsCharge is called before initialize_charge
43 <    ! we need to call it agin to make sure all of the precalculation
44 <    ! value is correct. for the time being, just a quick hack
45 <    call setCutoffsCharge(ecr, rt)
46 <    return
47 <  end subroutine initialize_charge
48 <  
49 <    
50 <  subroutine setCutoffsCharge(this_ecr, this_rt)
51 <    real(kind=dp), intent(in) :: this_ecr, this_rt
52 <    ecr = this_ecr
53 <    rt = this_rt
54 <
55 <    ! pre converts from fundamental charge to kcal/mol
56 <    pre = 332.06508_DP
57 <    
58 <    !if (.not.haveDielectric)then
59 <    !   write(default_error,*) 'Dielect constant in charge module is not set!'
60 <    !endif
61 <    
62 <    haveCutoffs = .true.
63 <    
64 <    return
65 <  end subroutine setCutoffsCharge
66 <
28 >      
29    subroutine createChargeMap(status)
30      integer :: nAtypes
31      integer :: status
# Line 98 | Line 60 | contains
60      haveChargeMap = .true.
61      
62    end subroutine createChargeMap
101 –  
102 –  
103 –  subroutine do_charge_pair(atom1, atom2, d, rij, r2, dc, rcij, rc2, mfact, &
104 –       pot, f, do_pot, do_stress, molecular_cutoffs)
63      
64 <    logical :: do_pot, do_stress, molecular_cutoffs
64 >  subroutine do_charge_pair(atom1, atom2, d, rij, r2, sw, vpair, &
65 >       pot, f, do_pot, do_stress)
66 >    
67 >    logical :: do_pot, do_stress
68  
69      integer atom1, atom2, me1, me2, id1, id2
70      integer :: localError
71 <    real(kind=dp) :: rij, r2, q1, q2, rcij, rc2
71 >    real(kind=dp) :: rij, r2, q1, q2, sw, vpair
72      real(kind=dp) :: drdx, drdy, drdz, dudr, fx, fy, fz
73 <    real(kind=dp) :: taper, dtdr, vterm
73 >    real(kind=dp) :: vterm
74  
75      real( kind = dp ) :: pot
76 <    real( kind = dp ), dimension(3) :: d, dc
76 >    real( kind = dp ), dimension(3) :: d
77      real( kind = dp ), dimension(3,nLocal) :: f
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
78      
121 –    if (.not. haveCutoffs) then
122 –       write(default_error,*) 'charge-charge does not have cutoffs set!'
123 –       return
124 –    endif
125 –
79      if (.not.haveChargeMap) then
80         localError = 0
81         call createChargeMap(localError)
# Line 143 | Line 96 | contains
96      q1 = ChargeMap(me1)%charge
97      q2 = ChargeMap(me2)%charge
98  
99 <    if (molecular_cutoffs) then
147 <       theR = rcij
148 <    else
149 <       theR = rij
150 <    endif
99 >    vterm = pre * q1 * q2 / rij
100  
101 <    if (theR.le.ecr) then
153 <      
154 <       if (theR.lt.rt) then
155 <          taper = 1.0d0
156 <          dtdr = 0.0d0
157 <       else
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
101 >    dudr =  -sw * vterm  / rij
102  
103 <       vterm = pre * q1 * q2 / rij
103 >    drdx = d(1) / rij
104 >    drdy = d(2) / rij
105 >    drdz = d(3) / rij
106 >    
107 >    fx = dudr * drdx
108 >    fy = dudr * drdy
109 >    fz = dudr * drdz          
110  
111 <       if (molecular_cutoffs) then
112 <
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 <          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 <
111 >    vpair = vpair + vterm
112 >    
113   #ifdef IS_MPI
114 <       if (do_pot) then
115 <          pot_Row(atom1) = pot_Row(atom1) + vterm*taper*0.5
116 <          pot_Col(atom2) = pot_Col(atom2) + vterm*taper*0.5
117 <       endif
118 <
119 <       f_Row(1,atom1) = f_Row(1,atom1) + fx
120 <       f_Row(2,atom1) = f_Row(2,atom1) + fy
121 <       f_Row(3,atom1) = f_Row(3,atom1) + fz
122 <
123 <       f_Col(1,atom2) = f_Col(1,atom2) - fx
124 <       f_Col(2,atom2) = f_Col(2,atom2) - fy
125 <       f_Col(3,atom2) = f_Col(3,atom2) - fz
126 <
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 <
114 >    if (do_pot) then
115 >       pot_Row(atom1) = pot_Row(atom1) + sw*vterm*0.5
116 >       pot_Col(atom2) = pot_Col(atom2) + sw*vterm*0.5
117 >    endif
118 >    
119 >    f_Row(1,atom1) = f_Row(1,atom1) + fx
120 >    f_Row(2,atom1) = f_Row(2,atom1) + fy
121 >    f_Row(3,atom1) = f_Row(3,atom1) + fz
122 >    
123 >    f_Col(1,atom2) = f_Col(1,atom2) - fx
124 >    f_Col(2,atom2) = f_Col(2,atom2) - fy
125 >    f_Col(3,atom2) = f_Col(3,atom2) - fz
126 >    
127   #else
128  
129 <       if (do_pot) pot = pot + vterm*taper
130 <
131 <       f(1,atom1) = f(1,atom1) + fx
132 <       f(2,atom1) = f(2,atom1) + fy
133 <       f(3,atom1) = f(3,atom1) + fz
134 <
135 <       f(1,atom2) = f(1,atom2) - fx
136 <       f(2,atom2) = f(2,atom2) - fy
137 <       f(3,atom2) = f(3,atom2) - fz
138 <
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 <
129 >    if (do_pot) pot = pot + sw*vterm
130 >    
131 >    f(1,atom1) = f(1,atom1) + fx
132 >    f(2,atom1) = f(2,atom1) + fy
133 >    f(3,atom1) = f(3,atom1) + fz
134 >    
135 >    f(1,atom2) = f(1,atom2) - fx
136 >    f(2,atom2) = f(2,atom2) - fy
137 >    f(3,atom2) = f(3,atom2) - fz
138 >        
139   #endif
140 <
141 <       if (do_stress) then
142 <
140 >    
141 >    if (do_stress) then
142 >      
143   #ifdef IS_MPI
144 <          id1 = tagRow(atom1)
145 <          id2 = tagColumn(atom2)
144 >       id1 = tagRow(atom1)
145 >       id2 = tagColumn(atom2)
146   #else
147 <          id1 = atom1
148 <          id2 = atom2
147 >       id1 = atom1
148 >       id2 = atom2
149   #endif
255 –
256 –
257 –          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
258 –
259 –             ! because the d vector is the rj - ri vector, and
260 –             ! because fx, fy, fz are the force on atom i, we need a
261 –             ! negative sign here:
262 –            
263 –             tau_Temp(1) = tau_Temp(1) - d(1) * fx
264 –             tau_Temp(2) = tau_Temp(2) - d(1) * fy
265 –             tau_Temp(3) = tau_Temp(3) - d(1) * fz
266 –             tau_Temp(4) = tau_Temp(4) - d(2) * fx
267 –             tau_Temp(5) = tau_Temp(5) - d(2) * fy
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))
286 –            
287 –          endif
288 –       endif
150        
151 +      
152 +       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
153 +          
154 +          ! because the d vector is the rj - ri vector, and
155 +          ! because fx, fy, fz are the force on atom i, we need a
156 +          ! negative sign here:
157 +          
158 +          tau_Temp(1) = tau_Temp(1) - d(1) * fx
159 +          tau_Temp(2) = tau_Temp(2) - d(1) * fy
160 +          tau_Temp(3) = tau_Temp(3) - d(1) * fz
161 +          tau_Temp(4) = tau_Temp(4) - d(2) * fx
162 +          tau_Temp(5) = tau_Temp(5) - d(2) * fy
163 +          tau_Temp(6) = tau_Temp(6) - d(2) * fz
164 +          tau_Temp(8) = tau_Temp(8) - d(3) * fy
165 +          tau_Temp(9) = tau_Temp(9) - d(3) * fz
166 +          
167 +          
168 +          virial_Temp = virial_Temp + &
169 +               (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
170 +          
171 +       endif
172      endif
173 +    
174 +    
175      return
176    end subroutine do_charge_pair
177    

Diff Legend

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