ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/sticky.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/sticky.F90 (file contents):
Revision 2220 by chrisfen, Thu May 5 14:47:35 2005 UTC vs.
Revision 2229 by chrisfen, Tue May 17 22:35:01 2005 UTC

# Line 50 | Line 50
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.10 2005-05-17 22:35:01 chrisfen Exp $, $Date: 2005-05-17 22:35:01 $, $Name: not supported by cvs2svn $, $Revision: 1.10 $
54  
55   module sticky
56  
# Line 530 | Line 530 | contains
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
# Line 548 | Line 549 | contains
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) :: frac1, frac2, 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
# Line 588 | Line 602 | contains
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 = rI4*rI3
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
# Line 626 | Line 645 | contains
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 +       frac1 = 0.5d0
668 +       frac2 = 0.5d0
669 +      
670 +       wi = 2.0d0*(xi2-yi2)*zi*rI3
671 +       wj = 2.0d0*(xj2-yj2)*zj*rI3
672 +      
673 + !       prodVal = zihat*zjhat
674 + !       if (prodVal .ge. 0.0d0) then
675 + !         wi = 0.0d0
676 + !         wj = 0.0d0
677 + !       endif
678  
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))
679         wi2 = wi*wi
680         wj2 = wj*wj
681  
682 <      
644 <       w = wi*wi2+wj*wj2
682 >       w = frac1*wi*wi2 + frac2*wi + wj*wj2
683  
684 <       zif = zi/rij - 0.6d0
685 <       zis = zi/rij + 0.8d0
684 >       zif = zihat - 0.6d0
685 >       zis = zihat + 0.8d0
686  
687 <       zjf = zj/rij - 0.6d0
688 <       zjs = zj/rij + 0.8d0
687 >       zjf = zjhat - 0.6d0
688 >       zjs = zjhat + 0.8d0
689  
690         wip = zif*zif*zis*zis - w0
691         wjp = zjf*zjf*zjs*zjs - w0
692         wp = wip + wjp
693 +        
694 +       !wip = zihat - 0.2d0
695 +       !wjp = zjhat - 0.2d0
696 +       !wip3 = wip*wip*wip
697 +       !wjp3 = wjp*wjp*wjp
698 +      
699 +       !wp = wip3*wip + wjp3*wjp
700  
701         vpair = vpair + 0.5d0*(v0*s*w + v0p*sp*wp)
702 +      
703         if (do_pot) then
704   #ifdef IS_MPI
705 <          pot_row(atom1) = pot_row(atom1) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw
706 <          pot_col(atom2) = pot_col(atom2) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw
705 >         pot_row(atom1) = pot_row(atom1) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw
706 >         pot_col(atom2) = pot_col(atom2) + 0.25d0*(v0*s*w + v0p*sp*wp)*sw
707   #else
708 <          pot = pot + 0.5d0*(v0*s*w + v0p*sp*wp)*sw
708 >         pot = pot + 0.5d0*(v0*s*w + v0p*sp*wp)*sw
709   #endif  
710         endif
711  
712 < !       dwidx = 1.5d0*rootwi*( 4.0d0*xi*zi/r3 - 6.0d0*xi*zi*(xi2-yi2)/r5 )
713 < !       dwidy = 1.5d0*rootwi*( -4.0d0*yi*zi/r3 - 6.0d0*yi*zi*(xi2-yi2)/r5 )
714 < !       dwidz = 1.5d0*rootwi*( 2.0d0*(xi2-yi2)/r3 - 6.0d0*zi2*(xi2-yi2)/r5 )
669 <
670 < !       dwjdx = 1.5d0*rootwj*( 4.0d0*xj*zj/r3  - 6.0d0*xj*zj*(xj2-yj2)/r5 )
671 < !       dwjdy = 1.5d0*rootwj*( -4.0d0*yj*zj/r3  - 6.0d0*yj*zj*(xj2-yj2)/r5 )
672 < !       dwjdz = 1.5d0*rootwj*( 2.0d0*(xj2-yj2)/r3  - 6.0d0*zj2*(xj2-yj2)/r5 )
712 >       dwidx = ( 4.0d0*xi*zi*rI3 - 6.0d0*xi*zi*(xi2-yi2)*rI5 )
713 >       dwidy = ( -4.0d0*yi*zi*rI3 - 6.0d0*yi*zi*(xi2-yi2)*rI5 )
714 >       dwidz = ( 2.0d0*(xi2-yi2)*rI3 - 6.0d0*zi2*(xi2-yi2)*rI5 )
715        
716 <       dwidx = 3.0d0*wi2*( 4.0d0*xi*zi/r3 - 6.0d0*xi*zi*(xi2-yi2)/r5 )
717 <       dwidy = 3.0d0*wi2*( -4.0d0*yi*zi/r3 - 6.0d0*yi*zi*(xi2-yi2)/r5 )
718 <       dwidz = 3.0d0*wi2*( 2.0d0*(xi2-yi2)/r3 - 6.0d0*zi2*(xi2-yi2)/r5 )
716 >       dwidx = frac1*3.0d0*wi2*dwidx + frac2*dwidx
717 >       dwidy = frac1*3.0d0*wi2*dwidy + frac2*dwidy
718 >       dwidz = frac1*3.0d0*wi2*dwidz + frac2*dwidz
719  
720 <       dwjdx = 3.0d0*wj2*( 4.0d0*xj*zj/r3  - 6.0d0*xj*zj*(xj2-yj2)/r5 )
721 <       dwjdy = 3.0d0*wj2*( -4.0d0*yj*zj/r3  - 6.0d0*yj*zj*(xj2-yj2)/r5 )
722 <       dwjdz = 3.0d0*wj2*( 2.0d0*(xj2-yj2)/r3  - 6.0d0*zj2*(xj2-yj2)/r5 )
720 >       dwjdx = ( 4.0d0*xj*zj*rI3  - 6.0d0*xj*zj*(xj2-yj2)*rI5 )
721 >       dwjdy = ( -4.0d0*yj*zj*rI3  - 6.0d0*yj*zj*(xj2-yj2)*rI5 )
722 >       dwjdz = ( 2.0d0*(xj2-yj2)*rI3  - 6.0d0*zj2*(xj2-yj2)*rI5 )
723  
724 +       dwjdx = frac1*3.0d0*wj2*dwjdx + frac2*dwjdx
725 +       dwjdy = frac1*3.0d0*wj2*dwjdy + frac2*dwjdy
726 +       dwjdz = frac1*3.0d0*wj2*dwjdz + frac2*dwjdz
727 +      
728         uglyi = zif*zif*zis + zif*zis*zis
729         uglyj = zjf*zjf*zjs + zjf*zjs*zjs
730  
731 <       dwipdx = -2.0d0*xi*zi*uglyi/r3
732 <       dwipdy = -2.0d0*yi*zi*uglyi/r3
733 <       dwipdz = 2.0d0*(1.0d0/rij - zi2/r3)*uglyi
731 >       dwipdx = -2.0d0*xi*zi*uglyi*rI3
732 >       dwipdy = -2.0d0*yi*zi*uglyi*rI3
733 >       dwipdz = 2.0d0*(rI - zi2*rI3)*uglyi
734  
735 <       dwjpdx = -2.0d0*xj*zj*uglyj/r3
736 <       dwjpdy = -2.0d0*yj*zj*uglyj/r3
737 <       dwjpdz = 2.0d0*(1.0d0/rij - zj2/r3)*uglyj
735 >       dwjpdx = -2.0d0*xj*zj*uglyj*rI3
736 >       dwjpdy = -2.0d0*yj*zj*uglyj*rI3
737 >       dwjpdz = 2.0d0*(rI - zj2*rI3)*uglyj
738  
739 < !       dwidux = 1.5d0*rootwi*( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))/r3 )
740 < !       dwiduy = 1.5d0*rootwi*( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))/r3 )
741 < !       dwiduz = 1.5d0*rootwi*( -8.0d0*xi*yi*zi/r3 )
739 >       !dwipdx = -4.0d0*wip3*zi*xihat
740 >       !dwipdy = -4.0d0*wip3*zi*yihat
741 >       !dwipdz = -4.0d0*wip3*(zi2 - 1.0d0)*rI
742  
743 < !       dwjdux = 1.5d0*rootwj*( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))/r3 )
744 < !       dwjduy = 1.5d0*rootwj*( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))/r3 )
745 < !       dwjduz = 1.5d0*rootwj*( -8.0d0*xj*yj*zj/r3 )
743 >       !dwjpdx = -4.0d0*wjp3*zj*xjhat
744 >       !dwjpdy = -4.0d0*wjp3*zj*yjhat
745 >       !dwjpdz = -4.0d0*wjp3*(zj2 - 1.0d0)*rI
746        
747 <       dwidux = 3.0d0*wi2*( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))/r3 )
748 <       dwiduy = 3.0d0*wi2*( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))/r3 )
749 <       dwiduz = 3.0d0*wi2*( -8.0d0*xi*yi*zi/r3 )
747 >       !dwipdx = 0.0d0
748 >       !dwipdy = 0.0d0
749 >       !dwipdz = 0.0d0
750  
751 <       dwjdux = 3.0d0*wj2*( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))/r3 )
752 <       dwjduy = 3.0d0*wj2*( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))/r3 )
753 <       dwjduz = 3.0d0*wj2*( -8.0d0*xj*yj*zj/r3 )
751 >       !dwjpdx = 0.0d0
752 >       !dwjpdy = 0.0d0
753 >       !dwjpdz = 0.0d0
754 >      
755 >       dwidux = ( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))*rI3 )
756 >       dwiduy = ( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))*rI3 )
757 >       dwiduz = ( -8.0d0*xi*yi*zi*rI3 )
758  
759 <       dwipdux =  2.0d0*yi*uglyi/rij
760 <       dwipduy = -2.0d0*xi*uglyi/rij
759 >       dwidux = frac1*3.0d0*wi2*dwidux + frac2*dwidux
760 >       dwiduy = frac1*3.0d0*wi2*dwiduy + frac2*dwiduy
761 >       dwiduz = frac1*3.0d0*wi2*dwiduz + frac2*dwiduz
762 >
763 >       dwjdux = ( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))*rI3 )
764 >       dwjduy = ( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))*rI3 )
765 >       dwjduz = ( -8.0d0*xj*yj*zj*rI3 )
766 >
767 >       dwjdux = frac1*3.0d0*wj2*dwjdux + frac2*dwjdux
768 >       dwjduy = frac1*3.0d0*wj2*dwjduy + frac2*dwjduy
769 >       dwjduz = frac1*3.0d0*wj2*dwjduz + frac2*dwjduz
770 >
771 >       dwipdux =  2.0d0*yi*uglyi*rI
772 >       dwipduy = -2.0d0*xi*uglyi*rI
773         dwipduz =  0.0d0
774  
775 <       dwjpdux =  2.0d0*yj*uglyj/rij
776 <       dwjpduy = -2.0d0*xj*uglyj/rij
775 >       dwjpdux =  2.0d0*yj*uglyj*rI
776 >       dwjpduy = -2.0d0*xj*uglyj*rI
777         dwjpduz =  0.0d0
778  
779 +       !dwipdux =  4.0d0*wip3*yihat
780 +       !dwipduy = -4.0d0*wip3*xihat
781 +       !dwipduz =  0.0d0
782 +
783 +       !dwjpdux =  4.0d0*wjp3*yjhat
784 +       !dwjpduy = -4.0d0*wjp3*xjhat
785 +       !dwjpduz =  0.0d0
786 +
787 +       !dwipdux = 0.0d0
788 +       !dwipduy = 0.0d0
789 +       !dwipduz = 0.0d0
790 +
791 +       !dwjpdux = 0.0d0
792 +       !dwjpduy = 0.0d0
793 +       !dwjpduz = 0.0d0
794 +      
795         ! do the torques first since they are easy:
796         ! remember that these are still in the body fixed axes
797  
# Line 814 | Line 892 | contains
892  
893         ! now assemble these with the radial-only terms:
894  
895 <       fxradial = 0.5d0*(v0*dsdr*drdx*w + v0p*dspdr*drdx*wp + fxii + fxji)
896 <       fyradial = 0.5d0*(v0*dsdr*drdy*w + v0p*dspdr*drdy*wp + fyii + fyji)
897 <       fzradial = 0.5d0*(v0*dsdr*drdz*w + v0p*dspdr*drdz*wp + fzii + fzji)
895 >       fxradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdx + fxii + fxji)
896 >       fyradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdy + fyii + fyji)
897 >       fzradial = 0.5d0*((v0*dsdr*w + v0p*dspdr*wp)*drdz + fzii + fzji)
898  
899   #ifdef IS_MPI
900         f_Row(1,atom1) = f_Row(1,atom1) + fxradial

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines