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 1948 by gezelter, Fri Jan 14 20:31:16 2005 UTC vs.
Revision 2277 by chrisfen, Fri Aug 26 21:30:41 2005 UTC

# Line 48 | Line 48
48   !! Corresponds to the force field defined in ssd_FF.cpp
49   !! @author Charles F. Vardeman II
50   !! @author Matthew Meineke
51 < !! @author Christopher Fennel
51 > !! @author Christopher Fennell
52   !! @author J. Daniel Gezelter
53 < !! @version $Id: sticky.F90,v 1.4 2005-01-14 20:31:16 gezelter Exp $, $Date: 2005-01-14 20:31:16 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
53 > !! @version $Id: sticky.F90,v 1.14 2005-08-26 21:30:41 chrisfen Exp $, $Date: 2005-08-26 21:30:41 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $
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 +  public :: getStickyCut
75 +  public :: getStickyPowerCut
76  
73 –
77    type :: StickyList
78       integer :: c_ident
79       real( kind = dp ) :: w0 = 0.0_dp
# Line 82 | Line 85 | module sticky
85       real( kind = dp ) :: rup = 0.0_dp
86       real( kind = dp ) :: rbig = 0.0_dp
87    end type StickyList
88 <  
88 >
89    type(StickyList), dimension(:),allocatable :: StickyMap
90  
91   contains
# Line 96 | Line 99 | contains
99      real( kind = dp ), intent(in) :: rlp, rup
100      integer :: nATypes, myATID
101  
102 <    
102 >
103      isError = 0
104      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
105 <    
105 >
106      !! Be simple-minded and assume that we need a StickyMap that
107      !! is the same size as the total number of atom types
108  
# Line 128 | Line 131 | contains
131      StickyMap(myATID)%c_ident = c_ident
132  
133      ! we could pass all 5 parameters if we felt like it...
134 <    
134 >
135      StickyMap(myATID)%w0 = w0
136      StickyMap(myATID)%v0 = v0
137      StickyMap(myATID)%v0p = v0p
# Line 142 | Line 145 | contains
145      else
146         StickyMap(myATID)%rbig = StickyMap(myATID)%rup
147      endif
148 <  
148 >
149      return
150    end subroutine newStickyType
151 +
152 +  function getStickyCut(atomID) result(cutValue)
153 +    integer, intent(in) :: atomID
154 +    real(kind=dp) :: cutValue
155 +
156 +    cutValue = StickyMap(atomID)%rbig
157 +  end function getStickyCut
158 +
159 +  function getStickyPowerCut(atomID) result(cutValue)
160 +    integer, intent(in) :: atomID
161 +    real(kind=dp) :: cutValue
162  
163 +    cutValue = StickyMap(atomID)%rbig
164 +  end function getStickyPowerCut
165 +
166    subroutine do_sticky_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
167         pot, A, f, t, do_pot)
168 <    
168 >
169      !! This routine does only the sticky portion of the SSD potential
170      !! [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)].
171      !! The Lennard-Jones and dipolar interaction must be handled separately.
172 <    
172 >
173      !! We assume that the rotation matrices have already been calculated
174      !! and placed in the A array.
175  
# Line 186 | Line 203 | contains
203      real (kind=dp) :: radcomxj, radcomyj, radcomzj
204      integer :: id1, id2
205      integer :: me1, me2
206 <   real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig
206 >    real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig
207  
208 < if (.not.allocated(StickyMap)) then
208 >    if (.not.allocated(StickyMap)) then
209         call handleError("sticky", "no StickyMap was present before first call of do_sticky_pair!")
210         return
211      end if
212 <    
212 >
213   #ifdef IS_MPI
214      me1 = atid_Row(atom1)
215      me2 = atid_Col(atom2)
# Line 460 | Line 477 | if (.not.allocated(StickyMap)) then
477         id1 = atom1
478         id2 = atom2
479   #endif
480 <      
480 >
481         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
482 <          
482 >
483            fpair(1) = fpair(1) + fxradial
484            fpair(2) = fpair(2) + fyradial
485            fpair(3) = fpair(3) + fzradial
486 <          
486 >
487         endif
488      endif
489    end subroutine do_sticky_pair
490  
491    !! calculates the switching functions and their derivatives for a given
492    subroutine calc_sw_fnc(r, rl, ru, rlp, rup, s, sp, dsdr, dspdr)
493 <    
493 >
494      real (kind=dp), intent(in) :: r, rl, ru, rlp, rup
495      real (kind=dp), intent(inout) :: s, sp, dsdr, dspdr
496 <    
496 >
497      ! distances must be in angstroms
498 <    
498 >
499      if (r.lt.rl) then
500         s = 1.0d0
501         dsdr = 0.0d0
# Line 502 | Line 519 | if (.not.allocated(StickyMap)) then
519              ((rup - rlp)**3)
520         dspdr = 6.0d0*(r-rup)*(r-rlp)/((rup - rlp)**3)      
521      endif
522 <    
522 >
523      return
524    end subroutine calc_sw_fnc
525 +
526 +  subroutine destroyStickyTypes()  
527 +    if(allocated(StickyMap)) deallocate(StickyMap)
528 +  end subroutine destroyStickyTypes
529 +  
530 +  subroutine do_sticky_power_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
531 +       pot, A, f, t, do_pot)
532 +    !! We assume that the rotation matrices have already been calculated
533 +    !! and placed in the A array.
534 +    
535 +    !! i and j are pointers to the two SSD atoms
536 +    
537 +    integer, intent(in) :: atom1, atom2
538 +    real (kind=dp), intent(inout) :: rij, r2
539 +    real (kind=dp), dimension(3), intent(in) :: d
540 +    real (kind=dp), dimension(3), intent(inout) :: fpair
541 +    real (kind=dp) :: pot, vpair, sw
542 +    real (kind=dp), dimension(9,nLocal) :: A
543 +    real (kind=dp), dimension(3,nLocal) :: f
544 +    real (kind=dp), dimension(3,nLocal) :: t
545 +    logical, intent(in) :: do_pot
546 +
547 +    real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
548 +    real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
549 +    real (kind=dp) :: rI, rI2, rI3, rI4, rI5, rI6, rI7, s, sp, dsdr, dspdr
550 +    real (kind=dp) :: wi, wj, w, wi2, wj2, eScale, v0scale
551 +    real (kind=dp) :: dwidx, dwidy, dwidz, dwjdx, dwjdy, dwjdz
552 +    real (kind=dp) :: dwidux, dwiduy, dwiduz, dwjdux, dwjduy, dwjduz
553 +    real (kind=dp) :: drdx, drdy, drdz
554 +    real (kind=dp) :: txi, tyi, tzi, txj, tyj, tzj
555 +    real (kind=dp) :: fxii, fyii, fzii, fxjj, fyjj, fzjj
556 +    real (kind=dp) :: fxij, fyij, fzij, fxji, fyji, fzji      
557 +    real (kind=dp) :: fxradial, fyradial, fzradial
558 +    real (kind=dp) :: rijtest, rjitest
559 +    real (kind=dp) :: radcomxi, radcomyi, radcomzi
560 +    real (kind=dp) :: radcomxj, radcomyj, radcomzj
561 +    integer :: id1, id2
562 +    integer :: me1, me2
563 +    real (kind=dp) :: w0, v0, v0p, rl, ru, rlp, rup, rbig
564 +    real (kind=dp) :: zi3, zi4, zi5, zj3, zj4, zj5
565 +    real (kind=dp) :: frac1, frac2
566 +          
567 +    if (.not.allocated(StickyMap)) then
568 +       call handleError("sticky", "no StickyMap was present before first call of do_sticky_power_pair!")
569 +       return
570 +    end if
571 +
572 + #ifdef IS_MPI
573 +    me1 = atid_Row(atom1)
574 +    me2 = atid_Col(atom2)
575 + #else
576 +    me1 = atid(atom1)
577 +    me2 = atid(atom2)
578 + #endif
579 +
580 +    if (me1.eq.me2) then
581 +       w0  = StickyMap(me1)%w0
582 +       v0  = StickyMap(me1)%v0
583 +       v0p = StickyMap(me1)%v0p
584 +       rl  = StickyMap(me1)%rl
585 +       ru  = StickyMap(me1)%ru
586 +       rlp = StickyMap(me1)%rlp
587 +       rup = StickyMap(me1)%rup
588 +       rbig = StickyMap(me1)%rbig
589 +    else
590 +       ! This is silly, but if you want 2 sticky types in your
591 +       ! simulation, we'll let you do it with the Lorentz-
592 +       ! Berthelot mixing rules.
593 +       ! (Warning: you'll be SLLLLLLLLLLLLLLLOOOOOOOOOOWWWWWWWWWWW)
594 +       rl   = 0.5_dp * ( StickyMap(me1)%rl + StickyMap(me2)%rl )
595 +       ru   = 0.5_dp * ( StickyMap(me1)%ru + StickyMap(me2)%ru )
596 +       rlp  = 0.5_dp * ( StickyMap(me1)%rlp + StickyMap(me2)%rlp )
597 +       rup  = 0.5_dp * ( StickyMap(me1)%rup + StickyMap(me2)%rup )
598 +       rbig = max(ru, rup)
599 +       w0  = sqrt( StickyMap(me1)%w0   * StickyMap(me2)%w0  )
600 +       v0  = sqrt( StickyMap(me1)%v0   * StickyMap(me2)%v0  )
601 +       v0p = sqrt( StickyMap(me1)%v0p  * StickyMap(me2)%v0p )
602 +    endif
603 +
604 +    if ( rij .LE. rbig ) then
605 +
606 +       rI = 1.0d0/rij
607 +       rI2 = rI*rI
608 +       rI3 = rI2*rI
609 +       rI4 = rI2*rI2
610 +       rI5 = rI3*rI2
611 +       rI6 = rI3*rI3
612 +       rI7 = rI4*rI3
613 +              
614 +       drdx = d(1) * rI
615 +       drdy = d(2) * rI
616 +       drdz = d(3) * rI
617 +
618 + #ifdef IS_MPI
619 +       ! rotate the inter-particle separation into the two different
620 +       ! body-fixed coordinate systems:
621 +
622 +       xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
623 +       yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
624 +       zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
625 +
626 +       ! negative sign because this is the vector from j to i:
627 +
628 +       xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
629 +       yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
630 +       zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
631 + #else
632 +       ! rotate the inter-particle separation into the two different
633 +       ! body-fixed coordinate systems:
634 +
635 +       xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
636 +       yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
637 +       zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
638 +
639 +       ! negative sign because this is the vector from j to i:
640 +
641 +       xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
642 +       yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
643 +       zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
644 + #endif
645 +
646 +       xi2 = xi*xi
647 +       yi2 = yi*yi
648 +       zi2 = zi*zi
649 +       zi3 = zi2*zi
650 +       zi4 = zi2*zi2
651 +       zi5 = zi3*zi2
652 +       xihat = xi*rI
653 +       yihat = yi*rI
654 +       zihat = zi*rI
655 +      
656 +       xj2 = xj*xj
657 +       yj2 = yj*yj
658 +       zj2 = zj*zj
659 +       zj3 = zj2*zj
660 +       zj4 = zj2*zj2
661 +       zj5 = zj3*zj2
662 +       xjhat = xj*rI
663 +       yjhat = yj*rI
664 +       zjhat = zj*rI
665 +      
666 +       call calc_sw_fnc(rij, rl, ru, rlp, rup, s, sp, dsdr, dspdr)
667 +          
668 +       frac1 = 0.25d0
669 +       frac2 = 0.75d0
670 +      
671 +       wi = 2.0d0*(xi2-yi2)*zi*rI3
672 +       wj = 2.0d0*(xj2-yj2)*zj*rI3
673 +      
674 +       wi2 = wi*wi
675 +       wj2 = wj*wj
676 +
677 +       w = frac1*wi*wi2 + frac2*wi + frac1*wj*wj2 + frac2*wj + v0p
678 +
679 +       vpair = vpair + 0.5d0*(v0*s*w)
680 +      
681 +       if (do_pot) then
682 + #ifdef IS_MPI
683 +         pot_row(atom1) = pot_row(atom1) + 0.25d0*(v0*s*w)*sw
684 +         pot_col(atom2) = pot_col(atom2) + 0.25d0*(v0*s*w)*sw
685 + #else
686 +         pot = pot + 0.5d0*(v0*s*w)*sw
687 + #endif  
688 +       endif
689 +
690 +       dwidx = ( 4.0d0*xi*zi*rI3 - 6.0d0*xi*zi*(xi2-yi2)*rI5 )
691 +       dwidy = ( -4.0d0*yi*zi*rI3 - 6.0d0*yi*zi*(xi2-yi2)*rI5 )
692 +       dwidz = ( 2.0d0*(xi2-yi2)*rI3 - 6.0d0*zi2*(xi2-yi2)*rI5 )
693 +      
694 +       dwidx = frac1*3.0d0*wi2*dwidx + frac2*dwidx
695 +       dwidy = frac1*3.0d0*wi2*dwidy + frac2*dwidy
696 +       dwidz = frac1*3.0d0*wi2*dwidz + frac2*dwidz
697 +
698 +       dwjdx = ( 4.0d0*xj*zj*rI3  - 6.0d0*xj*zj*(xj2-yj2)*rI5 )
699 +       dwjdy = ( -4.0d0*yj*zj*rI3  - 6.0d0*yj*zj*(xj2-yj2)*rI5 )
700 +       dwjdz = ( 2.0d0*(xj2-yj2)*rI3  - 6.0d0*zj2*(xj2-yj2)*rI5 )
701 +
702 +       dwjdx = frac1*3.0d0*wj2*dwjdx + frac2*dwjdx
703 +       dwjdy = frac1*3.0d0*wj2*dwjdy + frac2*dwjdy
704 +       dwjdz = frac1*3.0d0*wj2*dwjdz + frac2*dwjdz
705 +      
706 +       dwidux = ( 4.0d0*(yi*zi2 + 0.5d0*yi*(xi2-yi2))*rI3 )
707 +       dwiduy = ( 4.0d0*(xi*zi2 - 0.5d0*xi*(xi2-yi2))*rI3 )
708 +       dwiduz = ( -8.0d0*xi*yi*zi*rI3 )
709 +
710 +       dwidux = frac1*3.0d0*wi2*dwidux + frac2*dwidux
711 +       dwiduy = frac1*3.0d0*wi2*dwiduy + frac2*dwiduy
712 +       dwiduz = frac1*3.0d0*wi2*dwiduz + frac2*dwiduz
713 +
714 +       dwjdux = ( 4.0d0*(yj*zj2 + 0.5d0*yj*(xj2-yj2))*rI3 )
715 +       dwjduy = ( 4.0d0*(xj*zj2 - 0.5d0*xj*(xj2-yj2))*rI3 )
716 +       dwjduz = ( -8.0d0*xj*yj*zj*rI3 )
717 +
718 +       dwjdux = frac1*3.0d0*wj2*dwjdux + frac2*dwjdux
719 +       dwjduy = frac1*3.0d0*wj2*dwjduy + frac2*dwjduy
720 +       dwjduz = frac1*3.0d0*wj2*dwjduz + frac2*dwjduz
721 +
722 +       ! do the torques first since they are easy:
723 +       ! remember that these are still in the body fixed axes
724 +
725 +       txi = 0.5d0*(v0*s*dwidux)*sw
726 +       tyi = 0.5d0*(v0*s*dwiduy)*sw
727 +       tzi = 0.5d0*(v0*s*dwiduz)*sw
728 +
729 +       txj = 0.5d0*(v0*s*dwjdux)*sw
730 +       tyj = 0.5d0*(v0*s*dwjduy)*sw
731 +       tzj = 0.5d0*(v0*s*dwjduz)*sw
732 +
733 +       ! go back to lab frame using transpose of rotation matrix:
734 +
735 + #ifdef IS_MPI
736 +       t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
737 +            a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
738 +       t_Row(2,atom1) = t_Row(2,atom1) + a_Row(2,atom1)*txi + &
739 +            a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
740 +       t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
741 +            a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
742 +
743 +       t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
744 +            a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
745 +       t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
746 +            a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
747 +       t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
748 +            a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
749 + #else
750 +       t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
751 +       t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
752 +       t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
753 +
754 +       t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
755 +       t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
756 +       t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
757 + #endif    
758 +       ! Now, on to the forces:
759 +
760 +       ! first rotate the i terms back into the lab frame:
761 +
762 +       radcomxi = (v0*s*dwidx)*sw
763 +       radcomyi = (v0*s*dwidy)*sw
764 +       radcomzi = (v0*s*dwidz)*sw
765 +
766 +       radcomxj = (v0*s*dwjdx)*sw
767 +       radcomyj = (v0*s*dwjdy)*sw
768 +       radcomzj = (v0*s*dwjdz)*sw
769 +
770 + #ifdef IS_MPI    
771 +       fxii = a_Row(1,atom1)*(radcomxi) + &
772 +            a_Row(4,atom1)*(radcomyi) + &
773 +            a_Row(7,atom1)*(radcomzi)
774 +       fyii = a_Row(2,atom1)*(radcomxi) + &
775 +            a_Row(5,atom1)*(radcomyi) + &
776 +            a_Row(8,atom1)*(radcomzi)
777 +       fzii = a_Row(3,atom1)*(radcomxi) + &
778 +            a_Row(6,atom1)*(radcomyi) + &
779 +            a_Row(9,atom1)*(radcomzi)
780 +
781 +       fxjj = a_Col(1,atom2)*(radcomxj) + &
782 +            a_Col(4,atom2)*(radcomyj) + &
783 +            a_Col(7,atom2)*(radcomzj)
784 +       fyjj = a_Col(2,atom2)*(radcomxj) + &
785 +            a_Col(5,atom2)*(radcomyj) + &
786 +            a_Col(8,atom2)*(radcomzj)
787 +       fzjj = a_Col(3,atom2)*(radcomxj)+ &
788 +            a_Col(6,atom2)*(radcomyj) + &
789 +            a_Col(9,atom2)*(radcomzj)
790 + #else
791 +       fxii = a(1,atom1)*(radcomxi) + &
792 +            a(4,atom1)*(radcomyi) + &
793 +            a(7,atom1)*(radcomzi)
794 +       fyii = a(2,atom1)*(radcomxi) + &
795 +            a(5,atom1)*(radcomyi) + &
796 +            a(8,atom1)*(radcomzi)
797 +       fzii = a(3,atom1)*(radcomxi) + &
798 +            a(6,atom1)*(radcomyi) + &
799 +            a(9,atom1)*(radcomzi)
800 +
801 +       fxjj = a(1,atom2)*(radcomxj) + &
802 +            a(4,atom2)*(radcomyj) + &
803 +            a(7,atom2)*(radcomzj)
804 +       fyjj = a(2,atom2)*(radcomxj) + &
805 +            a(5,atom2)*(radcomyj) + &
806 +            a(8,atom2)*(radcomzj)
807 +       fzjj = a(3,atom2)*(radcomxj)+ &
808 +            a(6,atom2)*(radcomyj) + &
809 +            a(9,atom2)*(radcomzj)
810 + #endif
811 +
812 +       fxij = -fxii
813 +       fyij = -fyii
814 +       fzij = -fzii
815 +
816 +       fxji = -fxjj
817 +       fyji = -fyjj
818 +       fzji = -fzjj
819 +
820 +       ! now assemble these with the radial-only terms:
821 +
822 +       fxradial = 0.5d0*(v0*dsdr*w*drdx + fxii + fxji)
823 +       fyradial = 0.5d0*(v0*dsdr*w*drdy + fyii + fyji)
824 +       fzradial = 0.5d0*(v0*dsdr*w*drdz + fzii + fzji)
825 +
826 + #ifdef IS_MPI
827 +       f_Row(1,atom1) = f_Row(1,atom1) + fxradial
828 +       f_Row(2,atom1) = f_Row(2,atom1) + fyradial
829 +       f_Row(3,atom1) = f_Row(3,atom1) + fzradial
830 +
831 +       f_Col(1,atom2) = f_Col(1,atom2) - fxradial
832 +       f_Col(2,atom2) = f_Col(2,atom2) - fyradial
833 +       f_Col(3,atom2) = f_Col(3,atom2) - fzradial
834 + #else
835 +       f(1,atom1) = f(1,atom1) + fxradial
836 +       f(2,atom1) = f(2,atom1) + fyradial
837 +       f(3,atom1) = f(3,atom1) + fzradial
838 +
839 +       f(1,atom2) = f(1,atom2) - fxradial
840 +       f(2,atom2) = f(2,atom2) - fyradial
841 +       f(3,atom2) = f(3,atom2) - fzradial
842 + #endif
843 +
844 + #ifdef IS_MPI
845 +       id1 = AtomRowToGlobal(atom1)
846 +       id2 = AtomColToGlobal(atom2)
847 + #else
848 +       id1 = atom1
849 +       id2 = atom2
850 + #endif
851 +
852 +       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
853 +
854 +          fpair(1) = fpair(1) + fxradial
855 +          fpair(2) = fpair(2) + fyradial
856 +          fpair(3) = fpair(3) + fzradial
857 +
858 +       endif
859 +    endif
860 +  end subroutine do_sticky_power_pair
861 +
862   end module sticky

Diff Legend

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