50 |
|
!! @author Matthew Meineke |
51 |
|
!! @author Christopher Fennell |
52 |
|
!! @author J. Daniel Gezelter |
53 |
< |
!! @version $Id: sticky.F90,v 1.15 2005-10-12 18:59:16 chuckv Exp $, $Date: 2005-10-12 18:59:16 $, $Name: not supported by cvs2svn $, $Revision: 1.15 $ |
53 |
> |
!! @version $Id: sticky.F90,v 1.19 2006-04-20 21:02:00 chrisfen Exp $, $Date: 2006-04-20 21:02:00 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $ |
54 |
|
|
55 |
|
module sticky |
56 |
|
|
60 |
|
use vector_class |
61 |
|
use simulation |
62 |
|
use status |
63 |
+ |
use interpolation |
64 |
|
#ifdef IS_MPI |
65 |
|
use mpiSimulation |
66 |
|
#endif |
87 |
|
real( kind = dp ) :: rlp = 0.0_dp |
88 |
|
real( kind = dp ) :: rup = 0.0_dp |
89 |
|
real( kind = dp ) :: rbig = 0.0_dp |
90 |
+ |
type(cubicSpline) :: stickySpline |
91 |
+ |
type(cubicSpline) :: stickySplineP |
92 |
|
end type StickyList |
93 |
|
|
94 |
|
type(StickyList), dimension(:),allocatable :: StickyMap |
95 |
+ |
logical, save :: hasStickyMap = .false. |
96 |
|
|
97 |
|
contains |
98 |
|
|
103 |
|
real( kind = dp ), intent(in) :: w0, v0, v0p |
104 |
|
real( kind = dp ), intent(in) :: rl, ru |
105 |
|
real( kind = dp ), intent(in) :: rlp, rup |
106 |
+ |
real( kind = dp ), dimension(2) :: rCubVals, sCubVals, rpCubVals, spCubVals |
107 |
|
integer :: nATypes, myATID |
108 |
|
|
109 |
|
|
152 |
|
else |
153 |
|
StickyMap(myATID)%rbig = StickyMap(myATID)%rup |
154 |
|
endif |
155 |
+ |
|
156 |
+ |
! build the 2 cubic splines for the sticky switching functions |
157 |
+ |
|
158 |
+ |
rCubVals(1) = rl |
159 |
+ |
rCubVals(2) = ru |
160 |
+ |
sCubVals(1) = 1.0d0 |
161 |
+ |
sCubVals(2) = 0.0d0 |
162 |
+ |
call newSpline(StickyMap(myATID)%stickySpline, rCubVals, sCubVals, .true.) |
163 |
+ |
rpCubVals(1) = rlp |
164 |
+ |
rpCubVals(2) = rup |
165 |
+ |
spCubVals(1) = 1.0d0 |
166 |
+ |
spCubVals(2) = 0.0d0 |
167 |
+ |
call newSpline(StickyMap(myATID)%stickySplineP,rpCubVals,spCubVals,.true.) |
168 |
|
|
169 |
+ |
hasStickyMap = .true. |
170 |
+ |
|
171 |
|
return |
172 |
|
end subroutine newStickyType |
173 |
|
|
225 |
|
real (kind=dp) :: radcomxj, radcomyj, radcomzj |
226 |
|
integer :: id1, id2 |
227 |
|
integer :: me1, me2 |
228 |
< |
real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig |
228 |
> |
real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig, dx |
229 |
|
|
210 |
– |
if (.not.allocated(StickyMap)) then |
211 |
– |
call handleError("sticky", "no StickyMap was present before first call of do_sticky_pair!") |
212 |
– |
return |
213 |
– |
end if |
214 |
– |
|
230 |
|
#ifdef IS_MPI |
231 |
|
me1 = atid_Row(atom1) |
232 |
|
me2 = atid_Col(atom2) |
304 |
|
yj2 = yj*yj |
305 |
|
zj2 = zj*zj |
306 |
|
|
292 |
– |
call calc_sw_fnc(rij, rl, ru, rlp, rup, s, sp, dsdr, dspdr) |
307 |
|
|
308 |
+ |
! calculate the switching info. from the splines |
309 |
+ |
if (me1.eq.me2) then |
310 |
+ |
s = 0.0d0 |
311 |
+ |
dsdr = 0.0d0 |
312 |
+ |
sp = 0.0d0 |
313 |
+ |
dspdr = 0.0d0 |
314 |
+ |
|
315 |
+ |
if (rij.lt.ru) then |
316 |
+ |
if (rij.lt.rl) then |
317 |
+ |
s = 1.0d0 |
318 |
+ |
dsdr = 0.0d0 |
319 |
+ |
else |
320 |
+ |
! we are in the switching region |
321 |
+ |
dx = rij - rl |
322 |
+ |
s = StickyMap(me1)%stickySpline%y(1) + & |
323 |
+ |
dx*(dx*(StickyMap(me1)%stickySpline%c(1) + & |
324 |
+ |
dx*StickyMap(me1)%stickySpline%d(1))) |
325 |
+ |
dsdr = dx*(2.0d0 * StickyMap(me1)%stickySpline%c(1) + & |
326 |
+ |
3.0d0 * dx * StickyMap(me1)%stickySpline%d(1)) |
327 |
+ |
endif |
328 |
+ |
endif |
329 |
+ |
if (rij.lt.rup) then |
330 |
+ |
if (rij.lt.rlp) then |
331 |
+ |
sp = 1.0d0 |
332 |
+ |
dspdr = 0.0d0 |
333 |
+ |
else |
334 |
+ |
! we are in the switching region |
335 |
+ |
dx = rij - rlp |
336 |
+ |
sp = StickyMap(me1)%stickySplineP%y(1) + & |
337 |
+ |
dx*(dx*(StickyMap(me1)%stickySplineP%c(1) + & |
338 |
+ |
dx*StickyMap(me1)%stickySplineP%d(1))) |
339 |
+ |
dspdr = dx*(2.0d0 * StickyMap(me1)%stickySplineP%c(1) + & |
340 |
+ |
3.0d0 * dx * StickyMap(me1)%stickySplineP%d(1)) |
341 |
+ |
endif |
342 |
+ |
endif |
343 |
+ |
else |
344 |
+ |
! calculate the switching function explicitly rather than from |
345 |
+ |
! the splines with mixed sticky maps |
346 |
+ |
call calc_sw_fnc(rij, rl, ru, rlp, rup, s, sp, dsdr, dspdr) |
347 |
+ |
endif |
348 |
+ |
|
349 |
|
wi = 2.0d0*(xi2-yi2)*zi / r3 |
350 |
|
wj = 2.0d0*(xj2-yj2)*zj / r3 |
351 |
|
w = wi+wj |
363 |
|
vpair = vpair + 0.5d0*(v0*s*w + v0p*sp*wp) |
364 |
|
if (do_pot) then |
365 |
|
#ifdef IS_MPI |
366 |
< |
pot_row(STICKY_POT,atom1) = pot_row(STICKY_POT,atom1) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw |
367 |
< |
pot_col(STICKY_POT,atom2) = pot_col(STICKY_POT,atom2) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw |
366 |
> |
pot_row(HB_POT,atom1) = pot_row(HB_POT,atom1) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw |
367 |
> |
pot_col(HB_POT,atom2) = pot_col(HB_POT,atom2) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw |
368 |
|
#else |
369 |
|
pot = pot + 0.5d0*(v0*s*w + v0p*sp*wp)*sw |
370 |
|
#endif |
552 |
|
real (kind=dp), intent(inout) :: s, sp, dsdr, dspdr |
553 |
|
|
554 |
|
! distances must be in angstroms |
555 |
< |
|
556 |
< |
if (r.lt.rl) then |
557 |
< |
s = 1.0d0 |
558 |
< |
dsdr = 0.0d0 |
559 |
< |
elseif (r.gt.ru) then |
560 |
< |
s = 0.0d0 |
561 |
< |
dsdr = 0.0d0 |
562 |
< |
else |
563 |
< |
s = ((ru + 2.0d0*r - 3.0d0*rl) * (ru-r)**2) / & |
564 |
< |
((ru - rl)**3) |
565 |
< |
dsdr = 6.0d0*(r-ru)*(r-rl)/((ru - rl)**3) |
555 |
> |
s = 0.0d0 |
556 |
> |
dsdr = 0.0d0 |
557 |
> |
sp = 0.0d0 |
558 |
> |
dspdr = 0.0d0 |
559 |
> |
|
560 |
> |
if (r.lt.ru) then |
561 |
> |
if (r.lt.rl) then |
562 |
> |
s = 1.0d0 |
563 |
> |
dsdr = 0.0d0 |
564 |
> |
else |
565 |
> |
s = ((ru + 2.0d0*r - 3.0d0*rl) * (ru-r)**2) / & |
566 |
> |
((ru - rl)**3) |
567 |
> |
dsdr = 6.0d0*(r-ru)*(r-rl)/((ru - rl)**3) |
568 |
> |
endif |
569 |
|
endif |
570 |
|
|
571 |
< |
if (r.lt.rlp) then |
572 |
< |
sp = 1.0d0 |
573 |
< |
dspdr = 0.0d0 |
574 |
< |
elseif (r.gt.rup) then |
575 |
< |
sp = 0.0d0 |
576 |
< |
dspdr = 0.0d0 |
577 |
< |
else |
578 |
< |
sp = ((rup + 2.0d0*r - 3.0d0*rlp) * (rup-r)**2) / & |
579 |
< |
((rup - rlp)**3) |
522 |
< |
dspdr = 6.0d0*(r-rup)*(r-rlp)/((rup - rlp)**3) |
571 |
> |
if (r.lt.rup) then |
572 |
> |
if (r.lt.rlp) then |
573 |
> |
sp = 1.0d0 |
574 |
> |
dspdr = 0.0d0 |
575 |
> |
else |
576 |
> |
sp = ((rup + 2.0d0*r - 3.0d0*rlp) * (rup-r)**2) / & |
577 |
> |
((rup - rlp)**3) |
578 |
> |
dspdr = 6.0d0*(r-rup)*(r-rlp)/((rup - rlp)**3) |
579 |
> |
endif |
580 |
|
endif |
581 |
|
|
582 |
|
return |
739 |
|
|
740 |
|
if (do_pot) then |
741 |
|
#ifdef IS_MPI |
742 |
< |
pot_row(atom1) = pot_row(atom1) + 0.25d0*(v0*s*w)*sw |
743 |
< |
pot_col(atom2) = pot_col(atom2) + 0.25d0*(v0*s*w)*sw |
742 |
> |
pot_row(HB_POT,atom1) = pot_row(HB_POT,atom1) + 0.25d0*(v0*s*w)*sw |
743 |
> |
pot_col(HB_POT,atom2) = pot_col(HB_POT,atom2) + 0.25d0*(v0*s*w)*sw |
744 |
|
#else |
745 |
|
pot = pot + 0.5d0*(v0*s*w)*sw |
746 |
|
#endif |