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

Comparing trunk/OOPSE/libmdtools/calc_dipole_dipole.F90 (file contents):
Revision 483 by gezelter, Wed Apr 9 04:06:43 2003 UTC vs.
Revision 1150 by gezelter, Fri May 7 21:35:05 2004 UTC

# Line 5 | Line 5 | module dipole_dipole
5    use atype_module
6    use vector_class
7    use simulation
8 +  use status
9   #ifdef IS_MPI
10    use mpiSimulation
11   #endif
12    implicit none
13  
14 <  real(kind=dp), save :: rrf = 0.0
14 >  PRIVATE
15 >  real(kind=dp), save :: ecr = 0.0
16    real(kind=dp), save :: rt  = 0.0
17 <  real(kind=dp), save :: rrfsq = 0.0
18 <  real(kind=dp), save :: pre = 0.0
19 <  logical, save :: dipole_initialized = .false.
17 >   real(kind=dp), save :: pre = 0.0
18 >  logical, save :: haveCutoffs = .false.
19 >  logical, save :: haveMomentMap = .false.
20  
21 +  public::setCutoffsDipole
22 +  public::do_dipole_pair
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   contains
31      
32 <  subroutine initialize_dipole(this_rrf, this_rt)
33 <    real(kind=dp), intent(in) :: this_rrf, this_rt
34 <    rrf = this_rrf
32 >  subroutine setCutoffsDipole(this_ecr, this_rt)
33 >    real(kind=dp), intent(in) :: this_ecr, this_rt
34 >    ecr = this_ecr
35      rt = this_rt    
36 <    rrfsq = rrf * rrf
37 <    ! pre converts from mu in units of debye to kcal/mol
36 >
37 >      ! pre converts from mu in units of debye to kcal/mol
38      pre = 14.38362_dp
39  
40 <    dipole_initialized = .true.
40 >    haveCutoffs = .true.
41      
42      return
43 <  end subroutine initialize_dipole
43 >  end subroutine setCutoffsDipole
44  
45 +  subroutine createMomentMap(status)
46 +    integer :: nAtypes
47 +    integer :: status
48 +    integer :: i
49 +    real (kind=DP) :: thisDP
50 +    logical :: thisProperty
51  
52 <  subroutine do_dipole_pair(atom1, atom2, d, rij, r2, pot, u_l, f, t, &
53 <       do_pot, do_stress)
52 >    status = 0
53 >
54 >    nAtypes = getSize(atypes)
55      
56 +    if (nAtypes == 0) then
57 +       status = -1
58 +       return
59 +    end if
60 +    
61 +    if (.not. allocated(MomentMap)) then
62 +       allocate(MomentMap(nAtypes))
63 +    endif
64 +
65 +    do i = 1, nAtypes
66 +
67 +       call getElementProperty(atypes, i, "is_DP", thisProperty)
68 +
69 +       if (thisProperty) then
70 +          call getElementProperty(atypes, i, "dipole_moment", thisDP)
71 +          MomentMap(i)%dipole_moment = thisDP
72 +       endif
73 +      
74 +    end do
75 +    
76 +    haveMomentMap = .true.
77 +    
78 +  end subroutine createMomentMap
79 +  
80 +  subroutine do_dipole_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, &
81 +       u_l, f, t, do_pot, do_stress)
82 +    
83      logical :: do_pot, do_stress
84  
85 <    integer atom1, atom2, me1, me2
85 >    integer atom1, atom2, me1, me2, id1, id2
86 >    integer :: localError
87      real(kind=dp) :: rij, mu1, mu2
88      real(kind=dp) :: dfact1, dfact2, dip2, r2, r3, r5
89      real(kind=dp) :: dudx, dudy, dudz, dudu1x, dudu1y, dudu1z
90      real(kind=dp) :: dudu2x, dudu2y, dudu2z, rdotu1, rdotu2, u1dotu2
91 <    real(kind=dp) :: taper, dtdr, vterm
91 >    real(kind=dp) :: sw, vpair, vterm
92  
93      real( kind = dp ) :: pot
94      real( kind = dp ), dimension(3) :: d
95 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
96 <    real( kind = dp ), dimension(3,getNlocal()) :: f
97 <    real( kind = dp ), dimension(3,getNlocal()) :: t
95 >    real( kind = dp ), dimension(3,nLocal) :: u_l
96 >    real( kind = dp ), dimension(3,nLocal) :: f
97 >    real( kind = dp ), dimension(3,nLocal) :: t
98      
99      real (kind = dp), dimension(3) :: ul1
100      real (kind = dp), dimension(3) :: ul2
101  
102 <    if (.not.dipole_initialized) then
103 <       write(default_error,*) 'Dipole-dipole not initialized!'
102 >    if (.not. haveCutoffs) then
103 >       write(default_error,*) 'Dipole-dipole does not have cutoffs set!'
104         return
105      endif
106 <    
106 >
107 >    if (.not.haveMomentMap) then
108 >       localError = 0
109 >       call createMomentMap(localError)
110 >       if ( localError .ne. 0 ) then
111 >          call handleError("dipole-dipole", "MomentMap creation failed!")
112 >          return
113 >       end if
114 >    endif
115 >
116   #ifdef IS_MPI
117      me1 = atid_Row(atom1)
118      ul1(1) = u_l_Row(1,atom1)
# Line 80 | Line 135 | contains
135      ul2(3) = u_l(3,atom2)
136   #endif
137  
138 <    if( atom1 .eq. 2 )then
139 <       write (0,*) 'ul =', ul1(1), ul1(2), ul1(3)
140 <    endif
141 <
142 <    if( atom2 .eq. 2 )then
143 <       write (0,*) 'ul =', ul2(1), ul2(2), ul2(3)
144 <    endif
145 <
146 <
147 <    call getElementProperty(atypes, me1, "dipole_moment", mu1)
148 <    call getElementProperty(atypes, me2, "dipole_moment", mu2)
149 <
150 <    if (rij.le.rrf) then
151 <      
152 <       if (rij.lt.rt) then
153 <          taper = 1.0d0
154 <          dtdr = 0.0d0
155 <       else
156 <          taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3)
102 <          dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
103 <       endif
104 <
105 <       r3 = r2*rij
106 <       r5 = r3*r2
107 <      
108 <       rdotu1 = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
109 <       rdotu2 = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
110 <       u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
111 <
112 <       dip2 = pre * mu1 * mu2
113 <       dfact1 = 3.0d0*dip2 / r2
114 <       dfact2 = 3.0d0*dip2 / r5
115 <
116 <       vterm = dip2*((u1dotu2/r3) - 3.0d0*(rdotu1*rdotu2/r5))
117 <
118 <       if (do_pot) then
138 >    mu1 = MomentMap(me1)%dipole_moment
139 >    mu2 = MomentMap(me2)%dipole_moment
140 >    
141 >    r3 = r2*rij
142 >    r5 = r3*r2
143 >    
144 >    rdotu1 = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
145 >    rdotu2 = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
146 >    u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
147 >    
148 >    dip2 = pre * mu1 * mu2
149 >    dfact1 = 3.0d0*dip2 / r2
150 >    dfact2 = 3.0d0*dip2 / r5
151 >    
152 >    vterm = dip2*((u1dotu2/r3) - 3.0d0*(rdotu1*rdotu2/r5))
153 >    
154 >    vpair = vpair + vterm*sw
155 >    
156 >    if (do_pot) then
157   #ifdef IS_MPI
158 <          pot_row(atom1) = pot_row(atom1) + 0.5d0*vterm*taper
159 <          pot_col(atom2) = pot_col(atom2) + 0.5d0*vterm*taper
158 >       pot_row(atom1) = pot_row(atom1) + 0.5d0*vterm*sw
159 >       pot_col(atom2) = pot_col(atom2) + 0.5d0*vterm*sw
160   #else
161 <          pot = pot + vterm*taper
161 >       pot = pot + vterm*sw
162   #endif
163 <       endif
164 <      
165 <       dudx = (-dfact1 * d(1) * ((u1dotu2/r3) - &
166 <            (5.0d0*(rdotu1*rdotu2)/r5)) -  &
167 <            dfact2*(ul1(1)*rdotu2 + ul2(1)*rdotu1))*taper + &
168 <            vterm*dtdr*d(1)/rij
169 <      
170 <       dudy = (-dfact1 * d(2) * ((u1dotu2/r3) - &
171 <            (5.0d0*(rdotu1*rdotu2)/r5)) -  &
172 <            dfact2*(ul1(2)*rdotu2 + ul2(2)*rdotu1))*taper + &
173 <            vterm*dtdr*d(2)/rij
174 <      
175 <       dudz = (-dfact1 * d(3) * ((u1dotu2/r3) - &
176 <            (5.0d0*(rdotu1*rdotu2)/r5)) -  &
177 <            dfact2*(ul1(3)*rdotu2 + ul2(3)*rdotu1))*taper +  &
178 <            vterm*dtdr*d(3)/rij
179 <
180 <       dudu1x = (dip2*((ul2(1)/r3) - (3.0d0*d(1)*rdotu2/r5)))*taper
181 <       dudu1y = (dip2*((ul2(2)/r3) - (3.0d0*d(2)*rdotu2/r5)))*taper
182 <       dudu1z = (dip2*((ul2(3)/r3) - (3.0d0*d(3)*rdotu2/r5)))*taper
183 <
184 <       dudu2x = (dip2*((ul1(1)/r3) - (3.0d0*d(1)*rdotu1/r5)))*taper
185 <       dudu2y = (dip2*((ul1(2)/r3) - (3.0d0*d(2)*rdotu1/r5)))*taper
148 <       dudu2z = (dip2*((ul1(3)/r3) - (3.0d0*d(3)*rdotu1/r5)))*taper
149 <      
150 <
163 >    endif
164 >    
165 >    dudx = (-dfact1 * d(1) * ((u1dotu2/r3) - &
166 >         (5.0d0*(rdotu1*rdotu2)/r5)) -  &
167 >            dfact2*(ul1(1)*rdotu2 + ul2(1)*rdotu1))*sw
168 >    
169 >    dudy = (-dfact1 * d(2) * ((u1dotu2/r3) - &
170 >         (5.0d0*(rdotu1*rdotu2)/r5)) -  &
171 >            dfact2*(ul1(2)*rdotu2 + ul2(2)*rdotu1))*sw
172 >    
173 >    dudz = (-dfact1 * d(3) * ((u1dotu2/r3) - &
174 >         (5.0d0*(rdotu1*rdotu2)/r5)) -  &
175 >         dfact2*(ul1(3)*rdotu2 + ul2(3)*rdotu1))*sw
176 >    
177 >    dudu1x = (dip2*((ul2(1)/r3) - (3.0d0*d(1)*rdotu2/r5)))*sw
178 >    dudu1y = (dip2*((ul2(2)/r3) - (3.0d0*d(2)*rdotu2/r5)))*sw
179 >    dudu1z = (dip2*((ul2(3)/r3) - (3.0d0*d(3)*rdotu2/r5)))*sw
180 >    
181 >    dudu2x = (dip2*((ul1(1)/r3) - (3.0d0*d(1)*rdotu1/r5)))*sw
182 >    dudu2y = (dip2*((ul1(2)/r3) - (3.0d0*d(2)*rdotu1/r5)))*sw
183 >    dudu2z = (dip2*((ul1(3)/r3) - (3.0d0*d(3)*rdotu1/r5)))*sw
184 >    
185 >    
186   #ifdef IS_MPI
187 <       f_Row(1,atom1) = f_Row(1,atom1) + dudx
188 <       f_Row(2,atom1) = f_Row(2,atom1) + dudy
189 <       f_Row(3,atom1) = f_Row(3,atom1) + dudz
190 <
191 <       f_Col(1,atom2) = f_Col(1,atom2) - dudx
192 <       f_Col(2,atom2) = f_Col(2,atom2) - dudy
193 <       f_Col(3,atom2) = f_Col(3,atom2) - dudz
194 <      
195 <       t_Row(1,atom1) = t_Row(1,atom1) - ul1(2)*dudu1z + ul1(3)*dudu1y
196 <       t_Row(2,atom1) = t_Row(2,atom1) - ul1(3)*dudu1x + ul1(1)*dudu1z
197 <       t_Row(3,atom1) = t_Row(3,atom1) - ul1(1)*dudu1y + ul1(2)*dudu1x
198 <      
199 <       t_Col(1,atom2) = t_Col(1,atom2) - ul2(2)*dudu2z + ul2(3)*dudu2y
200 <       t_Col(2,atom2) = t_Col(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
201 <       t_Col(3,atom2) = t_Col(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
187 >    f_Row(1,atom1) = f_Row(1,atom1) + dudx
188 >    f_Row(2,atom1) = f_Row(2,atom1) + dudy
189 >    f_Row(3,atom1) = f_Row(3,atom1) + dudz
190 >    
191 >    f_Col(1,atom2) = f_Col(1,atom2) - dudx
192 >    f_Col(2,atom2) = f_Col(2,atom2) - dudy
193 >    f_Col(3,atom2) = f_Col(3,atom2) - dudz
194 >    
195 >    t_Row(1,atom1) = t_Row(1,atom1) - ul1(2)*dudu1z + ul1(3)*dudu1y
196 >    t_Row(2,atom1) = t_Row(2,atom1) - ul1(3)*dudu1x + ul1(1)*dudu1z
197 >    t_Row(3,atom1) = t_Row(3,atom1) - ul1(1)*dudu1y + ul1(2)*dudu1x
198 >    
199 >    t_Col(1,atom2) = t_Col(1,atom2) - ul2(2)*dudu2z + ul2(3)*dudu2y
200 >    t_Col(2,atom2) = t_Col(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
201 >    t_Col(3,atom2) = t_Col(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
202   #else
203 <       f(1,atom1) = f(1,atom1) + dudx
204 <       f(2,atom1) = f(2,atom1) + dudy
205 <       f(3,atom1) = f(3,atom1) + dudz
206 <      
207 <       f(1,atom2) = f(1,atom2) - dudx
208 <       f(2,atom2) = f(2,atom2) - dudy
209 <       f(3,atom2) = f(3,atom2) - dudz
210 <      
211 <       t(1,atom1) = t(1,atom1) - ul1(2)*dudu1z + ul1(3)*dudu1y
212 <       t(2,atom1) = t(2,atom1) - ul1(3)*dudu1x + ul1(1)*dudu1z
213 <       t(3,atom1) = t(3,atom1) - ul1(1)*dudu1y + ul1(2)*dudu1x
214 <      
215 <       t(1,atom2) = t(1,atom2) - ul2(2)*dudu2z + ul2(3)*dudu2y
216 <       t(2,atom2) = t(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
217 <       t(3,atom2) = t(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
203 >    f(1,atom1) = f(1,atom1) + dudx
204 >    f(2,atom1) = f(2,atom1) + dudy
205 >    f(3,atom1) = f(3,atom1) + dudz
206 >    
207 >    f(1,atom2) = f(1,atom2) - dudx
208 >    f(2,atom2) = f(2,atom2) - dudy
209 >    f(3,atom2) = f(3,atom2) - dudz
210 >    
211 >    t(1,atom1) = t(1,atom1) - ul1(2)*dudu1z + ul1(3)*dudu1y
212 >    t(2,atom1) = t(2,atom1) - ul1(3)*dudu1x + ul1(1)*dudu1z
213 >    t(3,atom1) = t(3,atom1) - ul1(1)*dudu1y + ul1(2)*dudu1x
214 >    
215 >    t(1,atom2) = t(1,atom2) - ul2(2)*dudu2z + ul2(3)*dudu2y
216 >    t(2,atom2) = t(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
217 >    t(3,atom2) = t(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
218   #endif
219 +    
220 +    if (do_stress) then          
221        
222 <       if (do_stress) then          
223 <
224 <          if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then
225 <
226 <             tau_Temp(1) = tau_Temp(1) + dudx * d(1)
227 <             tau_Temp(2) = tau_Temp(2) + dudx * d(2)
228 <             tau_Temp(3) = tau_Temp(3) + dudx * d(3)
229 <             tau_Temp(4) = tau_Temp(4) + dudy * d(1)
230 <             tau_Temp(5) = tau_Temp(5) + dudy * d(2)
231 <             tau_Temp(6) = tau_Temp(6) + dudy * d(3)
232 <             tau_Temp(7) = tau_Temp(7) + dudz * d(1)
233 <             tau_Temp(8) = tau_Temp(8) + dudz * d(2)
234 <             tau_Temp(9) = tau_Temp(9) + dudz * d(3)
235 <             virial_Temp = virial_Temp + &
236 <                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
237 <
238 <          endif          
222 > #ifdef IS_MPI
223 >       id1 = tagRow(atom1)
224 >       id2 = tagColumn(atom2)
225 > #else
226 >       id1 = atom1
227 >       id2 = atom2
228 > #endif                
229 >       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
230 >          
231 >          ! because the d vector is the rj - ri vector, and
232 >          ! because dudx, dudy, dudz are the (positive) force on
233 >          ! atom i (negative on atom j) we need a negative sign here:
234 >          
235 >          tau_Temp(1) = tau_Temp(1) - d(1) * dudx
236 >          tau_Temp(2) = tau_Temp(2) - d(1) * dudy
237 >          tau_Temp(3) = tau_Temp(3) - d(1) * dudz
238 >          tau_Temp(4) = tau_Temp(4) - d(2) * dudx
239 >          tau_Temp(5) = tau_Temp(5) - d(2) * dudy
240 >          tau_Temp(6) = tau_Temp(6) - d(2) * dudz
241 >          tau_Temp(7) = tau_Temp(7) - d(3) * dudx
242 >          tau_Temp(8) = tau_Temp(8) - d(3) * dudy
243 >          tau_Temp(9) = tau_Temp(9) - d(3) * dudz
244 >          
245 >          virial_Temp = virial_Temp + &
246 >               (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
247 >          
248         endif
203 –      
249      endif
250      
251      return

Diff Legend

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