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

Diff Legend

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