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 2390 by chrisfen, Wed Oct 19 19:24: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 >  logical, save :: preRFCalculated = .false.
101 >
102 > #ifdef __IFC
103 > ! error function for ifc version > 7.
104 >  double precision, external :: derfc
105 > #endif
106 >  
107    public :: setElectrostaticSummationMethod
108    public :: setElectrostaticCutoffRadius
109    public :: setDampedWolfAlpha
110    public :: setReactionFieldDielectric
111 +  public :: setReactionFieldPrefactor
112    public :: newElectrostaticType
113    public :: setCharge
114    public :: setDipoleMoment
# Line 98 | Line 119 | module electrostatic_module
119    public :: getDipoleMoment
120    public :: pre22
121    public :: destroyElectrostaticTypes
122 +  public :: accumulate_rf
123 +  public :: accumulate_self_rf
124 +  public :: reaction_field_final
125 +  public :: rf_correct_forces
126  
127    type :: Electrostatic
128       integer :: c_ident
# Line 117 | Line 142 | contains
142   contains
143  
144    subroutine setElectrostaticSummationMethod(the_ESM)
120 –
145      integer, intent(in) :: the_ESM    
146  
147      if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
148         call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
149      endif
150  
151 +    summationMethod = the_ESM
152 +
153    end subroutine setElectrostaticSummationMethod
154  
155 <  subroutine setElectrostaticCutoffRadius(thisRcut)
155 >  subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
156      real(kind=dp), intent(in) :: thisRcut
157 +    real(kind=dp), intent(in) :: thisRsw
158      defaultCutoff = thisRcut
159 +    rrf = defaultCutoff
160 +    rt = thisRsw
161      haveDefaultCutoff = .true.
162    end subroutine setElectrostaticCutoffRadius
163  
# Line 143 | Line 172 | contains
172      dielectric = thisDielectric
173      haveDielectric = .true.
174    end subroutine setReactionFieldDielectric
175 +
176 +  subroutine setReactionFieldPrefactor
177 +    if (haveDefaultCutoff .and. haveDielectric) then
178 +       defaultCutoff2 = defaultCutoff*defaultCutoff
179 +       preRF = pre22 * 2.0d0*(dielectric-1.0d0) / &
180 +            ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
181 +       preRFCalculated = .true.
182 +    else if (.not.haveDefaultCutoff) then
183 +       call handleError("setReactionFieldPrefactor", "Default cutoff not set")
184 +    else
185 +       call handleError("setReactionFieldPrefactor", "Dielectric not set")
186 +    endif
187 +  end subroutine setReactionFieldPrefactor
188  
189    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
190         is_SplitDipole, is_Quadrupole, is_Tap, status)
# Line 355 | Line 397 | contains
397    end function getDipoleMoment
398  
399    subroutine checkSummationMethod()
400 +
401 +    if (.not.haveDefaultCutoff) then
402 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
403 +    endif
404  
405 +    rcuti = 1.0d0 / defaultCutoff
406 +    rcuti2 = rcuti*rcuti
407 +    rcuti3 = rcuti2*rcuti
408 +    rcuti4 = rcuti2*rcuti2
409 +
410      if (summationMethod .eq. DAMPED_WOLF) then
411         if (.not.haveDWAconstants) then
412            
# Line 363 | Line 414 | contains
414               call handleError("checkSummationMethod", "no Damping Alpha set!")
415            endif
416            
417 <          if (.not.have....)
418 <          constEXP =
419 <          constERFC =
420 <          
417 >          if (.not.haveDefaultCutoff) then
418 >             call handleError("checkSummationMethod", "no Default Cutoff set!")
419 >          endif
420 >
421 >          constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
422 >          constERFC = derfc(dampingAlpha*defaultCutoff)
423 >          invRootPi = 0.56418958354775628695d0
424 >          alphaPi = 2*dampingAlpha*invRootPi
425 >  
426            haveDWAconstants = .true.
427         endif
428      endif
429  
430 +    if (summationMethod .eq. REACTION_FIELD) then
431 +       if (.not.haveDielectric) then
432 +          call handleError("checkSummationMethod", "no reaction field Dielectric set!")
433 +       endif
434 +    endif
435 +
436 +    summationMethodChecked = .true.
437    end subroutine checkSummationMethod
438  
439  
440  
441    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
442 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
442 >       vpair, fpair, pot, eFrame, f, t, do_pot)
443  
444      logical, intent(in) :: do_pot
445  
446      integer, intent(in) :: atom1, atom2
447      integer :: localError
385 –    integer, intent(in) :: corrMethod
448  
449 <    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
449 >    real(kind=dp), intent(in) :: rij, r2, sw
450      real(kind=dp), intent(in), dimension(3) :: d
451      real(kind=dp), intent(inout) :: vpair
452      real(kind=dp), intent(inout), dimension(3) :: fpair
453  
454 <    real( kind = dp ) :: pot, swi
454 >    real( kind = dp ) :: pot
455      real( kind = dp ), dimension(9,nLocal) :: eFrame
456      real( kind = dp ), dimension(3,nLocal) :: f
457      real( kind = dp ), dimension(3,nLocal) :: t
# Line 414 | Line 476 | contains
476      real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
477      real (kind=dp) :: xhat, yhat, zhat
478      real (kind=dp) :: dudx, dudy, dudz
479 <    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
480 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
479 >    real (kind=dp) :: scale, sc2, bigR
480 >    real (kind=dp) :: varERFC, varEXP
481 >    real (kind=dp) :: limScale
482  
483      if (.not.allocated(ElectrostaticMap)) then
484         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 424 | Line 487 | contains
487  
488      if (.not.summationMethodChecked) then
489         call checkSummationMethod()
490 +      
491      endif
492  
493  
# Line 438 | Line 502 | contains
502      !! some variables we'll need independent of electrostatic type:
503  
504      riji = 1.0d0 / rij
505 <
505 >  
506      xhat = d(1) * riji
507      yhat = d(2) * riji
508      zhat = d(3) * riji
509  
446 –    rcuti2 = rcuti*rcuti
447 –    rcuti3 = rcuti2*rcuti
448 –    rcuti4 = rcuti2*rcuti2
449 –
450 –    swi = 1.0d0 / sw
451 –
510      !! logicals
511      i_is_Charge = ElectrostaticMap(me1)%is_Charge
512      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 567 | Line 625 | contains
625         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
626      endif
627    
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 –
628      epot = 0.0_dp
629      dudx = 0.0_dp
630      dudy = 0.0_dp
# Line 592 | Line 642 | contains
642  
643         if (j_is_Charge) then
644  
645 <          if (corrMethod .eq. 1) then
645 >          if (summationMethod .eq. UNDAMPED_WOLF) then
646 >
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 +
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            else
608 –             vterm = pre11 * q_i * q_j * riji
675  
676 +             vterm = pre11 * q_i * q_j * riji
677               vpair = vpair + vterm
678 <             epot = epot + sw * vterm
678 >             epot = epot + sw*vterm
679              
680               dudr  = - sw * vterm * riji
681              
# Line 622 | Line 689 | contains
689  
690         if (j_is_Dipole) then
691  
692 <          pref = sw * pre12 * q_i * mu_j
692 >          pref = pre12 * q_i * mu_j
693  
694 <          if (corrMethod .eq. 1) then
694 >          if (summationMethod .eq. UNDAMPED_WOLF) then
695               ri2 = riji * riji
696               ri3 = ri2 * riji
697  
698 +             pref = pre12 * q_i * mu_j
699               vterm = - pref * ct_j * (ri2 - rcuti2)
700 <             vpair = vpair + swi*vterm
701 <             epot = epot + vterm
700 >             vpair = vpair + vterm
701 >             epot = epot + sw*vterm
702              
703               !! this has a + sign in the () because the rij vector is
704               !! r_j - r_i and the charge-dipole potential takes the origin
705               !! as the point dipole, which is atom j in this case.
706              
707 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
707 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
708                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
709 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
709 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
710                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
711 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
711 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
712                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
713              
714 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
715 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
716 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
714 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
715 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
716 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
717  
718            else
719               if (j_is_SplitDipole) then
# Line 660 | Line 728 | contains
728               ri2 = ri * ri
729               ri3 = ri2 * ri
730               sc2 = scale * scale
731 <            
731 >
732 >             pref = pre12 * q_i * mu_j
733               vterm = - pref * ct_j * ri2 * scale
734 <             vpair = vpair + swi * vterm
735 <             epot = epot + vterm
734 >             vpair = vpair + vterm
735 >             epot = epot + sw*vterm
736              
737               !! this has a + sign in the () because the rij vector is
738               !! r_j - r_i and the charge-dipole potential takes the origin
739               !! as the point dipole, which is atom j in this case.
740              
741 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
742 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
743 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
741 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
742 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
743 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
744              
745 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
746 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
747 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
745 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
746 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
747 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
748  
749            endif
750         endif
# Line 688 | Line 757 | contains
757            cy2 = cy_j * cy_j
758            cz2 = cz_j * cz_j
759  
760 <
761 <          pref =  sw * pre14 * q_i / 3.0_dp
693 <
694 <          if (corrMethod .eq. 1) then
760 >          if (summationMethod .eq. UNDAMPED_WOLF) then
761 >             pref =  pre14 * q_i / 3.0_dp
762               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
763                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
764                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
765               vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
766                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
767                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
768 <             vpair = vpair + swi*( vterm1 - vterm2 )
769 <             epot = epot + ( vterm1 - vterm2 )
768 >             vpair = vpair + ( vterm1 - vterm2 )
769 >             epot = epot + sw*( vterm1 - vterm2 )
770              
771               dudx = dudx - (5.0_dp * &
772 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
772 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
773                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
774                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
775                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 710 | Line 777 | contains
777                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
778                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
779               dudy = dudy - (5.0_dp * &
780 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
780 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
781                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
782                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
783                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 718 | Line 785 | contains
785                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
786                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
787               dudz = dudz - (5.0_dp * &
788 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
788 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
789                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
790                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
791                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 726 | Line 793 | contains
793                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
794                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
795              
796 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
796 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
797                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
798 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
798 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
799                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
800 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
800 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
801                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
802              
803 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
803 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
804                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
805 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
805 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
806                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
807 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
807 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
808                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
809              
810 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
810 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
811                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
812 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
812 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
813                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
814 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
814 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
815                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
816          
817            else
818 +             pref =  pre14 * q_i / 3.0_dp
819               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
820                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
821                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
822 <             vpair = vpair + swi * vterm
823 <             epot = epot + vterm
822 >             vpair = vpair + vterm
823 >             epot = epot + sw*vterm
824              
825 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
825 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
826                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
827                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
828                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
829 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
829 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
830                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
831                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
832                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
833 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
833 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
834                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
835                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
836                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
837              
838 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
839 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
840 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
838 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
839 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
840 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
841              
842 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
843 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
844 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
842 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
843 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
844 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
845              
846 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
847 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
848 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
846 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
847 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
848 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
849            
850            endif
851         endif
# Line 786 | Line 854 | contains
854      if (i_is_Dipole) then
855  
856         if (j_is_Charge) then
857 <
858 <          pref = sw * pre12 * q_j * mu_i
859 <
860 <          if (corrMethod .eq. 1) then
857 >          
858 >          pref = pre12 * q_j * mu_i
859 >          
860 >          if (summationMethod .eq. UNDAMPED_WOLF) then
861               ri2 = riji * riji
862               ri3 = ri2 * riji
863  
864 +             pref = pre12 * q_j * mu_i
865               vterm = pref * ct_i * (ri2 - rcuti2)
866 <             vpair = vpair + swi * vterm
867 <             epot = epot + vterm
866 >             vpair = vpair + vterm
867 >             epot = epot + sw*vterm
868              
869               !! this has a + sign in the () because the rij vector is
870               !! r_j - r_i and the charge-dipole potential takes the origin
871               !! as the point dipole, which is atom j in this case.
872              
873 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
873 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
874                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
875 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
875 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
876                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
877 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
877 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
878                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
879              
880 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
881 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
882 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
880 >             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
881 >             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
882 >             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
883  
884            else
885               if (i_is_SplitDipole) then
# Line 825 | Line 894 | contains
894               ri2 = ri * ri
895               ri3 = ri2 * ri
896               sc2 = scale * scale
897 <            
897 >
898 >             pref = pre12 * q_j * mu_i
899               vterm = pref * ct_i * ri2 * scale
900 <             vpair = vpair + swi * vterm
901 <             epot = epot + vterm
900 >             vpair = vpair + vterm
901 >             epot = epot + sw*vterm
902              
903 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
904 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
905 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
903 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
904 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
905 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
906              
907 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
908 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
909 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
907 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
908 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
909 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
910            endif
911         endif
912 <
912 >      
913         if (j_is_Dipole) then
914  
915 <          pref = sw * pre22 * mu_i * mu_j
846 <
847 <          if (corrMethod .eq. 1) then
915 >          if (summationMethod .eq. UNDAMPED_WOLF) then
916               ri2 = riji * riji
917               ri3 = ri2 * riji
918               ri4 = ri2 * ri2
919  
920 +             pref = pre22 * mu_i * mu_j
921               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
922 <             vpair = vpair + swi * vterm
923 <             epot = epot + vterm
922 >             vpair = vpair + vterm
923 >             epot = epot + sw*vterm
924              
925               a1 = 5.0d0 * ct_i * ct_j - ct_ij
926              
927 <             dudx = dudx + pref*3.0d0*ri4 &
928 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
929 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
930 <             dudy = dudy + pref*3.0d0*ri4 &
931 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
932 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
933 <             dudz = dudz + pref*3.0d0*ri4 &
934 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
935 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
927 >             dudx = dudx + sw*pref*3.0d0*ri4 &
928 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
929 >                         - sw*pref*3.0d0*rcuti4 &
930 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
931 >             dudy = dudy + sw*pref*3.0d0*ri4 &
932 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
933 >                         - sw*pref*3.0d0*rcuti4 &
934 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
935 >             dudz = dudz + sw*pref*3.0d0*ri4 &
936 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
937 >                         - sw*pref*3.0d0*rcuti4 &
938 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
939              
940 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
940 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
941                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
942 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
942 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
943                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
944 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
944 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
945                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
946 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
946 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
947                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
948 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
948 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
949                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
950 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
950 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
951                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
952 +
953            else
881 –            
954               if (i_is_SplitDipole) then
955                  if (j_is_SplitDipole) then
956                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 905 | Line 977 | contains
977               ri4 = ri2 * ri2
978               sc2 = scale * scale
979              
980 +             pref = pre22 * mu_i * mu_j
981               vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
982 <             vpair = vpair + swi * vterm
983 <             epot = epot + vterm
982 >             vpair = vpair + vterm
983 >             epot = epot + sw*vterm
984              
985               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
986              
987 <             dudx = dudx + pref*3.0d0*ri4*scale &
988 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
989 <             dudy = dudy + pref*3.0d0*ri4*scale &
990 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
991 <             dudz = dudz + pref*3.0d0*ri4*scale &
992 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
920 <            
921 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
922 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
923 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
924 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
925 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
926 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
987 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
988 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
989 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
990 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
991 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
992 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
993              
994 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
995 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
996 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
997 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
998 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
999 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
994 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
995 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
996 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
997 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
998 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
999 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1000 >            
1001 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1002 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1003 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1004 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1005 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1006 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1007            endif
1008         endif
1009      endif
# Line 945 | Line 1018 | contains
1018            cy2 = cy_i * cy_i
1019            cz2 = cz_i * cz_i
1020  
1021 <          pref = sw * pre14 * q_j / 3.0_dp
1022 <
950 <          if (corrMethod .eq. 1) then
1021 >          if (summationMethod .eq. UNDAMPED_WOLF) then
1022 >             pref = pre14 * q_j / 3.0_dp
1023               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1024                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1025                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1026               vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1027                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1028                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1029 <             vpair = vpair + swi * ( vterm1 - vterm2 )
1030 <             epot = epot + ( vterm1 - vterm2 )
1029 >             vpair = vpair + ( vterm1 - vterm2 )
1030 >             epot = epot + sw*( vterm1 - vterm2 )
1031              
1032 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
1033 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1032 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1033 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1034                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1035                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1036                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1037                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1038                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1039 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
1040 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1039 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1040 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1041                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1042                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1043                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1044                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1045                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1046 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
1047 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1046 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1047 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1048                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1049                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1050                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1051                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1052                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1053              
1054 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1054 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1055                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1056 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1056 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1057                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1058 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1058 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1059                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1060              
1061 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1061 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1062                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1063 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1063 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1064                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1065 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1065 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1066                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1067              
1068 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1068 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1069                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1070 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1070 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1071                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1072 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1072 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1073                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1074  
1075            else
1076 +             pref = pre14 * q_j / 3.0_dp
1077               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1078                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1079                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1080 <             vpair = vpair + swi * vterm
1081 <             epot = epot + vterm
1080 >             vpair = vpair + vterm
1081 >             epot = epot + sw*vterm
1082              
1083 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1083 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1084                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1085                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1086                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1087 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1087 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1088                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1089                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1090                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1091 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1091 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1092                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1093                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1094                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1095              
1096 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1097 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1098 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1096 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1097 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1098 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1099              
1100 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1101 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1102 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1100 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1101 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1102 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1103              
1104 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1105 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1106 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1104 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1105 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1106 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1107            endif
1108         endif
1109      endif
# Line 1038 | Line 1111 | contains
1111  
1112      if (do_pot) then
1113   #ifdef IS_MPI
1114 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1115 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1114 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1115 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1116   #else
1117         pot = pot + epot
1118   #endif
# Line 1144 | Line 1217 | contains
1217      return
1218    end subroutine doElectrostaticPair
1219  
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
1175 –
1220    subroutine destroyElectrostaticTypes()
1221  
1222      if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1223  
1224    end subroutine destroyElectrostaticTypes
1225 +
1226 +  subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper)
1227 +
1228 +    integer, intent(in) :: atom1, atom2
1229 +    real (kind = dp), intent(in) :: rij
1230 +    real (kind = dp), dimension(9,nLocal) :: eFrame
1231 +
1232 +    integer :: me1, me2
1233 +    real (kind = dp), intent(in) :: taper
1234 +    real (kind = dp):: mu1, mu2
1235 +    real (kind = dp), dimension(3) :: ul1
1236 +    real (kind = dp), dimension(3) :: ul2  
1237 +
1238 +    integer :: localError
1239 +
1240 + #ifdef IS_MPI
1241 +    me1 = atid_Row(atom1)
1242 +    ul1(1) = eFrame_Row(3,atom1)
1243 +    ul1(2) = eFrame_Row(6,atom1)
1244 +    ul1(3) = eFrame_Row(9,atom1)
1245 +
1246 +    me2 = atid_Col(atom2)
1247 +    ul2(1) = eFrame_Col(3,atom2)
1248 +    ul2(2) = eFrame_Col(6,atom2)
1249 +    ul2(3) = eFrame_Col(9,atom2)
1250 + #else
1251 +    me1 = atid(atom1)
1252 +    ul1(1) = eFrame(3,atom1)
1253 +    ul1(2) = eFrame(6,atom1)
1254 +    ul1(3) = eFrame(9,atom1)
1255 +
1256 +    me2 = atid(atom2)
1257 +    ul2(1) = eFrame(3,atom2)
1258 +    ul2(2) = eFrame(6,atom2)
1259 +    ul2(3) = eFrame(9,atom2)
1260 + #endif
1261 +
1262 +    mu1 = getDipoleMoment(me1)
1263 +    mu2 = getDipoleMoment(me2)
1264 +
1265 + #ifdef IS_MPI
1266 +    rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
1267 +    rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
1268 +    rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
1269 +
1270 +    rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
1271 +    rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
1272 +    rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
1273 + #else
1274 +    rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
1275 +    rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
1276 +    rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
1277  
1278 +    rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
1279 +    rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
1280 +    rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper    
1281 + #endif
1282 +    return  
1283 +  end subroutine accumulate_rf
1284 +
1285 +  subroutine accumulate_self_rf(atom1, mu1, eFrame)
1286 +
1287 +    integer, intent(in) :: atom1
1288 +    real(kind=dp), intent(in) :: mu1
1289 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1290 +
1291 +    !! should work for both MPI and non-MPI version since this is not pairwise.
1292 +    rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
1293 +    rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
1294 +    rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1
1295 +
1296 +    return
1297 +  end subroutine accumulate_self_rf
1298 +
1299 +  subroutine reaction_field_final(a1, mu1, eFrame, rfpot, t, do_pot)
1300 +
1301 +    integer, intent(in) :: a1
1302 +    real (kind=dp), intent(in) :: mu1
1303 +    real (kind=dp), intent(inout) :: rfpot
1304 +    logical, intent(in) :: do_pot
1305 +    real (kind = dp), dimension(9,nLocal) :: eFrame
1306 +    real (kind = dp), dimension(3,nLocal) :: t
1307 +
1308 +    integer :: localError
1309 +
1310 +    if (.not.preRFCalculated) then
1311 +       call setReactionFieldPrefactor()
1312 +    endif
1313 +
1314 +    ! compute torques on dipoles:
1315 +    ! pre converts from mu in units of debye to kcal/mol
1316 +
1317 +    ! The torque contribution is dipole cross reaction_field  
1318 +
1319 +    t(1,a1) = t(1,a1) + preRF*mu1*(eFrame(6,a1)*rf(3,a1) - &
1320 +                                   eFrame(9,a1)*rf(2,a1))
1321 +    t(2,a1) = t(2,a1) + preRF*mu1*(eFrame(9,a1)*rf(1,a1) - &
1322 +                                   eFrame(3,a1)*rf(3,a1))
1323 +    t(3,a1) = t(3,a1) + preRF*mu1*(eFrame(3,a1)*rf(2,a1) - &
1324 +                                   eFrame(6,a1)*rf(1,a1))
1325 +
1326 +    ! the potential contribution is -1/2 dipole dot reaction_field
1327 +
1328 +    if (do_pot) then
1329 +       rfpot = rfpot - 0.5d0 * preRF * mu1 * &
1330 +            (rf(1,a1)*eFrame(3,a1) + rf(2,a1)*eFrame(6,a1) + &
1331 +             rf(3,a1)*eFrame(9,a1))
1332 +    endif
1333 +
1334 +    return
1335 +  end subroutine reaction_field_final
1336 +
1337 +  subroutine rf_correct_forces(atom1, atom2, d, rij, eFrame, taper, f, fpair)
1338 +
1339 +    integer, intent(in) :: atom1, atom2
1340 +    real(kind=dp), dimension(3), intent(in) :: d
1341 +    real(kind=dp), intent(in) :: rij, taper
1342 +    real( kind = dp ), dimension(9,nLocal) :: eFrame
1343 +    real( kind = dp ), dimension(3,nLocal) :: f
1344 +    real( kind = dp ), dimension(3), intent(inout) :: fpair
1345 +
1346 +    real (kind = dp), dimension(3) :: ul1
1347 +    real (kind = dp), dimension(3) :: ul2
1348 +    real (kind = dp) :: dtdr
1349 +    real (kind = dp) :: dudx, dudy, dudz, u1dotu2
1350 +    integer :: me1, me2, id1, id2
1351 +    real (kind = dp) :: mu1, mu2
1352 +
1353 +    integer :: localError
1354 +
1355 +    if (.not.preRFCalculated) then
1356 +       call setReactionFieldPrefactor()
1357 +    endif
1358 +
1359 +    if (rij.le.rrf) then
1360 +
1361 +       if (rij.lt.rt) then
1362 +          dtdr = 0.0d0
1363 +       else
1364 +          !         write(*,*) 'rf correct in taper region'
1365 +          dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
1366 +       endif
1367 +
1368 + #ifdef IS_MPI
1369 +       me1 = atid_Row(atom1)
1370 +       ul1(1) = eFrame_Row(3,atom1)
1371 +       ul1(2) = eFrame_Row(6,atom1)
1372 +       ul1(3) = eFrame_Row(9,atom1)
1373 +
1374 +       me2 = atid_Col(atom2)
1375 +       ul2(1) = eFrame_Col(3,atom2)
1376 +       ul2(2) = eFrame_Col(6,atom2)
1377 +       ul2(3) = eFrame_Col(9,atom2)
1378 + #else
1379 +       me1 = atid(atom1)
1380 +       ul1(1) = eFrame(3,atom1)
1381 +       ul1(2) = eFrame(6,atom1)
1382 +       ul1(3) = eFrame(9,atom1)
1383 +
1384 +       me2 = atid(atom2)
1385 +       ul2(1) = eFrame(3,atom2)
1386 +       ul2(2) = eFrame(6,atom2)
1387 +       ul2(3) = eFrame(9,atom2)
1388 + #endif
1389 +
1390 +       mu1 = getDipoleMoment(me1)
1391 +       mu2 = getDipoleMoment(me2)
1392 +
1393 +       u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
1394 +
1395 +       dudx = - preRF*mu1*mu2*u1dotu2*dtdr*d(1)/rij
1396 +       dudy = - preRF*mu1*mu2*u1dotu2*dtdr*d(2)/rij
1397 +       dudz = - preRF*mu1*mu2*u1dotu2*dtdr*d(3)/rij
1398 +
1399 + #ifdef IS_MPI
1400 +       f_Row(1,atom1) = f_Row(1,atom1) + dudx
1401 +       f_Row(2,atom1) = f_Row(2,atom1) + dudy
1402 +       f_Row(3,atom1) = f_Row(3,atom1) + dudz
1403 +
1404 +       f_Col(1,atom2) = f_Col(1,atom2) - dudx
1405 +       f_Col(2,atom2) = f_Col(2,atom2) - dudy
1406 +       f_Col(3,atom2) = f_Col(3,atom2) - dudz
1407 + #else
1408 +       f(1,atom1) = f(1,atom1) + dudx
1409 +       f(2,atom1) = f(2,atom1) + dudy
1410 +       f(3,atom1) = f(3,atom1) + dudz
1411 +
1412 +       f(1,atom2) = f(1,atom2) - dudx
1413 +       f(2,atom2) = f(2,atom2) - dudy
1414 +       f(3,atom2) = f(3,atom2) - dudz
1415 + #endif
1416 +
1417 + #ifdef IS_MPI
1418 +       id1 = AtomRowToGlobal(atom1)
1419 +       id2 = AtomColToGlobal(atom2)
1420 + #else
1421 +       id1 = atom1
1422 +       id2 = atom2
1423 + #endif
1424 +
1425 +       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1426 +
1427 +          fpair(1) = fpair(1) + dudx
1428 +          fpair(2) = fpair(2) + dudy
1429 +          fpair(3) = fpair(3) + dudz
1430 +
1431 +       endif
1432 +
1433 +    end if
1434 +    return
1435 +  end subroutine rf_correct_forces
1436 +
1437   end module electrostatic_module

Diff Legend

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