ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/gb.F90
Revision: 1386
Committed: Fri Oct 23 18:41:09 2009 UTC (15 years, 10 months ago) by gezelter
File size: 14262 byte(s)
Log Message:
removing MPI responsibilities from the lowest level force routines.  This is
in preparation for migrating these routines (LJ, electrostatic, eam, suttonchen,
gay-berne, sticky) to C++

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. Acknowledgement of the program authors must be made in any
10 !! publication of scientific results based in part on use of the
11 !! program. An acceptable form of acknowledgement is citation of
12 !! the article in which the program was described (Matthew
13 !! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 !! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 !! Parallel Simulation Engine for Molecular Dynamics,"
16 !! J. Comput. Chem. 26, pp. 252-271 (2005))
17 !!
18 !! 2. Redistributions of source code must retain the above copyright
19 !! notice, this list of conditions and the following disclaimer.
20 !!
21 !! 3. Redistributions in binary form must reproduce the above copyright
22 !! notice, this list of conditions and the following disclaimer in the
23 !! documentation and/or other materials provided with the
24 !! distribution.
25 !!
26 !! This software is provided "AS IS," without a warranty of any
27 !! kind. All express or implied conditions, representations and
28 !! warranties, including any implied warranty of merchantability,
29 !! fitness for a particular purpose or non-infringement, are hereby
30 !! excluded. The University of Notre Dame and its licensors shall not
31 !! be liable for any damages suffered by licensee as a result of
32 !! using, modifying or distributing the software or its
33 !! derivatives. In no event will the University of Notre Dame or its
34 !! licensors be liable for any lost revenue, profit or data, or for
35 !! direct, indirect, special, consequential, incidental or punitive
36 !! damages, however caused and regardless of the theory of liability,
37 !! arising out of the use of or inability to use software, even if the
38 !! University of Notre Dame has been advised of the possibility of
39 !! such damages.
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(atom1, atom2, atid1, atid2, d, r, r2, sw, vdwMult, vpair, fpair, &
289 pot, A1, A2, f1, t1, t2, do_pot)
290
291 integer, intent(in) :: atom1, atom2, 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 logical, intent(in) :: do_pot
301 real (kind = dp), dimension(3) :: ul1, ul2, rxu1, rxu2, uxu, rhat
302
303 real (kind = dp) :: sigma0, dw, eps0, x2, xa2, xai2, xp2, xpap2, xpapi2
304 real (kind = dp) :: e1, e2, eps, sigma, s3, s03, au2, bu2, au, bu, a, b, g, g2
305 real (kind = dp) :: U, BigR, R3, R6, R7, R12, R13, H, Hp, fx, fy, fz
306 real (kind = dp) :: dUdr, dUda, dUdb, dUdg, pref1, pref2
307 logical :: i_is_lj, j_is_lj
308
309 if (.not.haveMixingMap) then
310 call createGBMixingMap()
311 endif
312
313 gbt1 = GBMap%atidToGBtype(atid1)
314 gbt2 = GBMap%atidToGBtype(atid2)
315
316 i_is_LJ = GBMap%GBTypes(gbt1)%isLJ
317 j_is_LJ = GBMap%GBTypes(gbt2)%isLJ
318
319 sigma0 = GBMixingMap(gbt1, gbt2)%sigma0
320 dw = GBMixingMap(gbt1, gbt2)%dw
321 eps0 = GBMixingMap(gbt1, gbt2)%eps0
322 x2 = GBMixingMap(gbt1, gbt2)%x2
323 xa2 = GBMixingMap(gbt1, gbt2)%xa2
324 xai2 = GBMixingMap(gbt1, gbt2)%xai2
325 xp2 = GBMixingMap(gbt1, gbt2)%xp2
326 xpap2 = GBMixingMap(gbt1, gbt2)%xpap2
327 xpapi2 = GBMixingMap(gbt1, gbt2)%xpapi2
328
329 ul1(1) = A1(7)
330 ul1(2) = A1(8)
331 ul1(3) = A1(9)
332
333 ul2(1) = A2(7)
334 ul2(2) = A2(8)
335 ul2(3) = A2(9)
336
337 if (i_is_LJ) then
338 a = 0.0_dp
339 ul1 = 0.0_dp
340 else
341 a = d(1)*ul1(1) + d(2)*ul1(2) + d(3)*ul1(3)
342 endif
343
344 if (j_is_LJ) then
345 b = 0.0_dp
346 ul2 = 0.0_dp
347 else
348 b = d(1)*ul2(1) + d(2)*ul2(2) + d(3)*ul2(3)
349 endif
350
351 if (i_is_LJ.or.j_is_LJ) then
352 g = 0.0_dp
353 else
354 g = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
355 endif
356
357 au = a / r
358 bu = b / r
359
360 au2 = au * au
361 bu2 = bu * bu
362 g2 = g * g
363
364 H = (xa2 * au2 + xai2 * bu2 - 2.0_dp*x2*au*bu*g) / (1.0_dp - x2*g2)
365 Hp = (xpap2*au2 + xpapi2*bu2 - 2.0_dp*xp2*au*bu*g) / (1.0_dp - xp2*g2)
366
367 sigma = sigma0 / sqrt(1.0_dp - H)
368 e1 = 1.0_dp / sqrt(1.0_dp - x2*g2)
369 e2 = 1.0_dp - Hp
370 eps = eps0 * (e1**nu) * (e2**mu)
371 BigR = dw*sigma0 / (r - sigma + dw*sigma0)
372
373 R3 = BigR*BigR*BigR
374 R6 = R3*R3
375 R7 = R6 * BigR
376 R12 = R6*R6
377 R13 = R6*R7
378
379 U = vdwMult * 4.0_dp * eps * (R12 - R6)
380
381 s3 = sigma*sigma*sigma
382 s03 = sigma0*sigma0*sigma0
383
384 pref1 = - vdwMult * 8.0_dp * eps * mu * (R12 - R6) / (e2 * r)
385
386 pref2 = vdwMult * 8.0_dp * eps * s3 * (6.0_dp*R13 - 3.0_dp*R7) / (dw*r*s03)
387
388 dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 + H))
389
390 dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0_dp - xp2 * g2) &
391 + pref2 * (xa2 * au - x2 *bu*g) / (1.0_dp - x2 * g2)
392
393 dUdb = pref1 * (xpapi2*bu - xp2*au*g) / (1.0_dp - xp2 * g2) &
394 + pref2 * (xai2 * bu - x2 *au*g) / (1.0_dp - x2 * g2)
395
396 dUdg = 4.0_dp * eps * nu * (R12 - R6) * x2 * g / (1.0_dp - x2*g2) &
397 + 8.0_dp * eps * mu * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) / &
398 (1.0_dp - xp2 * g2) / e2 &
399 + 8.0_dp * eps * s3 * (3.0_dp * R7 - 6.0_dp * R13) * &
400 (x2 * au * bu - H * x2 * g) / (1.0_dp - x2 * g2) / (dw * s03)
401
402
403 rhat = d / r
404
405 fx = dUdr * rhat(1) + dUda * ul1(1) + dUdb * ul2(1)
406 fy = dUdr * rhat(2) + dUda * ul1(2) + dUdb * ul2(2)
407 fz = dUdr * rhat(3) + dUda * ul1(3) + dUdb * ul2(3)
408
409 rxu1 = cross_product(d, ul1)
410 rxu2 = cross_product(d, ul2)
411 uxu = cross_product(ul1, ul2)
412
413 !!$ write(*,*) 'pref = ' , pref1, pref2
414 !!$ write(*,*) 'rxu1 = ' , rxu1(1), rxu1(2), rxu1(3)
415 !!$ write(*,*) 'rxu2 = ' , rxu2(1), rxu2(2), rxu2(3)
416 !!$ write(*,*) 'uxu = ' , uxu(1), uxu(2), uxu(3)
417 !!$ write(*,*) 'dUda = ', dUda, dudb, dudg, dudr
418 !!$ write(*,*) 'H = ', H,hp,sigma, e1, e2, BigR
419 !!$ write(*,*) 'chi = ', xa2, xai2, x2
420 !!$ write(*,*) 'chip = ', xpap2, xpapi2, xp2
421 !!$ write(*,*) 'eps = ', eps0, e1, e2, eps
422 !!$ write(*,*) 'U =', U, pref1, pref2
423 !!$ write(*,*) 'f =', fx, fy, fz
424 !!$ write(*,*) 'au =', au, bu, g
425 !!$
426
427 pot = pot + U*sw
428
429 f1(1) = f1(1) + fx
430 f1(2) = f1(2) + fy
431 f1(3) = f1(3) + fz
432
433 t1(1) = t1(1) + dUda*rxu1(1) - dUdg*uxu(1)
434 t1(2) = t1(2) + dUda*rxu1(2) - dUdg*uxu(2)
435 t1(3) = t1(3) + dUda*rxu1(3) - dUdg*uxu(3)
436
437 t2(1) = t2(1) + dUdb*rxu2(1) + dUdg*uxu(1)
438 t2(2) = t2(2) + dUdb*rxu2(2) + dUdg*uxu(2)
439 t2(3) = t2(3) + dUdb*rxu2(3) + dUdg*uxu(3)
440
441 vpair = vpair + U*sw
442
443 return
444 end subroutine do_gb_pair
445
446 subroutine destroyGBTypes()
447
448 GBMap%nGBtypes = 0
449 GBMap%currentGBtype = 0
450
451 if (associated(GBMap%GBtypes)) then
452 deallocate(GBMap%GBtypes)
453 GBMap%GBtypes => null()
454 end if
455
456 if (associated(GBMap%atidToGBtype)) then
457 deallocate(GBMap%atidToGBtype)
458 GBMap%atidToGBtype => null()
459 end if
460
461 haveMixingMap = .false.
462
463 end subroutine destroyGBTypes
464
465 end module gayberne
466