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 2715 by chrisfen, Sun Apr 16 02:51:16 2006 UTC vs.
Revision 3070 by chrisfen, Sun Oct 29 22:41:17 2006 UTC

# Line 76 | Line 76 | module electrostatic_module
76    !! This unit is also known affectionately as an esu centi-barn.
77    real(kind=dp), parameter :: pre14 = 69.13373_dp
78  
79 +  real(kind=dp), parameter :: zero = 0.0_dp
80 +  
81 +  !! conversions for the simulation box dipole moment
82 +  real(kind=dp), parameter :: chargeToC = 1.60217733e-19_dp
83 +  real(kind=dp), parameter :: angstromToM = 1.0e-10_dp
84 +  real(kind=dp), parameter :: debyeToCm = 3.33564095198e-30_dp
85 +
86 +  !! number of points for electrostatic splines
87 +  integer, parameter :: np = 100
88 +
89    !! variables to handle different summation methods for long-range
90    !! electrostatics:
91    integer, save :: summationMethod = NONE
# Line 85 | Line 95 | module electrostatic_module
95    real(kind=DP), save :: defaultCutoff2 = 0.0_DP
96    logical, save :: haveDefaultCutoff = .false.
97    real(kind=DP), save :: dampingAlpha = 0.0_DP
98 <  real(kind=DP), save :: alpha2 = 0.0_DP
98 >  real(kind=DP), save :: alpha2 = 0.0_DP
99 >  real(kind=DP), save :: alpha4 = 0.0_DP
100 >  real(kind=DP), save :: alpha6 = 0.0_DP
101 >  real(kind=DP), save :: alpha8 = 0.0_DP
102    logical, save :: haveDampingAlpha = .false.
103    real(kind=DP), save :: dielectric = 1.0_DP
104    logical, save :: haveDielectric = .false.
# Line 101 | Line 114 | module electrostatic_module
114    real(kind=dp), save :: rrfsq = 1.0_DP
115    real(kind=dp), save :: preRF = 0.0_DP
116    real(kind=dp), save :: preRF2 = 0.0_DP
117 <  real(kind=dp), save :: f0 = 1.0_DP
118 <  real(kind=dp), save :: f1 = 1.0_DP
119 <  real(kind=dp), save :: f2 = 0.0_DP
120 <  real(kind=dp), save :: f3 = 0.0_DP
121 <  real(kind=dp), save :: f4 = 0.0_DP
122 <  real(kind=dp), save :: f0c = 1.0_DP
123 <  real(kind=dp), save :: f1c = 1.0_DP
124 <  real(kind=dp), save :: f2c = 0.0_DP
125 <  real(kind=dp), save :: f3c = 0.0_DP
126 <  real(kind=dp), save :: f4c = 0.0_DP
117 >  real(kind=dp), save :: erfcVal = 1.0_DP
118 >  real(kind=dp), save :: derfcVal = 0.0_DP
119 >  type(cubicSpline), save :: erfcSpline
120 >  logical, save :: haveElectroSpline = .false.
121 >  real(kind=dp), save :: c1 = 1.0_DP
122 >  real(kind=dp), save :: c2 = 1.0_DP
123 >  real(kind=dp), save :: c3 = 0.0_DP
124 >  real(kind=dp), save :: c4 = 0.0_DP
125 >  real(kind=dp), save :: c5 = 0.0_DP
126 >  real(kind=dp), save :: c6 = 0.0_DP
127 >  real(kind=dp), save :: c1c = 1.0_DP
128 >  real(kind=dp), save :: c2c = 1.0_DP
129 >  real(kind=dp), save :: c3c = 0.0_DP
130 >  real(kind=dp), save :: c4c = 0.0_DP
131 >  real(kind=dp), save :: c5c = 0.0_DP
132 >  real(kind=dp), save :: c6c = 0.0_DP
133 >  real(kind=dp), save :: one_third = 1.0_DP / 3.0_DP
134  
135   #if defined(__IFC) || defined(__PGI)
136   ! error function for ifc version > 7.
137 <  double precision, external :: derfc
137 >  real(kind=dp), external :: erfc
138   #endif
139    
140    public :: setElectrostaticSummationMethod
# Line 122 | Line 142 | module electrostatic_module
142    public :: setElectrostaticCutoffRadius
143    public :: setDampingAlpha
144    public :: setReactionFieldDielectric
145 <  public :: buildElectroSplines
145 >  public :: buildElectroSpline
146    public :: newElectrostaticType
147    public :: setCharge
148    public :: setDipoleMoment
# Line 134 | Line 154 | module electrostatic_module
154    public :: destroyElectrostaticTypes
155    public :: self_self
156    public :: rf_self_excludes
157 +  public :: accumulate_box_dipole
158  
138 –
159    type :: Electrostatic
160       integer :: c_ident
161       logical :: is_Charge = .false.
# Line 151 | Line 171 | contains
171  
172    type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap
173  
174 +  logical, save :: hasElectrostaticMap
175 +
176   contains
177  
178    subroutine setElectrostaticSummationMethod(the_ESM)
# Line 183 | Line 205 | contains
205      real(kind=dp), intent(in) :: thisAlpha
206      dampingAlpha = thisAlpha
207      alpha2 = dampingAlpha*dampingAlpha
208 +    alpha4 = alpha2*alpha2
209 +    alpha6 = alpha4*alpha2
210 +    alpha8 = alpha4*alpha4
211      haveDampingAlpha = .true.
212    end subroutine setDampingAlpha
213    
# Line 192 | Line 217 | contains
217      haveDielectric = .true.
218    end subroutine setReactionFieldDielectric
219  
220 <  subroutine buildElectroSplines()
221 <  end subroutine buildElectroSplines
220 >  subroutine buildElectroSpline()
221 >    real( kind = dp ), dimension(np) :: xvals, yvals
222 >    real( kind = dp ) :: dx, rmin, rval
223 >    integer :: i
224  
225 +    rmin = 0.0_dp
226 +
227 +    dx = (defaultCutoff-rmin) / dble(np-1)
228 +    
229 +    do i = 1, np
230 +       rval = rmin + dble(i-1)*dx
231 +       xvals(i) = rval
232 +       yvals(i) = erfc(dampingAlpha*rval)
233 +    enddo
234 +
235 +    call newSpline(erfcSpline, xvals, yvals, .true.)
236 +
237 +    haveElectroSpline = .true.
238 +  end subroutine buildElectroSpline
239 +
240    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
241         is_SplitDipole, is_Quadrupole, is_Tap, status)
242  
# Line 222 | Line 264 | contains
264            return
265         end if
266  
267 <       if (.not. allocated(ElectrostaticMap)) then
226 <          allocate(ElectrostaticMap(nAtypes))
227 <       endif
267 >       allocate(ElectrostaticMap(nAtypes))
268  
269      end if
270  
# Line 242 | Line 282 | contains
282      ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
283      ElectrostaticMap(myATID)%is_Tap = is_Tap
284  
285 +    hasElectrostaticMap = .true.
286 +
287    end subroutine newElectrostaticType
288  
289    subroutine setCharge(c_ident, charge, status)
# Line 253 | Line 295 | contains
295      status = 0
296      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
297  
298 <    if (.not.allocated(ElectrostaticMap)) then
298 >    if (.not.hasElectrostaticMap) then
299         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setCharge!")
300         status = -1
301         return
# Line 283 | Line 325 | contains
325      status = 0
326      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
327  
328 <    if (.not.allocated(ElectrostaticMap)) then
328 >    if (.not.hasElectrostaticMap) then
329         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setDipoleMoment!")
330         status = -1
331         return
# Line 313 | Line 355 | contains
355      status = 0
356      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
357  
358 <    if (.not.allocated(ElectrostaticMap)) then
358 >    if (.not.hasElectrostaticMap) then
359         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setSplitDipoleDistance!")
360         status = -1
361         return
# Line 343 | Line 385 | contains
385      status = 0
386      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
387  
388 <    if (.not.allocated(ElectrostaticMap)) then
388 >    if (.not.hasElectrostaticMap) then
389         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setQuadrupoleMoments!")
390         status = -1
391         return
# Line 374 | Line 416 | contains
416      integer :: localError
417      real(kind=dp) :: c
418  
419 <    if (.not.allocated(ElectrostaticMap)) then
419 >    if (.not.hasElectrostaticMap) then
420         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
421         return
422      end if
# Line 392 | Line 434 | contains
434      integer :: localError
435      real(kind=dp) :: dm
436  
437 <    if (.not.allocated(ElectrostaticMap)) then
437 >    if (.not.hasElectrostaticMap) then
438         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
439         return
440      end if
# Line 411 | Line 453 | contains
453         call handleError("checkSummationMethod", "no Default Cutoff set!")
454      endif
455  
456 <    rcuti = 1.0d0 / defaultCutoff
456 >    rcuti = 1.0_dp / defaultCutoff
457      rcuti2 = rcuti*rcuti
458      rcuti3 = rcuti2*rcuti
459      rcuti4 = rcuti2*rcuti2
# Line 426 | Line 468 | contains
468         endif
469  
470         constEXP = exp(-alpha2*defaultCutoff2)
471 <       invRootPi = 0.56418958354775628695d0
472 <       alphaPi = 2.0d0*dampingAlpha*invRootPi
473 <       f0c = derfc(dampingAlpha*defaultCutoff)
474 <       f1c = alphaPi*defaultCutoff*constEXP + f0c
475 <       f2c = alphaPi*2.0d0*alpha2*constEXP
476 <       f3c = alphaPi*2.0d0*alpha2*constEXP*defaultCutoff2*defaultCutoff
471 >       invRootPi = 0.56418958354775628695_dp
472 >       alphaPi = 2.0_dp*dampingAlpha*invRootPi
473 >
474 >       c1c = erfc(dampingAlpha*defaultCutoff) * rcuti
475 >       c2c = alphaPi*constEXP*rcuti + c1c*rcuti
476 >       c3c = 2.0_dp*alphaPi*alpha2 + 3.0_dp*c2c*rcuti
477 >       c4c = 4.0_dp*alphaPi*alpha4 + 5.0_dp*c3c*rcuti2
478 >       c5c = 8.0_dp*alphaPi*alpha6 + 7.0_dp*c4c*rcuti2
479 >       c6c = 16.0_dp*alphaPi*alpha8 + 9.0_dp*c5c*rcuti2
480 >    else
481 >       c1c = rcuti
482 >       c2c = c1c*rcuti
483 >       c3c = 3.0_dp*c2c*rcuti
484 >       c4c = 5.0_dp*c3c*rcuti2
485 >       c5c = 7.0_dp*c4c*rcuti2
486 >       c6c = 9.0_dp*c5c*rcuti2
487      endif
488  
489      if (summationMethod .eq. REACTION_FIELD) then
490         if (haveDielectric) then
491            defaultCutoff2 = defaultCutoff*defaultCutoff
492 <          preRF = (dielectric-1.0d0) / &
493 <               ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
494 <          preRF2 = 2.0d0*preRF
492 >          preRF = (dielectric-1.0_dp) / &
493 >               ((2.0_dp*dielectric+1.0_dp)*defaultCutoff2*defaultCutoff)
494 >          preRF2 = 2.0_dp*preRF
495         else
496            call handleError("checkSummationMethod", "Dielectric not set")
497         endif
498        
499      endif
500  
501 +    if (.not.haveElectroSpline) then
502 +       call buildElectroSpline()
503 +    end if
504 +
505      summationMethodChecked = .true.
506    end subroutine checkSummationMethod
507  
# Line 493 | Line 549 | contains
549      real (kind=dp) :: varEXP
550      real (kind=dp) :: pot_term
551      real (kind=dp) :: preVal, rfVal
552 <    real (kind=dp) :: f13, f134
552 >    real (kind=dp) :: c2ri, c3ri, c4rij
553 >    real (kind=dp) :: cti3, ctj3, ctidotj
554 >    real (kind=dp) :: preSw, preSwSc
555 >    real (kind=dp) :: xhatdot2, yhatdot2, zhatdot2
556 >    real (kind=dp) :: xhatc4, yhatc4, zhatc4
557  
498 –    if (.not.allocated(ElectrostaticMap)) then
499 –       call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
500 –       return
501 –    end if
502 –
558      if (.not.summationMethodChecked) then
559         call checkSummationMethod()
560      endif
# Line 514 | Line 569 | contains
569  
570      !! some variables we'll need independent of electrostatic type:
571  
572 <    riji = 1.0d0 / rij
572 >    riji = 1.0_dp / rij
573    
574      xhat = d(1) * riji
575      yhat = d(2) * riji
# Line 553 | Line 608 | contains
608         if (i_is_SplitDipole) then
609            d_i = ElectrostaticMap(me1)%split_dipole_distance
610         endif
611 <
611 >       duduz_i = zero
612      endif
613  
614      if (i_is_Quadrupole) then
# Line 584 | Line 639 | contains
639         cx_i = ux_i(1)*xhat + ux_i(2)*yhat + ux_i(3)*zhat
640         cy_i = uy_i(1)*xhat + uy_i(2)*yhat + uy_i(3)*zhat
641         cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
642 +       dudux_i = zero
643 +       duduy_i = zero
644 +       duduz_i = zero
645      endif
646  
647      if (j_is_Charge) then
# Line 606 | Line 664 | contains
664         if (j_is_SplitDipole) then
665            d_j = ElectrostaticMap(me2)%split_dipole_distance
666         endif
667 +       duduz_j = zero
668      endif
669  
670      if (j_is_Quadrupole) then
# Line 636 | Line 695 | contains
695         cx_j = ux_j(1)*xhat + ux_j(2)*yhat + ux_j(3)*zhat
696         cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
697         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
698 +       dudux_j = zero
699 +       duduy_j = zero
700 +       duduz_j = zero
701      endif
702    
703 <    epot = 0.0_dp
704 <    dudx = 0.0_dp
705 <    dudy = 0.0_dp
706 <    dudz = 0.0_dp
703 >    epot = zero
704 >    dudx = zero
705 >    dudy = zero
706 >    dudz = zero  
707  
646 –    dudux_i = 0.0_dp
647 –    duduy_i = 0.0_dp
648 –    duduz_i = 0.0_dp
649 –
650 –    dudux_j = 0.0_dp
651 –    duduy_j = 0.0_dp
652 –    duduz_j = 0.0_dp
653 –
708      if (i_is_Charge) then
709  
710         if (j_is_Charge) then
711            if (screeningMethod .eq. DAMPED) then
712 <             f0 = derfc(dampingAlpha*rij)
713 <             varEXP = exp(-alpha2*rij*rij)
714 <             f1 = alphaPi*rij*varEXP + f0
712 >             ! assemble the damping variables
713 >             call lookupUniformSpline1d(erfcSpline, rij, erfcVal, derfcVal)
714 >             c1 = erfcVal*riji
715 >             c2 = (-derfcVal + c1)*riji
716 >          else
717 >             c1 = riji
718 >             c2 = c1*riji
719            endif
720  
721            preVal = pre11 * q_i * q_j
722  
723            if (summationMethod .eq. SHIFTED_POTENTIAL) then
724 <             vterm = preVal * (riji*f0 - rcuti*f0c)
724 >             vterm = preVal * (c1 - c1c)
725              
726 <             dudr  = -sw * preVal * riji * riji * f1
726 >             dudr  = -sw * preVal * c2
727    
728            elseif (summationMethod .eq. SHIFTED_FORCE) then
729 <             vterm = preVal * ( riji*f0 - rcuti*f0c + &
672 <                  f1c*rcuti2*(rij-defaultCutoff) )
729 >             vterm = preVal * ( c1 - c1c + c2c*(rij - defaultCutoff) )
730              
731 <             dudr  = -sw*preVal * (riji*riji*f1 - rcuti2*f1c)
731 >             dudr  = sw * preVal * (c2c - c2)
732    
733            elseif (summationMethod .eq. REACTION_FIELD) then
734               rfVal = preRF*rij*rij
735               vterm = preVal * ( riji + rfVal )
736              
737 <             dudr  = sw * preVal * ( 2.0d0*rfVal - riji )*riji
737 >             dudr  = sw * preVal * ( 2.0_dp*rfVal - riji )*riji
738    
739            else
740 <             vterm = preVal * riji*f0
740 >             vterm = preVal * riji*erfcVal
741              
742 <             dudr  = - sw * preVal * riji*riji*f1
742 >             dudr  = - sw * preVal * c2
743    
744            endif
745  
# Line 696 | Line 753 | contains
753         endif
754  
755         if (j_is_Dipole) then
756 <          if (screeningMethod .eq. DAMPED) then
700 <             f0 = derfc(dampingAlpha*rij)
701 <             varEXP = exp(-alpha2*rij*rij)
702 <             f1 = alphaPi*rij*varEXP + f0
703 <             f3 = alphaPi*2.0d0*alpha2*varEXP*rij*rij*rij
704 <          endif
705 <
756 >          ! pref is used by all the possible methods
757            pref = pre12 * q_i * mu_j
758 +          preSw = sw*pref
759  
760            if (summationMethod .eq. REACTION_FIELD) then
761               ri2 = riji * riji
# Line 713 | Line 765 | contains
765               vpair = vpair + vterm
766               epot = epot + sw*vterm
767              
768 <             !! this has a + sign in the () because the rij vector is
769 <             !! r_j - r_i and the charge-dipole potential takes the origin
770 <             !! as the point dipole, which is atom j in this case.
771 <            
772 <             dudx = dudx - sw*pref*( ri3*(uz_j(1) - 3.0d0*ct_j*xhat) - &
773 <                                     preRF2*uz_j(1) )
774 <             dudy = dudy - sw*pref*( ri3*(uz_j(2) - 3.0d0*ct_j*yhat) - &
775 <                                     preRF2*uz_j(2) )
776 <             dudz = dudz - sw*pref*( ri3*(uz_j(3) - 3.0d0*ct_j*zhat) - &
725 <                                     preRF2*uz_j(3) )        
726 <             duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
727 <             duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
728 <             duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
768 >             dudx = dudx - preSw*( ri3*(uz_j(1) - 3.0_dp*ct_j*xhat) - &
769 >                  preRF2*uz_j(1) )
770 >             dudy = dudy - preSw*( ri3*(uz_j(2) - 3.0_dp*ct_j*yhat) - &
771 >                  preRF2*uz_j(2) )
772 >             dudz = dudz - preSw*( ri3*(uz_j(3) - 3.0_dp*ct_j*zhat) - &
773 >                  preRF2*uz_j(3) )        
774 >             duduz_j(1) = duduz_j(1) - preSw * xhat * ( ri2 - preRF2*rij )
775 >             duduz_j(2) = duduz_j(2) - preSw * yhat * ( ri2 - preRF2*rij )
776 >             duduz_j(3) = duduz_j(3) - preSw * zhat * ( ri2 - preRF2*rij )
777  
778            else
779 +             ! determine the inverse r used if we have split dipoles
780               if (j_is_SplitDipole) then
781                  BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
782                  ri = 1.0_dp / BigR
# Line 736 | Line 785 | contains
785                  ri = riji
786                  scale = 1.0_dp
787               endif
788 <            
740 <             ri2 = ri * ri
741 <             ri3 = ri2 * ri
788 >
789               sc2 = scale * scale
790  
791 <             pot_term =  ri2 * scale * f1
792 <             vterm = - pref * ct_j * pot_term
791 >             if (screeningMethod .eq. DAMPED) then
792 >                ! assemble the damping variables
793 >                call lookupUniformSpline1d(erfcSpline, rij, erfcVal, derfcVal)
794 >                c1 = erfcVal*ri
795 >                c2 = (-derfcVal + c1)*ri
796 >                c3 = -2.0_dp*derfcVal*alpha2 + 3.0_dp*c2*ri
797 >             else
798 >                c1 = ri
799 >                c2 = c1*ri
800 >                c3 = 3.0_dp*c2*ri
801 >             endif
802 >            
803 >             c2ri = c2*ri
804 >
805 >             ! calculate the potential
806 >             pot_term =  scale * c2
807 >             vterm = -pref * ct_j * pot_term
808               vpair = vpair + vterm
809               epot = epot + sw*vterm
810              
811 <             !! this has a + sign in the () because the rij vector is
812 <             !! r_j - r_i and the charge-dipole potential takes the origin
813 <             !! as the point dipole, which is atom j in this case.
814 <            
753 <             dudx = dudx - sw*pref * ri3 * ( uz_j(1)*f1 - &
754 <                  ct_j*xhat*sc2*( 3.0d0*f1 + f3 ) )
755 <             dudy = dudy - sw*pref * ri3 * ( uz_j(2)*f1 - &
756 <                  ct_j*yhat*sc2*( 3.0d0*f1 + f3 ) )
757 <             dudz = dudz - sw*pref * ri3 * ( uz_j(3)*f1 - &
758 <                  ct_j*zhat*sc2*( 3.0d0*f1 + f3 ) )
811 >             ! calculate derivatives for forces and torques
812 >             dudx = dudx - preSw*( uz_j(1)*c2ri - ct_j*xhat*sc2*c3 )
813 >             dudy = dudy - preSw*( uz_j(2)*c2ri - ct_j*yhat*sc2*c3 )
814 >             dudz = dudz - preSw*( uz_j(3)*c2ri - ct_j*zhat*sc2*c3 )
815                            
816 <             duduz_j(1) = duduz_j(1) - sw*pref * pot_term * xhat
817 <             duduz_j(2) = duduz_j(2) - sw*pref * pot_term * yhat
818 <             duduz_j(3) = duduz_j(3) - sw*pref * pot_term * zhat
816 >             duduz_j(1) = duduz_j(1) - preSw * pot_term * xhat
817 >             duduz_j(2) = duduz_j(2) - preSw * pot_term * yhat
818 >             duduz_j(3) = duduz_j(3) - preSw * pot_term * zhat
819  
820            endif
821         endif
822  
823         if (j_is_Quadrupole) then
824 <          if (screeningMethod .eq. DAMPED) then
769 <             f0 = derfc(dampingAlpha*rij)
770 <             varEXP = exp(-alpha2*rij*rij)
771 <             f1 = alphaPi*rij*varEXP + f0
772 <             f2 = alphaPi*2.0d0*alpha2*varEXP
773 <             f3 = f2*rij*rij*rij
774 <             f4 = 2.0d0*alpha2*f2*rij
775 <          endif
776 <
777 <          ri2 = riji * riji
778 <          ri3 = ri2 * riji
779 <          ri4 = ri2 * ri2
824 >          ! first precalculate some necessary variables
825            cx2 = cx_j * cx_j
826            cy2 = cy_j * cy_j
827            cz2 = cz_j * cz_j
828 +          pref =  pre14 * q_i * one_third
829 +          
830 +          if (screeningMethod .eq. DAMPED) then
831 +             ! assemble the damping variables
832 +             call lookupUniformSpline1d(erfcSpline, rij, erfcVal, derfcVal)
833 +             c1 = erfcVal*riji
834 +             c2 = (-derfcVal + c1)*riji
835 +             c3 = -2.0_dp*derfcVal*alpha2 + 3.0_dp*c2*riji
836 +             c4 = -4.0_dp*derfcVal*alpha4 + 5.0_dp*c3*riji*riji
837 +          else
838 +             c1 = riji
839 +             c2 = c1*riji
840 +             c3 = 3.0_dp*c2*riji
841 +             c4 = 5.0_dp*c3*riji*riji
842 +          endif
843  
844 <          pref =  pre14 * q_i / 3.0_dp
845 <          pot_term = ri3*(qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
846 <               qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
847 <               qzz_j * (3.0_dp*cz2 - 1.0_dp))
848 <          vterm = pref * (pot_term*f1 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f2)
844 >          ! precompute variables for convenience
845 >          preSw = sw*pref
846 >          c2ri = c2*riji
847 >          c3ri = c3*riji
848 >          c4rij = c4*rij
849 >          xhatdot2 = 2.0_dp*xhat*c3
850 >          yhatdot2 = 2.0_dp*yhat*c3
851 >          zhatdot2 = 2.0_dp*zhat*c3
852 >          xhatc4 = xhat*c4rij
853 >          yhatc4 = yhat*c4rij
854 >          zhatc4 = zhat*c4rij
855 >
856 >          ! calculate the potential
857 >          pot_term = ( qxx_j*(cx2*c3 - c2ri) + qyy_j*(cy2*c3 - c2ri) + &
858 >               qzz_j*(cz2*c3 - c2ri) )
859 >          vterm = pref * pot_term
860            vpair = vpair + vterm
861            epot = epot + sw*vterm
862 +
863 +          ! calculate derivatives for the forces and torques
864 +          dudx = dudx - preSw * ( &
865 +               qxx_j*(cx2*xhatc4 - (2.0_dp*cx_j*ux_j(1) + xhat)*c3ri) + &
866 +               qyy_j*(cy2*xhatc4 - (2.0_dp*cy_j*uy_j(1) + xhat)*c3ri) + &
867 +               qzz_j*(cz2*xhatc4 - (2.0_dp*cz_j*uz_j(1) + xhat)*c3ri) )
868 +          dudy = dudy - preSw * ( &
869 +               qxx_j*(cx2*yhatc4 - (2.0_dp*cx_j*ux_j(2) + yhat)*c3ri) + &
870 +               qyy_j*(cy2*yhatc4 - (2.0_dp*cy_j*uy_j(2) + yhat)*c3ri) + &
871 +               qzz_j*(cz2*yhatc4 - (2.0_dp*cz_j*uz_j(2) + yhat)*c3ri) )
872 +          dudz = dudz - preSw * ( &
873 +               qxx_j*(cx2*zhatc4 - (2.0_dp*cx_j*ux_j(3) + zhat)*c3ri) + &
874 +               qyy_j*(cy2*zhatc4 - (2.0_dp*cy_j*uy_j(3) + zhat)*c3ri) + &
875 +               qzz_j*(cz2*zhatc4 - (2.0_dp*cz_j*uz_j(3) + zhat)*c3ri) )
876            
877 <          dudx = dudx - sw*pref*pot_term*riji*xhat*(5.0d0*f1 + f3) + &
878 <               sw*pref*ri4 * ( &
879 <               qxx_j*(2.0_dp*cx_j*ux_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
795 <               qyy_j*(2.0_dp*cy_j*uy_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
796 <               qzz_j*(2.0_dp*cz_j*uz_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) ) &
797 <               + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
798 <          dudy = dudy - sw*pref*pot_term*riji*yhat*(5.0d0*f1 + f3) + &
799 <               sw*pref*ri4 * ( &
800 <               qxx_j*(2.0_dp*cx_j*ux_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
801 <               qyy_j*(2.0_dp*cy_j*uy_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
802 <               qzz_j*(2.0_dp*cz_j*uz_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) ) &
803 <               + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
804 <          dudz = dudz - sw*pref*pot_term*riji*zhat*(5.0d0*f1 + f3) + &
805 <               sw*pref*ri4 * ( &
806 <               qxx_j*(2.0_dp*cx_j*ux_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
807 <               qyy_j*(2.0_dp*cy_j*uy_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
808 <               qzz_j*(2.0_dp*cz_j*uz_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) ) &
809 <               + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
877 >          dudux_j(1) = dudux_j(1) + preSw*(qxx_j*cx_j*xhatdot2)
878 >          dudux_j(2) = dudux_j(2) + preSw*(qxx_j*cx_j*yhatdot2)
879 >          dudux_j(3) = dudux_j(3) + preSw*(qxx_j*cx_j*zhatdot2)
880            
881 <          dudux_j(1) = dudux_j(1) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*xhat) &
882 <               * (3.0d0*f1 + f3) )
883 <          dudux_j(2) = dudux_j(2) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*yhat) &
814 <               * (3.0d0*f1 + f3) )
815 <          dudux_j(3) = dudux_j(3) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*zhat) &
816 <               * (3.0d0*f1 + f3) )
881 >          duduy_j(1) = duduy_j(1) + preSw*(qyy_j*cy_j*xhatdot2)
882 >          duduy_j(2) = duduy_j(2) + preSw*(qyy_j*cy_j*yhatdot2)
883 >          duduy_j(3) = duduy_j(3) + preSw*(qyy_j*cy_j*zhatdot2)
884            
885 <          duduy_j(1) = duduy_j(1) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*xhat) &
886 <               * (3.0d0*f1 + f3) )
887 <          duduy_j(2) = duduy_j(2) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*yhat) &
888 <               * (3.0d0*f1 + f3) )
822 <          duduy_j(3) = duduy_j(3) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*zhat) &
823 <               * (3.0d0*f1 + f3) )
824 <          
825 <          duduz_j(1) = duduz_j(1) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*xhat) &
826 <               * (3.0d0*f1 + f3) )
827 <          duduz_j(2) = duduz_j(2) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*yhat) &
828 <               * (3.0d0*f1 + f3) )
829 <          duduz_j(3) = duduz_j(3) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*zhat) &
830 <               * (3.0d0*f1 + f3) )
885 >          duduz_j(1) = duduz_j(1) + preSw*(qzz_j*cz_j*xhatdot2)
886 >          duduz_j(2) = duduz_j(2) + preSw*(qzz_j*cz_j*yhatdot2)
887 >          duduz_j(3) = duduz_j(3) + preSw*(qzz_j*cz_j*zhatdot2)
888 >
889            
890         endif
891      endif
# Line 835 | Line 893 | contains
893      if (i_is_Dipole) then
894  
895         if (j_is_Charge) then
896 <          if (screeningMethod .eq. DAMPED) then
897 <             f0 = derfc(dampingAlpha*rij)
898 <             varEXP = exp(-alpha2*rij*rij)
841 <             f1 = alphaPi*rij*varEXP + f0
842 <             f3 = alphaPi*2.0d0*alpha2*varEXP*rij*rij*rij
843 <          endif
844 <          
845 <          pref = pre12 * q_j * mu_i
846 <          
847 <          if (summationMethod .eq. SHIFTED_POTENTIAL) then
848 <             ri2 = riji * riji
849 <             ri3 = ri2 * riji
850 <            
851 <             pot_term = ri2*f1 - rcuti2*f1c
852 <             vterm = pref * ct_i * pot_term
853 <             vpair = vpair + vterm
854 <             epot = epot + sw*vterm
855 <            
856 <             dudx = dudx + sw*pref*( ri3*(uz_i(1)*f1-ct_i*xhat*(3.0d0*f1+f3)) )
857 <             dudy = dudy + sw*pref*( ri3*(uz_i(2)*f1-ct_i*yhat*(3.0d0*f1+f3)) )
858 <             dudz = dudz + sw*pref*( ri3*(uz_i(3)*f1-ct_i*zhat*(3.0d0*f1+f3)) )
859 <            
860 <             duduz_i(1) = duduz_i(1) + sw*pref * xhat * pot_term
861 <             duduz_i(2) = duduz_i(2) + sw*pref * yhat * pot_term
862 <             duduz_i(3) = duduz_i(3) + sw*pref * zhat * pot_term
896 >          ! variables used by all the methods
897 >          pref = pre12 * q_j * mu_i                      
898 >          preSw = sw*pref
899  
900 <          elseif (summationMethod .eq. SHIFTED_FORCE) then
865 <             ri2 = riji * riji
866 <             ri3 = ri2 * riji
900 >          if (summationMethod .eq. REACTION_FIELD) then
901  
868 –             !! might need a -(f1c-f0c) or dct_i/dr in the derivative term...
869 –             pot_term = ri2*f1 - rcuti2*f1c + &
870 –                  (2.0d0*rcuti3*f1c + f2c)*( rij - defaultCutoff )
871 –             vterm = pref * ct_i * pot_term
872 –             vpair = vpair + vterm
873 –             epot = epot + sw*vterm
874 –            
875 –             dudx = dudx + sw*pref*( ri3*(uz_i(1)*f1-ct_i*xhat*(3.0d0*f1+f3)) &
876 –                  - rcuti3*(uz_i(1)*f1c-ct_i*xhat*(3.0d0*f1c+f3c)) )
877 –             dudy = dudy + sw*pref*( ri3*(uz_i(2)*f1-ct_i*yhat*(3.0d0*f1+f3)) &
878 –                  - rcuti3*(uz_i(1)*f1c-ct_i*xhat*(3.0d0*f1c+f3c)) )
879 –             dudz = dudz + sw*pref*( ri3*(uz_i(3)*f1-ct_i*zhat*(3.0d0*f1+f3)) &
880 –                  - rcuti3*(uz_i(1)*f1c-ct_i*xhat*(3.0d0*f1c+f3c)) )
881 –            
882 –             duduz_i(1) = duduz_i(1) + sw*pref * xhat * pot_term
883 –             duduz_i(2) = duduz_i(2) + sw*pref * yhat * pot_term
884 –             duduz_i(3) = duduz_i(3) + sw*pref * zhat * pot_term
885 –
886 –          elseif (summationMethod .eq. REACTION_FIELD) then
902               ri2 = riji * riji
903               ri3 = ri2 * riji
904  
# Line 891 | Line 906 | contains
906               vpair = vpair + vterm
907               epot = epot + sw*vterm
908              
909 <             dudx = dudx + sw*pref * ( ri3*(uz_i(1) - 3.0d0*ct_i*xhat) - &
909 >             dudx = dudx + preSw * ( ri3*(uz_i(1) - 3.0_dp*ct_i*xhat) - &
910                    preRF2*uz_i(1) )
911 <             dudy = dudy + sw*pref * ( ri3*(uz_i(2) - 3.0d0*ct_i*yhat) - &
911 >             dudy = dudy + preSw * ( ri3*(uz_i(2) - 3.0_dp*ct_i*yhat) - &
912                    preRF2*uz_i(2) )
913 <             dudz = dudz + sw*pref * ( ri3*(uz_i(3) - 3.0d0*ct_i*zhat) - &
913 >             dudz = dudz + preSw * ( ri3*(uz_i(3) - 3.0_dp*ct_i*zhat) - &
914                    preRF2*uz_i(3) )
915              
916 <             duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF2*rij )
917 <             duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF2*rij )
918 <             duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF2*rij )
916 >             duduz_i(1) = duduz_i(1) + preSw * xhat * ( ri2 - preRF2*rij )
917 >             duduz_i(2) = duduz_i(2) + preSw * yhat * ( ri2 - preRF2*rij )
918 >             duduz_i(3) = duduz_i(3) + preSw * zhat * ( ri2 - preRF2*rij )
919  
920            else
921 +             ! determine inverse r if we are using split dipoles
922               if (i_is_SplitDipole) then
923                  BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
924                  ri = 1.0_dp / BigR
# Line 911 | Line 927 | contains
927                  ri = riji
928                  scale = 1.0_dp
929               endif
930 <            
915 <             ri2 = ri * ri
916 <             ri3 = ri2 * ri
930 >
931               sc2 = scale * scale
932 +              
933 +             if (screeningMethod .eq. DAMPED) then
934 +                ! assemble the damping variables
935 +                call lookupUniformSpline1d(erfcSpline, rij, erfcVal, derfcVal)
936 +                c1 = erfcVal*ri
937 +                c2 = (-derfcVal + c1)*ri
938 +                c3 = -2.0_dp*derfcVal*alpha2 + 3.0_dp*c2*ri
939 +             else
940 +                c1 = ri
941 +                c2 = c1*ri
942 +                c3 = 3.0_dp*c2*ri
943 +             endif
944 +            
945 +             c2ri = c2*ri
946  
947 <             pot_term = ri2 * f1 * scale
947 >             ! calculate the potential
948 >             pot_term = c2 * scale
949               vterm = pref * ct_i * pot_term
950               vpair = vpair + vterm
951               epot = epot + sw*vterm
952 +
953 +             ! calculate derivatives for the forces and torques
954 +             dudx = dudx + preSw * ( uz_i(1)*c2ri - ct_i*xhat*sc2*c3 )
955 +             dudy = dudy + preSw * ( uz_i(2)*c2ri - ct_i*yhat*sc2*c3 )
956 +             dudz = dudz + preSw * ( uz_i(3)*c2ri - ct_i*zhat*sc2*c3 )
957 +
958 +             duduz_i(1) = duduz_i(1) + preSw * pot_term * xhat
959 +             duduz_i(2) = duduz_i(2) + preSw * pot_term * yhat
960 +             duduz_i(3) = duduz_i(3) + preSw * pot_term * zhat
961              
924 –             dudx = dudx + sw*pref * ri3 * ( uz_i(1)*f1 - &
925 –                  ct_i*xhat*sc2*( 3.0d0*f1 + f3 ) )
926 –             dudy = dudy + sw*pref * ri3 * ( uz_i(2)*f1 - &
927 –                  ct_i*yhat*sc2*( 3.0d0*f1 + f3 ) )
928 –             dudz = dudz + sw*pref * ri3 * ( uz_i(3)*f1 - &
929 –                  ct_i*zhat*sc2*( 3.0d0*f1 + f3 ) )
930 –            
931 –             duduz_i(1) = duduz_i(1) + sw*pref * pot_term * xhat
932 –             duduz_i(2) = duduz_i(2) + sw*pref * pot_term * yhat
933 –             duduz_i(3) = duduz_i(3) + sw*pref * pot_term * zhat
962            endif
963         endif
964        
965         if (j_is_Dipole) then
966 <          if (screeningMethod .eq. DAMPED) then
939 <             f0 = derfc(dampingAlpha*rij)
940 <             varEXP = exp(-alpha2*rij*rij)
941 <             f1 = alphaPi*rij*varEXP + f0
942 <             f2 = alphaPi*2.0d0*alpha2*varEXP
943 <             f3 = f2*rij*rij*rij
944 <             f4 = 2.0d0*alpha2*f3*rij*rij
945 <          endif
946 <
966 >          ! variables used by all methods
967            ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
948 –          
949 –          ri2 = riji * riji
950 –          ri3 = ri2 * riji
951 –          ri4 = ri2 * ri2
952 –          
968            pref = pre22 * mu_i * mu_j
969 +          preSw = sw*pref
970  
971            if (summationMethod .eq. REACTION_FIELD) then
972 <             vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
972 >             ri2 = riji * riji
973 >             ri3 = ri2 * riji
974 >             ri4 = ri2 * ri2
975 >
976 >             vterm = pref*( ri3*(ct_ij - 3.0_dp * ct_i * ct_j) - &
977                    preRF2*ct_ij )
978               vpair = vpair + vterm
979               epot = epot + sw*vterm
980              
981 <             a1 = 5.0d0 * ct_i * ct_j - ct_ij
981 >             a1 = 5.0_dp * ct_i * ct_j - ct_ij
982              
983 <             dudx = dudx + sw*pref*3.0d0*ri4 &
984 <                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
985 <             dudy = dudy + sw*pref*3.0d0*ri4 &
966 <                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
967 <             dudz = dudz + sw*pref*3.0d0*ri4 &
968 <                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
983 >             dudx = dudx + preSw*3.0_dp*ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
984 >             dudy = dudy + preSw*3.0_dp*ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
985 >             dudz = dudz + preSw*3.0_dp*ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
986              
987 <             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
987 >             duduz_i(1) = duduz_i(1) + preSw*(ri3*(uz_j(1)-3.0_dp*ct_j*xhat) &
988                    - preRF2*uz_j(1))
989 <             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
989 >             duduz_i(2) = duduz_i(2) + preSw*(ri3*(uz_j(2)-3.0_dp*ct_j*yhat) &
990                    - preRF2*uz_j(2))
991 <             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
991 >             duduz_i(3) = duduz_i(3) + preSw*(ri3*(uz_j(3)-3.0_dp*ct_j*zhat) &
992                    - preRF2*uz_j(3))
993 <             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
993 >             duduz_j(1) = duduz_j(1) + preSw*(ri3*(uz_i(1)-3.0_dp*ct_i*xhat) &
994                    - preRF2*uz_i(1))
995 <             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
995 >             duduz_j(2) = duduz_j(2) + preSw*(ri3*(uz_i(2)-3.0_dp*ct_i*yhat) &
996                    - preRF2*uz_i(2))
997 <             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
997 >             duduz_j(3) = duduz_j(3) + preSw*(ri3*(uz_i(3)-3.0_dp*ct_i*zhat) &
998                    - preRF2*uz_i(3))
999  
1000            else
# Line 999 | Line 1016 | contains
1016                     scale = 1.0_dp
1017                  endif
1018               endif
1019 <            
1019 >
1020 >             if (screeningMethod .eq. DAMPED) then
1021 >                ! assemble the damping variables
1022 >                call lookupUniformSpline1d(erfcSpline, rij, erfcVal, derfcVal)
1023 >                c1 = erfcVal*ri
1024 >                c2 = (-derfcVal + c1)*ri
1025 >                c3 = -2.0_dp*derfcVal*alpha2 + 3.0_dp*c2*ri
1026 >                c4 = -4.0_dp*derfcVal*alpha4 + 5.0_dp*c3*ri*ri
1027 >             else
1028 >                c1 = ri
1029 >                c2 = c1*ri
1030 >                c3 = 3.0_dp*c2*ri
1031 >                c4 = 5.0_dp*c3*ri*ri
1032 >             endif
1033 >
1034 >             ! precompute variables for convenience
1035               sc2 = scale * scale
1036 +             cti3 = ct_i*sc2*c3
1037 +             ctj3 = ct_j*sc2*c3
1038 +             ctidotj = ct_i * ct_j * sc2        
1039 +             preSwSc = preSw*scale
1040 +             c2ri = c2*ri
1041 +             c3ri = c3*ri
1042 +             c4rij = c4*rij
1043  
1044 <             pot_term = (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1045 <             vterm = pref * ( ri3*pot_term*f1 + (ct_i * ct_j)*f2 )
1044 >
1045 >             ! calculate the potential
1046 >             pot_term = (ct_ij*c2ri - ctidotj*c3)
1047 >             vterm = pref * pot_term
1048               vpair = vpair + vterm
1049               epot = epot + sw*vterm
1009 –            
1010 –             f13 = f1+f3
1011 –             f134 = f13 + f4
1012 –            
1013 – !!$             dudx = dudx + sw*pref * ( ri4*scale*( &
1014 – !!$                  3.0d0*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))*f1 &
1015 – !!$                  - pot_term*f3) &
1016 – !!$                  + 2.0d0*ct_i*ct_j*xhat*(ct_i*uz_j(1)+ct_j*uz_i(1))*f3 &
1017 – !!$                  + (ct_i * ct_j)*f4 )
1018 – !!$             dudy = dudy + sw*pref * ( ri4*scale*( &
1019 – !!$                  3.0d0*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))*f1 &
1020 – !!$                  - pot_term*f3) &
1021 – !!$                  + 2.0d0*ct_i*ct_j*yhat*(ct_i*uz_j(2)+ct_j*uz_i(2))*f3 &
1022 – !!$                  + (ct_i * ct_j)*f4 )
1023 – !!$             dudz = dudz + sw*pref * ( ri4*scale*( &
1024 – !!$                  3.0d0*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))*f1 &
1025 – !!$                  - pot_term*f3) &
1026 – !!$                  + 2.0d0*ct_i*ct_j*zhat*(ct_i*uz_j(3)+ct_j*uz_i(3))*f3 &
1027 – !!$                  + (ct_i * ct_j)*f4 )
1050  
1051 <             dudx = dudx + sw*pref * ( ri4*scale*( &
1052 <                  15.0d0*(ct_i * ct_j * sc2)*xhat*f134 - &
1053 <                  3.0d0*(ct_i*uz_j(1) + ct_j*uz_i(1) + ct_ij*xhat)*f134) )
1054 <             dudy = dudy + sw*pref * ( ri4*scale*( &
1055 <                  15.0d0*(ct_i * ct_j * sc2)*yhat*f134 - &
1056 <                  3.0d0*(ct_i*uz_j(2) + ct_j*uz_i(2) + ct_ij*yhat)*f134) )
1057 <             dudz = dudz + sw*pref * ( ri4*scale*( &
1058 <                  15.0d0*(ct_i * ct_j * sc2)*zhat*f134 - &
1059 <                  3.0d0*(ct_i*uz_j(3) + ct_j*uz_i(3) + ct_ij*zhat)*f134) )
1051 >             ! calculate derivatives for the forces and torques
1052 >             dudx = dudx + preSwSc * ( ctidotj*xhat*c4rij - &
1053 >                  (ct_i*uz_j(1) + ct_j*uz_i(1) + ct_ij*xhat)*c3ri )
1054 >             dudy = dudy + preSwSc * ( ctidotj*yhat*c4rij - &
1055 >                  (ct_i*uz_j(2) + ct_j*uz_i(2) + ct_ij*yhat)*c3ri )
1056 >             dudz = dudz + preSwSc * ( ctidotj*zhat*c4rij - &
1057 >                  (ct_i*uz_j(3) + ct_j*uz_i(3) + ct_ij*zhat)*c3ri )
1058 >
1059 >             duduz_i(1) = duduz_i(1) + preSw * ( uz_j(1)*c2ri - ctj3*xhat )
1060 >             duduz_i(2) = duduz_i(2) + preSw * ( uz_j(2)*c2ri - ctj3*yhat )
1061 >             duduz_i(3) = duduz_i(3) + preSw * ( uz_j(3)*c2ri - ctj3*zhat )
1062              
1063 <             duduz_i(1) = duduz_i(1) + sw*pref * &
1064 <                  ( ri3*(uz_j(1) - 3.0d0*ct_j*xhat*sc2)*f1 + (ct_j*xhat)*f2 )
1065 <             duduz_i(2) = duduz_i(2) + sw*pref * &
1066 <                  ( ri3*(uz_j(2) - 3.0d0*ct_j*yhat*sc2)*f1 + (ct_j*yhat)*f2 )
1043 <             duduz_i(3) = duduz_i(3) + sw*pref * &
1044 <                  ( ri3*(uz_j(3) - 3.0d0*ct_j*zhat*sc2)*f1 + (ct_j*zhat)*f2 )
1045 <            
1046 <             duduz_j(1) = duduz_j(1) + sw*pref * &
1047 <                  ( ri3*(uz_i(1) - 3.0d0*ct_i*xhat*sc2)*f1 + (ct_i*xhat)*f2 )
1048 <             duduz_j(2) = duduz_j(2) + sw*pref * &
1049 <                  ( ri3*(uz_i(2) - 3.0d0*ct_i*yhat*sc2)*f1 + (ct_i*yhat)*f2 )
1050 <             duduz_j(3) = duduz_j(3) + sw*pref * &
1051 <                  ( ri3*(uz_i(3) - 3.0d0*ct_i*zhat*sc2)*f1 + (ct_i*zhat)*f2 )
1063 >             duduz_j(1) = duduz_j(1) + preSw * ( uz_i(1)*c2ri - cti3*xhat )
1064 >             duduz_j(2) = duduz_j(2) + preSw * ( uz_i(2)*c2ri - cti3*yhat )
1065 >             duduz_j(3) = duduz_j(3) + preSw * ( uz_i(3)*c2ri - cti3*zhat )
1066 >
1067            endif
1068         endif
1069      endif
1070  
1071      if (i_is_Quadrupole) then
1072         if (j_is_Charge) then
1073 <          if (screeningMethod .eq. DAMPED) then
1059 <             f0 = derfc(dampingAlpha*rij)
1060 <             varEXP = exp(-alpha2*rij*rij)
1061 <             f1 = alphaPi*rij*varEXP + f0
1062 <             f2 = alphaPi*2.0d0*alpha2*varEXP
1063 <             f3 = f2*rij*rij*rij
1064 <             f4 = 2.0d0*alpha2*f2*rij
1065 <          endif
1066 <
1067 <          ri2 = riji * riji
1068 <          ri3 = ri2 * riji
1069 <          ri4 = ri2 * ri2
1073 >          ! precompute some necessary variables
1074            cx2 = cx_i * cx_i
1075            cy2 = cy_i * cy_i
1076            cz2 = cz_i * cz_i
1077 +          pref = pre14 * q_j * one_third
1078  
1079 <          pref = pre14 * q_j / 3.0_dp
1080 <          pot_term = ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1081 <                            qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1082 <                            qzz_i * (3.0_dp*cz2 - 1.0_dp))
1083 <          vterm = pref * (pot_term*f1 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f2)
1079 >          if (screeningMethod .eq. DAMPED) then
1080 >             ! assemble the damping variables
1081 >             call lookupUniformSpline1d(erfcSpline, rij, erfcVal, derfcVal)
1082 >             c1 = erfcVal*riji
1083 >             c2 = (-derfcVal + c1)*riji
1084 >             c3 = -2.0_dp*derfcVal*alpha2 + 3.0_dp*c2*riji
1085 >             c4 = -4.0_dp*derfcVal*alpha4 + 5.0_dp*c3*riji*riji
1086 >          else
1087 >             c1 = riji
1088 >             c2 = c1*riji
1089 >             c3 = 3.0_dp*c2*riji
1090 >             c4 = 5.0_dp*c3*riji*riji
1091 >          endif
1092 >          
1093 >          ! precompute some variables for convenience
1094 >          preSw = sw*pref
1095 >          c2ri = c2*riji
1096 >          c3ri = c3*riji
1097 >          c4rij = c4*rij
1098 >          xhatdot2 = 2.0_dp*xhat*c3
1099 >          yhatdot2 = 2.0_dp*yhat*c3
1100 >          zhatdot2 = 2.0_dp*zhat*c3
1101 >          xhatc4 = xhat*c4rij
1102 >          yhatc4 = yhat*c4rij
1103 >          zhatc4 = zhat*c4rij
1104 >
1105 >          ! calculate the potential
1106 >          pot_term = ( qxx_i * (cx2*c3 - c2ri) + qyy_i * (cy2*c3 - c2ri) + &
1107 >               qzz_i * (cz2*c3 - c2ri) )
1108 >
1109 >          vterm = pref * pot_term
1110            vpair = vpair + vterm
1111            epot = epot + sw*vterm
1112 <          
1113 <          dudx = dudx - sw*pref*pot_term*riji*xhat*(5.0d0*f1 + f3) + &
1114 <               sw*pref*ri4 * ( &
1115 <               qxx_i*(2.0_dp*cx_i*ux_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
1116 <               qyy_i*(2.0_dp*cy_i*uy_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
1117 <               qzz_i*(2.0_dp*cz_i*uz_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) ) &
1118 <               + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1119 <          dudy = dudy - sw*pref*pot_term*riji*yhat*(5.0d0*f1 + f3) + &
1120 <               sw*pref*ri4 * ( &
1121 <               qxx_i*(2.0_dp*cx_i*ux_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
1122 <               qyy_i*(2.0_dp*cy_i*uy_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
1123 <               qzz_i*(2.0_dp*cz_i*uz_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) ) &
1124 <               + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1125 <          dudz = dudz - sw*pref*pot_term*riji*zhat*(5.0d0*f1 + f3) + &
1095 <               sw*pref*ri4 * ( &
1096 <               qxx_i*(2.0_dp*cx_i*ux_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
1097 <               qyy_i*(2.0_dp*cy_i*uy_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
1098 <               qzz_i*(2.0_dp*cz_i*uz_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) ) &
1099 <               + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1100 <          
1101 <          dudux_i(1) = dudux_i(1) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*xhat) &
1102 <               * (3.0d0*f1 + f3) )
1103 <          dudux_i(2) = dudux_i(2) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*yhat) &
1104 <               * (3.0d0*f1 + f3) )
1105 <          dudux_i(3) = dudux_i(3) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*zhat) &
1106 <               * (3.0d0*f1 + f3) )
1112 >
1113 >          ! calculate the derivatives for the forces and torques
1114 >          dudx = dudx - preSw * ( &
1115 >               qxx_i*(cx2*xhatc4 - (2.0_dp*cx_i*ux_i(1) + xhat)*c3ri) + &
1116 >               qyy_i*(cy2*xhatc4 - (2.0_dp*cy_i*uy_i(1) + xhat)*c3ri) + &
1117 >               qzz_i*(cz2*xhatc4 - (2.0_dp*cz_i*uz_i(1) + xhat)*c3ri) )
1118 >          dudy = dudy - preSw * ( &
1119 >               qxx_i*(cx2*yhatc4 - (2.0_dp*cx_i*ux_i(2) + yhat)*c3ri) + &
1120 >               qyy_i*(cy2*yhatc4 - (2.0_dp*cy_i*uy_i(2) + yhat)*c3ri) + &
1121 >               qzz_i*(cz2*yhatc4 - (2.0_dp*cz_i*uz_i(2) + yhat)*c3ri) )
1122 >          dudz = dudz - preSw * ( &
1123 >               qxx_i*(cx2*zhatc4 - (2.0_dp*cx_i*ux_i(3) + zhat)*c3ri) + &
1124 >               qyy_i*(cy2*zhatc4 - (2.0_dp*cy_i*uy_i(3) + zhat)*c3ri) + &
1125 >               qzz_i*(cz2*zhatc4 - (2.0_dp*cz_i*uz_i(3) + zhat)*c3ri) )
1126            
1127 <          duduy_i(1) = duduy_i(1) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*xhat) &
1128 <               * (3.0d0*f1 + f3) )
1129 <          duduy_i(2) = duduy_i(2) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*yhat) &
1111 <               * (3.0d0*f1 + f3) )
1112 <          duduy_i(3) = duduy_i(3) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*zhat) &
1113 <               * (3.0d0*f1 + f3) )
1127 >          dudux_i(1) = dudux_i(1) + preSw*(qxx_i*cx_i*xhatdot2)
1128 >          dudux_i(2) = dudux_i(2) + preSw*(qxx_i*cx_i*yhatdot2)
1129 >          dudux_i(3) = dudux_i(3) + preSw*(qxx_i*cx_i*zhatdot2)
1130            
1131 <          duduz_i(1) = duduz_i(1) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*xhat) &
1132 <               * (3.0d0*f1 + f3) )
1133 <          duduz_i(2) = duduz_i(2) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*yhat) &
1134 <               * (3.0d0*f1 + f3) )
1135 <          duduz_i(3) = duduz_i(3) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*zhat) &
1136 <               * (3.0d0*f1 + f3) )
1137 <
1131 >          duduy_i(1) = duduy_i(1) + preSw*(qyy_i*cy_i*xhatdot2)
1132 >          duduy_i(2) = duduy_i(2) + preSw*(qyy_i*cy_i*yhatdot2)
1133 >          duduy_i(3) = duduy_i(3) + preSw*(qyy_i*cy_i*zhatdot2)
1134 >          
1135 >          duduz_i(1) = duduz_i(1) + preSw*(qzz_i*cz_i*xhatdot2)
1136 >          duduz_i(2) = duduz_i(2) + preSw*(qzz_i*cz_i*yhatdot2)
1137 >          duduz_i(3) = duduz_i(3) + preSw*(qzz_i*cz_i*zhatdot2)
1138         endif
1139      endif
1140  
1141  
1142      if (do_pot) then
1143   #ifdef IS_MPI
1144 <       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1145 <       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1144 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5_dp*epot
1145 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5_dp*epot
1146   #else
1147         pot = pot + epot
1148   #endif
# Line 1243 | Line 1259 | contains
1259      integer :: atid1
1260      real(kind=dp), dimension(9,nLocal) :: eFrame
1261      real(kind=dp), dimension(3,nLocal) :: t
1262 <    real(kind=dp) :: mu1, c1
1262 >    real(kind=dp) :: mu1, chg1
1263      real(kind=dp) :: preVal, epot, mypot
1264      real(kind=dp) :: eix, eiy, eiz
1265  
# Line 1259 | Line 1275 | contains
1275            mu1 = getDipoleMoment(atid1)
1276            
1277            preVal = pre22 * preRF2 * mu1*mu1
1278 <          mypot = mypot - 0.5d0*preVal
1278 >          mypot = mypot - 0.5_dp*preVal
1279            
1280            ! The self-correction term adds into the reaction field vector
1281            
# Line 1282 | Line 1298 | contains
1298      elseif ( (summationMethod .eq. SHIFTED_FORCE) .or. &
1299           (summationMethod .eq. SHIFTED_POTENTIAL) ) then
1300         if (ElectrostaticMap(atid1)%is_Charge) then
1301 <          c1 = getCharge(atid1)
1301 >          chg1 = getCharge(atid1)
1302            
1303            if (screeningMethod .eq. DAMPED) then
1304 <             mypot = mypot - (f0c * rcuti * 0.5_dp + &
1289 <                  dampingAlpha*invRootPi) * c1 * c1    
1304 >             mypot = mypot - 0.5_dp*(c1c + alphaPi) * chg1 * chg1    
1305              
1306            else            
1307 <             mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1307 >             mypot = mypot - 0.5_dp*(rcuti * chg1 * chg1)
1308              
1309            endif
1310         endif
# Line 1329 | Line 1344 | contains
1344         call checkSummationMethod()
1345      endif
1346  
1347 <    dudx = 0.0d0
1348 <    dudy = 0.0d0
1349 <    dudz = 0.0d0
1347 >    dudx = zero
1348 >    dudy = zero
1349 >    dudz = zero
1350  
1351 <    riji = 1.0d0/rij
1351 >    riji = 1.0_dp/rij
1352  
1353      xhat = d(1) * riji
1354      yhat = d(2) * riji
# Line 1357 | Line 1372 | contains
1372        
1373         myPot = myPot + sw*vterm
1374        
1375 <       dudr  = sw*preVal * 2.0d0*rfVal*riji
1375 >       dudr  = sw*preVal * 2.0_dp*rfVal*riji
1376        
1377         dudx = dudx + dudr * xhat
1378         dudy = dudy + dudr * yhat
# Line 1378 | Line 1393 | contains
1393         vterm = - pref * ct_j * ( ri2 - preRF2*rij )
1394         myPot = myPot + sw*vterm
1395        
1396 <       dudx = dudx - sw*pref*( ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1396 >       dudx = dudx - sw*pref*( ri3*(uz_j(1)-3.0_dp*ct_j*xhat) &
1397              - preRF2*uz_j(1) )
1398 <       dudy = dudy - sw*pref*( ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1398 >       dudy = dudy - sw*pref*( ri3*(uz_j(2)-3.0_dp*ct_j*yhat) &
1399              - preRF2*uz_j(2) )
1400 <       dudz = dudz - sw*pref*( ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1400 >       dudz = dudz - sw*pref*( ri3*(uz_j(3)-3.0_dp*ct_j*zhat) &
1401              - preRF2*uz_j(3) )
1402        
1403         duduz_j(1) = duduz_j(1) - sw * pref * xhat * ( ri2 - preRF2*rij )
# Line 1404 | Line 1419 | contains
1419         vterm = pref * ct_i * ( ri2 - preRF2*rij )
1420         myPot = myPot + sw*vterm
1421        
1422 <       dudx = dudx + sw*pref*( ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1422 >       dudx = dudx + sw*pref*( ri3*(uz_i(1)-3.0_dp*ct_i*xhat) &
1423              - preRF2*uz_i(1) )
1424 <       dudy = dudy + sw*pref*( ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1424 >       dudy = dudy + sw*pref*( ri3*(uz_i(2)-3.0_dp*ct_i*yhat) &
1425              - preRF2*uz_i(2) )
1426 <       dudz = dudz + sw*pref*( ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1426 >       dudz = dudz + sw*pref*( ri3*(uz_i(3)-3.0_dp*ct_i*zhat) &
1427              - preRF2*uz_i(3) )
1428        
1429         duduz_i(1) = duduz_i(1) + sw * pref * xhat * ( ri2 - preRF2*rij )
# Line 1440 | Line 1455 | end module electrostatic_module
1455      return
1456    end subroutine rf_self_excludes
1457  
1458 +  subroutine accumulate_box_dipole(atom1, eFrame, d, pChg, nChg, pChgPos, &
1459 +       nChgPos, dipVec, pChgCount, nChgCount)
1460 +    integer, intent(in) :: atom1
1461 +    logical :: i_is_Charge
1462 +    logical :: i_is_Dipole
1463 +    integer :: atid1
1464 +    integer :: pChgCount
1465 +    integer :: nChgCount
1466 +    real(kind=dp), intent(in), dimension(3) :: d
1467 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1468 +    real(kind=dp) :: pChg
1469 +    real(kind=dp) :: nChg
1470 +    real(kind=dp), dimension(3) :: pChgPos
1471 +    real(kind=dp), dimension(3) :: nChgPos
1472 +    real(kind=dp), dimension(3) :: dipVec
1473 +    real(kind=dp), dimension(3) :: uz_i
1474 +    real(kind=dp), dimension(3) :: pos
1475 +    real(kind=dp) :: q_i, mu_i
1476 +    real(kind=dp) :: pref, preVal
1477 +
1478 +    if (.not.summationMethodChecked) then
1479 +       call checkSummationMethod()
1480 +    endif
1481 +
1482 +    ! this is a local only array, so we use the local atom type id's:
1483 +    atid1 = atid(atom1)
1484 +    i_is_Charge = ElectrostaticMap(atid1)%is_Charge
1485 +    i_is_Dipole = ElectrostaticMap(atid1)%is_Dipole
1486 +    
1487 +    if (i_is_Charge) then
1488 +       q_i = ElectrostaticMap(atid1)%charge
1489 +       ! convert to the proper units
1490 +       q_i = q_i * chargeToC
1491 +       pos = d * angstromToM
1492 +
1493 +       if (q_i.le.0.0_dp) then
1494 +          nChg = nChg - q_i
1495 +          nChgPos(1) = nChgPos(1) + pos(1)
1496 +          nChgPos(2) = nChgPos(2) + pos(2)
1497 +          nChgPos(3) = nChgPos(3) + pos(3)
1498 +          nChgCount = nChgCount + 1
1499 +
1500 +       else
1501 +          pChg = pChg + q_i
1502 +          pChgPos(1) = pChgPos(1) + pos(1)
1503 +          pChgPos(2) = pChgPos(2) + pos(2)
1504 +          pChgPos(3) = pChgPos(3) + pos(3)
1505 +          pChgCount = pChgCount + 1
1506 +
1507 +       endif
1508 +
1509 +    endif
1510 +    
1511 +    if (i_is_Dipole) then
1512 +       mu_i = ElectrostaticMap(atid1)%dipole_moment
1513 +       uz_i(1) = eFrame(3,atom1)
1514 +       uz_i(2) = eFrame(6,atom1)
1515 +       uz_i(3) = eFrame(9,atom1)
1516 +       ! convert to the proper units
1517 +       mu_i = mu_i * debyeToCm
1518 +
1519 +       dipVec(1) = dipVec(1) + uz_i(1)*mu_i
1520 +       dipVec(2) = dipVec(2) + uz_i(2)*mu_i
1521 +       dipVec(3) = dipVec(3) + uz_i(3)*mu_i
1522 +
1523 +    endif
1524 +  
1525 +    return
1526 +  end subroutine accumulate_box_dipole
1527 +
1528   end module electrostatic_module

Diff Legend

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