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 2395 by chrisfen, Mon Oct 24 14:06:36 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  
124    type :: Electrostatic
125       integer :: c_ident
# Line 117 | Line 139 | contains
139   contains
140  
141    subroutine setElectrostaticSummationMethod(the_ESM)
120 –
142      integer, intent(in) :: the_ESM    
143  
144      if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
145         call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
146      endif
147  
148 +    summationMethod = the_ESM
149 +
150    end subroutine setElectrostaticSummationMethod
151  
152 <  subroutine setElectrostaticCutoffRadius(thisRcut)
152 >  subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
153      real(kind=dp), intent(in) :: thisRcut
154 +    real(kind=dp), intent(in) :: thisRsw
155      defaultCutoff = thisRcut
156 +    rrf = defaultCutoff
157 +    rt = thisRsw
158      haveDefaultCutoff = .true.
159    end subroutine setElectrostaticCutoffRadius
160  
# Line 143 | Line 169 | contains
169      dielectric = thisDielectric
170      haveDielectric = .true.
171    end subroutine setReactionFieldDielectric
172 +
173 +  subroutine setReactionFieldPrefactor
174 +    if (haveDefaultCutoff .and. haveDielectric) then
175 +       defaultCutoff2 = defaultCutoff*defaultCutoff
176 +       preRF = (dielectric-1.0d0) / &
177 +            ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
178 +       preRF2 = 2.0d0*preRF
179 +       preRFCalculated = .true.
180 +    else if (.not.haveDefaultCutoff) then
181 +       call handleError("setReactionFieldPrefactor", "Default cutoff not set")
182 +    else
183 +       call handleError("setReactionFieldPrefactor", "Dielectric not set")
184 +    endif
185 +  end subroutine setReactionFieldPrefactor
186  
187    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
188         is_SplitDipole, is_Quadrupole, is_Tap, status)
# Line 356 | Line 396 | contains
396  
397    subroutine checkSummationMethod()
398  
399 +    if (.not.haveDefaultCutoff) then
400 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
401 +    endif
402 +
403 +    rcuti = 1.0d0 / defaultCutoff
404 +    rcuti2 = rcuti*rcuti
405 +    rcuti3 = rcuti2*rcuti
406 +    rcuti4 = rcuti2*rcuti2
407 +
408      if (summationMethod .eq. DAMPED_WOLF) then
409         if (.not.haveDWAconstants) then
410            
# Line 363 | Line 412 | contains
412               call handleError("checkSummationMethod", "no Damping Alpha set!")
413            endif
414            
415 <          if (.not.have....)
416 <          constEXP =
417 <          constERFC =
418 <          
415 >          if (.not.haveDefaultCutoff) then
416 >             call handleError("checkSummationMethod", "no Default Cutoff set!")
417 >          endif
418 >
419 >          constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
420 >          constERFC = derfc(dampingAlpha*defaultCutoff)
421 >          invRootPi = 0.56418958354775628695d0
422 >          alphaPi = 2*dampingAlpha*invRootPi
423 >  
424            haveDWAconstants = .true.
425         endif
426      endif
427  
428 +    if (summationMethod .eq. REACTION_FIELD) then
429 +       if (.not.haveDielectric) then
430 +          call handleError("checkSummationMethod", "no reaction field Dielectric set!")
431 +       endif
432 +    endif
433 +
434 +    summationMethodChecked = .true.
435    end subroutine checkSummationMethod
436  
437  
438  
439    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
440 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
440 >       vpair, fpair, pot, eFrame, f, t, do_pot)
441  
442      logical, intent(in) :: do_pot
443  
444      integer, intent(in) :: atom1, atom2
445      integer :: localError
385 –    integer, intent(in) :: corrMethod
446  
447 <    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
447 >    real(kind=dp), intent(in) :: rij, r2, sw
448      real(kind=dp), intent(in), dimension(3) :: d
449      real(kind=dp), intent(inout) :: vpair
450      real(kind=dp), intent(inout), dimension(3) :: fpair
451  
452 <    real( kind = dp ) :: pot, swi
452 >    real( kind = dp ) :: pot
453      real( kind = dp ), dimension(9,nLocal) :: eFrame
454      real( kind = dp ), dimension(3,nLocal) :: f
455      real( kind = dp ), dimension(3,nLocal) :: t
# Line 414 | Line 474 | contains
474      real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
475      real (kind=dp) :: xhat, yhat, zhat
476      real (kind=dp) :: dudx, dudy, dudz
477 <    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
478 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
477 >    real (kind=dp) :: scale, sc2, bigR
478 >    real (kind=dp) :: varERFC, varEXP
479 >    real (kind=dp) :: limScale
480 >    real (kind=dp) :: preVal, rfVal
481  
482      if (.not.allocated(ElectrostaticMap)) then
483         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 426 | Line 488 | contains
488         call checkSummationMethod()
489      endif
490  
491 +    if (.not.preRFCalculated) then
492 +       call setReactionFieldPrefactor()
493 +    endif
494  
495   #ifdef IS_MPI
496      me1 = atid_Row(atom1)
# Line 438 | Line 503 | contains
503      !! some variables we'll need independent of electrostatic type:
504  
505      riji = 1.0d0 / rij
506 <
506 >  
507      xhat = d(1) * riji
508      yhat = d(2) * riji
509      zhat = d(3) * riji
510  
446 –    rcuti2 = rcuti*rcuti
447 –    rcuti3 = rcuti2*rcuti
448 –    rcuti4 = rcuti2*rcuti2
449 –
450 –    swi = 1.0d0 / sw
451 –
511      !! logicals
512      i_is_Charge = ElectrostaticMap(me1)%is_Charge
513      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 567 | Line 626 | contains
626         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
627      endif
628    
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 –
629      epot = 0.0_dp
630      dudx = 0.0_dp
631      dudy = 0.0_dp
# Line 592 | Line 643 | contains
643  
644         if (j_is_Charge) then
645  
646 <          if (corrMethod .eq. 1) then
646 >          if (summationMethod .eq. UNDAMPED_WOLF) then
647               vterm = pre11 * q_i * q_j * (riji - rcuti)
648 +             vpair = vpair + vterm
649 +             epot = epot + sw*vterm
650 +            
651 +             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
652 +            
653 +             dudx = dudx + dudr * d(1)
654 +             dudy = dudy + dudr * d(2)
655 +             dudz = dudz + dudr * d(3)
656  
657 +          elseif (summationMethod .eq. DAMPED_WOLF) then
658 +             varERFC = derfc(dampingAlpha*rij)
659 +             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
660 +             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
661               vpair = vpair + vterm
662 <             epot = epot + sw * vterm
662 >             epot = epot + sw*vterm
663              
664 <             dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
664 >             dudr  = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
665 >                                                  + alphaPi*varEXP) &
666 >                                                 - (constERFC*rcuti2 &
667 >                                                    + alphaPi*constEXP)) )
668              
669               dudx = dudx + dudr * d(1)
670               dudy = dudy + dudr * d(2)
671               dudz = dudz + dudr * d(3)
672  
673 +          elseif (summationMethod .eq. REACTION_FIELD) then
674 +             preVal = pre11 * q_i * q_j
675 +             rfVal = preRF*rij*rij
676 +             vterm = preVal * ( riji + rfVal )
677 +             vpair = vpair + vterm
678 +             epot = epot + sw*vterm
679 +            
680 +             dudr  = sw * preVal * ( 2.0d0*rfVal - riji )*riji
681 +            
682 +             dudx = dudx + dudr * xhat
683 +             dudy = dudy + dudr * yhat
684 +             dudz = dudz + dudr * zhat
685 +
686            else
687               vterm = pre11 * q_i * q_j * riji
609 –
688               vpair = vpair + vterm
689 <             epot = epot + sw * vterm
689 >             epot = epot + sw*vterm
690              
691               dudr  = - sw * vterm * riji
692              
# Line 622 | Line 700 | contains
700  
701         if (j_is_Dipole) then
702  
703 <          pref = sw * pre12 * q_i * mu_j
703 >          pref = pre12 * q_i * mu_j
704  
705 <          if (corrMethod .eq. 1) then
705 >          if (summationMethod .eq. UNDAMPED_WOLF) then
706               ri2 = riji * riji
707               ri3 = ri2 * riji
708  
709 +             pref = pre12 * q_i * mu_j
710               vterm = - pref * ct_j * (ri2 - rcuti2)
711 <             vpair = vpair + swi*vterm
712 <             epot = epot + vterm
711 >             vpair = vpair + vterm
712 >             epot = epot + sw*vterm
713              
714               !! this has a + sign in the () because the rij vector is
715               !! r_j - r_i and the charge-dipole potential takes the origin
716               !! as the point dipole, which is atom j in this case.
717              
718 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
718 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
719                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
720 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
720 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
721                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
722 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
722 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
723                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
724              
725 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
726 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
727 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
725 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
726 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
727 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
728  
729 +          elseif (summationMethod .eq. REACTION_FIELD) then
730 +             ri2 = ri * ri
731 +             ri3 = ri2 * ri
732 +    
733 +             pref = pre12 * q_i * mu_j
734 +             vterm = - pref * ct_j * ( ri2 - preRF2*rij )
735 +             vpair = vpair + vterm
736 +             epot = epot + sw*vterm
737 +            
738 +             !! this has a + sign in the () because the rij vector is
739 +             !! r_j - r_i and the charge-dipole potential takes the origin
740 +             !! as the point dipole, which is atom j in this case.
741 +            
742 +             dudx = dudx - sw*pref*( ri3*(uz_j(1) - 3.0d0*ct_j*xhat) - &
743 +                                     preRF2*uz_j(1) )
744 +             dudy = dudy - sw*pref*( ri3*(uz_j(2) - 3.0d0*ct_j*yhat) - &
745 +                                     preRF2*uz_j(2) )
746 +             dudz = dudz - sw*pref*( ri3*(uz_j(3) - 3.0d0*ct_j*zhat) - &
747 +                                     preRF2*uz_j(3) )        
748 +             duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
749 +             duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
750 +             duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
751 +
752            else
753               if (j_is_SplitDipole) then
754                  BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
# Line 660 | Line 762 | contains
762               ri2 = ri * ri
763               ri3 = ri2 * ri
764               sc2 = scale * scale
765 <            
765 >
766 >             pref = pre12 * q_i * mu_j
767               vterm = - pref * ct_j * ri2 * scale
768 <             vpair = vpair + swi * vterm
769 <             epot = epot + vterm
768 >             vpair = vpair + vterm
769 >             epot = epot + sw*vterm
770              
771               !! this has a + sign in the () because the rij vector is
772               !! r_j - r_i and the charge-dipole potential takes the origin
773               !! as the point dipole, which is atom j in this case.
774              
775 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
776 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
777 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
775 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
776 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
777 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
778              
779 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
780 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
781 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
779 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
780 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
781 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
782  
783            endif
784         endif
# Line 688 | Line 791 | contains
791            cy2 = cy_j * cy_j
792            cz2 = cz_j * cz_j
793  
794 <
795 <          pref =  sw * pre14 * q_i / 3.0_dp
693 <
694 <          if (corrMethod .eq. 1) then
794 >          if (summationMethod .eq. UNDAMPED_WOLF) then
795 >             pref =  pre14 * q_i / 3.0_dp
796               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
797                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
798                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
799               vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
800                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
801                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
802 <             vpair = vpair + swi*( vterm1 - vterm2 )
803 <             epot = epot + ( vterm1 - vterm2 )
802 >             vpair = vpair + ( vterm1 - vterm2 )
803 >             epot = epot + sw*( vterm1 - vterm2 )
804              
805               dudx = dudx - (5.0_dp * &
806 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
806 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
807                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
808                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
809                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 710 | Line 811 | contains
811                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
812                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
813               dudy = dudy - (5.0_dp * &
814 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
814 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
815                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
816                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
817                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 718 | Line 819 | contains
819                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
820                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
821               dudz = dudz - (5.0_dp * &
822 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
822 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
823                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
824                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
825                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 726 | Line 827 | contains
827                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
828                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
829              
830 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
830 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
831                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
832 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
832 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
833                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
834 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
834 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
835                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
836              
837 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
837 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
838                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
839 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
839 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
840                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
841 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
841 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
842                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
843              
844 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
844 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
845                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
846 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
846 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
847                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
848 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
848 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
849                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
850          
851            else
852 +             pref =  pre14 * q_i / 3.0_dp
853               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
854                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
855                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
856 <             vpair = vpair + swi * vterm
857 <             epot = epot + vterm
856 >             vpair = vpair + vterm
857 >             epot = epot + sw*vterm
858              
859 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
859 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
860                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
861                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
862                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
863 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
863 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
864                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
865                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
866                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
867 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
867 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
868                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
869                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
870                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
871              
872 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
873 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
874 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
872 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
873 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
874 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
875              
876 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
877 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
878 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
876 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
877 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
878 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
879              
880 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
881 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
882 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
880 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
881 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
882 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
883            
884            endif
885         endif
# Line 786 | Line 888 | contains
888      if (i_is_Dipole) then
889  
890         if (j_is_Charge) then
891 <
892 <          pref = sw * pre12 * q_j * mu_i
893 <
894 <          if (corrMethod .eq. 1) then
891 >          
892 >          pref = pre12 * q_j * mu_i
893 >          
894 >          if (summationMethod .eq. UNDAMPED_WOLF) then
895               ri2 = riji * riji
896               ri3 = ri2 * riji
897  
898 +             pref = pre12 * q_j * mu_i
899               vterm = pref * ct_i * (ri2 - rcuti2)
900 <             vpair = vpair + swi * vterm
901 <             epot = epot + vterm
900 >             vpair = vpair + vterm
901 >             epot = epot + sw*vterm
902              
903               !! this has a + sign in the () because the rij vector is
904               !! r_j - r_i and the charge-dipole potential takes the origin
905               !! as the point dipole, which is atom j in this case.
906              
907 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
907 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
908                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
909 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
909 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
910                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
911 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
911 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
912                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
913              
914 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
915 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
916 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
914 >             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
915 >             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
916 >             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
917  
918 +          elseif (summationMethod .eq. REACTION_FIELD) then
919 +             ri2 = ri * ri
920 +             ri3 = ri2 * ri
921 +
922 +             pref = pre12 * q_j * mu_i
923 +             vterm = pref * ct_i * ( ri2 - preRF*rij )
924 +             vpair = vpair + vterm
925 +             epot = epot + sw*vterm
926 +            
927 +             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0*ct_i*xhat - &
928 +                                             preRF*uz_i(1) )
929 +             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0*ct_i*yhat - &
930 +                                             preRF*uz_i(2) )
931 +             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0*ct_i*zhat - &
932 +                                             preRF*uz_i(3) )
933 +            
934 +             duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF*rij )
935 +             duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF*rij )
936 +             duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF*rij )
937 +
938            else
939               if (i_is_SplitDipole) then
940                  BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
# Line 825 | Line 948 | contains
948               ri2 = ri * ri
949               ri3 = ri2 * ri
950               sc2 = scale * scale
951 <            
951 >
952 >             pref = pre12 * q_j * mu_i
953               vterm = pref * ct_i * ri2 * scale
954 <             vpair = vpair + swi * vterm
955 <             epot = epot + vterm
954 >             vpair = vpair + vterm
955 >             epot = epot + sw*vterm
956              
957 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
958 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
959 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
957 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
958 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
959 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
960              
961 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
962 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
963 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
961 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
962 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
963 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
964            endif
965         endif
966 <
966 >      
967         if (j_is_Dipole) then
968  
969 <          pref = sw * pre22 * mu_i * mu_j
846 <
847 <          if (corrMethod .eq. 1) then
969 >          if (summationMethod .eq. UNDAMPED_WOLF) then
970               ri2 = riji * riji
971               ri3 = ri2 * riji
972               ri4 = ri2 * ri2
973  
974 +             pref = pre22 * mu_i * mu_j
975               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
976 <             vpair = vpair + swi * vterm
977 <             epot = epot + vterm
976 >             vpair = vpair + vterm
977 >             epot = epot + sw*vterm
978              
979               a1 = 5.0d0 * ct_i * ct_j - ct_ij
980              
981 <             dudx = dudx + pref*3.0d0*ri4 &
982 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
983 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
984 <             dudy = dudy + pref*3.0d0*ri4 &
985 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
986 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
987 <             dudz = dudz + pref*3.0d0*ri4 &
988 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
989 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
981 >             dudx = dudx + sw*pref*3.0d0*ri4 &
982 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
983 >                         - sw*pref*3.0d0*rcuti4 &
984 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
985 >             dudy = dudy + sw*pref*3.0d0*ri4 &
986 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
987 >                         - sw*pref*3.0d0*rcuti4 &
988 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
989 >             dudz = dudz + sw*pref*3.0d0*ri4 &
990 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
991 >                         - sw*pref*3.0d0*rcuti4 &
992 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
993              
994 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
994 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
995                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
996 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
996 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
997                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
998 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
998 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
999                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
1000 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
1000 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1001                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
1002 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
1002 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1003                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
1004 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
1004 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1005                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
1006 <          else
1006 >
1007 >         elseif (summationMethod .eq. REACTION_FIELD) then
1008 >             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1009 >
1010 >             ri2 = riji * riji
1011 >             ri3 = ri2 * riji
1012 >             ri4 = ri2 * ri2
1013 >
1014 >             pref = pre22 * mu_i * mu_j
1015 >              
1016 >             vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
1017 >                  preRF2*ct_ij )
1018 >             vpair = vpair + vterm
1019 >             epot = epot + sw*vterm
1020              
1021 +             a1 = 5.0d0 * ct_i * ct_j - ct_ij
1022 +            
1023 +             dudx = dudx + sw*pref*3.0d0*ri4 &
1024 +                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1025 +             dudy = dudy + sw*pref*3.0d0*ri4 &
1026 +                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1027 +             dudz = dudz + sw*pref*3.0d0*ri4 &
1028 +                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1029 +            
1030 +             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1031 +                  - preRF2*uz_j(1))
1032 +             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1033 +                  - preRF2*uz_j(2))
1034 +             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1035 +                  - preRF2*uz_j(3))
1036 +             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1037 +                  - preRF2*uz_i(1))
1038 +             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1039 +                  - preRF2*uz_i(2))
1040 +             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1041 +                  - preRF2*uz_i(3))
1042 +
1043 +          else
1044               if (i_is_SplitDipole) then
1045                  if (j_is_SplitDipole) then
1046                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 905 | Line 1067 | contains
1067               ri4 = ri2 * ri2
1068               sc2 = scale * scale
1069              
1070 +             pref = pre22 * mu_i * mu_j
1071               vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1072 <             vpair = vpair + swi * vterm
1073 <             epot = epot + vterm
1072 >             vpair = vpair + vterm
1073 >             epot = epot + sw*vterm
1074              
1075               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1076              
1077 <             dudx = dudx + pref*3.0d0*ri4*scale &
1078 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1079 <             dudy = dudy + pref*3.0d0*ri4*scale &
1080 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1081 <             dudz = dudz + pref*3.0d0*ri4*scale &
1082 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1077 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
1078 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1079 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
1080 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1081 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
1082 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1083              
1084 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
1085 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1086 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
1087 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1088 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
1089 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1084 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
1085 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1086 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
1087 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1088 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
1089 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1090              
1091 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
1092 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1093 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
1094 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1095 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
1096 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1091 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1092 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1093 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1094 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1095 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1096 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1097            endif
1098         endif
1099      endif
# Line 945 | Line 1108 | contains
1108            cy2 = cy_i * cy_i
1109            cz2 = cz_i * cz_i
1110  
1111 <          pref = sw * pre14 * q_j / 3.0_dp
1112 <
950 <          if (corrMethod .eq. 1) then
1111 >          if (summationMethod .eq. UNDAMPED_WOLF) then
1112 >             pref = pre14 * q_j / 3.0_dp
1113               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1114                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1115                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1116               vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1117                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1118                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1119 <             vpair = vpair + swi * ( vterm1 - vterm2 )
1120 <             epot = epot + ( vterm1 - vterm2 )
1119 >             vpair = vpair + ( vterm1 - vterm2 )
1120 >             epot = epot + sw*( vterm1 - vterm2 )
1121              
1122 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
1123 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1122 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1123 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1124                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1125                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1126                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1127                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1128                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1129 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
1130 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1129 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1130 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1131                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1132                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1133                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1134                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1135                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1136 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
1137 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1136 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1137 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1138                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1139                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1140                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1141                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1142                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1143              
1144 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1144 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1145                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1146 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1146 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1147                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1148 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1148 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1149                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1150              
1151 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1151 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1152                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1153 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1153 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1154                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1155 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1155 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1156                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1157              
1158 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1158 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1159                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1160 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1160 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1161                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1162 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1162 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1163                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1164  
1165            else
1166 +             pref = pre14 * q_j / 3.0_dp
1167               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1168                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1169                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1170 <             vpair = vpair + swi * vterm
1171 <             epot = epot + vterm
1170 >             vpair = vpair + vterm
1171 >             epot = epot + sw*vterm
1172              
1173 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1173 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1174                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1175                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1176                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1177 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1177 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1178                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1179                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1180                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1181 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1181 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1182                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1183                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1184                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1185              
1186 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1187 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1188 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1186 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1187 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1188 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1189              
1190 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1191 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1192 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1190 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1191 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1192 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1193              
1194 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1195 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1196 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1194 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1195 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1196 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1197            endif
1198         endif
1199      endif
# Line 1038 | Line 1201 | contains
1201  
1202      if (do_pot) then
1203   #ifdef IS_MPI
1204 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1205 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1204 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1205 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1206   #else
1207         pot = pot + epot
1208   #endif
# Line 1144 | Line 1307 | contains
1307      return
1308    end subroutine doElectrostaticPair
1309  
1310 <  !! calculates the switching functions and their derivatives for a given
1148 <  subroutine calc_switch(r, mu, scale, dscale)
1310 >  subroutine destroyElectrostaticTypes()
1311  
1312 <    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
1312 >    if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1313  
1314 <    ! distances must be in angstroms
1315 <    rl = 2.75d0
1316 <    ru = 3.75d0
1317 <    mulow = 0.0d0 !3.3856d0 ! 1.84 * 1.84
1318 <    minRatio = mulow / (mu*mu)
1319 <    scaleVal = 1.0d0 - minRatio
1314 >  end subroutine destroyElectrostaticTypes
1315 >
1316 >  subroutine rf_self_self(atom1, eFrame, rfpot, t, do_pot)
1317 >    logical, intent(in) :: do_pot
1318 >    integer, intent(in) :: atom1
1319 >    integer :: atid1
1320 >    real(kind=dp), dimension(9,nLocal) :: eFrame
1321 >    real(kind=dp), dimension(3,nLocal) :: t
1322 >    real(kind=dp) :: mu1
1323 >    real(kind=dp) :: preVal, epot, rfpot
1324 >    real(kind=dp) :: eix, eiy, eiz
1325 >
1326 >    ! this is a local only array, so we use the local atom type id's:
1327 >    atid1 = atid(atom1)
1328      
1329 <    if (r.lt.rl) then
1330 <       scale = minRatio
1331 <       dscale = 0.0d0
1332 <    elseif (r.gt.ru) then
1333 <       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
1329 >    if (ElectrostaticMap(atid1)%is_Dipole) then
1330 >       mu1 = getDipoleMoment(atid1)
1331 >      
1332 >       preVal = pre22 * preRF2 * mu1*mu1
1333 >       rfpot = rfpot - 0.5d0*preVal
1334  
1335 <  subroutine destroyElectrostaticTypes()
1335 >       ! The self-correction term adds into the reaction field vector
1336 >      
1337 >       eix = preVal * eFrame(3,atom1)
1338 >       eiy = preVal * eFrame(6,atom1)
1339 >       eiz = preVal * eFrame(9,atom1)
1340  
1341 <    if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1341 >       ! once again, this is self-self, so only the local arrays are needed
1342 >       ! even for MPI jobs:
1343 >      
1344 >       t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1345 >            eFrame(9,atom1)*eiy
1346 >       t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1347 >            eFrame(3,atom1)*eiz
1348 >       t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1349 >            eFrame(6,atom1)*eix
1350  
1351 <  end subroutine destroyElectrostaticTypes
1351 >    endif
1352 >    
1353 >    return
1354 >  end subroutine rf_self_self
1355  
1356   end module electrostatic_module

Diff Legend

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