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: 329
Committed: Wed Mar 12 22:27:59 2003 UTC (21 years, 5 months ago) by gezelter
File size: 5679 byte(s)
Log Message: 
Stick a fork in it.  It's rare.

File Contents

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