ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/calc_dipole_dipole.F90
Revision: 328
Committed: Wed Mar 12 20:00:58 2003 UTC (21 years, 5 months ago) by gezelter
File size: 5698 byte(s)
Log Message: 
bye bye atypeGlobals (part2)

File Contents

# User Rev Content
1 chuckv 307 module dipole_dipole
2    
3     use simulation
4     use definitions
5     use forceGlobals
6 gezelter 328 use atype_module
7 gezelter 317 use vector_class
8 gezelter 309 #ifdef IS_MPI
9     use mpiSimulation
10     #endif
11     implicit none
12 chuckv 307
13     contains
14    
15 gezelter 326 subroutine do_dipole_pair(atom1, atom2, d, rij, pot, u_l, f, t, &
16     do_pot, do_stress)
17 chuckv 307
18 gezelter 326 logical :: do_pot, do_stress
19    
20 gezelter 317 integer atom1, atom2, me1, me2
21     double precision rij, mu1, mu2
22 chuckv 307 double precision dfact1, dfact2, dip2, r2, r3, r5, pre
23     double precision dudx, dudy, dudz, dudu1x, dudu1y, dudu1z
24     double precision dudu2x, dudu2y, dudu2z, rdotu1, rdotu2, u1dotu2
25     double precision taper, dtdr, vterm
26 gezelter 309
27     real( kind = dp ) :: pot, rt, rrf
28 gezelter 317 real( kind = dp ), dimension(3) :: d
29 gezelter 309 real( kind = dp ), dimension(3,getNlocal()) :: u_l
30     real( kind = dp ), dimension(3,getNlocal()) :: f
31     real( kind = dp ), dimension(3,getNlocal()) :: t
32 chuckv 307
33     real (kind = dp), dimension(3) :: ul1
34     real (kind = dp), dimension(3) :: ul2
35 gezelter 317
36 chuckv 307
37 gezelter 309 #ifdef IS_MPI
38 gezelter 317 me1 = atid_Row(atom1)
39 gezelter 309 ul1(1) = u_l_Row(1,atom1)
40     ul1(2) = u_l_Row(2,atom1)
41     ul1(3) = u_l_Row(3,atom1)
42    
43 gezelter 317 me2 = atid_Col(atom2)
44 gezelter 309 ul2(1) = u_l_Col(1,atom2)
45     ul2(2) = u_l_Col(2,atom2)
46     ul2(3) = u_l_Col(3,atom2)
47     #else
48 gezelter 317 me1 = atid(atom1)
49 gezelter 309 ul1(1) = u_l(1,atom1)
50     ul1(2) = u_l(2,atom1)
51     ul1(3) = u_l(3,atom1)
52    
53 gezelter 317 me2 = atid(atom2)
54     ul2(1) = u_l(1,atom2)
55     ul2(2) = u_l(2,atom2)
56     ul2(3) = u_l(3,atom2)
57 gezelter 309 #endif
58    
59 gezelter 317 call getElementProperty(atypes, me1, "dipole_moment", mu1)
60     call getElementProperty(atypes, me2, "dipole_moment", mu2)
61    
62 gezelter 312 rrf = getRrf()
63     rt = getRt()
64    
65 chuckv 307 ! pre converts from mu in units of debye to kcal/mol
66     pre = 14.38362d0
67    
68     if (rij.le.rrf) then
69    
70     if (rij.lt.rt) then
71     taper = 1.0d0
72     dtdr = 0.0d0
73     else
74     taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3)
75     dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
76     endif
77    
78     r2 = rij*rij
79     r3 = r2*rij
80     r5 = r3*r2
81    
82 gezelter 317 rdotu1 = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
83     rdotu2 = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
84 chuckv 307 u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
85    
86 gezelter 317 dip2 = pre * mu1 * mu2
87 chuckv 307
88     dfact1 = 3.0d0*dip2 / r2
89     dfact2 = 3.0d0*dip2 / r5
90    
91     vterm = dip2*((u1dotu2/r3) - 3.0d0*(rdotu1*rdotu2/r5))
92    
93 gezelter 326 if (do_pot) then
94 gezelter 309 #ifdef IS_MPI
95 gezelter 326 pot_row(atom1) = pot_row(atom1) + 0.5d0*vterm*taper
96     pot_col(atom2) = pot_col(atom2) + 0.5d0*vterm*taper
97 gezelter 309 #else
98 gezelter 326 pot = pot + vterm*taper
99 gezelter 309 #endif
100 gezelter 326 endif
101 chuckv 307
102 gezelter 317 dudx = (-dfact1 * d(1) * ((u1dotu2/r3) - &
103 chuckv 307 (5.0d0*(rdotu1*rdotu2)/r5)) - &
104     dfact2*(ul1(1)*rdotu2 + ul2(1)*rdotu1))*taper + &
105 gezelter 317 vterm*dtdr*d(1)/rij
106 chuckv 307
107 gezelter 317 dudy = (-dfact1 * d(2) * ((u1dotu2/r3) - &
108 chuckv 307 (5.0d0*(rdotu1*rdotu2)/r5)) - &
109     dfact2*(ul1(2)*rdotu2 + ul2(2)*rdotu1))*taper + &
110 gezelter 317 vterm*dtdr*d(2)/rij
111 chuckv 307
112 gezelter 317 dudz = (-dfact1 * d(3) * ((u1dotu2/r3) - &
113 chuckv 307 (5.0d0*(rdotu1*rdotu2)/r5)) - &
114     dfact2*(ul1(3)*rdotu2 + ul2(3)*rdotu1))*taper + &
115 gezelter 317 vterm*dtdr*d(3)/rij
116 chuckv 307
117 gezelter 317 dudu1x = (dip2*((ul2(1)/r3) - (3.0d0*d(1)*rdotu2/r5)))*taper
118     dudu1y = (dip2*((ul2(2)/r3) - (3.0d0*d(2)*rdotu2/r5)))*taper
119     dudu1z = (dip2*((ul2(3)/r3) - (3.0d0*d(3)*rdotu2/r5)))*taper
120 chuckv 307
121 gezelter 317 dudu2x = (dip2*((ul1(1)/r3) - (3.0d0*d(1)*rdotu1/r5)))*taper
122     dudu2y = (dip2*((ul1(2)/r3) - (3.0d0*d(2)*rdotu1/r5)))*taper
123     dudu2z = (dip2*((ul1(3)/r3) - (3.0d0*d(3)*rdotu1/r5)))*taper
124 chuckv 307
125    
126     #ifdef IS_MPI
127 gezelter 309 f_Row(1,atom1) = f_Row(1,atom1) + dudx
128     f_Row(2,atom1) = f_Row(2,atom1) + dudy
129     f_Row(3,atom1) = f_Row(3,atom1) + dudz
130 chuckv 307
131 gezelter 309 f_Col(1,atom2) = f_Col(1,atom2) - dudx
132     f_Col(2,atom2) = f_Col(2,atom2) - dudy
133     f_Col(3,atom2) = f_Col(3,atom2) - dudz
134 chuckv 307
135 gezelter 309 t_Row(1,atom1) = t_Row(1,atom1) - ul1(2)*dudu1z + ul1(3)*dudu1y
136     t_Row(2,atom1) = t_Row(2,atom1) - ul1(3)*dudu1x + ul1(1)*dudu1z
137     t_Row(3,atom1) = t_Row(3,atom1) - ul1(1)*dudu1y + ul1(2)*dudu1x
138 chuckv 307
139 gezelter 309 t_Col(1,atom2) = t_Col(1,atom2) - ul2(2)*dudu2z + ul2(3)*dudu2y
140     t_Col(2,atom2) = t_Col(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
141     t_Col(3,atom2) = t_Col(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
142 chuckv 307 #else
143 gezelter 309 f(1,atom1) = f(1,atom1) + dudx
144     f(2,atom1) = f(2,atom1) + dudy
145     f(3,atom1) = f(3,atom1) + dudz
146 chuckv 307
147 gezelter 309 f(1,atom2) = f(1,atom2) - dudx
148     f(2,atom2) = f(2,atom2) - dudy
149     f(3,atom2) = f(3,atom2) - dudz
150 chuckv 307
151 gezelter 309 t(1,atom1) = t(1,atom1) - ul1(2)*dudu1z + ul1(3)*dudu1y
152     t(2,atom1) = t(2,atom1) - ul1(3)*dudu1x + ul1(1)*dudu1z
153     t(3,atom1) = t(3,atom1) - ul1(1)*dudu1y + ul1(2)*dudu1x
154 chuckv 307
155 gezelter 309 t(1,atom2) = t(1,atom2) - ul2(2)*dudu2z + ul2(3)*dudu2y
156     t(2,atom2) = t(2,atom2) - ul2(3)*dudu2x + ul2(1)*dudu2z
157     t(3,atom2) = t(3,atom2) - ul2(1)*dudu2y + ul2(2)*dudu2x
158 chuckv 307 #endif
159    
160 gezelter 326 if (do_stress) then
161 gezelter 317 tau_Temp(1) = tau_Temp(1) + dudx * d(1)
162     tau_Temp(2) = tau_Temp(2) + dudx * d(2)
163     tau_Temp(3) = tau_Temp(3) + dudx * d(3)
164     tau_Temp(4) = tau_Temp(4) + dudy * d(1)
165     tau_Temp(5) = tau_Temp(5) + dudy * d(2)
166     tau_Temp(6) = tau_Temp(6) + dudy * d(3)
167     tau_Temp(7) = tau_Temp(7) + dudz * d(1)
168     tau_Temp(8) = tau_Temp(8) + dudz * d(2)
169     tau_Temp(9) = tau_Temp(9) + dudz * d(3)
170 gezelter 309 virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
171 chuckv 307 endif
172    
173     endif
174    
175     return
176 gezelter 309 end subroutine do_dipole_pair
177    
178 chuckv 307 end module dipole_dipole