| 50 |
|
!! @author Matthew Meineke |
| 51 |
|
!! @author Christopher Fennell |
| 52 |
|
!! @author J. Daniel Gezelter |
| 53 |
< |
!! @version $Id: sticky.F90,v 1.5 2005-03-12 19:05:16 chrisfen Exp $, $Date: 2005-03-12 19:05:16 $, $Name: not supported by cvs2svn $, $Revision: 1.5 $ |
| 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 |
