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 2162 by chrisfen, Mon Apr 11 20:19:22 2005 UTC vs.
Revision 2405 by chrisfen, Tue Nov 1 19:24:57 2005 UTC

# Line 40 | Line 40 | module electrostatic_module
40   !!
41  
42   module electrostatic_module
43 <  
43 >
44    use force_globals
45    use definitions
46    use atype_module
# 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 68 | Line 74 | module electrostatic_module
74    !! This unit is also known affectionately as an esu centi-barn.
75    real(kind=dp), parameter :: pre14 = 69.13373_dp
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 = 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 +  real(kind=dp), save :: rcuti = 0.0_DP
90 +  real(kind=dp), save :: rcuti2 = 0.0_DP
91 +  real(kind=dp), save :: rcuti3 = 0.0_DP
92 +  real(kind=dp), save :: rcuti4 = 0.0_DP
93 +  real(kind=dp), save :: alphaPi = 0.0_DP
94 +  real(kind=dp), save :: invRootPi = 0.0_DP
95 +  real(kind=dp), save :: rrf = 1.0_DP
96 +  real(kind=dp), save :: rt = 1.0_DP
97 +  real(kind=dp), save :: rrfsq = 1.0_DP
98 +  real(kind=dp), save :: preRF = 0.0_DP
99 +  real(kind=dp), save :: preRF2 = 0.0_DP
100 +
101 + #ifdef __IFC
102 + ! error function for ifc version > 7.
103 +  double precision, external :: derfc
104 + #endif
105 +  
106 +  public :: setElectrostaticSummationMethod
107 +  public :: setElectrostaticCutoffRadius
108 +  public :: setDampedWolfAlpha
109 +  public :: setReactionFieldDielectric
110    public :: newElectrostaticType
111    public :: setCharge
112    public :: setDipoleMoment
# Line 76 | Line 115 | module electrostatic_module
115    public :: doElectrostaticPair
116    public :: getCharge
117    public :: getDipoleMoment
118 <  public :: pre22
118 >  public :: destroyElectrostaticTypes
119 >  public :: self_self
120 >  public :: rf_self_excludes
121  
122    type :: Electrostatic
123       integer :: c_ident
# Line 84 | Line 125 | module electrostatic_module
125       logical :: is_Dipole = .false.
126       logical :: is_SplitDipole = .false.
127       logical :: is_Quadrupole = .false.
128 +     logical :: is_Tap = .false.
129       real(kind=DP) :: charge = 0.0_DP
130       real(kind=DP) :: dipole_moment = 0.0_DP
131       real(kind=DP) :: split_dipole_distance = 0.0_DP
# Line 93 | Line 135 | contains
135    type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap
136  
137   contains
138 +
139 +  subroutine setElectrostaticSummationMethod(the_ESM)
140 +    integer, intent(in) :: the_ESM    
141 +
142 +    if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
143 +       call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
144 +    endif
145 +
146 +    summationMethod = the_ESM
147 +
148 +  end subroutine setElectrostaticSummationMethod
149 +
150 +  subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
151 +    real(kind=dp), intent(in) :: thisRcut
152 +    real(kind=dp), intent(in) :: thisRsw
153 +    defaultCutoff = thisRcut
154 +    rrf = defaultCutoff
155 +    rt = thisRsw
156 +    haveDefaultCutoff = .true.
157 +  end subroutine setElectrostaticCutoffRadius
158 +
159 +  subroutine setDampedWolfAlpha(thisAlpha)
160 +    real(kind=dp), intent(in) :: thisAlpha
161 +    dampingAlpha = thisAlpha
162 +    haveDampingAlpha = .true.
163 +  end subroutine setDampedWolfAlpha
164 +  
165 +  subroutine setReactionFieldDielectric(thisDielectric)
166 +    real(kind=dp), intent(in) :: thisDielectric
167 +    dielectric = thisDielectric
168 +    haveDielectric = .true.
169 +  end subroutine setReactionFieldDielectric
170  
171    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
172 <       is_SplitDipole, is_Quadrupole, status)
173 <    
172 >       is_SplitDipole, is_Quadrupole, is_Tap, status)
173 >
174      integer, intent(in) :: c_ident
175      logical, intent(in) :: is_Charge
176      logical, intent(in) :: is_Dipole
177      logical, intent(in) :: is_SplitDipole
178      logical, intent(in) :: is_Quadrupole
179 +    logical, intent(in) :: is_Tap
180      integer, intent(out) :: status
181      integer :: nAtypes, myATID, i, j
182  
183      status = 0
184      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
185 <    
185 >
186      !! Be simple-minded and assume that we need an ElectrostaticMap that
187      !! is the same size as the total number of atom types
188  
189      if (.not.allocated(ElectrostaticMap)) then
190 <      
190 >
191         nAtypes = getSize(atypes)
192 <    
192 >
193         if (nAtypes == 0) then
194            status = -1
195            return
196         end if
197 <      
197 >
198         if (.not. allocated(ElectrostaticMap)) then
199            allocate(ElectrostaticMap(nAtypes))
200         endif
201 <      
201 >
202      end if
203  
204      if (myATID .gt. size(ElectrostaticMap)) then
205         status = -1
206         return
207      endif
208 <    
208 >
209      ! set the values for ElectrostaticMap for this atom type:
210  
211      ElectrostaticMap(myATID)%c_ident = c_ident
# Line 138 | Line 213 | contains
213      ElectrostaticMap(myATID)%is_Dipole = is_Dipole
214      ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole
215      ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
216 <    
216 >    ElectrostaticMap(myATID)%is_Tap = is_Tap
217 >
218    end subroutine newElectrostaticType
219  
220    subroutine setCharge(c_ident, charge, status)
# Line 166 | Line 242 | contains
242         call handleError("electrostatic", "Attempt to setCharge of an atom type that is not a charge!")
243         status = -1
244         return
245 <    endif      
245 >    endif
246  
247      ElectrostaticMap(myATID)%charge = charge
248    end subroutine setCharge
# Line 257 | Line 333 | contains
333         status = -1
334         return
335      endif
336 <    
336 >
337      do i = 1, 3
338 <          ElectrostaticMap(myATID)%quadrupole_moments(i) = &
339 <               quadrupole_moments(i)
340 <       enddo
338 >       ElectrostaticMap(myATID)%quadrupole_moments(i) = &
339 >            quadrupole_moments(i)
340 >    enddo
341  
342    end subroutine setQuadrupoleMoments
343  
344 <  
344 >
345    function getCharge(atid) result (c)
346      integer, intent(in) :: atid
347      integer :: localError
348      real(kind=dp) :: c
349 <    
349 >
350      if (.not.allocated(ElectrostaticMap)) then
351         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
352         return
353      end if
354 <    
354 >
355      if (.not.ElectrostaticMap(atid)%is_Charge) then
356         call handleError("electrostatic", "getCharge was called for an atom type that isn't a charge!")
357         return
358      endif
359 <    
359 >
360      c = ElectrostaticMap(atid)%charge
361    end function getCharge
362  
# Line 288 | Line 364 | contains
364      integer, intent(in) :: atid
365      integer :: localError
366      real(kind=dp) :: dm
367 <    
367 >
368      if (.not.allocated(ElectrostaticMap)) then
369         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
370         return
371      end if
372 <    
372 >
373      if (.not.ElectrostaticMap(atid)%is_Dipole) then
374         call handleError("electrostatic", "getDipoleMoment was called for an atom type that isn't a dipole!")
375         return
376      endif
377 <    
377 >
378      dm = ElectrostaticMap(atid)%dipole_moment
379    end function getDipoleMoment
380 +
381 +  subroutine checkSummationMethod()
382 +
383 +    if (.not.haveDefaultCutoff) then
384 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
385 +    endif
386 +
387 +    rcuti = 1.0d0 / defaultCutoff
388 +    rcuti2 = rcuti*rcuti
389 +    rcuti3 = rcuti2*rcuti
390 +    rcuti4 = rcuti2*rcuti2
391 +
392 +    if (summationMethod .eq. DAMPED_WOLF) then
393 +       if (.not.haveDampingAlpha) then
394 +          call handleError("checkSummationMethod", "no Damping Alpha set!")
395 +       endif
396 +      
397 +       if (.not.haveDefaultCutoff) then
398 +          call handleError("checkSummationMethod", "no Default Cutoff set!")
399 +       endif
400 +
401 +       constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
402 +       constERFC = derfc(dampingAlpha*defaultCutoff)
403 +       invRootPi = 0.56418958354775628695d0
404 +       alphaPi = 2*dampingAlpha*invRootPi
405 +      
406 +    endif
407  
408 +    if (summationMethod .eq. REACTION_FIELD) then
409 +       if (haveDielectric) then
410 +          defaultCutoff2 = defaultCutoff*defaultCutoff
411 +          preRF = (dielectric-1.0d0) / &
412 +               ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
413 +          preRF2 = 2.0d0*preRF
414 +       else
415 +          call handleError("checkSummationMethod", "Dielectric not set")
416 +       endif
417 +      
418 +    endif
419 +
420 +    summationMethodChecked = .true.
421 +  end subroutine checkSummationMethod
422 +
423 +
424    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
425         vpair, fpair, pot, eFrame, f, t, do_pot)
426 <    
426 >
427      logical, intent(in) :: do_pot
428 <    
428 >
429      integer, intent(in) :: atom1, atom2
430      integer :: localError
431  
432      real(kind=dp), intent(in) :: rij, r2, sw
433      real(kind=dp), intent(in), dimension(3) :: d
434      real(kind=dp), intent(inout) :: vpair
435 <    real(kind=dp), intent(inout), dimension(3) :: fpair
435 >    real(kind=dp), intent(inout), dimension(3) :: fpair    
436  
437      real( kind = dp ) :: pot
438      real( kind = dp ), dimension(9,nLocal) :: eFrame
439      real( kind = dp ), dimension(3,nLocal) :: f
440      real( kind = dp ), dimension(3,nLocal) :: t
441 <    
441 >
442      real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
443      real (kind = dp), dimension(3) :: ux_j, uy_j, uz_j
444      real (kind = dp), dimension(3) :: dudux_i, duduy_i, duduz_i
# Line 327 | Line 446 | contains
446  
447      logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole
448      logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole
449 +    logical :: i_is_Tap, j_is_Tap
450      integer :: me1, me2, id1, id2
451      real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
452      real (kind=dp) :: qxx_i, qyy_i, qzz_i
# Line 336 | Line 456 | contains
456      real (kind=dp) :: cx2, cy2, cz2
457      real (kind=dp) :: ct_i, ct_j, ct_ij, a1
458      real (kind=dp) :: riji, ri, ri2, ri3, ri4
459 <    real (kind=dp) :: pref, vterm, epot, dudr    
459 >    real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
460      real (kind=dp) :: xhat, yhat, zhat
461      real (kind=dp) :: dudx, dudy, dudz
462      real (kind=dp) :: scale, sc2, bigR
463 +    real (kind=dp) :: varERFC, varEXP
464 +    real (kind=dp) :: limScale
465 +    real (kind=dp) :: preVal, rfVal
466  
467      if (.not.allocated(ElectrostaticMap)) then
468         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
469         return
470      end if
471  
472 +    if (.not.summationMethodChecked) then
473 +       call checkSummationMethod()
474 +    endif
475 +
476   #ifdef IS_MPI
477      me1 = atid_Row(atom1)
478      me2 = atid_Col(atom2)
# Line 357 | Line 484 | contains
484      !! some variables we'll need independent of electrostatic type:
485  
486      riji = 1.0d0 / rij
487 <
487 >  
488      xhat = d(1) * riji
489      yhat = d(2) * riji
490      zhat = d(3) * riji
491  
492      !! logicals
366
493      i_is_Charge = ElectrostaticMap(me1)%is_Charge
494      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
495      i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole
496      i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole
497 +    i_is_Tap = ElectrostaticMap(me1)%is_Tap
498  
499      j_is_Charge = ElectrostaticMap(me2)%is_Charge
500      j_is_Dipole = ElectrostaticMap(me2)%is_Dipole
501      j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole
502      j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole
503 +    j_is_Tap = ElectrostaticMap(me2)%is_Tap
504  
505      if (i_is_Charge) then
506         q_i = ElectrostaticMap(me1)%charge      
507      endif
508 <    
508 >
509      if (i_is_Dipole) then
510         mu_i = ElectrostaticMap(me1)%dipole_moment
511   #ifdef IS_MPI
# Line 394 | Line 522 | contains
522         if (i_is_SplitDipole) then
523            d_i = ElectrostaticMap(me1)%split_dipole_distance
524         endif
525 <      
525 >
526      endif
527  
528      if (i_is_Quadrupole) then
# Line 427 | Line 555 | contains
555         cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
556      endif
557  
430
558      if (j_is_Charge) then
559         q_j = ElectrostaticMap(me2)%charge      
560      endif
561 <    
561 >
562      if (j_is_Dipole) then
563         mu_j = ElectrostaticMap(me2)%dipole_moment
564   #ifdef IS_MPI
# Line 479 | Line 606 | contains
606         cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
607         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
608      endif
609 <
609 >  
610      epot = 0.0_dp
611      dudx = 0.0_dp
612      dudy = 0.0_dp
# Line 496 | Line 623 | contains
623      if (i_is_Charge) then
624  
625         if (j_is_Charge) then
499          
500          vterm = pre11 * q_i * q_j * riji
501          vpair = vpair + vterm
502          epot = epot + sw*vterm
626  
627 <          dudr  = - sw * vterm * riji
627 >          if (summationMethod .eq. UNDAMPED_WOLF) then
628 >             vterm = pre11 * q_i * q_j * (riji - rcuti)
629 >             vpair = vpair + vterm
630 >             epot = epot + sw*vterm
631 >            
632 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)
633 >            
634 >             dudx = dudx + dudr * xhat
635 >             dudy = dudy + dudr * yhat
636 >             dudz = dudz + dudr * zhat
637  
638 <          dudx = dudx + dudr * xhat
639 <          dudy = dudy + dudr * yhat
640 <          dudz = dudz + dudr * zhat
641 <      
642 <       endif
638 >          elseif (summationMethod .eq. DAMPED_WOLF) then
639 >             varERFC = derfc(dampingAlpha*rij)
640 >             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
641 >             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
642 >             vpair = vpair + vterm
643 >             epot = epot + sw*vterm
644 >            
645 >             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
646 >                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
647 >                  + alphaPi*constEXP*rcuti)) )
648 >            
649 >             dudx = dudx + dudr * xhat
650 >             dudy = dudy + dudr * yhat
651 >             dudz = dudz + dudr * zhat
652  
653 <       if (j_is_Dipole) then
653 >          elseif (summationMethod .eq. REACTION_FIELD) then
654 >             preVal = pre11 * q_i * q_j
655 >             rfVal = preRF*rij*rij
656 >             vterm = preVal * ( riji + rfVal )
657 >            
658 >             vpair = vpair + vterm
659 >             epot = epot + sw*vterm
660 >            
661 >             dudr  = sw * preVal * ( 2.0d0*rfVal - riji )*riji
662 >            
663 >             dudx = dudx + dudr * xhat
664 >             dudy = dudy + dudr * yhat
665 >             dudz = dudz + dudr * zhat
666  
514          if (j_is_SplitDipole) then
515             BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
516             ri = 1.0_dp / BigR
517             scale = rij * ri
667            else
668 <             ri = riji
669 <             scale = 1.0_dp
668 >             vterm = pre11 * q_i * q_j * riji
669 >             vpair = vpair + vterm
670 >             epot = epot + sw*vterm
671 >            
672 >             dudr  = - sw * vterm * riji
673 >            
674 >             dudx = dudx + dudr * xhat
675 >             dudy = dudy + dudr * yhat
676 >             dudz = dudz + dudr * zhat
677 >
678            endif
679  
680 <          ri2 = ri * ri
681 <          ri3 = ri2 * ri
682 <          sc2 = scale * scale
683 <            
680 >       endif
681 >
682 >       if (j_is_Dipole) then
683 >
684            pref = pre12 * q_i * mu_j
528          vterm = - pref * ct_j * ri2 * scale
529          vpair = vpair + vterm
530          epot = epot + sw * vterm
685  
686 <          !! this has a + sign in the () because the rij vector is
687 <          !! r_j - r_i and the charge-dipole potential takes the origin
688 <          !! as the point dipole, which is atom j in this case.
686 >          if (summationMethod .eq. UNDAMPED_WOLF) then
687 >             ri2 = riji * riji
688 >             ri3 = ri2 * riji
689  
690 <          dudx = dudx - pref * sw * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
691 <          dudy = dudy - pref * sw * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
692 <          dudz = dudz - pref * sw * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
690 >             pref = pre12 * q_i * mu_j
691 >             vterm = - pref * ct_j * (ri2 - rcuti2)
692 >             vpair = vpair + vterm
693 >             epot = epot + sw*vterm
694 >            
695 >             !! this has a + sign in the () because the rij vector is
696 >             !! r_j - r_i and the charge-dipole potential takes the origin
697 >             !! as the point dipole, which is atom j in this case.
698 >            
699 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
700 >                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
701 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
702 >                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
703 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
704 >                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
705 >            
706 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
707 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
708 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
709  
710 <          duduz_j(1) = duduz_j(1) - pref * sw * ri2 * xhat * scale
711 <          duduz_j(2) = duduz_j(2) - pref * sw * ri2 * yhat * scale
712 <          duduz_j(3) = duduz_j(3) - pref * sw * ri2 * zhat * scale
713 <          
710 >          elseif (summationMethod .eq. REACTION_FIELD) then
711 >             ri2 = riji * riji
712 >             ri3 = ri2 * riji
713 >    
714 >             pref = pre12 * q_i * mu_j
715 >             vterm = - pref * ct_j * ( ri2 - preRF2*rij )
716 >             vpair = vpair + vterm
717 >             epot = epot + sw*vterm
718 >            
719 >             !! this has a + sign in the () because the rij vector is
720 >             !! r_j - r_i and the charge-dipole potential takes the origin
721 >             !! as the point dipole, which is atom j in this case.
722 >            
723 >             dudx = dudx - sw*pref*( ri3*(uz_j(1) - 3.0d0*ct_j*xhat) - &
724 >                                     preRF2*uz_j(1) )
725 >             dudy = dudy - sw*pref*( ri3*(uz_j(2) - 3.0d0*ct_j*yhat) - &
726 >                                     preRF2*uz_j(2) )
727 >             dudz = dudz - sw*pref*( ri3*(uz_j(3) - 3.0d0*ct_j*zhat) - &
728 >                                     preRF2*uz_j(3) )        
729 >             duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
730 >             duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
731 >             duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
732 >
733 >          else
734 >             if (j_is_SplitDipole) then
735 >                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
736 >                ri = 1.0_dp / BigR
737 >                scale = rij * ri
738 >             else
739 >                ri = riji
740 >                scale = 1.0_dp
741 >             endif
742 >            
743 >             ri2 = ri * ri
744 >             ri3 = ri2 * ri
745 >             sc2 = scale * scale
746 >
747 >             pref = pre12 * q_i * mu_j
748 >             vterm = - pref * ct_j * ri2 * scale
749 >             vpair = vpair + vterm
750 >             epot = epot + sw*vterm
751 >            
752 >             !! this has a + sign in the () because the rij vector is
753 >             !! r_j - r_i and the charge-dipole potential takes the origin
754 >             !! as the point dipole, which is atom j in this case.
755 >            
756 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
757 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
758 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
759 >            
760 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
761 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
762 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
763 >
764 >          endif
765         endif
766  
767         if (j_is_Quadrupole) then
# Line 551 | Line 772 | contains
772            cy2 = cy_j * cy_j
773            cz2 = cz_j * cz_j
774  
775 <
776 <          pref =  pre14 * q_i / 3.0_dp
777 <          vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
778 <               qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
779 <               qzz_j * (3.0_dp*cz2 - 1.0_dp))
780 <          vpair = vpair + vterm
781 <          epot = epot + sw * vterm
782 <
783 <          dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( &
784 <               qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
785 <               qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
786 <               qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
787 <          dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + pref * sw * ri4 * ( &
788 <               qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
789 <               qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
790 <               qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
791 <          dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + pref * sw * ri4 * ( &
792 <               qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
793 <               qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
794 <               qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
775 >          if (summationMethod .eq. UNDAMPED_WOLF) then
776 >             pref =  pre14 * q_i / 3.0_dp
777 >             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
778 >                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
779 >                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
780 >             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
781 >                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
782 >                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
783 >             vpair = vpair + ( vterm1 - vterm2 )
784 >             epot = epot + sw*( vterm1 - vterm2 )
785 >            
786 >             dudx = dudx - (5.0_dp * &
787 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
788 >                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
789 >                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
790 >                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
791 >                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
792 >                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
793 >                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
794 >             dudy = dudy - (5.0_dp * &
795 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
796 >                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
797 >                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
798 >                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
799 >                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
800 >                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
801 >                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
802 >             dudz = dudz - (5.0_dp * &
803 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
804 >                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
805 >                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
806 >                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
807 >                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
808 >                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
809 >                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
810 >            
811 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
812 >                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
813 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
814 >                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
815 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
816 >                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
817 >            
818 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
819 >                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
820 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
821 >                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
822 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
823 >                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
824 >            
825 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
826 >                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
827 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
828 >                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
829 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
830 >                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
831 >        
832 >          else
833 >             pref =  pre14 * q_i / 3.0_dp
834 >             vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
835 >                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
836 >                  qzz_j * (3.0_dp*cz2 - 1.0_dp))
837 >             vpair = vpair + vterm
838 >             epot = epot + sw*vterm
839 >            
840 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
841 >                  qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
842 >                  qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
843 >                  qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
844 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
845 >                  qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
846 >                  qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
847 >                  qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
848 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
849 >                  qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
850 >                  qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
851 >                  qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
852 >            
853 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
854 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
855 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
856 >            
857 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
858 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
859 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
860 >            
861 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
862 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
863 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
864            
865 <          dudux_j(1) = dudux_j(1) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*xhat)
576 <          dudux_j(2) = dudux_j(2) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*yhat)
577 <          dudux_j(3) = dudux_j(3) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*zhat)
578 <
579 <          duduy_j(1) = duduy_j(1) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*xhat)
580 <          duduy_j(2) = duduy_j(2) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*yhat)
581 <          duduy_j(3) = duduy_j(3) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*zhat)
582 <
583 <          duduz_j(1) = duduz_j(1) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*xhat)
584 <          duduz_j(2) = duduz_j(2) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*yhat)
585 <          duduz_j(3) = duduz_j(3) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*zhat)
865 >          endif
866         endif
587
867      endif
868 <  
868 >
869      if (i_is_Dipole) then
870 <      
870 >
871         if (j_is_Charge) then
872 +          
873 +          pref = pre12 * q_j * mu_i
874 +          
875 +          if (summationMethod .eq. UNDAMPED_WOLF) then
876 +             ri2 = riji * riji
877 +             ri3 = ri2 * riji
878  
879 <          if (i_is_SplitDipole) then
880 <             BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
881 <             ri = 1.0_dp / BigR
882 <             scale = rij * ri
598 <          else
599 <             ri = riji
600 <             scale = 1.0_dp
601 <          endif
602 <
603 <          ri2 = ri * ri
604 <          ri3 = ri2 * ri
605 <          sc2 = scale * scale
879 >             pref = pre12 * q_j * mu_i
880 >             vterm = pref * ct_i * (ri2 - rcuti2)
881 >             vpair = vpair + vterm
882 >             epot = epot + sw*vterm
883              
884 <          pref = pre12 * q_j * mu_i
885 <          vterm = pref * ct_i * ri2 * scale
886 <          vpair = vpair + vterm
887 <          epot = epot + sw * vterm
884 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
885 >                  - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
886 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
887 >                  - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
888 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
889 >                  - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
890 >            
891 >             duduz_i(1) = duduz_i(1) + sw*pref*( ri2*xhat - d(1)*rcuti3 )
892 >             duduz_i(2) = duduz_i(2) + sw*pref*( ri2*yhat - d(2)*rcuti3 )
893 >             duduz_i(3) = duduz_i(3) + sw*pref*( ri2*zhat - d(3)*rcuti3 )
894  
895 <          dudx = dudx + pref * sw * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
896 <          dudy = dudy + pref * sw * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
897 <          dudz = dudz + pref * sw * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
895 >          elseif (summationMethod .eq. REACTION_FIELD) then
896 >             ri2 = riji * riji
897 >             ri3 = ri2 * riji
898  
899 <          duduz_i(1) = duduz_i(1) + pref * sw * ri2 * xhat * scale
900 <          duduz_i(2) = duduz_i(2) + pref * sw * ri2 * yhat * scale
901 <          duduz_i(3) = duduz_i(3) + pref * sw * ri2 * zhat * scale
902 <       endif
899 >             pref = pre12 * q_j * mu_i
900 >             vterm = pref * ct_i * ( ri2 - preRF2*rij )
901 >             vpair = vpair + vterm
902 >             epot = epot + sw*vterm
903 >            
904 >             dudx = dudx + sw*pref * ( ri3*(uz_i(1) - 3.0d0*ct_i*xhat) - &
905 >                  preRF2*uz_i(1) )
906 >             dudy = dudy + sw*pref * ( ri3*(uz_i(2) - 3.0d0*ct_i*yhat) - &
907 >                  preRF2*uz_i(2) )
908 >             dudz = dudz + sw*pref * ( ri3*(uz_i(3) - 3.0d0*ct_i*zhat) - &
909 >                  preRF2*uz_i(3) )
910 >            
911 >             duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF2*rij )
912 >             duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF2*rij )
913 >             duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF2*rij )
914  
621       if (j_is_Dipole) then
622
623          if (i_is_SplitDipole) then
624             if (j_is_SplitDipole) then
625                BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
626             else
627                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
628             endif
629             ri = 1.0_dp / BigR
630             scale = rij * ri                
915            else
916 <             if (j_is_SplitDipole) then
917 <                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
916 >             if (i_is_SplitDipole) then
917 >                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
918                  ri = 1.0_dp / BigR
919 <                scale = rij * ri                            
920 <             else                
919 >                scale = rij * ri
920 >             else
921                  ri = riji
922                  scale = 1.0_dp
923               endif
924 +            
925 +             ri2 = ri * ri
926 +             ri3 = ri2 * ri
927 +             sc2 = scale * scale
928 +
929 +             pref = pre12 * q_j * mu_i
930 +             vterm = pref * ct_i * ri2 * scale
931 +             vpair = vpair + vterm
932 +             epot = epot + sw*vterm
933 +            
934 +             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
935 +             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
936 +             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
937 +            
938 +             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
939 +             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
940 +             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
941            endif
942 +       endif
943 +      
944 +       if (j_is_Dipole) then
945  
946 <          ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
947 <
948 <          ri2 = ri * ri
949 <          ri3 = ri2 * ri
950 <          ri4 = ri2 * ri2
951 <          sc2 = scale * scale
952 <
953 <          pref = pre22 * mu_i * mu_j
954 <          vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
955 <          vpair = vpair + vterm
956 <          epot = epot + sw * vterm
946 >          if (summationMethod .eq. UNDAMPED_WOLF) then
947 > !!$             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 > !!$
953 > !!$             pref = pre22 * mu_i * mu_j
954 > !!$             vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
955 > !!$             vpair = vpair + vterm
956 > !!$             epot = epot + sw*vterm
957 > !!$            
958 > !!$             a1 = 5.0d0 * ct_i * ct_j - ct_ij
959 > !!$            
960 > !!$             dudx = dudx + sw*pref*3.0d0*( &
961 > !!$                  ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
962 > !!$                  - rcuti4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) )
963 > !!$             dudy = dudy + sw*pref*3.0d0*( &
964 > !!$                  ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
965 > !!$                  - rcuti4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) )
966 > !!$             dudz = dudz + sw*pref*3.0d0*( &
967 > !!$                  ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
968 > !!$                  - rcuti4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) )
969 > !!$            
970 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
971 > !!$                  - rcuti3*(uz_j(1) - 3.0d0*ct_j*xhat))
972 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
973 > !!$                  - rcuti3*(uz_j(2) - 3.0d0*ct_j*yhat))
974 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
975 > !!$                  - rcuti3*(uz_j(3) - 3.0d0*ct_j*zhat))
976 > !!$             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
977 > !!$                  - rcuti3*(uz_i(1) - 3.0d0*ct_i*xhat))
978 > !!$             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
979 > !!$                  - rcuti3*(uz_i(2) - 3.0d0*ct_i*yhat))
980 > !!$             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
981 > !!$                  - rcuti3*(uz_i(3) - 3.0d0*ct_i*zhat))
982            
983 <          a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
983 >          elseif (summationMethod .eq. DAMPED_WOLF) then
984 >             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
985 >            
986 >             ri2 = riji * riji
987 >             ri3 = ri2 * riji
988 >             ri4 = ri2 * ri2
989 >             sc2 = scale * scale
990 >            
991 >             pref = pre22 * mu_i * mu_j
992 >             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j)
993 >             vpair = vpair + vterm
994 >             epot = epot + sw*vterm
995 >            
996 >             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
997 >            
998 >             dudx = dudx + sw*pref*3.0d0*ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
999 >             dudy = dudy + sw*pref*3.0d0*ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1000 >             dudz = dudz + sw*pref*3.0d0*ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1001 >            
1002 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 *(uz_j(1) - 3.0d0*ct_j*xhat)
1003 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 *(uz_j(2) - 3.0d0*ct_j*yhat)
1004 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 *(uz_j(3) - 3.0d0*ct_j*zhat)
1005 >            
1006 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 *(uz_i(1) - 3.0d0*ct_i*xhat)
1007 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 *(uz_i(2) - 3.0d0*ct_i*yhat)
1008 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 *(uz_i(3) - 3.0d0*ct_i*zhat)
1009 >            
1010 >          elseif (summationMethod .eq. REACTION_FIELD) then
1011 >             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1012  
1013 <          dudx=dudx+pref*sw*3.0d0*ri4*scale*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1014 <          dudy=dudy+pref*sw*3.0d0*ri4*scale*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1015 <          dudz=dudz+pref*sw*3.0d0*ri4*scale*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1013 >             ri2 = riji * riji
1014 >             ri3 = ri2 * riji
1015 >             ri4 = ri2 * ri2
1016  
1017 <          duduz_i(1) = duduz_i(1) + pref*sw*ri3*(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1018 <          duduz_i(2) = duduz_i(2) + pref*sw*ri3*(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1019 <          duduz_i(3) = duduz_i(3) + pref*sw*ri3*(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1017 >             pref = pre22 * mu_i * mu_j
1018 >              
1019 >             vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
1020 >                  preRF2*ct_ij )
1021 >             vpair = vpair + vterm
1022 >             epot = epot + sw*vterm
1023 >            
1024 >             a1 = 5.0d0 * ct_i * ct_j - ct_ij
1025 >            
1026 >             dudx = dudx + sw*pref*3.0d0*ri4 &
1027 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1028 >             dudy = dudy + sw*pref*3.0d0*ri4 &
1029 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1030 >             dudz = dudz + sw*pref*3.0d0*ri4 &
1031 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1032 >            
1033 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1034 >                  - preRF2*uz_j(1))
1035 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1036 >                  - preRF2*uz_j(2))
1037 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1038 >                  - preRF2*uz_j(3))
1039 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1040 >                  - preRF2*uz_i(1))
1041 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1042 >                  - preRF2*uz_i(2))
1043 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1044 >                  - preRF2*uz_i(3))
1045  
1046 <          duduz_j(1) = duduz_j(1) + pref*sw*ri3*(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1047 <          duduz_j(2) = duduz_j(2) + pref*sw*ri3*(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1048 <          duduz_j(3) = duduz_j(3) + pref*sw*ri3*(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1046 >          else
1047 >             if (i_is_SplitDipole) then
1048 >                if (j_is_SplitDipole) then
1049 >                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
1050 >                else
1051 >                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
1052 >                endif
1053 >                ri = 1.0_dp / BigR
1054 >                scale = rij * ri                
1055 >             else
1056 >                if (j_is_SplitDipole) then
1057 >                   BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
1058 >                   ri = 1.0_dp / BigR
1059 >                   scale = rij * ri                            
1060 >                else                
1061 >                   ri = riji
1062 >                   scale = 1.0_dp
1063 >                endif
1064 >             endif
1065 >            
1066 >             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1067 >            
1068 >             ri2 = ri * ri
1069 >             ri3 = ri2 * ri
1070 >             ri4 = ri2 * ri2
1071 >             sc2 = scale * scale
1072 >            
1073 >             pref = pre22 * mu_i * mu_j
1074 >             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1075 >             vpair = vpair + vterm
1076 >             epot = epot + sw*vterm
1077 >            
1078 >             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1079 >            
1080 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
1081 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1082 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
1083 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1084 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
1085 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1086 >            
1087 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
1088 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1089 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
1090 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1091 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
1092 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1093 >            
1094 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1095 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1096 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1097 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1098 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1099 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1100 >          endif
1101         endif
668
1102      endif
1103  
1104      if (i_is_Quadrupole) then
1105         if (j_is_Charge) then
1106 <          
1106 >
1107            ri2 = riji * riji
1108            ri3 = ri2 * riji
1109            ri4 = ri2 * ri2
1110            cx2 = cx_i * cx_i
1111            cy2 = cy_i * cy_i
1112            cz2 = cz_i * cz_i
1113 <          
1114 <          pref = pre14 * q_j / 3.0_dp
1115 <          vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1116 <               qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1117 <               qzz_i * (3.0_dp*cz2 - 1.0_dp))
1118 <          vpair = vpair + vterm
1119 <          epot = epot + sw * vterm
1120 <          
1121 <          dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( &
1122 <               qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1123 <               qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1124 <               qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1125 <          dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + pref * sw * ri4 * ( &
1126 <               qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1127 <               qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1128 <               qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1129 <          dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + pref * sw * ri4 * ( &
1130 <               qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1131 <               qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1132 <               qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1133 <          
1134 <          dudux_i(1) = dudux_i(1) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*xhat)
1135 <          dudux_i(2) = dudux_i(2) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*yhat)
1136 <          dudux_i(3) = dudux_i(3) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*zhat)
1137 <          
1138 <          duduy_i(1) = duduy_i(1) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*xhat)
1139 <          duduy_i(2) = duduy_i(2) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*yhat)
1140 <          duduy_i(3) = duduy_i(3) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*zhat)
1141 <          
1142 <          duduz_i(1) = duduz_i(1) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*xhat)
1143 <          duduz_i(2) = duduz_i(2) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*yhat)
1144 <          duduz_i(3) = duduz_i(3) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*zhat)
1113 >
1114 >          if (summationMethod .eq. UNDAMPED_WOLF) then
1115 >             pref = pre14 * q_j / 3.0_dp
1116 >             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1117 >                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1118 >                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1119 >             vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1120 >                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1121 >                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1122 >             vpair = vpair + ( vterm1 - vterm2 )
1123 >             epot = epot + sw*( vterm1 - vterm2 )
1124 >            
1125 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1126 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1127 >                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1128 >                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1129 >                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1130 >                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1131 >                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1132 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1133 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1134 >                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1135 >                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1136 >                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1137 >                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1138 >                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1139 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1140 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1141 >                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1142 >                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1143 >                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1144 >                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1145 >                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1146 >            
1147 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1148 >                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1149 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1150 >                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1151 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1152 >                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1153 >            
1154 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1155 >                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1156 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1157 >                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1158 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1159 >                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1160 >            
1161 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1162 >                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1163 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1164 >                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1165 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1166 >                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1167 >
1168 >          else
1169 >             pref = pre14 * q_j / 3.0_dp
1170 >             vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1171 >                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1172 >                  qzz_i * (3.0_dp*cz2 - 1.0_dp))
1173 >             vpair = vpair + vterm
1174 >             epot = epot + sw*vterm
1175 >            
1176 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1177 >                  qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1178 >                  qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1179 >                  qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1180 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1181 >                  qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1182 >                  qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1183 >                  qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1184 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1185 >                  qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1186 >                  qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1187 >                  qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1188 >            
1189 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1190 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1191 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1192 >            
1193 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1194 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1195 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1196 >            
1197 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1198 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1199 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1200 >          endif
1201         endif
1202      endif
1203 <      
1204 <    
1203 >
1204 >
1205      if (do_pot) then
1206   #ifdef IS_MPI
1207 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1208 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1207 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1208 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1209   #else
1210         pot = pot + epot
1211   #endif
1212      endif
1213 <        
1213 >
1214   #ifdef IS_MPI
1215      f_Row(1,atom1) = f_Row(1,atom1) + dudx
1216      f_Row(2,atom1) = f_Row(2,atom1) + dudy
1217      f_Row(3,atom1) = f_Row(3,atom1) + dudz
1218 <    
1218 >
1219      f_Col(1,atom2) = f_Col(1,atom2) - dudx
1220      f_Col(2,atom2) = f_Col(2,atom2) - dudy
1221      f_Col(3,atom2) = f_Col(3,atom2) - dudz
1222 <    
1222 >
1223      if (i_is_Dipole .or. i_is_Quadrupole) then
1224         t_Row(1,atom1)=t_Row(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1225         t_Row(2,atom1)=t_Row(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
# Line 765 | Line 1254 | contains
1254      f(1,atom1) = f(1,atom1) + dudx
1255      f(2,atom1) = f(2,atom1) + dudy
1256      f(3,atom1) = f(3,atom1) + dudz
1257 <    
1257 >
1258      f(1,atom2) = f(1,atom2) - dudx
1259      f(2,atom2) = f(2,atom2) - dudy
1260      f(3,atom2) = f(3,atom2) - dudz
1261 <    
1261 >
1262      if (i_is_Dipole .or. i_is_Quadrupole) then
1263         t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1264         t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
# Line 801 | Line 1290 | contains
1290      endif
1291  
1292   #endif
1293 <    
1293 >
1294   #ifdef IS_MPI
1295      id1 = AtomRowToGlobal(atom1)
1296      id2 = AtomColToGlobal(atom2)
# Line 811 | Line 1300 | contains
1300   #endif
1301  
1302      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1303 <      
1303 >
1304         fpair(1) = fpair(1) + dudx
1305         fpair(2) = fpair(2) + dudy
1306         fpair(3) = fpair(3) + dudz
# Line 820 | Line 1309 | contains
1309  
1310      return
1311    end subroutine doElectrostaticPair
1312 <  
1312 >
1313 >  subroutine destroyElectrostaticTypes()
1314 >
1315 >    if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1316 >
1317 >  end subroutine destroyElectrostaticTypes
1318 >
1319 >  subroutine self_self(atom1, eFrame, mypot, t, do_pot)
1320 >    logical, intent(in) :: do_pot
1321 >    integer, intent(in) :: atom1
1322 >    integer :: atid1
1323 >    real(kind=dp), dimension(9,nLocal) :: eFrame
1324 >    real(kind=dp), dimension(3,nLocal) :: t
1325 >    real(kind=dp) :: mu1, c1
1326 >    real(kind=dp) :: preVal, epot, mypot
1327 >    real(kind=dp) :: eix, eiy, eiz
1328 >
1329 >    ! this is a local only array, so we use the local atom type id's:
1330 >    atid1 = atid(atom1)
1331 >
1332 >    if (.not.summationMethodChecked) then
1333 >       call checkSummationMethod()
1334 >    endif
1335 >    
1336 >    if (summationMethod .eq. REACTION_FIELD) then
1337 >       if (ElectrostaticMap(atid1)%is_Dipole) then
1338 >          mu1 = getDipoleMoment(atid1)
1339 >          
1340 >          preVal = pre22 * preRF2 * mu1*mu1
1341 >          mypot = mypot - 0.5d0*preVal
1342 >          
1343 >          ! The self-correction term adds into the reaction field vector
1344 >          
1345 >          eix = preVal * eFrame(3,atom1)
1346 >          eiy = preVal * eFrame(6,atom1)
1347 >          eiz = preVal * eFrame(9,atom1)
1348 >          
1349 >          ! once again, this is self-self, so only the local arrays are needed
1350 >          ! even for MPI jobs:
1351 >          
1352 >          t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1353 >               eFrame(9,atom1)*eiy
1354 >          t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1355 >               eFrame(3,atom1)*eiz
1356 >          t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1357 >               eFrame(6,atom1)*eix
1358 >          
1359 >       endif
1360 >
1361 >    elseif (summationMethod .eq. UNDAMPED_WOLF) then
1362 >       if (ElectrostaticMap(atid1)%is_Charge) then
1363 >          c1 = getCharge(atid1)
1364 >          
1365 >          mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1366 >       endif
1367 >      
1368 >    elseif (summationMethod .eq. DAMPED_WOLF) then
1369 >       if (ElectrostaticMap(atid1)%is_Charge) then
1370 >          c1 = getCharge(atid1)
1371 >          
1372 >          mypot = mypot - (constERFC * rcuti * 0.5_dp + &
1373 >               dampingAlpha*invRootPi) * c1 * c1      
1374 >       endif
1375 >    endif
1376 >    
1377 >    return
1378 >  end subroutine self_self
1379 >
1380 >  subroutine rf_self_excludes(atom1, atom2, sw, eFrame, d, rij, vpair, myPot, &
1381 >       f, t, do_pot)
1382 >    logical, intent(in) :: do_pot
1383 >    integer, intent(in) :: atom1
1384 >    integer, intent(in) :: atom2
1385 >    logical :: i_is_Charge, j_is_Charge
1386 >    logical :: i_is_Dipole, j_is_Dipole
1387 >    integer :: atid1
1388 >    integer :: atid2
1389 >    real(kind=dp), intent(in) :: rij
1390 >    real(kind=dp), intent(in) :: sw
1391 >    real(kind=dp), intent(in), dimension(3) :: d
1392 >    real(kind=dp), intent(inout) :: vpair
1393 >    real(kind=dp), dimension(9,nLocal) :: eFrame
1394 >    real(kind=dp), dimension(3,nLocal) :: f
1395 >    real(kind=dp), dimension(3,nLocal) :: t
1396 >    real (kind = dp), dimension(3) :: duduz_i
1397 >    real (kind = dp), dimension(3) :: duduz_j
1398 >    real (kind = dp), dimension(3) :: uz_i
1399 >    real (kind = dp), dimension(3) :: uz_j
1400 >    real(kind=dp) :: q_i, q_j, mu_i, mu_j
1401 >    real(kind=dp) :: xhat, yhat, zhat
1402 >    real(kind=dp) :: ct_i, ct_j
1403 >    real(kind=dp) :: ri2, ri3, riji, vterm
1404 >    real(kind=dp) :: pref, preVal, rfVal, myPot
1405 >    real(kind=dp) :: dudx, dudy, dudz, dudr
1406 >
1407 >    if (.not.summationMethodChecked) then
1408 >       call checkSummationMethod()
1409 >    endif
1410 >
1411 >    dudx = 0.0d0
1412 >    dudy = 0.0d0
1413 >    dudz = 0.0d0
1414 >
1415 >    riji = 1.0d0/rij
1416 >
1417 >    xhat = d(1) * riji
1418 >    yhat = d(2) * riji
1419 >    zhat = d(3) * riji
1420 >
1421 >    ! this is a local only array, so we use the local atom type id's:
1422 >    atid1 = atid(atom1)
1423 >    atid2 = atid(atom2)
1424 >    i_is_Charge = ElectrostaticMap(atid1)%is_Charge
1425 >    j_is_Charge = ElectrostaticMap(atid2)%is_Charge
1426 >    i_is_Dipole = ElectrostaticMap(atid1)%is_Dipole
1427 >    j_is_Dipole = ElectrostaticMap(atid2)%is_Dipole
1428 >
1429 >    if (i_is_Charge.and.j_is_Charge) then
1430 >       q_i = ElectrostaticMap(atid1)%charge
1431 >       q_j = ElectrostaticMap(atid2)%charge
1432 >      
1433 >       preVal = pre11 * q_i * q_j
1434 >       rfVal = preRF*rij*rij
1435 >       vterm = preVal * rfVal
1436 >      
1437 >       myPot = myPot + sw*vterm
1438 >      
1439 >       dudr  = sw*preVal * 2.0d0*rfVal*riji
1440 >      
1441 >       dudx = dudx + dudr * xhat
1442 >       dudy = dudy + dudr * yhat
1443 >       dudz = dudz + dudr * zhat
1444 >      
1445 >    elseif (i_is_Charge.and.j_is_Dipole) then
1446 >       q_i = ElectrostaticMap(atid1)%charge
1447 >       mu_j = ElectrostaticMap(atid2)%dipole_moment
1448 >       uz_j(1) = eFrame(3,atom2)
1449 >       uz_j(2) = eFrame(6,atom2)
1450 >       uz_j(3) = eFrame(9,atom2)
1451 >       ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
1452 >      
1453 >       ri2 = riji * riji
1454 >       ri3 = ri2 * riji
1455 >      
1456 >       pref = pre12 * q_i * mu_j
1457 >       vterm = - pref * ct_j * ( ri2 - preRF2*rij )
1458 >       myPot = myPot + sw*vterm
1459 >      
1460 >       dudx = dudx - sw*pref*( ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1461 >            - preRF2*uz_j(1) )
1462 >       dudy = dudy - sw*pref*( ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1463 >            - preRF2*uz_j(2) )
1464 >       dudz = dudz - sw*pref*( ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1465 >            - preRF2*uz_j(3) )
1466 >      
1467 >       duduz_j(1) = duduz_j(1) - sw * pref * xhat * ( ri2 - preRF2*rij )
1468 >       duduz_j(2) = duduz_j(2) - sw * pref * yhat * ( ri2 - preRF2*rij )
1469 >       duduz_j(3) = duduz_j(3) - sw * pref * zhat * ( ri2 - preRF2*rij )
1470 >      
1471 >    elseif (i_is_Dipole.and.j_is_Charge) then
1472 >       mu_i = ElectrostaticMap(atid1)%dipole_moment
1473 >       q_j = ElectrostaticMap(atid2)%charge
1474 >       uz_i(1) = eFrame(3,atom1)
1475 >       uz_i(2) = eFrame(6,atom1)
1476 >       uz_i(3) = eFrame(9,atom1)
1477 >       ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
1478 >      
1479 >       ri2 = riji * riji
1480 >       ri3 = ri2 * riji
1481 >      
1482 >       pref = pre12 * q_j * mu_i
1483 >       vterm = pref * ct_i * ( ri2 - preRF2*rij )
1484 >       myPot = myPot + sw*vterm
1485 >      
1486 >       dudx = dudx + sw*pref*( ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1487 >            - preRF2*uz_i(1) )
1488 >       dudy = dudy + sw*pref*( ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1489 >            - preRF2*uz_i(2) )
1490 >       dudz = dudz + sw*pref*( ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1491 >            - preRF2*uz_i(3) )
1492 >      
1493 >       duduz_i(1) = duduz_i(1) + sw * pref * xhat * ( ri2 - preRF2*rij )
1494 >       duduz_i(2) = duduz_i(2) + sw * pref * yhat * ( ri2 - preRF2*rij )
1495 >       duduz_i(3) = duduz_i(3) + sw * pref * zhat * ( ri2 - preRF2*rij )
1496 >      
1497 >    endif
1498 >      
1499 >
1500 >    ! accumulate the forces and torques resulting from the self term
1501 >    f(1,atom1) = f(1,atom1) + dudx
1502 >    f(2,atom1) = f(2,atom1) + dudy
1503 >    f(3,atom1) = f(3,atom1) + dudz
1504 >    
1505 >    f(1,atom2) = f(1,atom2) - dudx
1506 >    f(2,atom2) = f(2,atom2) - dudy
1507 >    f(3,atom2) = f(3,atom2) - dudz
1508 >    
1509 >    if (i_is_Dipole) then
1510 >       t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1511 >       t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1512 >       t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1513 >    elseif (j_is_Dipole) then
1514 >       t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1515 >       t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1516 >       t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1517 >    endif
1518 >
1519 >    return
1520 >  end subroutine rf_self_excludes
1521 >
1522   end module electrostatic_module

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines