ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/UseTheForce/DarkSide/gb.F90
Revision: 1467
Committed: Sat Jul 17 15:33:03 2010 UTC (15 years ago) by gezelter
File size: 14215 byte(s)
Log Message:
well, it compiles, but still segfaults

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 fForceOptions
52
53 implicit none
54
55 private
56 real(kind=dp), external :: getSigma
57 real(kind=dp), external :: getEpsilon
58
59 #define __FORTRAN90
60 #include "UseTheForce/DarkSide/fInteractionMap.h"
61
62 logical, save :: haveGBMap = .false.
63 logical, save :: haveMixingMap = .false.
64 real(kind=dp), save :: mu = 2.0_dp
65 real(kind=dp), save :: nu = 1.0_dp
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, &
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) :: pot, sw, vpair
296 real (kind=dp), dimension(9) :: A1, A2
297 real (kind=dp), dimension(3) :: f1
298 real (kind=dp), dimension(3) :: t1, t2
299 real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
300
301 real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
302 real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
303 real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
304 real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
305 logical :: i_is_lj, j_is_lj
306
307 if (.not.haveMixingMap) then
308 call createGBMixingMap()
309 endif
310
311 gbt1 = GBMap%atidToGBtype(atid1)
312 gbt2 = GBMap%atidToGBtype(atid2)
313
314 i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
315 j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
316
317 sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
318 dw = GBMixingMap(gbt1, gbt2)%dw
319 eps0 = GBMixingMap(gbt1, gbt2)%eps0
320 x2 = GBMixingMap(gbt1, gbt2)%x2
321 xa2 = GBMixingMap(gbt1, gbt2)%xa2
322 xai2 = GBMixingMap(gbt1, gbt2)%xai2
323 xp2 = GBMixingMap(gbt1, gbt2)%xp2
324 xpap2 = GBMixingMap(gbt1, gbt2)%xpap2
325 xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
326
327 ul1(1) = A1(7)
328 ul1(2) = A1(8)
329 ul1(3) = A1(9)
330
331 ul2(1) = A2(7)
332 ul2(2) = A2(8)
333 ul2(3) = A2(9)
334
335 if (i_is_LJ) then
336 a = 0.0_dp
337 ul1 = 0.0_dp
338 else
339 a = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
340 endif
341
342 if (j_is_LJ) then
343 b = 0.0_dp
344 ul2 = 0.0_dp
345 else
346 b = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
347 endif
348
349 if (i_is_LJ.or.j_is_LJ) then
350 g = 0.0_dp
351 else
352 g = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
353 endif
354
355 au = a / r
356 bu = b / r
357
358 au2 = au * au
359 bu2 = bu * bu
360 g2 = g * g
361
362 H = (xa2 * au2 + xai2 * bu2 - 2.0_dp*x2*au*bu*g) / (1.0_dp - x2*g2)
363 Hp = (xpap2*au2 + xpapi2*bu2 - 2.0_dp*xp2*au*bu*g) / (1.0_dp - xp2*g2)
364
365 sigma = sigma0 / sqrt(1.0_dp - H)
366 e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
367 e2 = 1.0_dp - Hp
368 eps = eps0 * (e1**nu) * (e2**mu)
369 BigR = dw*sigma0 / (r - sigma + dw*sigma0)
370
371 R3 = BigR*BigR*BigR
372 R6 = R3*R3
373 R7 = R6 * BigR
374 R12 = R6*R6
375 R13 = R6*R7
376
377 U = vdwMult * 4.0_dp * eps * (R12 - R6)
378
379 s3 = sigma*sigma*sigma
380 s03 = sigma0*sigma0*sigma0
381
382 pref1 = - vdwMult * 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
383
384 pref2 = vdwMult * 8.0_dp * eps * s3 * (6.0_dp*R13 - 3.0_dp*R7) / (dw*r*s03)
385
386 dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 + H))
387
388 dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0_dp - xp2 * g2) &
389 + pref2 * (xa2 * au - x2 *bu*g) / (1.0_dp - x2 * g2)
390
391 dUdb = pref1 * (xpapi2*bu - xp2*au*g) / (1.0_dp - xp2 * g2) &
392 + pref2 * (xai2 * bu - x2 *au*g) / (1.0_dp - x2 * g2)
393
394 dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
395 + 8.0_dp * eps * mu * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) / &
396 (1.0_dp - xp2 * g2) / e2 &
397 + 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
398 (x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
399
400
401 rhat = d / r
402
403 fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
404 fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
405 fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)
406
407 rxu1 = cross_product(d, ul1)
408 rxu2 = cross_product(d, ul2)
409 uxu = cross_product(ul1, ul2)
410
411 !!$ write(*,*) 'pref = ' , pref1, pref2
412 !!$ write(*,*) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3)
413 !!$ write(*,*) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3)
414 !!$ write(*,*) 'uxu = ' , uxu(1), uxu(2), uxu(3)
415 !!$ write(*,*) 'dUda = ', dUda, dudb, dudg, dudr
416 !!$ write(*,*) 'H = ', H,hp,sigma, e1, e2, BigR
417 !!$ write(*,*) 'chi = ', xa2, xai2, x2
418 !!$ write(*,*) 'chip = ', xpap2, xpapi2, xp2
419 !!$ write(*,*) 'eps = ', eps0, e1, e2, eps
420 !!$ write(*,*) 'U =', U, pref1, pref2
421 !!$ write(*,*) 'f =', fx, fy, fz
422 !!$ write(*,*) 'au =', au, bu, g
423 !!$
424
425 pot = pot + U*sw
426
427 f1(1) = f1(1) + fx
428 f1(2) = f1(2) + fy
429 f1(3) = f1(3) + fz
430
431 t1(1) = t1(1) + dUda*rxu1(1) - dUdg*uxu(1)
432 t1(2) = t1(2) + dUda*rxu1(2) - dUdg*uxu(2)
433 t1(3) = t1(3) + dUda*rxu1(3) - dUdg*uxu(3)
434
435 t2(1) = t2(1) + dUdb*rxu2(1) + dUdg*uxu(1)
436 t2(2) = t2(2) + dUdb*rxu2(2) + dUdg*uxu(2)
437 t2(3) = t2(3) + dUdb*rxu2(3) + dUdg*uxu(3)
438
439 vpair = vpair + U*sw
440
441 return
442 end subroutine do_gb_pair
443
444 subroutine destroyGBTypes()
445
446 GBMap%nGBtypes = 0
447 GBMap%currentGBtype = 0
448
449 if (associated(GBMap%GBtypes)) then
450 deallocate(GBMap%GBtypes)
451 GBMap%GBtypes => null()
452 end if
453
454 if (associated(GBMap%atidToGBtype)) then
455 deallocate(GBMap%atidToGBtype)
456 GBMap%atidToGBtype => null()
457 end if
458
459 haveMixingMap = .false.
460
461 end subroutine destroyGBTypes
462
463 end module gayberne
464

Properties

Name Value
svn:keywords Author Id Revision Date