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

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90 (file contents):
Revision 2301 by gezelter, Thu Sep 15 22:05:21 2005 UTC vs.
Revision 2399 by chrisfen, Wed Oct 26 23:31:40 2005 UTC

# Line 54 | Line 54 | module electrostatic_module
54  
55    PRIVATE
56  
57 +
58   #define __FORTRAN90
59 + #include "UseTheForce/DarkSide/fInteractionMap.h"
60   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61  
62 +
63    !! these prefactors convert the multipole interactions into kcal / mol
64    !! all were computed assuming distances are measured in angstroms
65    !! Charge-Charge, assuming charges are measured in electrons
# Line 73 | Line 76 | module electrostatic_module
76  
77    !! variables to handle different summation methods for long-range electrostatics:
78    integer, save :: summationMethod = NONE
79 +  logical, save :: summationMethodChecked = .false.
80    real(kind=DP), save :: defaultCutoff = 0.0_DP
81 +  real(kind=DP), save :: defaultCutoff2 = 0.0_DP
82    logical, save :: haveDefaultCutoff = .false.
83    real(kind=DP), save :: dampingAlpha = 0.0_DP
84    logical, save :: haveDampingAlpha = .false.
85 <  real(kind=DP), save :: dielectric = 0.0_DP
85 >  real(kind=DP), save :: dielectric = 1.0_DP
86    logical, save :: haveDielectric = .false.
87    real(kind=DP), save :: constERFC = 0.0_DP
88    real(kind=DP), save :: constEXP = 0.0_DP
89    logical, save :: haveDWAconstants = .false.
90 <
91 <
90 >  real(kind=dp), save :: rcuti = 0.0_DP
91 >  real(kind=dp), save :: rcuti2 = 0.0_DP
92 >  real(kind=dp), save :: rcuti3 = 0.0_DP
93 >  real(kind=dp), save :: rcuti4 = 0.0_DP
94 >  real(kind=dp), save :: alphaPi = 0.0_DP
95 >  real(kind=dp), save :: invRootPi = 0.0_DP
96 >  real(kind=dp), save :: rrf = 1.0_DP
97 >  real(kind=dp), save :: rt = 1.0_DP
98 >  real(kind=dp), save :: rrfsq = 1.0_DP
99 >  real(kind=dp), save :: preRF = 0.0_DP
100 >  real(kind=dp), save :: preRF2 = 0.0_DP
101 >  logical, save :: preRFCalculated = .false.
102 >
103 > #ifdef __IFC
104 > ! error function for ifc version > 7.
105 >  double precision, external :: derfc
106 > #endif
107 >  
108    public :: setElectrostaticSummationMethod
109    public :: setElectrostaticCutoffRadius
110    public :: setDampedWolfAlpha
111    public :: setReactionFieldDielectric
112 +  public :: setReactionFieldPrefactor
113    public :: newElectrostaticType
114    public :: setCharge
115    public :: setDipoleMoment
# Line 96 | Line 118 | module electrostatic_module
118    public :: doElectrostaticPair
119    public :: getCharge
120    public :: getDipoleMoment
99 –  public :: pre22
121    public :: destroyElectrostaticTypes
122 +  public :: rf_self_self
123 +  public :: rf_self_excludes
124  
125    type :: Electrostatic
126       integer :: c_ident
# Line 117 | Line 140 | contains
140   contains
141  
142    subroutine setElectrostaticSummationMethod(the_ESM)
120 –
143      integer, intent(in) :: the_ESM    
144  
145      if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
146         call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
147      endif
148  
149 +    summationMethod = the_ESM
150 +
151    end subroutine setElectrostaticSummationMethod
152  
153 <  subroutine setElectrostaticCutoffRadius(thisRcut)
153 >  subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
154      real(kind=dp), intent(in) :: thisRcut
155 +    real(kind=dp), intent(in) :: thisRsw
156      defaultCutoff = thisRcut
157 +    rrf = defaultCutoff
158 +    rt = thisRsw
159      haveDefaultCutoff = .true.
160    end subroutine setElectrostaticCutoffRadius
161  
# Line 144 | Line 171 | contains
171      haveDielectric = .true.
172    end subroutine setReactionFieldDielectric
173  
174 +  subroutine setReactionFieldPrefactor
175 +    if (haveDefaultCutoff .and. haveDielectric) then
176 +       defaultCutoff2 = defaultCutoff*defaultCutoff
177 +       preRF = (dielectric-1.0d0) / &
178 +            ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
179 +       preRF2 = 2.0d0*preRF
180 +       preRFCalculated = .true.
181 +    else if (.not.haveDefaultCutoff) then
182 +       call handleError("setReactionFieldPrefactor", "Default cutoff not set")
183 +    else
184 +       call handleError("setReactionFieldPrefactor", "Dielectric not set")
185 +    endif
186 +  end subroutine setReactionFieldPrefactor
187 +
188    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
189         is_SplitDipole, is_Quadrupole, is_Tap, status)
190  
# Line 355 | Line 396 | contains
396    end function getDipoleMoment
397  
398    subroutine checkSummationMethod()
399 +
400 +    if (.not.haveDefaultCutoff) then
401 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
402 +    endif
403 +
404 +    rcuti = 1.0d0 / defaultCutoff
405 +    rcuti2 = rcuti*rcuti
406 +    rcuti3 = rcuti2*rcuti
407 +    rcuti4 = rcuti2*rcuti2
408  
409      if (summationMethod .eq. DAMPED_WOLF) then
410         if (.not.haveDWAconstants) then
# Line 363 | Line 413 | contains
413               call handleError("checkSummationMethod", "no Damping Alpha set!")
414            endif
415            
416 <          if (.not.have....)
417 <          constEXP =
418 <          constERFC =
419 <          
416 >          if (.not.haveDefaultCutoff) then
417 >             call handleError("checkSummationMethod", "no Default Cutoff set!")
418 >          endif
419 >
420 >          constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
421 >          constERFC = derfc(dampingAlpha*defaultCutoff)
422 >          invRootPi = 0.56418958354775628695d0
423 >          alphaPi = 2*dampingAlpha*invRootPi
424 >  
425            haveDWAconstants = .true.
426         endif
427      endif
428  
429 +    if (summationMethod .eq. REACTION_FIELD) then
430 +       if (.not.haveDielectric) then
431 +          call handleError("checkSummationMethod", "no reaction field Dielectric set!")
432 +       endif
433 +    endif
434 +
435 +    summationMethodChecked = .true.
436    end subroutine checkSummationMethod
437  
438  
439  
440    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
441 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
441 >       vpair, fpair, pot, eFrame, f, t, do_pot, indirect_only)
442  
443 <    logical, intent(in) :: do_pot
443 >    logical, intent(in) :: do_pot, indirect_only
444  
445      integer, intent(in) :: atom1, atom2
446      integer :: localError
385 –    integer, intent(in) :: corrMethod
447  
448 <    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
448 >    real(kind=dp), intent(in) :: rij, r2, sw
449      real(kind=dp), intent(in), dimension(3) :: d
450      real(kind=dp), intent(inout) :: vpair
451      real(kind=dp), intent(inout), dimension(3) :: fpair
452  
453 <    real( kind = dp ) :: pot, swi
453 >    real( kind = dp ) :: pot
454      real( kind = dp ), dimension(9,nLocal) :: eFrame
455      real( kind = dp ), dimension(3,nLocal) :: f
456      real( kind = dp ), dimension(3,nLocal) :: t
# Line 414 | Line 475 | contains
475      real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
476      real (kind=dp) :: xhat, yhat, zhat
477      real (kind=dp) :: dudx, dudy, dudz
478 <    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
479 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
478 >    real (kind=dp) :: scale, sc2, bigR
479 >    real (kind=dp) :: varERFC, varEXP
480 >    real (kind=dp) :: limScale
481 >    real (kind=dp) :: preVal, rfVal
482  
483      if (.not.allocated(ElectrostaticMap)) then
484         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 426 | Line 489 | contains
489         call checkSummationMethod()
490      endif
491  
492 +    if (.not.preRFCalculated) then
493 +       call setReactionFieldPrefactor()
494 +    endif
495  
496   #ifdef IS_MPI
497      me1 = atid_Row(atom1)
# Line 438 | Line 504 | contains
504      !! some variables we'll need independent of electrostatic type:
505  
506      riji = 1.0d0 / rij
507 <
507 >  
508      xhat = d(1) * riji
509      yhat = d(2) * riji
510      zhat = d(3) * riji
511  
446 –    rcuti2 = rcuti*rcuti
447 –    rcuti3 = rcuti2*rcuti
448 –    rcuti4 = rcuti2*rcuti2
449 –
450 –    swi = 1.0d0 / sw
451 –
512      !! logicals
513      i_is_Charge = ElectrostaticMap(me1)%is_Charge
514      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 567 | Line 627 | contains
627         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
628      endif
629    
570 – !!$    switcher = 1.0d0
571 – !!$    dswitcher = 0.0d0
572 – !!$    ebalance = 0.0d0
573 – !!$    ! weaken the dipole interaction at close range for TAP water
574 – !!$    if (j_is_Tap .and. i_is_Tap) then
575 – !!$      call calc_switch(rij, mu_i, switcher, dswitcher)
576 – !!$    endif
577 –
630      epot = 0.0_dp
631      dudx = 0.0_dp
632      dudy = 0.0_dp
# Line 592 | Line 644 | contains
644  
645         if (j_is_Charge) then
646  
647 <          if (corrMethod .eq. 1) then
647 >          if (summationMethod .eq. UNDAMPED_WOLF) then
648               vterm = pre11 * q_i * q_j * (riji - rcuti)
649 +             vpair = vpair + vterm
650 +             epot = epot + sw*vterm
651 +            
652 +             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
653 +            
654 +             dudx = dudx + dudr * d(1)
655 +             dudy = dudy + dudr * d(2)
656 +             dudz = dudz + dudr * d(3)
657  
658 +          elseif (summationMethod .eq. DAMPED_WOLF) then
659 +             varERFC = derfc(dampingAlpha*rij)
660 +             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
661 +             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
662               vpair = vpair + vterm
663 <             epot = epot + sw * vterm
663 >             epot = epot + sw*vterm
664              
665 <             dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
665 >             dudr  = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
666 >                                                  + alphaPi*varEXP) &
667 >                                                 - (constERFC*rcuti2 &
668 >                                                    + alphaPi*constEXP)) )
669              
670               dudx = dudx + dudr * d(1)
671               dudy = dudy + dudr * d(2)
672               dudz = dudz + dudr * d(3)
673  
674 +          elseif (summationMethod .eq. REACTION_FIELD) then
675 +             preVal = pre11 * q_i * q_j
676 +             rfVal = preRF*rij*rij
677 +             vterm = preVal * ( riji + rfVal )
678 +            
679 +             vpair = vpair + vterm
680 +             epot = epot + sw*vterm
681 +            
682 +             dudr  = sw * preVal * ( 2.0d0*rfVal - riji )*riji
683 +            
684 +             dudx = dudx + dudr * xhat
685 +             dudy = dudy + dudr * yhat
686 +             dudz = dudz + dudr * zhat
687 +
688            else
689               vterm = pre11 * q_i * q_j * riji
609 –
690               vpair = vpair + vterm
691 <             epot = epot + sw * vterm
691 >             epot = epot + sw*vterm
692              
693               dudr  = - sw * vterm * riji
694              
# Line 622 | Line 702 | contains
702  
703         if (j_is_Dipole) then
704  
705 <          pref = sw * pre12 * q_i * mu_j
705 >          pref = pre12 * q_i * mu_j
706  
707 <          if (corrMethod .eq. 1) then
707 >          if (summationMethod .eq. UNDAMPED_WOLF) then
708               ri2 = riji * riji
709               ri3 = ri2 * riji
710  
711 +             pref = pre12 * q_i * mu_j
712               vterm = - pref * ct_j * (ri2 - rcuti2)
713 <             vpair = vpair + swi*vterm
714 <             epot = epot + vterm
713 >             vpair = vpair + vterm
714 >             epot = epot + sw*vterm
715              
716               !! this has a + sign in the () because the rij vector is
717               !! r_j - r_i and the charge-dipole potential takes the origin
718               !! as the point dipole, which is atom j in this case.
719              
720 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
720 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
721                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
722 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
722 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
723                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
724 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
724 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
725                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
726              
727 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
728 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
729 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
727 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
728 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
729 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
730  
731 +          elseif (summationMethod .eq. REACTION_FIELD) then
732 +             ri2 = riji * riji
733 +             ri3 = ri2 * riji
734 +    
735 +             pref = pre12 * q_i * mu_j
736 +             vterm = - pref * ct_j * ( ri2 - preRF2*rij )
737 +             vpair = vpair + vterm
738 +             epot = epot + sw*vterm
739 +            
740 +             !! this has a + sign in the () because the rij vector is
741 +             !! r_j - r_i and the charge-dipole potential takes the origin
742 +             !! as the point dipole, which is atom j in this case.
743 +            
744 +             dudx = dudx - sw*pref*( ri3*(uz_j(1) - 3.0d0*ct_j*xhat) - &
745 +                                     preRF2*uz_j(1) )
746 +             dudy = dudy - sw*pref*( ri3*(uz_j(2) - 3.0d0*ct_j*yhat) - &
747 +                                     preRF2*uz_j(2) )
748 +             dudz = dudz - sw*pref*( ri3*(uz_j(3) - 3.0d0*ct_j*zhat) - &
749 +                                     preRF2*uz_j(3) )        
750 +             duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
751 +             duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
752 +             duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
753 +
754            else
755               if (j_is_SplitDipole) then
756                  BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
# Line 660 | Line 764 | contains
764               ri2 = ri * ri
765               ri3 = ri2 * ri
766               sc2 = scale * scale
767 <            
767 >
768 >             pref = pre12 * q_i * mu_j
769               vterm = - pref * ct_j * ri2 * scale
770 <             vpair = vpair + swi * vterm
771 <             epot = epot + vterm
770 >             vpair = vpair + vterm
771 >             epot = epot + sw*vterm
772              
773               !! this has a + sign in the () because the rij vector is
774               !! r_j - r_i and the charge-dipole potential takes the origin
775               !! as the point dipole, which is atom j in this case.
776              
777 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
778 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
779 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
777 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
778 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
779 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
780              
781 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
782 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
783 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
781 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
782 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
783 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
784  
785            endif
786         endif
# Line 688 | Line 793 | contains
793            cy2 = cy_j * cy_j
794            cz2 = cz_j * cz_j
795  
796 <
797 <          pref =  sw * pre14 * q_i / 3.0_dp
693 <
694 <          if (corrMethod .eq. 1) then
796 >          if (summationMethod .eq. UNDAMPED_WOLF) then
797 >             pref =  pre14 * q_i / 3.0_dp
798               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
799                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
800                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
801               vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
802                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
803                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
804 <             vpair = vpair + swi*( vterm1 - vterm2 )
805 <             epot = epot + ( vterm1 - vterm2 )
804 >             vpair = vpair + ( vterm1 - vterm2 )
805 >             epot = epot + sw*( vterm1 - vterm2 )
806              
807               dudx = dudx - (5.0_dp * &
808 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
808 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
809                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
810                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
811                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 710 | Line 813 | contains
813                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
814                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
815               dudy = dudy - (5.0_dp * &
816 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
816 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
817                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
818                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
819                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 718 | Line 821 | contains
821                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
822                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
823               dudz = dudz - (5.0_dp * &
824 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
824 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
825                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
826                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
827                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 726 | Line 829 | contains
829                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
830                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
831              
832 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
832 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
833                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
834 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
834 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
835                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
836 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
836 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
837                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
838              
839 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
839 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
840                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
841 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
841 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
842                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
843 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
843 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
844                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
845              
846 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
846 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
847                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
848 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
848 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
849                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
850 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
850 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
851                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
852          
853            else
854 +             pref =  pre14 * q_i / 3.0_dp
855               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
856                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
857                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
858 <             vpair = vpair + swi * vterm
859 <             epot = epot + vterm
858 >             vpair = vpair + vterm
859 >             epot = epot + sw*vterm
860              
861 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
861 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
862                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
863                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
864                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
865 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
865 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
866                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
867                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
868                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
869 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
869 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
870                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
871                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
872                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
873              
874 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
875 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
876 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
874 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
875 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
876 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
877              
878 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
879 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
880 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
878 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
879 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
880 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
881              
882 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
883 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
884 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
882 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
883 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
884 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
885            
886            endif
887         endif
# Line 786 | Line 890 | contains
890      if (i_is_Dipole) then
891  
892         if (j_is_Charge) then
893 <
894 <          pref = sw * pre12 * q_j * mu_i
895 <
896 <          if (corrMethod .eq. 1) then
893 >          
894 >          pref = pre12 * q_j * mu_i
895 >          
896 >          if (summationMethod .eq. UNDAMPED_WOLF) then
897               ri2 = riji * riji
898               ri3 = ri2 * riji
899  
900 +             pref = pre12 * q_j * mu_i
901               vterm = pref * ct_i * (ri2 - rcuti2)
902 <             vpair = vpair + swi * vterm
903 <             epot = epot + vterm
902 >             vpair = vpair + vterm
903 >             epot = epot + sw*vterm
904              
905 <             !! this has a + sign in the () because the rij vector is
801 <             !! r_j - r_i and the charge-dipole potential takes the origin
802 <             !! as the point dipole, which is atom j in this case.
803 <            
804 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
905 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
906                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
907 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
907 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
908                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
909 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
909 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
910                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
911              
912 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
913 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
914 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
912 >             duduz_i(1) = duduz_i(1) + sw*pref*( ri2*xhat - d(1)*rcuti3 )
913 >             duduz_i(2) = duduz_i(2) + sw*pref*( ri2*yhat - d(2)*rcuti3 )
914 >             duduz_i(3) = duduz_i(3) + sw*pref*( ri2*zhat - d(3)*rcuti3 )
915  
916 +          elseif (summationMethod .eq. REACTION_FIELD) then
917 +             ri2 = riji * riji
918 +             ri3 = ri2 * riji
919 +
920 +             pref = pre12 * q_j * mu_i
921 +             vterm = pref * ct_i * ( ri2 - preRF2*rij )
922 +             vpair = vpair + vterm
923 +             epot = epot + sw*vterm
924 +            
925 +             dudx = dudx + sw*pref * ( ri3*(uz_i(1) - 3.0d0*ct_i*xhat) - &
926 +                  preRF2*uz_i(1) )
927 +             dudy = dudy + sw*pref * ( ri3*(uz_i(2) - 3.0d0*ct_i*yhat) - &
928 +                  preRF2*uz_i(2) )
929 +             dudz = dudz + sw*pref * ( ri3*(uz_i(3) - 3.0d0*ct_i*zhat) - &
930 +                  preRF2*uz_i(3) )
931 +            
932 +             duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF2*rij )
933 +             duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF2*rij )
934 +             duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF2*rij )
935 +
936            else
937               if (i_is_SplitDipole) then
938                  BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
# Line 825 | Line 946 | contains
946               ri2 = ri * ri
947               ri3 = ri2 * ri
948               sc2 = scale * scale
949 <            
949 >
950 >             pref = pre12 * q_j * mu_i
951               vterm = pref * ct_i * ri2 * scale
952 <             vpair = vpair + swi * vterm
953 <             epot = epot + vterm
952 >             vpair = vpair + vterm
953 >             epot = epot + sw*vterm
954              
955 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
956 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
957 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
955 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
956 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
957 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
958              
959 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
960 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
961 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
959 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
960 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
961 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
962            endif
963         endif
964 <
964 >      
965         if (j_is_Dipole) then
966  
967 <          pref = sw * pre22 * mu_i * mu_j
846 <
847 <          if (corrMethod .eq. 1) then
967 >          if (summationMethod .eq. UNDAMPED_WOLF) then
968               ri2 = riji * riji
969               ri3 = ri2 * riji
970               ri4 = ri2 * ri2
971  
972 +             pref = pre22 * mu_i * mu_j
973               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
974 <             vpair = vpair + swi * vterm
975 <             epot = epot + vterm
974 >             vpair = vpair + vterm
975 >             epot = epot + sw*vterm
976              
977               a1 = 5.0d0 * ct_i * ct_j - ct_ij
978              
979 <             dudx = dudx + pref*3.0d0*ri4 &
980 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
981 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
982 <             dudy = dudy + pref*3.0d0*ri4 &
983 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
984 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
985 <             dudz = dudz + pref*3.0d0*ri4 &
986 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
987 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
979 >             dudx = dudx + sw*pref*3.0d0*ri4 &
980 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
981 >                         - sw*pref*3.0d0*rcuti4 &
982 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
983 >             dudy = dudy + sw*pref*3.0d0*ri4 &
984 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
985 >                         - sw*pref*3.0d0*rcuti4 &
986 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
987 >             dudz = dudz + sw*pref*3.0d0*ri4 &
988 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
989 >                         - sw*pref*3.0d0*rcuti4 &
990 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
991              
992 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
992 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
993                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
994 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
994 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
995                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
996 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
996 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
997                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
998 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
998 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
999                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
1000 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
1000 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1001                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
1002 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
1002 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1003                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
1004 <          else
1004 >
1005 >         elseif (summationMethod .eq. REACTION_FIELD) then
1006 >             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1007 >
1008 >             ri2 = riji * riji
1009 >             ri3 = ri2 * riji
1010 >             ri4 = ri2 * ri2
1011 >
1012 >             pref = pre22 * mu_i * mu_j
1013 >              
1014 >             vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
1015 >                  preRF2*ct_ij )
1016 >             vpair = vpair + vterm
1017 >             epot = epot + sw*vterm
1018              
1019 +             a1 = 5.0d0 * ct_i * ct_j - ct_ij
1020 +            
1021 +             dudx = dudx + sw*pref*3.0d0*ri4 &
1022 +                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1023 +             dudy = dudy + sw*pref*3.0d0*ri4 &
1024 +                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1025 +             dudz = dudz + sw*pref*3.0d0*ri4 &
1026 +                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1027 +            
1028 +             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1029 +                  - preRF2*uz_j(1))
1030 +             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1031 +                  - preRF2*uz_j(2))
1032 +             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1033 +                  - preRF2*uz_j(3))
1034 +             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1035 +                  - preRF2*uz_i(1))
1036 +             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1037 +                  - preRF2*uz_i(2))
1038 +             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1039 +                  - preRF2*uz_i(3))
1040 +
1041 +          else
1042               if (i_is_SplitDipole) then
1043                  if (j_is_SplitDipole) then
1044                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 905 | Line 1065 | contains
1065               ri4 = ri2 * ri2
1066               sc2 = scale * scale
1067              
1068 +             pref = pre22 * mu_i * mu_j
1069               vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1070 <             vpair = vpair + swi * vterm
1071 <             epot = epot + vterm
1070 >             vpair = vpair + vterm
1071 >             epot = epot + sw*vterm
1072              
1073               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1074              
1075 <             dudx = dudx + pref*3.0d0*ri4*scale &
1076 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1077 <             dudy = dudy + pref*3.0d0*ri4*scale &
1078 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1079 <             dudz = dudz + pref*3.0d0*ri4*scale &
1080 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1075 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
1076 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1077 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
1078 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1079 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
1080 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1081              
1082 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
1083 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1084 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
1085 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1086 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
1087 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1082 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
1083 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1084 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
1085 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1086 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
1087 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1088              
1089 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
1090 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1091 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
1092 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1093 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
1094 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1089 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1090 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1091 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1092 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1093 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1094 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1095            endif
1096         endif
1097      endif
# Line 945 | Line 1106 | contains
1106            cy2 = cy_i * cy_i
1107            cz2 = cz_i * cz_i
1108  
1109 <          pref = sw * pre14 * q_j / 3.0_dp
1110 <
950 <          if (corrMethod .eq. 1) then
1109 >          if (summationMethod .eq. UNDAMPED_WOLF) then
1110 >             pref = pre14 * q_j / 3.0_dp
1111               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1112                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1113                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1114               vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1115                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1116                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1117 <             vpair = vpair + swi * ( vterm1 - vterm2 )
1118 <             epot = epot + ( vterm1 - vterm2 )
1117 >             vpair = vpair + ( vterm1 - vterm2 )
1118 >             epot = epot + sw*( vterm1 - vterm2 )
1119              
1120 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
1121 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1120 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1121 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1122                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1123                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1124                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1125                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1126                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1127 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
1128 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1127 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1128 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1129                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1130                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1131                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1132                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1133                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1134 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
1135 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1134 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1135 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1136                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1137                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1138                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1139                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1140                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1141              
1142 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1142 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1143                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1144 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1144 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1145                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1146 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1146 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1147                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1148              
1149 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1149 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1150                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1151 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1151 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1152                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1153 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1153 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1154                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1155              
1156 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1156 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1157                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1158 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1158 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1159                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1160 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1160 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1161                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1162  
1163            else
1164 +             pref = pre14 * q_j / 3.0_dp
1165               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1166                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1167                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1168 <             vpair = vpair + swi * vterm
1169 <             epot = epot + vterm
1168 >             vpair = vpair + vterm
1169 >             epot = epot + sw*vterm
1170              
1171 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1171 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1172                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1173                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1174                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1175 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1175 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1176                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1177                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1178                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1179 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1179 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1180                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1181                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1182                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1183              
1184 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1185 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1186 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1184 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1185 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1186 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1187              
1188 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1189 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1190 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1188 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1189 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1190 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1191              
1192 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1193 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1194 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1192 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1193 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1194 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1195            endif
1196         endif
1197      endif
# Line 1038 | Line 1199 | contains
1199  
1200      if (do_pot) then
1201   #ifdef IS_MPI
1202 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1203 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1202 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1203 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1204   #else
1205         pot = pot + epot
1206   #endif
# Line 1143 | Line 1304 | contains
1304  
1305      return
1306    end subroutine doElectrostaticPair
1146 –
1147 –  !! calculates the switching functions and their derivatives for a given
1148 –  subroutine calc_switch(r, mu, scale, dscale)
1149 –
1150 –    real (kind=dp), intent(in) :: r, mu
1151 –    real (kind=dp), intent(inout) :: scale, dscale
1152 –    real (kind=dp) :: rl, ru, mulow, minRatio, temp, scaleVal
1153 –
1154 –    ! distances must be in angstroms
1155 –    rl = 2.75d0
1156 –    ru = 3.75d0
1157 –    mulow = 0.0d0 !3.3856d0 ! 1.84 * 1.84
1158 –    minRatio = mulow / (mu*mu)
1159 –    scaleVal = 1.0d0 - minRatio
1160 –    
1161 –    if (r.lt.rl) then
1162 –       scale = minRatio
1163 –       dscale = 0.0d0
1164 –    elseif (r.gt.ru) then
1165 –       scale = 1.0d0
1166 –       dscale = 0.0d0
1167 –    else
1168 –       scale = 1.0d0 - scaleVal*((ru + 2.0d0*r - 3.0d0*rl) * (ru-r)**2) &
1169 –                        / ((ru - rl)**3)
1170 –       dscale = -scaleVal * 6.0d0 * (r-ru)*(r-rl)/((ru - rl)**3)    
1171 –    endif
1172 –        
1173 –    return
1174 –  end subroutine calc_switch
1307  
1308    subroutine destroyElectrostaticTypes()
1309  
# Line 1179 | Line 1311 | end module electrostatic_module
1311  
1312    end subroutine destroyElectrostaticTypes
1313  
1314 +  subroutine rf_self_self(atom1, eFrame, rfpot, t, do_pot)
1315 +    logical, intent(in) :: do_pot
1316 +    integer, intent(in) :: atom1
1317 +    integer :: atid1
1318 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1319 +    real(kind=dp), dimension(3,nLocal) :: t
1320 +    real(kind=dp) :: mu1
1321 +    real(kind=dp) :: preVal, epot, rfpot
1322 +    real(kind=dp) :: eix, eiy, eiz
1323 +
1324 +    if (.not.preRFCalculated) then
1325 +       call setReactionFieldPrefactor()
1326 +    endif
1327 +
1328 +    ! this is a local only array, so we use the local atom type id's:
1329 +    atid1 = atid(atom1)
1330 +    
1331 +    if (ElectrostaticMap(atid1)%is_Dipole) then
1332 +       mu1 = getDipoleMoment(atid1)
1333 +      
1334 +       preVal = pre22 * preRF2 * mu1*mu1
1335 +       rfpot = rfpot - 0.5d0*preVal
1336 +      
1337 +       ! The self-correction term adds into the reaction field vector
1338 +      
1339 +       eix = preVal * eFrame(3,atom1)
1340 +       eiy = preVal * eFrame(6,atom1)
1341 +       eiz = preVal * eFrame(9,atom1)
1342 +      
1343 +       ! once again, this is self-self, so only the local arrays are needed
1344 +       ! even for MPI jobs:
1345 +      
1346 +       t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1347 +            eFrame(9,atom1)*eiy
1348 +       t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1349 +            eFrame(3,atom1)*eiz
1350 +       t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1351 +            eFrame(6,atom1)*eix
1352 +      
1353 +    endif
1354 +    
1355 +    return
1356 +  end subroutine rf_self_self
1357 +
1358 +  subroutine rf_self_excludes(atom1, atom2, sw, eFrame, d, rij, vpair, rfpot, &
1359 +       f, t, do_pot)
1360 +    logical, intent(in) :: do_pot
1361 +    integer, intent(in) :: atom1
1362 +    integer, intent(in) :: atom2
1363 +    logical :: i_is_Charge, j_is_Charge
1364 +    logical :: i_is_Dipole, j_is_Dipole
1365 +    integer :: atid1
1366 +    integer :: atid2
1367 +    real(kind=dp), intent(in) :: rij
1368 +    real(kind=dp), intent(in) :: sw
1369 +    real(kind=dp), intent(in), dimension(3) :: d
1370 +    real(kind=dp), intent(inout) :: vpair
1371 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1372 +    real(kind=dp), dimension(3,nLocal) :: f
1373 +    real(kind=dp), dimension(3,nLocal) :: t
1374 +    real (kind = dp), dimension(3) :: duduz_i
1375 +    real (kind = dp), dimension(3) :: duduz_j
1376 +    real (kind = dp), dimension(3) :: uz_i
1377 +    real (kind = dp), dimension(3) :: uz_j
1378 +    real(kind=dp) :: q_i, q_j, mu_i, mu_j
1379 +    real(kind=dp) :: xhat, yhat, zhat
1380 +    real(kind=dp) :: ct_i, ct_j
1381 +    real(kind=dp) :: ri2, ri3, riji, vterm
1382 +    real(kind=dp) :: pref, preVal, rfVal, rfpot
1383 +    real(kind=dp) :: dudx, dudy, dudz, dudr
1384 +
1385 +    if (.not.preRFCalculated) then
1386 +       call setReactionFieldPrefactor()
1387 +    endif
1388 +
1389 +    dudx = 0.0d0
1390 +    dudy = 0.0d0
1391 +    dudz = 0.0d0
1392 +
1393 +    riji = 1.0d0/rij
1394 +
1395 +    xhat = d(1) * riji
1396 +    yhat = d(2) * riji
1397 +    zhat = d(3) * riji
1398 +
1399 +    ! this is a local only array, so we use the local atom type id's:
1400 +    atid1 = atid(atom1)
1401 +    atid2 = atid(atom2)
1402 +    i_is_Charge = ElectrostaticMap(atid1)%is_Charge
1403 +    j_is_Charge = ElectrostaticMap(atid2)%is_Charge
1404 +    i_is_Dipole = ElectrostaticMap(atid1)%is_Dipole
1405 +    j_is_Dipole = ElectrostaticMap(atid2)%is_Dipole
1406 +
1407 +    if (i_is_Charge.and.j_is_Charge) then
1408 +       q_i = ElectrostaticMap(atid1)%charge
1409 +       q_j = ElectrostaticMap(atid2)%charge
1410 +      
1411 +       preVal = pre11 * q_i * q_j
1412 +       rfVal = preRF*rij*rij
1413 +       vterm = preVal * rfVal
1414 +
1415 +       rfpot = rfpot + sw*vterm
1416 +      
1417 +       dudr  = sw*preVal * 2.0d0*rfVal*riji
1418 +              
1419 +       dudx = dudx + dudr * xhat
1420 +       dudy = dudy + dudr * yhat
1421 +       dudz = dudz + dudr * zhat
1422 +
1423 +    elseif (i_is_Charge.and.j_is_Dipole) then
1424 +       q_i = ElectrostaticMap(atid1)%charge
1425 +       mu_j = ElectrostaticMap(atid2)%dipole_moment
1426 +       uz_j(1) = eFrame(3,atom2)
1427 +       uz_j(2) = eFrame(6,atom2)
1428 +       uz_j(3) = eFrame(9,atom2)
1429 +       ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
1430 +
1431 +       ri2 = riji * riji
1432 +       ri3 = ri2 * riji
1433 +      
1434 +       pref = pre12 * q_i * mu_j
1435 +       vterm = - pref * ct_j * ( ri2 - preRF2*rij )
1436 +       rfpot = rfpot + sw*vterm
1437 +
1438 +       dudx = dudx - sw*pref*( ri3*(uz_j(1)-3.0d0*ct_j*xhat) - preRF2*uz_j(1) )
1439 +       dudy = dudy - sw*pref*( ri3*(uz_j(2)-3.0d0*ct_j*yhat) - preRF2*uz_j(2) )
1440 +       dudz = dudz - sw*pref*( ri3*(uz_j(3)-3.0d0*ct_j*zhat) - preRF2*uz_j(3) )
1441 +
1442 +       duduz_j(1) = duduz_j(1) - sw * pref * xhat * ( ri2 - preRF2*rij )
1443 +       duduz_j(2) = duduz_j(2) - sw * pref * yhat * ( ri2 - preRF2*rij )
1444 +       duduz_j(3) = duduz_j(3) - sw * pref * zhat * ( ri2 - preRF2*rij )
1445 +              
1446 +    elseif (i_is_Dipole.and.j_is_Charge) then
1447 +       mu_i = ElectrostaticMap(atid1)%dipole_moment
1448 +       q_j = ElectrostaticMap(atid2)%charge
1449 +       uz_i(1) = eFrame(3,atom1)
1450 +       uz_i(2) = eFrame(6,atom1)
1451 +       uz_i(3) = eFrame(9,atom1)
1452 +       ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
1453 +
1454 +       ri2 = riji * riji
1455 +       ri3 = ri2 * riji
1456 +      
1457 +       pref = pre12 * q_j * mu_i
1458 +       vterm = pref * ct_i * ( ri2 - preRF2*rij )
1459 +       rfpot = rfpot + sw*vterm
1460 +      
1461 +       dudx = dudx + sw*pref*( ri3*(uz_i(1)-3.0d0*ct_i*xhat) - preRF2*uz_i(1) )
1462 +       dudy = dudy + sw*pref*( ri3*(uz_i(2)-3.0d0*ct_i*yhat) - preRF2*uz_i(2) )
1463 +       dudz = dudz + sw*pref*( ri3*(uz_i(3)-3.0d0*ct_i*zhat) - preRF2*uz_i(3) )
1464 +      
1465 +       duduz_i(1) = duduz_i(1) + sw * pref * xhat * ( ri2 - preRF2*rij )
1466 +       duduz_i(2) = duduz_i(2) + sw * pref * yhat * ( ri2 - preRF2*rij )
1467 +       duduz_i(3) = duduz_i(3) + sw * pref * zhat * ( ri2 - preRF2*rij )
1468 +      
1469 +    endif
1470 +    
1471 +    ! accumulate the forces and torques resulting from the RF self term
1472 +    f(1,atom1) = f(1,atom1) + dudx
1473 +    f(2,atom1) = f(2,atom1) + dudy
1474 +    f(3,atom1) = f(3,atom1) + dudz
1475 +    
1476 +    f(1,atom2) = f(1,atom2) - dudx
1477 +    f(2,atom2) = f(2,atom2) - dudy
1478 +    f(3,atom2) = f(3,atom2) - dudz
1479 +    
1480 +    if (i_is_Dipole) then
1481 +       t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1482 +       t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1483 +       t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1484 +    elseif (j_is_Dipole) then
1485 +       t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1486 +       t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1487 +       t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1488 +    endif
1489 +
1490 +    return
1491 +  end subroutine rf_self_excludes
1492 +
1493   end module electrostatic_module

Diff Legend

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