ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/gb.F90
Revision: 1608
Committed: Wed Oct 20 04:02:48 2004 UTC (19 years, 8 months ago) by gezelter
File size: 12268 byte(s)
Log Message:
name sanity on the fortran side

File Contents

# Content
1 module gb_pair
2 use force_globals
3 use definitions
4 use simulation
5 #ifdef IS_MPI
6 use mpiSimulation
7 #endif
8
9 implicit none
10
11 PRIVATE
12
13 logical, save :: gb_pair_initialized = .false.
14 real(kind=dp), save :: gb_sigma
15 real(kind=dp), save :: gb_l2b_ratio
16 real(kind=dp), save :: gb_eps
17 real(kind=dp), save :: gb_eps_ratio
18 real(kind=dp), save :: gb_mu
19 real(kind=dp), save :: gb_nu
20
21 public :: check_gb_pair_FF
22 public :: set_gb_pair_params
23 public :: do_gb_pair
24
25 contains
26
27 subroutine check_gb_pair_FF(status)
28 integer :: status
29 status = -1
30 if (gb_pair_initialized) status = 0
31 return
32 end subroutine check_gb_pair_FF
33
34 subroutine set_gb_pair_params(sigma, l2b_ratio, eps, eps_ratio, mu, nu)
35 real( kind = dp ), intent(in) :: sigma, l2b_ratio, eps, eps_ratio
36 real( kind = dp ), intent(in) :: mu, nu
37
38 gb_sigma = sigma
39 gb_l2b_ratio = l2b_ratio
40 gb_eps = eps
41 gb_eps_ratio = eps_ratio
42 gb_mu = mu
43 gb_nu = nu
44
45 gb_pair_initialized = .true.
46 return
47 end subroutine set_gb_pair_params
48
49
50 subroutine do_gb_pair(atom1, atom2, d, r, r2, sw, vpair, fpair, &
51 pot, u_l, f, t, do_pot)
52
53 integer, intent(in) :: atom1, atom2
54 integer :: id1, id2
55 real (kind=dp), intent(inout) :: r, r2
56 real (kind=dp), dimension(3), intent(in) :: d
57 real (kind=dp), dimension(3), intent(inout) :: fpair
58 real (kind=dp) :: pot, sw, vpair
59 real (kind=dp), dimension(3,nLocal) :: u_l
60 real (kind=dp), dimension(3,nLocal) :: f
61 real (kind=dp), dimension(3,nLocal) :: t
62 logical, intent(in) :: do_pot
63 real (kind = dp), dimension(3) :: ul1
64 real (kind = dp), dimension(3) :: ul2
65
66 real(kind=dp) :: chi, chiprime, emu, s2
67 real(kind=dp) :: r4, rdotu1, rdotu2, u1dotu2, g, gp, gpi, gmu, gmum
68 real(kind=dp) :: curlyE, enu, enum, eps, dotsum, dotdiff, ds2, dd2
69 real(kind=dp) :: opXdot, omXdot, opXpdot, omXpdot, pref, gfact
70 real(kind=dp) :: BigR, Ri, Ri2, Ri6, Ri7, Ri12, Ri13, R126, R137
71 real(kind=dp) :: dru1dx, dru1dy, dru1dz
72 real(kind=dp) :: dru2dx, dru2dy, dru2dz
73 real(kind=dp) :: dBigRdx, dBigRdy, dBigRdz
74 real(kind=dp) :: dBigRdu1x, dBigRdu1y, dBigRdu1z
75 real(kind=dp) :: dBigRdu2x, dBigRdu2y, dBigRdu2z
76 real(kind=dp) :: dUdx, dUdy, dUdz
77 real(kind=dp) :: dUdu1x, dUdu1y, dUdu1z, dUdu2x, dUdu2y, dUdu2z
78 real(kind=dp) :: dcE, dcEdu1x, dcEdu1y, dcEdu1z, dcEdu2x, dcEdu2y, dcEdu2z
79 real(kind=dp) :: depsdu1x, depsdu1y, depsdu1z, depsdu2x, depsdu2y, depsdu2z
80 real(kind=dp) :: drdx, drdy, drdz
81 real(kind=dp) :: dgdx, dgdy, dgdz
82 real(kind=dp) :: dgdu1x, dgdu1y, dgdu1z, dgdu2x, dgdu2y, dgdu2z
83 real(kind=dp) :: dgpdx, dgpdy, dgpdz
84 real(kind=dp) :: dgpdu1x, dgpdu1y, dgpdu1z, dgpdu2x, dgpdu2y, dgpdu2z
85 real(kind=dp) :: line1a, line1bx, line1by, line1bz
86 real(kind=dp) :: line2a, line2bx, line2by, line2bz
87 real(kind=dp) :: line3a, line3b, line3, line3x, line3y, line3z
88 real(kind=dp) :: term1x, term1y, term1z, term1u1x, term1u1y, term1u1z
89 real(kind=dp) :: term1u2x, term1u2y, term1u2z
90 real(kind=dp) :: term2a, term2b, term2u1x, term2u1y, term2u1z
91 real(kind=dp) :: term2u2x, term2u2y, term2u2z
92 real(kind=dp) :: yick1, yick2, mess1, mess2
93
94 s2 = (gb_l2b_ratio)**2
95 emu = (gb_eps_ratio)**(1.0d0/gb_mu)
96
97 chi = (s2 - 1.0d0)/(s2 + 1.0d0)
98 chiprime = (1.0d0 - emu)/(1.0d0 + emu)
99
100 r4 = r2*r2
101
102 #ifdef IS_MPI
103 ul1(1) = u_l_Row(1,atom1)
104 ul1(2) = u_l_Row(2,atom1)
105 ul1(3) = u_l_Row(3,atom1)
106
107 ul2(1) = u_l_Col(1,atom2)
108 ul2(2) = u_l_Col(2,atom2)
109 ul2(3) = u_l_Col(3,atom2)
110 #else
111 ul1(1) = u_l(1,atom1)
112 ul1(2) = u_l(2,atom1)
113 ul1(3) = u_l(3,atom1)
114
115 ul2(1) = u_l(1,atom2)
116 ul2(2) = u_l(2,atom2)
117 ul2(3) = u_l(3,atom2)
118 #endif
119
120 dru1dx = ul1(1)
121 dru2dx = ul2(1)
122 dru1dy = ul1(2)
123 dru2dy = ul2(2)
124 dru1dz = ul1(3)
125 dru2dz = ul2(3)
126
127 drdx = d(1) / r
128 drdy = d(2) / r
129 drdz = d(3) / r
130
131 ! do some dot products:
132 ! NB the r in these dot products is the actual intermolecular vector,
133 ! and is not the unit vector in that direction.
134
135 rdotu1 = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
136 rdotu2 = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
137 u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
138
139 ! This stuff is all for the calculation of g(Chi) and dgdx
140 ! Line numbers roughly follow the lines in equation A25 of Luckhurst
141 ! et al. Liquid Crystals 8, 451-464 (1990).
142 ! We note however, that there are some major typos in that Appendix
143 ! of the Luckhurst paper, particularly in equations A23, A29 and A31
144 ! We have attempted to correct them below.
145
146 dotsum = rdotu1+rdotu2
147 dotdiff = rdotu1-rdotu2
148 ds2 = dotsum*dotsum
149 dd2 = dotdiff*dotdiff
150
151 opXdot = 1.0d0 + Chi*u1dotu2
152 omXdot = 1.0d0 - Chi*u1dotu2
153 opXpdot = 1.0d0 + ChiPrime*u1dotu2
154 omXpdot = 1.0d0 - ChiPrime*u1dotu2
155
156 line1a = dotsum/opXdot
157 line1bx = dru1dx + dru2dx
158 line1by = dru1dy + dru2dy
159 line1bz = dru1dz + dru2dz
160
161 line2a = dotdiff/omXdot
162 line2bx = dru1dx - dru2dx
163 line2by = dru1dy - dru2dy
164 line2bz = dru1dz - dru2dz
165
166 term1x = -Chi*(line1a*line1bx + line2a*line2bx)/r2
167 term1y = -Chi*(line1a*line1by + line2a*line2by)/r2
168 term1z = -Chi*(line1a*line1bz + line2a*line2bz)/r2
169
170 line3a = ds2/opXdot
171 line3b = dd2/omXdot
172 line3 = Chi*(line3a + line3b)/r4
173 line3x = d(1)*line3
174 line3y = d(2)*line3
175 line3z = d(3)*line3
176
177 dgdx = term1x + line3x
178 dgdy = term1y + line3y
179 dgdz = term1z + line3z
180
181 term1u1x = 2.0d0*(line1a+line2a)*d(1)
182 term1u1y = 2.0d0*(line1a+line2a)*d(2)
183 term1u1z = 2.0d0*(line1a+line2a)*d(3)
184 term1u2x = 2.0d0*(line1a-line2a)*d(1)
185 term1u2y = 2.0d0*(line1a-line2a)*d(2)
186 term1u2z = 2.0d0*(line1a-line2a)*d(3)
187
188 term2a = -line3a/opXdot
189 term2b = line3b/omXdot
190
191 term2u1x = Chi*ul2(1)*(term2a + term2b)
192 term2u1y = Chi*ul2(2)*(term2a + term2b)
193 term2u1z = Chi*ul2(3)*(term2a + term2b)
194 term2u2x = Chi*ul1(1)*(term2a + term2b)
195 term2u2y = Chi*ul1(2)*(term2a + term2b)
196 term2u2z = Chi*ul1(3)*(term2a + term2b)
197
198 pref = -Chi*0.5d0/r2
199
200 dgdu1x = pref*(term1u1x+term2u1x)
201 dgdu1y = pref*(term1u1y+term2u1y)
202 dgdu1z = pref*(term1u1z+term2u1z)
203 dgdu2x = pref*(term1u2x+term2u2x)
204 dgdu2y = pref*(term1u2y+term2u2y)
205 dgdu2z = pref*(term1u2z+term2u2z)
206
207 g = 1.0d0 - Chi*(line3a + line3b)/(2.0d0*r2)
208
209 BigR = (r - gb_sigma*(g**(-0.5d0)) + gb_sigma)/gb_sigma
210 Ri = 1.0d0/BigR
211 Ri2 = Ri*Ri
212 Ri6 = Ri2*Ri2*Ri2
213 Ri7 = Ri6*Ri
214 Ri12 = Ri6*Ri6
215 Ri13 = Ri6*Ri7
216
217 gfact = (g**(-1.5d0))*0.5d0
218
219 dBigRdx = drdx/gb_sigma + dgdx*gfact
220 dBigRdy = drdy/gb_sigma + dgdy*gfact
221 dBigRdz = drdz/gb_sigma + dgdz*gfact
222 dBigRdu1x = dgdu1x*gfact
223 dBigRdu1y = dgdu1y*gfact
224 dBigRdu1z = dgdu1z*gfact
225 dBigRdu2x = dgdu2x*gfact
226 dBigRdu2y = dgdu2y*gfact
227 dBigRdu2z = dgdu2z*gfact
228
229 ! Now, we must do it again for g(ChiPrime) and dgpdx
230
231 line1a = dotsum/opXpdot
232 line2a = dotdiff/omXpdot
233 term1x = -ChiPrime*(line1a*line1bx + line2a*line2bx)/r2
234 term1y = -ChiPrime*(line1a*line1by + line2a*line2by)/r2
235 term1z = -ChiPrime*(line1a*line1bz + line2a*line2bz)/r2
236 line3a = ds2/opXpdot
237 line3b = dd2/omXpdot
238 line3 = ChiPrime*(line3a + line3b)/r4
239 line3x = d(1)*line3
240 line3y = d(2)*line3
241 line3z = d(3)*line3
242
243 dgpdx = term1x + line3x
244 dgpdy = term1y + line3y
245 dgpdz = term1z + line3z
246
247 term1u1x = 2.0d0*(line1a+line2a)*d(1)
248 term1u1y = 2.0d0*(line1a+line2a)*d(2)
249 term1u1z = 2.0d0*(line1a+line2a)*d(3)
250 term1u2x = 2.0d0*(line1a-line2a)*d(1)
251 term1u2y = 2.0d0*(line1a-line2a)*d(2)
252 term1u2z = 2.0d0*(line1a-line2a)*d(3)
253
254 term2a = -line3a/opXpdot
255 term2b = line3b/omXpdot
256
257 term2u1x = ChiPrime*ul2(1)*(term2a + term2b)
258 term2u1y = ChiPrime*ul2(2)*(term2a + term2b)
259 term2u1z = ChiPrime*ul2(3)*(term2a + term2b)
260 term2u2x = ChiPrime*ul1(1)*(term2a + term2b)
261 term2u2y = ChiPrime*ul1(2)*(term2a + term2b)
262 term2u2z = ChiPrime*ul1(3)*(term2a + term2b)
263
264 pref = -ChiPrime*0.5d0/r2
265
266 dgpdu1x = pref*(term1u1x+term2u1x)
267 dgpdu1y = pref*(term1u1y+term2u1y)
268 dgpdu1z = pref*(term1u1z+term2u1z)
269 dgpdu2x = pref*(term1u2x+term2u2x)
270 dgpdu2y = pref*(term1u2y+term2u2y)
271 dgpdu2z = pref*(term1u2z+term2u2z)
272
273 gp = 1.0d0 - ChiPrime*(line3a + line3b)/(2.0d0*r2)
274 gmu = gp**gb_mu
275 gpi = 1.0d0 / gp
276 gmum = gmu*gpi
277
278 ! write(*,*) atom1, atom2, Chi, u1dotu2
279 curlyE = 1.0d0/dsqrt(1.0d0 - Chi*Chi*u1dotu2*u1dotu2)
280
281 dcE = (curlyE**3)*Chi*Chi*u1dotu2
282
283 dcEdu1x = dcE*ul2(1)
284 dcEdu1y = dcE*ul2(2)
285 dcEdu1z = dcE*ul2(3)
286 dcEdu2x = dcE*ul1(1)
287 dcEdu2y = dcE*ul1(2)
288 dcEdu2z = dcE*ul1(3)
289
290 enu = curlyE**gb_nu
291 enum = enu/curlyE
292
293 eps = gb_eps*enu*gmu
294
295 yick1 = gb_eps*enu*gb_mu*gmum
296 yick2 = gb_eps*gmu*gb_nu*enum
297
298 depsdu1x = yick1*dgpdu1x + yick2*dcEdu1x
299 depsdu1y = yick1*dgpdu1y + yick2*dcEdu1y
300 depsdu1z = yick1*dgpdu1z + yick2*dcEdu1z
301 depsdu2x = yick1*dgpdu2x + yick2*dcEdu2x
302 depsdu2y = yick1*dgpdu2y + yick2*dcEdu2y
303 depsdu2z = yick1*dgpdu2z + yick2*dcEdu2z
304
305 R126 = Ri12 - Ri6
306 R137 = 6.0d0*Ri7 - 12.0d0*Ri13
307
308 mess1 = gmu*R137
309 mess2 = R126*gb_mu*gmum
310
311 dUdx = 4.0d0*gb_eps*enu*(mess1*dBigRdx + mess2*dgpdx)*sw
312 dUdy = 4.0d0*gb_eps*enu*(mess1*dBigRdy + mess2*dgpdy)*sw
313 dUdz = 4.0d0*gb_eps*enu*(mess1*dBigRdz + mess2*dgpdz)*sw
314
315 dUdu1x = 4.0d0*(R126*depsdu1x + eps*R137*dBigRdu1x)*sw
316 dUdu1y = 4.0d0*(R126*depsdu1y + eps*R137*dBigRdu1y)*sw
317 dUdu1z = 4.0d0*(R126*depsdu1z + eps*R137*dBigRdu1z)*sw
318 dUdu2x = 4.0d0*(R126*depsdu2x + eps*R137*dBigRdu2x)*sw
319 dUdu2y = 4.0d0*(R126*depsdu2y + eps*R137*dBigRdu2y)*sw
320 dUdu2z = 4.0d0*(R126*depsdu2z + eps*R137*dBigRdu2z)*sw
321
322 #ifdef IS_MPI
323 f_Row(1,atom1) = f_Row(1,atom1) + dUdx
324 f_Row(2,atom1) = f_Row(2,atom1) + dUdy
325 f_Row(3,atom1) = f_Row(3,atom1) + dUdz
326
327 f_Col(1,atom2) = f_Col(1,atom2) - dUdx
328 f_Col(2,atom2) = f_Col(2,atom2) - dUdy
329 f_Col(3,atom2) = f_Col(3,atom2) - dUdz
330
331 t_Row(1,atom1) = t_Row(1,atom1) - ul1(2)*dUdu1z + ul1(3)*dUdu1y
332 t_Row(2,atom1) = t_Row(2,atom1) - ul1(3)*dUdu1x + ul1(1)*dUdu1z
333 t_Row(3,atom1) = t_Row(3,atom1) - ul1(1)*dUdu1y + ul1(2)*dUdu1x
334
335 t_Col(1,atom2) = t_Col(1,atom2) - ul2(2)*dUdu2z + ul2(3)*dUdu2y
336 t_Col(2,atom2) = t_Col(2,atom2) - ul2(3)*dUdu2x + ul2(1)*dUdu2z
337 t_Col(3,atom2) = t_Col(3,atom2) - ul2(1)*dUdu2y + ul2(2)*dUdu2x
338 #else
339 f(1,atom1) = f(1,atom1) + dUdx
340 f(2,atom1) = f(2,atom1) + dUdy
341 f(3,atom1) = f(3,atom1) + dUdz
342
343 f(1,atom2) = f(1,atom2) - dUdx
344 f(2,atom2) = f(2,atom2) - dUdy
345 f(3,atom2) = f(3,atom2) - dUdz
346
347 t(1,atom1) = t(1,atom1) - ul1(2)*dUdu1z + ul1(3)*dUdu1y
348 t(2,atom1) = t(2,atom1) - ul1(3)*dUdu1x + ul1(1)*dUdu1z
349 t(3,atom1) = t(3,atom1) - ul1(1)*dUdu1y + ul1(2)*dUdu1x
350
351 t(1,atom2) = t(1,atom2) - ul2(2)*dUdu2z + ul2(3)*dUdu2y
352 t(2,atom2) = t(2,atom2) - ul2(3)*dUdu2x + ul2(1)*dUdu2z
353 t(3,atom2) = t(3,atom2) - ul2(1)*dUdu2y + ul2(2)*dUdu2x
354 #endif
355
356 if (do_pot) then
357 #ifdef IS_MPI
358 pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*R126*sw
359 pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*R126*sw
360 #else
361 pot = pot + 4.0*eps*R126*sw
362 #endif
363 endif
364
365 vpair = vpair + 4.0*eps*R126
366 #ifdef IS_MPI
367 id1 = AtomRowToGlobal(atom1)
368 id2 = AtomColToGlobal(atom2)
369 #else
370 id1 = atom1
371 id2 = atom2
372 #endif
373
374 if (molMembershipList(id1) .ne. molMembershipList(id2)) then
375
376 fpair(1) = fpair(1) + dUdx
377 fpair(2) = fpair(2) + dUdy
378 fpair(3) = fpair(3) + dUdz
379
380 endif
381
382 return
383 end subroutine do_gb_pair
384
385 end module gb_pair
386
387
388 subroutine set_gb_pair_params(sigma, l2b_ratio, eps, eps_ratio, mu, nu)
389 use definitions, ONLY : dp
390 use gb_pair, ONLY : module_set_gb_pair_params => set_gb_pair_params
391 real( kind = dp ), intent(inout) :: sigma, l2b_ratio, eps, eps_ratio
392 real( kind = dp ), intent(inout) :: mu, nu
393 call module_set_gb_pair_params(sigma, l2b_ratio, eps, eps_ratio, mu, nu)
394 end subroutine set_gb_pair_params