50 |
|
!! @author Matthew Meineke |
51 |
|
!! @author Christopher Fennell |
52 |
|
!! @author J. Daniel Gezelter |
53 |
< |
!! @version $Id: sticky.F90,v 1.8 2005-05-05 14:47:35 chrisfen Exp $, $Date: 2005-05-05 14:47:35 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $ |
53 |
> |
!! @version $Id: sticky.F90,v 1.9 2005-05-12 19:43:48 chrisfen Exp $, $Date: 2005-05-12 19:43:48 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $ |
54 |
|
|
55 |
|
module sticky |
56 |
|
|
530 |
|
logical, intent(in) :: do_pot |
531 |
|
|
532 |
|
real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2 |
533 |
< |
real (kind=dp) :: r3, r5, r6, s, sp, dsdr, dspdr |
534 |
< |
real (kind=dp) :: wi, wj, w, wip, wjp, wp, wi2, wj2 |
533 |
> |
real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat |
534 |
> |
real (kind=dp) :: rI, rI2, rI3, rI4, rI5, rI6, rI7, s, sp, dsdr, dspdr |
535 |
> |
real (kind=dp) :: wi, wj, w, wip, wjp, wp, wi2, wj2, wip3, wjp3 |
536 |
|
real (kind=dp) :: dwidx, dwidy, dwidz, dwjdx, dwjdy, dwjdz |
537 |
|
real (kind=dp) :: dwipdx, dwipdy, dwipdz, dwjpdx, dwjpdy, dwjpdz |
538 |
|
real (kind=dp) :: dwidux, dwiduy, dwiduz, dwjdux, dwjduy, dwjduz |
549 |
|
integer :: id1, id2 |
550 |
|
integer :: me1, me2 |
551 |
|
real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig |
552 |
< |
|
552 |
> |
real (kind=dp) :: zi3, zi4, zi5, zj3, zj4, zj5 |
553 |
> |
real (kind=dp) :: oSw1, oSw2, prodVal |
554 |
> |
real (kind=dp) :: prei1, prei2, prei, prej1, prej2, prej |
555 |
> |
real (kind=dp) :: walt, walti, waltj, dwaltidx, dwaltidy, dwaltidz |
556 |
> |
real (kind=dp) :: dwaltjdx, dwaltjdy, dwaltjdz |
557 |
> |
real (kind=dp) :: dwaltidux, dwaltiduy, dwaltiduz |
558 |
> |
real (kind=dp) :: dwaltjdux, dwaltjduy, dwaltjduz |
559 |
> |
real (kind=dp) :: doSw1idx, doSw1idy, doSw1idz, doSw1jdx, doSw1jdy, doSw1jdz |
560 |
> |
real (kind=dp) :: doSw1idux, doSw1iduy, doSw1iduz |
561 |
> |
real (kind=dp) :: doSw1jdux, doSw1jduy, doSw1jduz |
562 |
> |
real (kind=dp) :: doSw2idx, doSw2idy, doSw2idz, doSw2jdx, doSw2jdy, doSw2jdz |
563 |
> |
real (kind=dp) :: doSw2idux, doSw2iduy, doSw2iduz |
564 |
> |
real (kind=dp) :: doSw2jdux, doSw2jduy, doSw2jduz |
565 |
> |
|
566 |
|
if (.not.allocated(StickyMap)) then |
567 |
|
call handleError("sticky", "no StickyMap was present before first call of do_sticky_power_pair!") |
568 |
|
return |
602 |
|
|
603 |
|
if ( rij .LE. rbig ) then |
604 |
|
|
605 |
< |
r3 = r2*rij |
606 |
< |
r5 = r3*r2 |
607 |
< |
|
608 |
< |
drdx = d(1) / rij |
609 |
< |
drdy = d(2) / rij |
610 |
< |
drdz = d(3) / rij |
605 |
> |
rI = 1.0d0/rij |
606 |
> |
rI2 = rI*rI |
607 |
> |
rI3 = rI2*rI |
608 |
> |
rI4 = rI2*rI2 |
609 |
> |
rI5 = rI3*rI2 |
610 |
> |
rI6 = rI3*rI3 |
611 |
> |
rI7 = rI5*rI2 |
612 |
> |
|
613 |
> |
drdx = d(1) * rI |
614 |
> |
drdy = d(2) * rI |
615 |
> |
drdz = d(3) * rI |
616 |
|
|
617 |
|
#ifdef IS_MPI |
618 |
|
! rotate the inter-particle separation into the two different |
645 |
|
xi2 = xi*xi |
646 |
|
yi2 = yi*yi |
647 |
|
zi2 = zi*zi |
648 |
< |
|
648 |
> |
zi3 = zi2*zi |
649 |
> |
zi4 = zi2*zi2 |
650 |
> |
zi5 = zi4*zi |
651 |
> |
xihat = xi*rI |
652 |
> |
yihat = yi*rI |
653 |
> |
zihat = zi*rI |
654 |
> |
|
655 |
|
xj2 = xj*xj |
656 |
|
yj2 = yj*yj |
657 |
|
zj2 = zj*zj |
658 |
< |
|
658 |
> |
zj3 = zj2*zj |
659 |
> |
zj4 = zj2*zj2 |
660 |
> |
zj5 = zj4*zj |
661 |
> |
xjhat = xj*rI |
662 |
> |
yjhat = yj*rI |
663 |
> |
zjhat = zj*rI |
664 |
> |
|
665 |
|
call calc_sw_fnc(rij, rl, ru, rlp, rup, s, sp, dsdr, dspdr) |
666 |
+ |
|
667 |
+ |
wi = 2.0d0*(xi2-yi2)*zi * rI3 |
668 |
+ |
wj = 2.0d0*(xj2-yj2)*zj * rI3 |
669 |
+ |
|
670 |
+ |
! prodVal = zihat*zjhat |
671 |
+ |
! if (prodVal .ge. 0.0d0) then |
672 |
+ |
! wi = 0.0d0 |
673 |
+ |
! wj = 0.0d0 |
674 |
+ |
! endif |
675 |
|
|
636 |
– |
wi = 2.0d0*(xi2-yi2)*zi / r3 |
637 |
– |
wj = 2.0d0*(xj2-yj2)*zj / r3 |
638 |
– |
!rootwi = sqrt(abs(wi)) |
639 |
– |
!rootwj = sqrt(abs(wj)) |
676 |
|
wi2 = wi*wi |
677 |
|
wj2 = wj*wj |
678 |
|
|
643 |
– |
|
679 |
|
w = wi*wi2+wj*wj2 |
680 |
|
|
681 |
< |
zif = zi/rij - 0.6d0 |
682 |
< |
zis = zi/rij + 0.8d0 |
681 |
> |
zif = zihat - 0.6d0 |
682 |
> |
zis = zihat + 0.8d0 |
683 |
|
|
684 |
< |
zjf = zj/rij - 0.6d0 |
685 |
< |
zjs = zj/rij + 0.8d0 |
684 |
> |
zjf = zjhat - 0.6d0 |
685 |
> |
zjs = zjhat + 0.8d0 |
686 |
|
|
687 |
|
wip = zif*zif*zis*zis - w0 |
688 |
|
wjp = zjf*zjf*zjs*zjs - w0 |
689 |
|
wp = wip + wjp |
690 |
+ |
|
691 |
+ |
!wip = zihat - 0.2d0 |
692 |
+ |
!wjp = zjhat - 0.2d0 |
693 |
+ |
!wip3 = wip*wip*wip |
694 |
+ |
!wjp3 = wjp*wjp*wjp |
695 |
+ |
|
696 |
+ |
!wp = wip3*wip + wjp3*wjp |
697 |
|
|
698 |
|
vpair = vpair + 0.5d0*(v0*s*w + v0p*sp*wp) |
699 |
+ |
|
700 |
|
if (do_pot) then |
701 |
|
#ifdef IS_MPI |
702 |
< |
pot_row(atom1) = pot_row(atom1) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw |
703 |
< |
pot_col(atom2) = pot_col(atom2) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw |
702 |
> |
pot_row(atom1) = pot_row(atom1) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw |
703 |
> |
pot_col(atom2) = pot_col(atom2) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw |
704 |
|
#else |
705 |
< |
pot = pot + 0.5d0*(v0*s*w + v0p*sp*wp)*sw |
705 |
> |
pot = pot + 0.5d0*(v0*s*w + v0p*sp*wp)*sw |
706 |
|
#endif |
707 |
|
endif |
708 |
|
|
709 |
< |
! dwidx = 1.5d0*rootwi*( 4.0d0*xi*zi/r3 - 6.0d0*xi*zi*(xi2-yi2)/r5 ) |
710 |
< |
! dwidy = 1.5d0*rootwi*( -4.0d0*yi*zi/r3 - 6.0d0*yi*zi*(xi2-yi2)/r5 ) |
711 |
< |
! dwidz = 1.5d0*rootwi*( 2.0d0*(xi2-yi2)/r3 - 6.0d0*zi2*(xi2-yi2)/r5 ) |
709 |
> |
dwidx = 3.0d0*wi2*( 4.0d0*xi*zi*rI3 - 6.0d0*xi*zi*(xi2-yi2)*rI5 ) |
710 |
> |
dwidy = 3.0d0*wi2*( -4.0d0*yi*zi*rI3 - 6.0d0*yi*zi*(xi2-yi2)*rI5 ) |
711 |
> |
dwidz = 3.0d0*wi2*( 2.0d0*(xi2-yi2)*rI3 - 6.0d0*zi2*(xi2-yi2)*rI5 ) |
712 |
|
|
713 |
< |
! dwjdx = 1.5d0*rootwj*( 4.0d0*xj*zj/r3 - 6.0d0*xj*zj*(xj2-yj2)/r5 ) |
714 |
< |
! dwjdy = 1.5d0*rootwj*( -4.0d0*yj*zj/r3 - 6.0d0*yj*zj*(xj2-yj2)/r5 ) |
715 |
< |
! dwjdz = 1.5d0*rootwj*( 2.0d0*(xj2-yj2)/r3 - 6.0d0*zj2*(xj2-yj2)/r5 ) |
673 |
< |
|
674 |
< |
dwidx = 3.0d0*wi2*( 4.0d0*xi*zi/r3 - 6.0d0*xi*zi*(xi2-yi2)/r5 ) |
675 |
< |
dwidy = 3.0d0*wi2*( -4.0d0*yi*zi/r3 - 6.0d0*yi*zi*(xi2-yi2)/r5 ) |
676 |
< |
dwidz = 3.0d0*wi2*( 2.0d0*(xi2-yi2)/r3 - 6.0d0*zi2*(xi2-yi2)/r5 ) |
713 |
> |
dwjdx = 3.0d0*wj2*( 4.0d0*xj*zj*rI3 - 6.0d0*xj*zj*(xj2-yj2)*rI5 ) |
714 |
> |
dwjdy = 3.0d0*wj2*( -4.0d0*yj*zj*rI3 - 6.0d0*yj*zj*(xj2-yj2)*rI5 ) |
715 |
> |
dwjdz = 3.0d0*wj2*( 2.0d0*(xj2-yj2)*rI3 - 6.0d0*zj2*(xj2-yj2)*rI5 ) |
716 |
|
|
678 |
– |
dwjdx = 3.0d0*wj2*( 4.0d0*xj*zj/r3 - 6.0d0*xj*zj*(xj2-yj2)/r5 ) |
679 |
– |
dwjdy = 3.0d0*wj2*( -4.0d0*yj*zj/r3 - 6.0d0*yj*zj*(xj2-yj2)/r5 ) |
680 |
– |
dwjdz = 3.0d0*wj2*( 2.0d0*(xj2-yj2)/r3 - 6.0d0*zj2*(xj2-yj2)/r5 ) |
681 |
– |
|
717 |
|
uglyi = zif*zif*zis + zif*zis*zis |
718 |
|
uglyj = zjf*zjf*zjs + zjf*zjs*zjs |
719 |
|
|
720 |
< |
dwipdx = -2.0d0*xi*zi*uglyi/r3 |
721 |
< |
dwipdy = -2.0d0*yi*zi*uglyi/r3 |
722 |
< |
dwipdz = 2.0d0*(1.0d0/rij - zi2/r3)*uglyi |
720 |
> |
dwipdx = -2.0d0*xi*zi*uglyi*rI3 |
721 |
> |
dwipdy = -2.0d0*yi*zi*uglyi*rI3 |
722 |
> |
dwipdz = 2.0d0*(rI - zi2*rI3)*uglyi |
723 |
|
|
724 |
< |
dwjpdx = -2.0d0*xj*zj*uglyj/r3 |
725 |
< |
dwjpdy = -2.0d0*yj*zj*uglyj/r3 |
726 |
< |
dwjpdz = 2.0d0*(1.0d0/rij - zj2/r3)*uglyj |
724 |
> |
dwjpdx = -2.0d0*xj*zj*uglyj*rI3 |
725 |
> |
dwjpdy = -2.0d0*yj*zj*uglyj*rI3 |
726 |
> |
dwjpdz = 2.0d0*(rI - zj2*rI3)*uglyj |
727 |
|
|
728 |
< |
! dwidux = 1.5d0*rootwi*( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))/r3 ) |
729 |
< |
! dwiduy = 1.5d0*rootwi*( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))/r3 ) |
730 |
< |
! dwiduz = 1.5d0*rootwi*( -8.0d0*xi*yi*zi/r3 ) |
728 |
> |
!dwipdx = -4.0d0*wip3*zi*xihat |
729 |
> |
!dwipdy = -4.0d0*wip3*zi*yihat |
730 |
> |
!dwipdz = -4.0d0*wip3*(zi2 - 1.0d0)*rI |
731 |
|
|
732 |
< |
! dwjdux = 1.5d0*rootwj*( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))/r3 ) |
733 |
< |
! dwjduy = 1.5d0*rootwj*( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))/r3 ) |
734 |
< |
! dwjduz = 1.5d0*rootwj*( -8.0d0*xj*yj*zj/r3 ) |
732 |
> |
!dwjpdx = -4.0d0*wjp3*zj*xjhat |
733 |
> |
!dwjpdy = -4.0d0*wjp3*zj*yjhat |
734 |
> |
!dwjpdz = -4.0d0*wjp3*(zj2 - 1.0d0)*rI |
735 |
|
|
736 |
< |
dwidux = 3.0d0*wi2*( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))/r3 ) |
737 |
< |
dwiduy = 3.0d0*wi2*( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))/r3 ) |
738 |
< |
dwiduz = 3.0d0*wi2*( -8.0d0*xi*yi*zi/r3 ) |
736 |
> |
!dwipdx = 0.0d0 |
737 |
> |
!dwipdy = 0.0d0 |
738 |
> |
!dwipdz = 0.0d0 |
739 |
|
|
740 |
< |
dwjdux = 3.0d0*wj2*( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))/r3 ) |
741 |
< |
dwjduy = 3.0d0*wj2*( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))/r3 ) |
742 |
< |
dwjduz = 3.0d0*wj2*( -8.0d0*xj*yj*zj/r3 ) |
740 |
> |
!dwjpdx = 0.0d0 |
741 |
> |
!dwjpdy = 0.0d0 |
742 |
> |
!dwjpdz = 0.0d0 |
743 |
> |
|
744 |
> |
dwidux = 3.0d0*wi2*( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))*rI3 ) |
745 |
> |
dwiduy = 3.0d0*wi2*( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))*rI3 ) |
746 |
> |
dwiduz = 3.0d0*wi2*( -8.0d0*xi*yi*zi*rI3 ) |
747 |
|
|
748 |
< |
dwipdux = 2.0d0*yi*uglyi/rij |
749 |
< |
dwipduy = -2.0d0*xi*uglyi/rij |
748 |
> |
dwjdux = 3.0d0*wj2*( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))*rI3 ) |
749 |
> |
dwjduy = 3.0d0*wj2*( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))*rI3 ) |
750 |
> |
dwjduz = 3.0d0*wj2*( -8.0d0*xj*yj*zj*rI3 ) |
751 |
> |
|
752 |
> |
dwipdux = 2.0d0*yi*uglyi*rI |
753 |
> |
dwipduy = -2.0d0*xi*uglyi*rI |
754 |
|
dwipduz = 0.0d0 |
755 |
|
|
756 |
< |
dwjpdux = 2.0d0*yj*uglyj/rij |
757 |
< |
dwjpduy = -2.0d0*xj*uglyj/rij |
756 |
> |
dwjpdux = 2.0d0*yj*uglyj*rI |
757 |
> |
dwjpduy = -2.0d0*xj*uglyj*rI |
758 |
|
dwjpduz = 0.0d0 |
759 |
|
|
760 |
+ |
!dwipdux = 4.0d0*wip3*yihat |
761 |
+ |
!dwipduy = -4.0d0*wip3*xihat |
762 |
+ |
!dwipduz = 0.0d0 |
763 |
+ |
|
764 |
+ |
!dwjpdux = 4.0d0*wjp3*yjhat |
765 |
+ |
!dwjpduy = -4.0d0*wjp3*xjhat |
766 |
+ |
!dwjpduz = 0.0d0 |
767 |
+ |
|
768 |
+ |
!dwipdux = 0.0d0 |
769 |
+ |
!dwipduy = 0.0d0 |
770 |
+ |
!dwipduz = 0.0d0 |
771 |
+ |
|
772 |
+ |
!dwjpdux = 0.0d0 |
773 |
+ |
!dwjpduy = 0.0d0 |
774 |
+ |
!dwjpduz = 0.0d0 |
775 |
+ |
|
776 |
|
! do the torques first since they are easy: |
777 |
|
! remember that these are still in the body fixed axes |
778 |
|
|
873 |
|
|
874 |
|
! now assemble these with the radial-only terms: |
875 |
|
|
876 |
< |
fxradial = 0.5d0*(v0*dsdr*drdx*w + v0p*dspdr*drdx*wp + fxii + fxji) |
877 |
< |
fyradial = 0.5d0*(v0*dsdr*drdy*w + v0p*dspdr*drdy*wp + fyii + fyji) |
878 |
< |
fzradial = 0.5d0*(v0*dsdr*drdz*w + v0p*dspdr*drdz*wp + fzii + fzji) |
876 |
> |
fxradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdx + fxii + fxji) |
877 |
> |
fyradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdy + fyii + fyji) |
878 |
> |
fzradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdz + fzii + fzji) |
879 |
|
|
880 |
|
#ifdef IS_MPI |
881 |
|
f_Row(1,atom1) = f_Row(1,atom1) + fxradial |