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 388 by chuckv, Fri Mar 21 22:11:50 2003 UTC vs.
Revision 1160 by gezelter, Tue May 11 21:31:15 2004 UTC

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

Diff Legend

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