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 2223 by chrisfen, Thu May 5 14:47:35 2005 UTC vs.
Revision 2224 by chrisfen, Thu May 12 19:43:48 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.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  
# 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) :: 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
# 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 = 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
# 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 +       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  
# Line 814 | Line 873 | contains
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

Diff Legend

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