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

Comparing trunk/OOPSE-2.0/src/UseTheForce/DarkSide/sticky.F90 (file contents):
Revision 2121 by chrisfen, Sat Mar 12 19:05:16 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.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.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 69 | Line 69 | module sticky
69  
70    public :: newStickyType
71    public :: do_sticky_pair
72 +  public :: destroyStickyTypes
73 +  public :: do_sticky_power_pair
74  
75  
76    type :: StickyList
# Line 82 | Line 84 | module sticky
84       real( kind = dp ) :: rup = 0.0_dp
85       real( kind = dp ) :: rbig = 0.0_dp
86    end type StickyList
87 <  
87 >
88    type(StickyList), dimension(:),allocatable :: StickyMap
89  
90   contains
# Line 96 | Line 98 | contains
98      real( kind = dp ), intent(in) :: rlp, rup
99      integer :: nATypes, myATID
100  
101 <    
101 >
102      isError = 0
103      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
104 <    
104 >
105      !! Be simple-minded and assume that we need a StickyMap that
106      !! is the same size as the total number of atom types
107  
# Line 128 | Line 130 | contains
130      StickyMap(myATID)%c_ident = c_ident
131  
132      ! we could pass all 5 parameters if we felt like it...
133 <    
133 >
134      StickyMap(myATID)%w0 = w0
135      StickyMap(myATID)%v0 = v0
136      StickyMap(myATID)%v0p = v0p
# Line 142 | Line 144 | contains
144      else
145         StickyMap(myATID)%rbig = StickyMap(myATID)%rup
146      endif
147 <  
147 >
148      return
149    end subroutine newStickyType
150  
151    subroutine do_sticky_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
152         pot, A, f, t, do_pot)
153 <    
153 >
154      !! This routine does only the sticky portion of the SSD potential
155      !! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
156      !! The Lennard-Jones and dipolar interaction must be handled separately.
157 <    
157 >
158      !! We assume that the rotation matrices have already been calculated
159      !! and placed in the A array.
160  
# Line 186 | Line 188 | contains
188      real (kind=dp) :: radcomxj, radcomyj, radcomzj
189      integer :: id1, id2
190      integer :: me1, me2
191 <   real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig
191 >    real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig
192  
193 < if (.not.allocated(StickyMap)) then
193 >    if (.not.allocated(StickyMap)) then
194         call handleError("sticky", "no StickyMap was present before first call of do_sticky_pair!")
195         return
196      end if
197 <    
197 >
198   #ifdef IS_MPI
199      me1 = atid_Row(atom1)
200      me2 = atid_Col(atom2)
# Line 460 | Line 462 | if (.not.allocated(StickyMap)) then
462         id1 = atom1
463         id2 = atom2
464   #endif
465 <      
465 >
466         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
467 <          
467 >
468            fpair(1) = fpair(1) + fxradial
469            fpair(2) = fpair(2) + fyradial
470            fpair(3) = fpair(3) + fzradial
471 <          
471 >
472         endif
473      endif
474    end subroutine do_sticky_pair
475  
476    !! calculates the switching functions and their derivatives for a given
477    subroutine calc_sw_fnc(r, rl, ru, rlp, rup, s, sp, dsdr, dspdr)
478 <    
478 >
479      real (kind=dp), intent(in) :: r, rl, ru, rlp, rup
480      real (kind=dp), intent(inout) :: s, sp, dsdr, dspdr
481 <    
481 >
482      ! distances must be in angstroms
483 <    
483 >
484      if (r.lt.rl) then
485         s = 1.0d0
486         dsdr = 0.0d0
# Line 502 | Line 504 | if (.not.allocated(StickyMap)) then
504              ((rup - rlp)**3)
505         dspdr = 6.0d0*(r-rup)*(r-rlp)/((rup - rlp)**3)      
506      endif
507 <    
507 >
508      return
509    end subroutine calc_sw_fnc
510 +
511 +  subroutine destroyStickyTypes()  
512 +    if(allocated(StickyMap)) deallocate(StickyMap)
513 +  end subroutine destroyStickyTypes
514 +  
515 +    subroutine do_sticky_power_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
516 +       pot, A, f, t, do_pot)
517 +    !! We assume that the rotation matrices have already been calculated
518 +    !! and placed in the A array.
519 +
520 +    !! i and j are pointers to the two SSD atoms
521 +
522 +    integer, intent(in) :: atom1, atom2
523 +    real (kind=dp), intent(inout) :: rij, r2
524 +    real (kind=dp), dimension(3), intent(in) :: d
525 +    real (kind=dp), dimension(3), intent(inout) :: fpair
526 +    real (kind=dp) :: pot, vpair, sw
527 +    real (kind=dp), dimension(9,nLocal) :: A
528 +    real (kind=dp), dimension(3,nLocal) :: f
529 +    real (kind=dp), dimension(3,nLocal) :: t
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) :: 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
539 +    real (kind=dp) :: dwipdux, dwipduy, dwipduz, dwjpdux, dwjpduy, dwjpduz
540 +    real (kind=dp) :: zif, zis, zjf, zjs, uglyi, uglyj
541 +    real (kind=dp) :: drdx, drdy, drdz
542 +    real (kind=dp) :: txi, tyi, tzi, txj, tyj, tzj
543 +    real (kind=dp) :: fxii, fyii, fzii, fxjj, fyjj, fzjj
544 +    real (kind=dp) :: fxij, fyij, fzij, fxji, fyji, fzji      
545 +    real (kind=dp) :: fxradial, fyradial, fzradial
546 +    real (kind=dp) :: rijtest, rjitest
547 +    real (kind=dp) :: radcomxi, radcomyi, radcomzi
548 +    real (kind=dp) :: radcomxj, radcomyj, radcomzj
549 +    integer :: id1, id2
550 +    integer :: me1, me2
551 +    real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig
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
569 +    end if
570 +
571 + #ifdef IS_MPI
572 +    me1 = atid_Row(atom1)
573 +    me2 = atid_Col(atom2)
574 + #else
575 +    me1 = atid(atom1)
576 +    me2 = atid(atom2)
577 + #endif
578 +
579 +    if (me1.eq.me2) then
580 +       w0  = StickyMap(me1)%w0
581 +       v0  = StickyMap(me1)%v0
582 +       v0p = StickyMap(me1)%v0p
583 +       rl  = StickyMap(me1)%rl
584 +       ru  = StickyMap(me1)%ru
585 +       rlp = StickyMap(me1)%rlp
586 +       rup = StickyMap(me1)%rup
587 +       rbig = StickyMap(me1)%rbig
588 +    else
589 +       ! This is silly, but if you want 2 sticky types in your
590 +       ! simulation, we'll let you do it with the Lorentz-
591 +       ! Berthelot mixing rules.
592 +       ! (Warning: you'll be SLLLLLLLLLLLLLLLOOOOOOOOOOWWWWWWWWWWW)
593 +       rl   = 0.5_dp * ( StickyMap(me1)%rl + StickyMap(me2)%rl )
594 +       ru   = 0.5_dp * ( StickyMap(me1)%ru + StickyMap(me2)%ru )
595 +       rlp  = 0.5_dp * ( StickyMap(me1)%rlp + StickyMap(me2)%rlp )
596 +       rup  = 0.5_dp * ( StickyMap(me1)%rup + StickyMap(me2)%rup )
597 +       rbig = max(ru, rup)
598 +       w0  = sqrt( StickyMap(me1)%w0   * StickyMap(me2)%w0  )
599 +       v0  = sqrt( StickyMap(me1)%v0   * StickyMap(me2)%v0  )
600 +       v0p = sqrt( StickyMap(me1)%v0p  * StickyMap(me2)%v0p )
601 +    endif
602 +
603 +    if ( rij .LE. rbig ) then
604 +
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
619 +       ! body-fixed coordinate systems:
620 +
621 +       xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
622 +       yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
623 +       zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
624 +
625 +       ! negative sign because this is the vector from j to i:
626 +
627 +       xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
628 +       yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
629 +       zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
630 + #else
631 +       ! rotate the inter-particle separation into the two different
632 +       ! body-fixed coordinate systems:
633 +
634 +       xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
635 +       yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
636 +       zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
637 +
638 +       ! negative sign because this is the vector from j to i:
639 +
640 +       xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
641 +       yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
642 +       zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
643 + #endif
644 +
645 +       xi2 = xi*xi
646 +       yi2 = yi*yi
647 +       zi2 = zi*zi
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 +       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 +
676 +       wi2 = wi*wi
677 +       wj2 = wj*wj
678 +
679 +       w = wi*wi2+wj*wj2
680 +
681 +       zif = zihat - 0.6d0
682 +       zis = zihat + 0.8d0
683 +
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
704 + #else
705 +         pot = pot + 0.5d0*(v0*s*w + v0p*sp*wp)*sw
706 + #endif  
707 +       endif
708 +
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 = 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 +
717 +       uglyi = zif*zif*zis + zif*zis*zis
718 +       uglyj = zjf*zjf*zjs + zjf*zjs*zjs
719 +
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*rI3
725 +       dwjpdy = -2.0d0*yj*zj*uglyj*rI3
726 +       dwjpdz = 2.0d0*(rI - zj2*rI3)*uglyj
727 +
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 +       !dwjpdx = -4.0d0*wjp3*zj*xjhat
733 +       !dwjpdy = -4.0d0*wjp3*zj*yjhat
734 +       !dwjpdz = -4.0d0*wjp3*(zj2 - 1.0d0)*rI
735 +      
736 +       !dwipdx = 0.0d0
737 +       !dwipdy = 0.0d0
738 +       !dwipdz = 0.0d0
739 +
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 +       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*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 +
779 +       txi = 0.5d0*(v0*s*dwidux + v0p*sp*dwipdux)*sw
780 +       tyi = 0.5d0*(v0*s*dwiduy + v0p*sp*dwipduy)*sw
781 +       tzi = 0.5d0*(v0*s*dwiduz + v0p*sp*dwipduz)*sw
782 +
783 +       txj = 0.5d0*(v0*s*dwjdux + v0p*sp*dwjpdux)*sw
784 +       tyj = 0.5d0*(v0*s*dwjduy + v0p*sp*dwjpduy)*sw
785 +       tzj = 0.5d0*(v0*s*dwjduz + v0p*sp*dwjpduz)*sw
786 +
787 +       ! go back to lab frame using transpose of rotation matrix:
788 +
789 + #ifdef IS_MPI
790 +       t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
791 +            a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
792 +       t_Row(2,atom1) = t_Row(2,atom1) + a_Row(2,atom1)*txi + &
793 +            a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
794 +       t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
795 +            a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
796 +
797 +       t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
798 +            a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
799 +       t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
800 +            a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
801 +       t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
802 +            a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
803 + #else
804 +       t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
805 +       t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
806 +       t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
807 +
808 +       t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
809 +       t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
810 +       t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
811 + #endif    
812 +       ! Now, on to the forces:
813 +
814 +       ! first rotate the i terms back into the lab frame:
815 +
816 +       radcomxi = (v0*s*dwidx+v0p*sp*dwipdx)*sw
817 +       radcomyi = (v0*s*dwidy+v0p*sp*dwipdy)*sw
818 +       radcomzi = (v0*s*dwidz+v0p*sp*dwipdz)*sw
819 +
820 +       radcomxj = (v0*s*dwjdx+v0p*sp*dwjpdx)*sw
821 +       radcomyj = (v0*s*dwjdy+v0p*sp*dwjpdy)*sw
822 +       radcomzj = (v0*s*dwjdz+v0p*sp*dwjpdz)*sw
823 +
824 + #ifdef IS_MPI    
825 +       fxii = a_Row(1,atom1)*(radcomxi) + &
826 +            a_Row(4,atom1)*(radcomyi) + &
827 +            a_Row(7,atom1)*(radcomzi)
828 +       fyii = a_Row(2,atom1)*(radcomxi) + &
829 +            a_Row(5,atom1)*(radcomyi) + &
830 +            a_Row(8,atom1)*(radcomzi)
831 +       fzii = a_Row(3,atom1)*(radcomxi) + &
832 +            a_Row(6,atom1)*(radcomyi) + &
833 +            a_Row(9,atom1)*(radcomzi)
834 +
835 +       fxjj = a_Col(1,atom2)*(radcomxj) + &
836 +            a_Col(4,atom2)*(radcomyj) + &
837 +            a_Col(7,atom2)*(radcomzj)
838 +       fyjj = a_Col(2,atom2)*(radcomxj) + &
839 +            a_Col(5,atom2)*(radcomyj) + &
840 +            a_Col(8,atom2)*(radcomzj)
841 +       fzjj = a_Col(3,atom2)*(radcomxj)+ &
842 +            a_Col(6,atom2)*(radcomyj) + &
843 +            a_Col(9,atom2)*(radcomzj)
844 + #else
845 +       fxii = a(1,atom1)*(radcomxi) + &
846 +            a(4,atom1)*(radcomyi) + &
847 +            a(7,atom1)*(radcomzi)
848 +       fyii = a(2,atom1)*(radcomxi) + &
849 +            a(5,atom1)*(radcomyi) + &
850 +            a(8,atom1)*(radcomzi)
851 +       fzii = a(3,atom1)*(radcomxi) + &
852 +            a(6,atom1)*(radcomyi) + &
853 +            a(9,atom1)*(radcomzi)
854 +
855 +       fxjj = a(1,atom2)*(radcomxj) + &
856 +            a(4,atom2)*(radcomyj) + &
857 +            a(7,atom2)*(radcomzj)
858 +       fyjj = a(2,atom2)*(radcomxj) + &
859 +            a(5,atom2)*(radcomyj) + &
860 +            a(8,atom2)*(radcomzj)
861 +       fzjj = a(3,atom2)*(radcomxj)+ &
862 +            a(6,atom2)*(radcomyj) + &
863 +            a(9,atom2)*(radcomzj)
864 + #endif
865 +
866 +       fxij = -fxii
867 +       fyij = -fyii
868 +       fzij = -fzii
869 +
870 +       fxji = -fxjj
871 +       fyji = -fyjj
872 +       fzji = -fzjj
873 +
874 +       ! now assemble these with the radial-only terms:
875 +
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
882 +       f_Row(2,atom1) = f_Row(2,atom1) + fyradial
883 +       f_Row(3,atom1) = f_Row(3,atom1) + fzradial
884 +
885 +       f_Col(1,atom2) = f_Col(1,atom2) - fxradial
886 +       f_Col(2,atom2) = f_Col(2,atom2) - fyradial
887 +       f_Col(3,atom2) = f_Col(3,atom2) - fzradial
888 + #else
889 +       f(1,atom1) = f(1,atom1) + fxradial
890 +       f(2,atom1) = f(2,atom1) + fyradial
891 +       f(3,atom1) = f(3,atom1) + fzradial
892 +
893 +       f(1,atom2) = f(1,atom2) - fxradial
894 +       f(2,atom2) = f(2,atom2) - fyradial
895 +       f(3,atom2) = f(3,atom2) - fzradial
896 + #endif
897 +
898 + #ifdef IS_MPI
899 +       id1 = AtomRowToGlobal(atom1)
900 +       id2 = AtomColToGlobal(atom2)
901 + #else
902 +       id1 = atom1
903 +       id2 = atom2
904 + #endif
905 +
906 +       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
907 +
908 +          fpair(1) = fpair(1) + fxradial
909 +          fpair(2) = fpair(2) + fyradial
910 +          fpair(3) = fpair(3) + fzradial
911 +
912 +       endif
913 +    endif
914 +  end subroutine do_sticky_power_pair
915 +
916   end module sticky

Diff Legend

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