48 |
|
!! Corresponds to the force field defined in ssd_FF.cpp |
49 |
|
!! @author Charles F. Vardeman II |
50 |
|
!! @author Matthew Meineke |
51 |
< |
!! @author Christopher Fennel |
51 |
> |
!! @author Christopher Fennell |
52 |
|
!! @author J. Daniel Gezelter |
53 |
< |
!! @version $Id: sticky.F90,v 1.4 2005-01-14 20:31:16 gezelter Exp $, $Date: 2005-01-14 20:31:16 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $ |
53 |
> |
!! @version $Id: sticky.F90,v 1.7 2005-04-15 22:03:49 gezelter Exp $, $Date: 2005-04-15 22:03:49 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $ |
54 |
|
|
55 |
|
module sticky |
56 |
|
|
69 |
|
|
70 |
|
public :: newStickyType |
71 |
|
public :: do_sticky_pair |
72 |
+ |
public :: destroyStickyTypes |
73 |
|
|
74 |
|
|
75 |
|
type :: StickyList |
83 |
|
real( kind = dp ) :: rup = 0.0_dp |
84 |
|
real( kind = dp ) :: rbig = 0.0_dp |
85 |
|
end type StickyList |
86 |
< |
|
86 |
> |
|
87 |
|
type(StickyList), dimension(:),allocatable :: StickyMap |
88 |
|
|
89 |
|
contains |
97 |
|
real( kind = dp ), intent(in) :: rlp, rup |
98 |
|
integer :: nATypes, myATID |
99 |
|
|
100 |
< |
|
100 |
> |
|
101 |
|
isError = 0 |
102 |
|
myATID = getFirstMatchingElement(atypes, "c_ident", c_ident) |
103 |
< |
|
103 |
> |
|
104 |
|
!! Be simple-minded and assume that we need a StickyMap that |
105 |
|
!! is the same size as the total number of atom types |
106 |
|
|
129 |
|
StickyMap(myATID)%c_ident = c_ident |
130 |
|
|
131 |
|
! we could pass all 5 parameters if we felt like it... |
132 |
< |
|
132 |
> |
|
133 |
|
StickyMap(myATID)%w0 = w0 |
134 |
|
StickyMap(myATID)%v0 = v0 |
135 |
|
StickyMap(myATID)%v0p = v0p |
143 |
|
else |
144 |
|
StickyMap(myATID)%rbig = StickyMap(myATID)%rup |
145 |
|
endif |
146 |
< |
|
146 |
> |
|
147 |
|
return |
148 |
|
end subroutine newStickyType |
149 |
|
|
150 |
|
subroutine do_sticky_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, & |
151 |
|
pot, A, f, t, do_pot) |
152 |
< |
|
152 |
> |
|
153 |
|
!! This routine does only the sticky portion of the SSD potential |
154 |
|
!! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)]. |
155 |
|
!! The Lennard-Jones and dipolar interaction must be handled separately. |
156 |
< |
|
156 |
> |
|
157 |
|
!! We assume that the rotation matrices have already been calculated |
158 |
|
!! and placed in the A array. |
159 |
|
|
187 |
|
real (kind=dp) :: radcomxj, radcomyj, radcomzj |
188 |
|
integer :: id1, id2 |
189 |
|
integer :: me1, me2 |
190 |
< |
real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig |
190 |
> |
real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig |
191 |
|
|
192 |
< |
if (.not.allocated(StickyMap)) then |
192 |
> |
if (.not.allocated(StickyMap)) then |
193 |
|
call handleError("sticky", "no StickyMap was present before first call of do_sticky_pair!") |
194 |
|
return |
195 |
|
end if |
196 |
< |
|
196 |
> |
|
197 |
|
#ifdef IS_MPI |
198 |
|
me1 = atid_Row(atom1) |
199 |
|
me2 = atid_Col(atom2) |
461 |
|
id1 = atom1 |
462 |
|
id2 = atom2 |
463 |
|
#endif |
464 |
< |
|
464 |
> |
|
465 |
|
if (molMembershipList(id1) .ne. molMembershipList(id2)) then |
466 |
< |
|
466 |
> |
|
467 |
|
fpair(1) = fpair(1) + fxradial |
468 |
|
fpair(2) = fpair(2) + fyradial |
469 |
|
fpair(3) = fpair(3) + fzradial |
470 |
< |
|
470 |
> |
|
471 |
|
endif |
472 |
|
endif |
473 |
|
end subroutine do_sticky_pair |
474 |
|
|
475 |
|
!! calculates the switching functions and their derivatives for a given |
476 |
|
subroutine calc_sw_fnc(r, rl, ru, rlp, rup, s, sp, dsdr, dspdr) |
477 |
< |
|
477 |
> |
|
478 |
|
real (kind=dp), intent(in) :: r, rl, ru, rlp, rup |
479 |
|
real (kind=dp), intent(inout) :: s, sp, dsdr, dspdr |
480 |
< |
|
480 |
> |
|
481 |
|
! distances must be in angstroms |
482 |
< |
|
482 |
> |
|
483 |
|
if (r.lt.rl) then |
484 |
|
s = 1.0d0 |
485 |
|
dsdr = 0.0d0 |
503 |
|
((rup - rlp)**3) |
504 |
|
dspdr = 6.0d0*(r-rup)*(r-rlp)/((rup - rlp)**3) |
505 |
|
endif |
506 |
< |
|
506 |
> |
|
507 |
|
return |
508 |
|
end subroutine calc_sw_fnc |
509 |
+ |
|
510 |
+ |
subroutine destroyStickyTypes() |
511 |
+ |
if(allocated(StickyMap)) deallocate(StickyMap) |
512 |
+ |
end subroutine destroyStickyTypes |
513 |
|
end module sticky |