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 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 2361 by gezelter, Wed Oct 12 21:00:50 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.3 2005-01-12 22:40:45 gezelter Exp $, $Date: 2005-01-12 22:40:45 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $
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
508 – end module sticky
527  
528 <  subroutine newStickyType(c_ident, w0, v0, v0p, rl, ru, rlp, rup, isError)
529 <
530 <    use definitions, ONLY : dp  
531 <    use sticky, ONLY : module_newStickyType => newStickyType
532 <
533 <    integer, intent(inout) :: c_ident, isError
534 <    real( kind = dp ), intent(inout) :: w0, v0, v0p, rl, ru, rlp, rup
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 <    call module_newStickyType(c_ident, w0, v0, v0p, rl, ru, rlp, rup, &
519 <         isError)
537 >    !! i and j are pointers to the two SSD atoms
538      
539 <  end subroutine newStickyType
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