ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/gb.F90
Revision: 1464
Committed: Fri Jul 9 19:29:05 2010 UTC (15 years ago) by gezelter
File size: 14209 byte(s)
Log Message:
removing cruft (atom numbers, do_pot, do_stress) from many modules and
force managers

File Contents

# Content
1 !!
2 !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 !!
4 !! The University of Notre Dame grants you ("Licensee") a
5 !! non-exclusive, royalty free, license to use, modify and
6 !! redistribute this software in source and binary code form, provided
7 !! that the following conditions are met:
8 !!
9 !! 1. Redistributions of source code must retain the above copyright
10 !! notice, this list of conditions and the following disclaimer.
11 !!
12 !! 2. Redistributions in binary form must reproduce the above copyright
13 !! notice, this list of conditions and the following disclaimer in the
14 !! documentation and/or other materials provided with the
15 !! distribution.
16 !!
17 !! This software is provided "AS IS," without a warranty of any
18 !! kind. All express or implied conditions, representations and
19 !! warranties, including any implied warranty of merchantability,
20 !! fitness for a particular purpose or non-infringement, are hereby
21 !! excluded. The University of Notre Dame and its licensors shall not
22 !! be liable for any damages suffered by licensee as a result of
23 !! using, modifying or distributing the software or its
24 !! derivatives. In no event will the University of Notre Dame or its
25 !! licensors be liable for any lost revenue, profit or data, or for
26 !! direct, indirect, special, consequential, incidental or punitive
27 !! damages, however caused and regardless of the theory of liability,
28 !! arising out of the use of or inability to use software, even if the
29 !! University of Notre Dame has been advised of the possibility of
30 !! such damages.
31 !!
32 !! SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 !! research, please cite the appropriate papers when you publish your
34 !! work. Good starting points are:
35 !!
36 !! [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 !! [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 !! [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39 !! [4] Vardeman & Gezelter, in progress (2009).
40 !!
41
42
43 module gayberne
44 use force_globals
45 use definitions
46 use simulation
47 use atype_module
48 use vector_class
49 use linearalgebra
50 use status
51 use lj
52 use fForceOptions
53
54 implicit none
55
56 private
57
58 #define __FORTRAN90
59 #include "UseTheForce/DarkSide/fInteractionMap.h"
60
61 logical, save :: haveGBMap = .false.
62 logical, save :: haveMixingMap = .false.
63 real(kind=dp), save :: mu = 2.0_dp
64 real(kind=dp), save :: nu = 1.0_dp
65
66
67 public :: newGBtype
68 public :: complete_GB_FF
69 public :: do_gb_pair
70 public :: getGayBerneCut
71 public :: destroyGBtypes
72
73 type :: GBtype
74 integer :: atid
75 real(kind = dp ) :: d
76 real(kind = dp ) :: l
77 real(kind = dp ) :: eps
78 real(kind = dp ) :: eps_ratio
79 real(kind = dp ) :: dw
80 logical :: isLJ
81 end type GBtype
82
83 type, private :: GBList
84 integer :: nGBtypes = 0
85 integer :: currentGBtype = 0
86 type(GBtype), pointer :: GBtypes(:) => null()
87 integer, pointer :: atidToGBtype(:) => null()
88 end type GBList
89
90 type(GBList), save :: GBMap
91
92 type :: GBMixParameters
93 real(kind=DP) :: sigma0
94 real(kind=DP) :: eps0
95 real(kind=DP) :: dw
96 real(kind=DP) :: x2
97 real(kind=DP) :: xa2
98 real(kind=DP) :: xai2
99 real(kind=DP) :: xp2
100 real(kind=DP) :: xpap2
101 real(kind=DP) :: xpapi2
102 end type GBMixParameters
103
104 type(GBMixParameters), dimension(:,:), allocatable :: GBMixingMap
105
106 contains
107
108 subroutine newGBtype(c_ident, d, l, eps, eps_ratio, dw, status)
109
110 integer, intent(in) :: c_ident
111 real( kind = dp ), intent(in) :: d, l, eps, eps_ratio, dw
112 integer, intent(out) :: status
113
114 integer :: nGBTypes, nLJTypes, ntypes, myATID
115 integer, pointer :: MatchList(:) => null()
116 integer :: current, i
117 status = 0
118
119 if (.not.associated(GBMap%GBtypes)) then
120
121 call getMatchingElementList(atypes, "is_GayBerne", .true., &
122 nGBtypes, MatchList)
123
124 call getMatchingElementList(atypes, "is_LennardJones", .true., &
125 nLJTypes, MatchList)
126
127 GBMap%nGBtypes = nGBtypes + nLJTypes
128
129 allocate(GBMap%GBtypes(nGBtypes + nLJTypes))
130
131 ntypes = getSize(atypes)
132
133 allocate(GBMap%atidToGBtype(ntypes))
134 endif
135
136 GBMap%currentGBtype = GBMap%currentGBtype + 1
137 current = GBMap%currentGBtype
138
139 myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
140
141 GBMap%atidToGBtype(myATID) = current
142 GBMap%GBtypes(current)%atid = myATID
143 GBMap%GBtypes(current)%d = d
144 GBMap%GBtypes(current)%l = l
145 GBMap%GBtypes(current)%eps = eps
146 GBMap%GBtypes(current)%eps_ratio = eps_ratio
147 GBMap%GBtypes(current)%dw = dw
148 GBMap%GBtypes(current)%isLJ = .false.
149
150 return
151 end subroutine newGBtype
152
153 subroutine complete_GB_FF(status)
154 integer :: status
155 integer :: i, j, l, m, lm, function_type
156 real(kind=dp) :: thisDP, sigma
157 integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
158 logical :: thisProperty
159
160 status = 0
161 if (GBMap%currentGBtype == 0) then
162 call handleError("complete_GB_FF", "No members in GBMap")
163 status = -1
164 return
165 end if
166
167 nAtypes = getSize(atypes)
168
169 if (nAtypes == 0) then
170 status = -1
171 return
172 end if
173
174 ! atypes comes from c side
175 do i = 1, nAtypes
176
177 myATID = getFirstMatchingElement(atypes, 'c_ident', i)
178 call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
179
180 if (thisProperty) then
181 GBMap%currentGBtype = GBMap%currentGBtype + 1
182 current = GBMap%currentGBtype
183
184 GBMap%atidToGBtype(myATID) = current
185 GBMap%GBtypes(current)%atid = myATID
186 GBMap%GBtypes(current)%isLJ = .true.
187 GBMap%GBtypes(current)%d = getSigma(myATID) / sqrt(2.0_dp)
188 GBMap%GBtypes(current)%l = GBMap%GBtypes(current)%d
189 GBMap%GBtypes(current)%eps = getEpsilon(myATID)
190 GBMap%GBtypes(current)%eps_ratio = 1.0_dp
191 GBMap%GBtypes(current)%dw = 1.0_dp
192
193 endif
194
195 end do
196
197 haveGBMap = .true.
198
199
200 end subroutine complete_GB_FF
201
202 subroutine createGBMixingMap()
203 integer :: nGBtypes, i, j
204 real (kind = dp) :: d1, l1, e1, er1, dw1
205 real (kind = dp) :: d2, l2, e2, er2, dw2
206 real (kind = dp) :: er, ermu, xp, ap2
207
208 if (GBMap%currentGBtype == 0) then
209 call handleError("GB", "No members in GBMap")
210 return
211 end if
212
213 nGBtypes = GBMap%nGBtypes
214
215 if (.not. allocated(GBMixingMap)) then
216 allocate(GBMixingMap(nGBtypes, nGBtypes))
217 endif
218
219 do i = 1, nGBtypes
220
221 d1 = GBMap%GBtypes(i)%d
222 l1 = GBMap%GBtypes(i)%l
223 e1 = GBMap%GBtypes(i)%eps
224 er1 = GBMap%GBtypes(i)%eps_ratio
225 dw1 = GBMap%GBtypes(i)%dw
226
227 do j = 1, nGBtypes
228
229 d2 = GBMap%GBtypes(j)%d
230 l2 = GBMap%GBtypes(j)%l
231 e2 = GBMap%GBtypes(j)%eps
232 er2 = GBMap%GBtypes(j)%eps_ratio
233 dw2 = GBMap%GBtypes(j)%dw
234
235 ! Cleaver paper uses sqrt of squares to get sigma0 for
236 ! mixed interactions.
237
238 GBMixingMap(i,j)%sigma0 = sqrt(d1*d1 + d2*d2)
239 GBMixingMap(i,j)%xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2)
240 GBMixingMap(i,j)%xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1)
241 GBMixingMap(i,j)%x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) / &
242 ((l2*l2 + d1*d1) * (l1*l1 + d2*d2))
243
244 ! assumed LB mixing rules for now:
245
246 GBMixingMap(i,j)%dw = 0.5_dp * (dw1 + dw2)
247 GBMixingMap(i,j)%eps0 = sqrt(e1 * e2)
248
249 er = sqrt(er1 * er2)
250 ermu = er**(1.0_dp / mu)
251 xp = (1.0_dp - ermu) / (1.0_dp + ermu)
252 ap2 = 1.0_dp / (1.0_dp + ermu)
253
254 GBMixingMap(i,j)%xp2 = xp*xp
255 GBMixingMap(i,j)%xpap2 = xp*ap2
256 GBMixingMap(i,j)%xpapi2 = xp/ap2
257 enddo
258 enddo
259 haveMixingMap = .true.
260 mu = getGayBerneMu()
261 nu = getGayBerneNu()
262 end subroutine createGBMixingMap
263
264
265 !! gay berne cutoff should be a parameter in globals, this is a temporary
266 !! work around - this should be fixed when gay berne is up and running
267
268 function getGayBerneCut(atomID) result(cutValue)
269 integer, intent(in) :: atomID
270 integer :: gbt1
271 real(kind=dp) :: cutValue, l, d
272
273 if (GBMap%currentGBtype == 0) then
274 call handleError("GB", "No members in GBMap")
275 return
276 end if
277
278 gbt1 = GBMap%atidToGBtype(atomID)
279 l = GBMap%GBtypes(gbt1)%l
280 d = GBMap%GBtypes(gbt1)%d
281
282 ! sigma is actually sqrt(2)*l for prolate ellipsoids
283
284 cutValue = 2.5_dp*sqrt(2.0_dp)*max(l,d)
285
286 end function getGayBerneCut
287
288 subroutine do_gb_pair(atid1, atid2, d, r, r2, sw, vdwMult, vpair, fpair, &
289 pot, A1, A2, f1, t1, t2)
290
291 integer, intent(in) :: atid1, atid2
292 integer :: gbt1, gbt2, id1, id2
293 real (kind=dp), intent(inout) :: r, r2, vdwMult
294 real (kind=dp), dimension(3), intent(in) :: d
295 real (kind=dp), dimension(3), intent(inout) :: fpair
296 real (kind=dp) :: pot, sw, vpair
297 real (kind=dp), dimension(9) :: A1, A2
298 real (kind=dp), dimension(3) :: f1
299 real (kind=dp), dimension(3) :: t1, t2
300 real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
301
302 real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
303 real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
304 real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
305 real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
306 logical :: i_is_lj, j_is_lj
307
308 if (.not.haveMixingMap) then
309 call createGBMixingMap()
310 endif
311
312 gbt1 = GBMap%atidToGBtype(atid1)
313 gbt2 = GBMap%atidToGBtype(atid2)
314
315 i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
316 j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
317
318 sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
319 dw = GBMixingMap(gbt1, gbt2)%dw
320 eps0 = GBMixingMap(gbt1, gbt2)%eps0
321 x2 = GBMixingMap(gbt1, gbt2)%x2
322 xa2 = GBMixingMap(gbt1, gbt2)%xa2
323 xai2 = GBMixingMap(gbt1, gbt2)%xai2
324 xp2 = GBMixingMap(gbt1, gbt2)%xp2
325 xpap2 = GBMixingMap(gbt1, gbt2)%xpap2
326 xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
327
328 ul1(1) = A1(7)
329 ul1(2) = A1(8)
330 ul1(3) = A1(9)
331
332 ul2(1) = A2(7)
333 ul2(2) = A2(8)
334 ul2(3) = A2(9)
335
336 if (i_is_LJ) then
337 a = 0.0_dp
338 ul1 = 0.0_dp
339 else
340 a = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
341 endif
342
343 if (j_is_LJ) then
344 b = 0.0_dp
345 ul2 = 0.0_dp
346 else
347 b = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
348 endif
349
350 if (i_is_LJ.or.j_is_LJ) then
351 g = 0.0_dp
352 else
353 g = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
354 endif
355
356 au = a / r
357 bu = b / r
358
359 au2 = au * au
360 bu2 = bu * bu
361 g2 = g * g
362
363 H = (xa2 * au2 + xai2 * bu2 - 2.0_dp*x2*au*bu*g) / (1.0_dp - x2*g2)
364 Hp = (xpap2*au2 + xpapi2*bu2 - 2.0_dp*xp2*au*bu*g) / (1.0_dp - xp2*g2)
365
366 sigma = sigma0 / sqrt(1.0_dp - H)
367 e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
368 e2 = 1.0_dp - Hp
369 eps = eps0 * (e1**nu) * (e2**mu)
370 BigR = dw*sigma0 / (r - sigma + dw*sigma0)
371
372 R3 = BigR*BigR*BigR
373 R6 = R3*R3
374 R7 = R6 * BigR
375 R12 = R6*R6
376 R13 = R6*R7
377
378 U = vdwMult * 4.0_dp * eps * (R12 - R6)
379
380 s3 = sigma*sigma*sigma
381 s03 = sigma0*sigma0*sigma0
382
383 pref1 = - vdwMult * 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
384
385 pref2 = vdwMult * 8.0_dp * eps * s3 * (6.0_dp*R13 - 3.0_dp*R7) / (dw*r*s03)
386
387 dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 + H))
388
389 dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0_dp - xp2 * g2) &
390 + pref2 * (xa2 * au - x2 *bu*g) / (1.0_dp - x2 * g2)
391
392 dUdb = pref1 * (xpapi2*bu - xp2*au*g) / (1.0_dp - xp2 * g2) &
393 + pref2 * (xai2 * bu - x2 *au*g) / (1.0_dp - x2 * g2)
394
395 dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
396 + 8.0_dp * eps * mu * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) / &
397 (1.0_dp - xp2 * g2) / e2 &
398 + 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
399 (x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
400
401
402 rhat = d / r
403
404 fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
405 fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
406 fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)
407
408 rxu1 = cross_product(d, ul1)
409 rxu2 = cross_product(d, ul2)
410 uxu = cross_product(ul1, ul2)
411
412 !!$ write(*,*) 'pref = ' , pref1, pref2
413 !!$ write(*,*) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3)
414 !!$ write(*,*) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3)
415 !!$ write(*,*) 'uxu = ' , uxu(1), uxu(2), uxu(3)
416 !!$ write(*,*) 'dUda = ', dUda, dudb, dudg, dudr
417 !!$ write(*,*) 'H = ', H,hp,sigma, e1, e2, BigR
418 !!$ write(*,*) 'chi = ', xa2, xai2, x2
419 !!$ write(*,*) 'chip = ', xpap2, xpapi2, xp2
420 !!$ write(*,*) 'eps = ', eps0, e1, e2, eps
421 !!$ write(*,*) 'U =', U, pref1, pref2
422 !!$ write(*,*) 'f =', fx, fy, fz
423 !!$ write(*,*) 'au =', au, bu, g
424 !!$
425
426 pot = pot + U*sw
427
428 f1(1) = f1(1) + fx
429 f1(2) = f1(2) + fy
430 f1(3) = f1(3) + fz
431
432 t1(1) = t1(1) + dUda*rxu1(1) - dUdg*uxu(1)
433 t1(2) = t1(2) + dUda*rxu1(2) - dUdg*uxu(2)
434 t1(3) = t1(3) + dUda*rxu1(3) - dUdg*uxu(3)
435
436 t2(1) = t2(1) + dUdb*rxu2(1) + dUdg*uxu(1)
437 t2(2) = t2(2) + dUdb*rxu2(2) + dUdg*uxu(2)
438 t2(3) = t2(3) + dUdb*rxu2(3) + dUdg*uxu(3)
439
440 vpair = vpair + U*sw
441
442 return
443 end subroutine do_gb_pair
444
445 subroutine destroyGBTypes()
446
447 GBMap%nGBtypes = 0
448 GBMap%currentGBtype = 0
449
450 if (associated(GBMap%GBtypes)) then
451 deallocate(GBMap%GBtypes)
452 GBMap%GBtypes => null()
453 end if
454
455 if (associated(GBMap%atidToGBtype)) then
456 deallocate(GBMap%atidToGBtype)
457 GBMap%atidToGBtype => null()
458 end if
459
460 haveMixingMap = .false.
461
462 end subroutine destroyGBTypes
463
464 end module gayberne
465

Properties

Name Value
svn:keywords Author Id Revision Date