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 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2409 by chrisfen, Wed Nov 2 20:36:25 2005 UTC

# Line 54 | Line 54 | module electrostatic_module
54  
55    PRIVATE
56  
57 +
58 + #define __FORTRAN90
59 + #include "UseTheForce/DarkSide/fInteractionMap.h"
60 + #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61 + #include "UseTheForce/DarkSide/fScreeningMethod.h"
62 +
63 +
64    !! these prefactors convert the multipole interactions into kcal / mol
65    !! all were computed assuming distances are measured in angstroms
66    !! Charge-Charge, assuming charges are measured in electrons
# Line 68 | Line 75 | module electrostatic_module
75    !! This unit is also known affectionately as an esu centi-barn.
76    real(kind=dp), parameter :: pre14 = 69.13373_dp
77  
78 +  !! variables to handle different summation methods for long-range electrostatics:
79 +  integer, save :: summationMethod = NONE
80 +  integer, save :: screeningMethod = UNDAMPED
81 +  logical, save :: summationMethodChecked = .false.
82 +  real(kind=DP), save :: defaultCutoff = 0.0_DP
83 +  real(kind=DP), save :: defaultCutoff2 = 0.0_DP
84 +  logical, save :: haveDefaultCutoff = .false.
85 +  real(kind=DP), save :: dampingAlpha = 0.0_DP
86 +  logical, save :: haveDampingAlpha = .false.
87 +  real(kind=DP), save :: dielectric = 1.0_DP
88 +  logical, save :: haveDielectric = .false.
89 +  real(kind=DP), save :: constERFC = 0.0_DP
90 +  real(kind=DP), save :: constEXP = 0.0_DP
91 +  real(kind=dp), save :: rcuti = 0.0_DP
92 +  real(kind=dp), save :: rcuti2 = 0.0_DP
93 +  real(kind=dp), save :: rcuti3 = 0.0_DP
94 +  real(kind=dp), save :: rcuti4 = 0.0_DP
95 +  real(kind=dp), save :: alphaPi = 0.0_DP
96 +  real(kind=dp), save :: invRootPi = 0.0_DP
97 +  real(kind=dp), save :: rrf = 1.0_DP
98 +  real(kind=dp), save :: rt = 1.0_DP
99 +  real(kind=dp), save :: rrfsq = 1.0_DP
100 +  real(kind=dp), save :: preRF = 0.0_DP
101 +  real(kind=dp), save :: preRF2 = 0.0_DP
102 +
103 + #ifdef __IFC
104 + ! error function for ifc version > 7.
105 +  double precision, external :: derfc
106 + #endif
107 +  
108 +  public :: setElectrostaticSummationMethod
109 +  public :: setScreeningMethod
110 +  public :: setElectrostaticCutoffRadius
111 +  public :: setDampingAlpha
112 +  public :: setReactionFieldDielectric
113    public :: newElectrostaticType
114    public :: setCharge
115    public :: setDipoleMoment
# Line 76 | Line 118 | module electrostatic_module
118    public :: doElectrostaticPair
119    public :: getCharge
120    public :: getDipoleMoment
79 –  public :: pre22
121    public :: destroyElectrostaticTypes
122 +  public :: self_self
123 +  public :: rf_self_excludes
124  
125    type :: Electrostatic
126       integer :: c_ident
# Line 85 | Line 128 | module electrostatic_module
128       logical :: is_Dipole = .false.
129       logical :: is_SplitDipole = .false.
130       logical :: is_Quadrupole = .false.
131 +     logical :: is_Tap = .false.
132       real(kind=DP) :: charge = 0.0_DP
133       real(kind=DP) :: dipole_moment = 0.0_DP
134       real(kind=DP) :: split_dipole_distance = 0.0_DP
# Line 95 | Line 139 | contains
139  
140   contains
141  
142 +  subroutine setElectrostaticSummationMethod(the_ESM)
143 +    integer, intent(in) :: the_ESM    
144 +
145 +    if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
146 +       call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
147 +    endif
148 +
149 +    summationMethod = the_ESM
150 +
151 +  end subroutine setElectrostaticSummationMethod
152 +
153 +  subroutine setScreeningMethod(the_SM)
154 +    integer, intent(in) :: the_SM    
155 +    screeningMethod = the_SM
156 +  end subroutine setScreeningMethod
157 +
158 +  subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
159 +    real(kind=dp), intent(in) :: thisRcut
160 +    real(kind=dp), intent(in) :: thisRsw
161 +    defaultCutoff = thisRcut
162 +    rrf = defaultCutoff
163 +    rt = thisRsw
164 +    haveDefaultCutoff = .true.
165 +  end subroutine setElectrostaticCutoffRadius
166 +
167 +  subroutine setDampingAlpha(thisAlpha)
168 +    real(kind=dp), intent(in) :: thisAlpha
169 +    dampingAlpha = thisAlpha
170 +    haveDampingAlpha = .true.
171 +  end subroutine setDampingAlpha
172 +  
173 +  subroutine setReactionFieldDielectric(thisDielectric)
174 +    real(kind=dp), intent(in) :: thisDielectric
175 +    dielectric = thisDielectric
176 +    haveDielectric = .true.
177 +  end subroutine setReactionFieldDielectric
178 +
179    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
180 <       is_SplitDipole, is_Quadrupole, status)
180 >       is_SplitDipole, is_Quadrupole, is_Tap, status)
181  
182      integer, intent(in) :: c_ident
183      logical, intent(in) :: is_Charge
184      logical, intent(in) :: is_Dipole
185      logical, intent(in) :: is_SplitDipole
186      logical, intent(in) :: is_Quadrupole
187 +    logical, intent(in) :: is_Tap
188      integer, intent(out) :: status
189      integer :: nAtypes, myATID, i, j
190  
# Line 139 | Line 221 | contains
221      ElectrostaticMap(myATID)%is_Dipole = is_Dipole
222      ElectrostaticMap(myATID)%is_SplitDipole = is_SplitDipole
223      ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
224 +    ElectrostaticMap(myATID)%is_Tap = is_Tap
225  
226    end subroutine newElectrostaticType
227  
# Line 303 | Line 386 | contains
386      dm = ElectrostaticMap(atid)%dipole_moment
387    end function getDipoleMoment
388  
389 +  subroutine checkSummationMethod()
390 +
391 +    if (.not.haveDefaultCutoff) then
392 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
393 +    endif
394 +
395 +    rcuti = 1.0d0 / defaultCutoff
396 +    rcuti2 = rcuti*rcuti
397 +    rcuti3 = rcuti2*rcuti
398 +    rcuti4 = rcuti2*rcuti2
399 +
400 +    if (screeningMethod .eq. DAMPED) then
401 +       if (.not.haveDampingAlpha) then
402 +          call handleError("checkSummationMethod", "no Damping Alpha set!")
403 +       endif
404 +      
405 +       if (.not.haveDefaultCutoff) then
406 +          call handleError("checkSummationMethod", "no Default Cutoff set!")
407 +       endif
408 +
409 +       constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
410 +       constERFC = derfc(dampingAlpha*defaultCutoff)
411 +       invRootPi = 0.56418958354775628695d0
412 +       alphaPi = 2*dampingAlpha*invRootPi
413 +      
414 +    endif
415 +
416 +    if (summationMethod .eq. REACTION_FIELD) then
417 +       if (haveDielectric) then
418 +          defaultCutoff2 = defaultCutoff*defaultCutoff
419 +          preRF = (dielectric-1.0d0) / &
420 +               ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
421 +          preRF2 = 2.0d0*preRF
422 +       else
423 +          call handleError("checkSummationMethod", "Dielectric not set")
424 +       endif
425 +      
426 +    endif
427 +
428 +    summationMethodChecked = .true.
429 +  end subroutine checkSummationMethod
430 +
431 + !!$
432 + !!$  subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
433 + !!$       vpair, fpair, pot, eFrame, f, t, do_pot)
434    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
435 <       vpair, fpair, pot, eFrame, f, t, do_pot)
435 >       vpair, fpair, pot, eFrame, f, t, do_pot, fstrs)
436  
437      logical, intent(in) :: do_pot
438  
# Line 314 | Line 442 | contains
442      real(kind=dp), intent(in) :: rij, r2, sw
443      real(kind=dp), intent(in), dimension(3) :: d
444      real(kind=dp), intent(inout) :: vpair
445 <    real(kind=dp), intent(inout), dimension(3) :: fpair
445 >    real(kind=dp), intent(inout), dimension(3) :: fpair    
446 >    real(kind=dp), intent(inout), dimension(3) :: fstrs
447  
448      real( kind = dp ) :: pot
449      real( kind = dp ), dimension(9,nLocal) :: eFrame
450      real( kind = dp ), dimension(3,nLocal) :: f
451 +    real( kind = dp ), dimension(3,nLocal) :: felec
452      real( kind = dp ), dimension(3,nLocal) :: t
453  
454      real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
# Line 328 | Line 458 | contains
458  
459      logical :: i_is_Charge, i_is_Dipole, i_is_SplitDipole, i_is_Quadrupole
460      logical :: j_is_Charge, j_is_Dipole, j_is_SplitDipole, j_is_Quadrupole
461 +    logical :: i_is_Tap, j_is_Tap
462      integer :: me1, me2, id1, id2
463      real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
464      real (kind=dp) :: qxx_i, qyy_i, qzz_i
# Line 337 | Line 468 | contains
468      real (kind=dp) :: cx2, cy2, cz2
469      real (kind=dp) :: ct_i, ct_j, ct_ij, a1
470      real (kind=dp) :: riji, ri, ri2, ri3, ri4
471 <    real (kind=dp) :: pref, vterm, epot, dudr    
471 >    real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
472      real (kind=dp) :: xhat, yhat, zhat
473      real (kind=dp) :: dudx, dudy, dudz
474      real (kind=dp) :: scale, sc2, bigR
475 +    real (kind=dp) :: varERFC, varEXP
476 +    real (kind=dp) :: limScale
477 +    real (kind=dp) :: preVal, rfVal
478  
479      if (.not.allocated(ElectrostaticMap)) then
480         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
481         return
482      end if
483  
484 +    if (.not.summationMethodChecked) then
485 +       call checkSummationMethod()
486 +    endif
487 +
488   #ifdef IS_MPI
489      me1 = atid_Row(atom1)
490      me2 = atid_Col(atom2)
# Line 355 | Line 493 | contains
493      me2 = atid(atom2)
494   #endif
495  
496 + !!$    if (rij .ge. defaultCutoff) then
497 + !!$       write(*,*) 'warning: rij = ', rij, ' rcut = ', defaultCutoff, ' sw = ', sw
498 + !!$    endif
499 +
500      !! some variables we'll need independent of electrostatic type:
501  
502      riji = 1.0d0 / rij
503 <
503 >  
504      xhat = d(1) * riji
505      yhat = d(2) * riji
506      zhat = d(3) * riji
507  
508      !! logicals
367 –
509      i_is_Charge = ElectrostaticMap(me1)%is_Charge
510      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
511      i_is_SplitDipole = ElectrostaticMap(me1)%is_SplitDipole
512      i_is_Quadrupole = ElectrostaticMap(me1)%is_Quadrupole
513 +    i_is_Tap = ElectrostaticMap(me1)%is_Tap
514  
515      j_is_Charge = ElectrostaticMap(me2)%is_Charge
516      j_is_Dipole = ElectrostaticMap(me2)%is_Dipole
517      j_is_SplitDipole = ElectrostaticMap(me2)%is_SplitDipole
518      j_is_Quadrupole = ElectrostaticMap(me2)%is_Quadrupole
519 +    j_is_Tap = ElectrostaticMap(me2)%is_Tap
520  
521      if (i_is_Charge) then
522         q_i = ElectrostaticMap(me1)%charge      
# Line 428 | Line 571 | contains
571         cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
572      endif
573  
431 –
574      if (j_is_Charge) then
575         q_j = ElectrostaticMap(me2)%charge      
576      endif
# Line 480 | Line 622 | contains
622         cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
623         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
624      endif
625 <
625 >  
626      epot = 0.0_dp
627      dudx = 0.0_dp
628      dudy = 0.0_dp
# Line 498 | Line 640 | contains
640  
641         if (j_is_Charge) then
642  
643 <          vterm = pre11 * q_i * q_j * riji
502 <          vpair = vpair + vterm
503 <          epot = epot + sw*vterm
643 >          if (summationMethod .eq. SHIFTED_POTENTIAL) then
644  
645 <          dudr  = - sw * vterm * riji
645 >             vterm = pre11 * q_i * q_j * (riji - rcuti)
646 >             vpair = vpair + vterm
647 >             epot = epot + sw*vterm
648 >            
649 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji)
650 >            
651 >             dudx = dudx + dudr * xhat
652 >             dudy = dudy + dudr * yhat
653 >             dudz = dudz + dudr * zhat
654  
655 <          dudx = dudx + dudr * xhat
656 <          dudy = dudy + dudr * yhat
657 <          dudz = dudz + dudr * zhat
655 > !!$          elseif (summationMethod .eq. DAMPED_WOLF) then
656 > !!$             varERFC = derfc(dampingAlpha*rij)
657 > !!$             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
658 > !!$             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
659 > !!$             vpair = vpair + vterm
660 > !!$             epot = epot + sw*vterm
661 > !!$            
662 > !!$!             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
663 > !!$!                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
664 > !!$!                  + alphaPi*constEXP*rcuti)) )
665 > !!$             dudr  = -sw*pre11*q_i*q_j * (varERFC*riji*riji &
666 > !!$                  + alphaPi*varEXP*riji)
667 > !!$            
668 > !!$             dudx = dudx + dudr * xhat
669 > !!$             dudy = dudy + dudr * yhat
670 > !!$             dudz = dudz + dudr * zhat
671  
672 <       endif
672 >          elseif (summationMethod .eq. SHIFTED_FORCE) then
673 >             vterm = pre11 * q_i * q_j * (riji + rij*rcuti2 - 2.0d0*rcuti)
674 >             vpair = vpair + vterm
675 >             epot = epot + sw*vterm
676 >            
677 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)
678 >                          
679 >             dudx = dudx + dudr * xhat
680 >             dudy = dudy + dudr * yhat
681 >             dudz = dudz + dudr * zhat
682  
683 <       if (j_is_Dipole) then
683 > !!$          elseif (summationMethod .eq. DAMPED_SHIFTED_FORCE) then
684 > !!$             varERFC = derfc(dampingAlpha*rij)
685 > !!$             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
686 > !!$             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
687 > !!$             vpair = vpair + vterm
688 > !!$             epot = epot + sw*vterm
689 > !!$            
690 > !!$!             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
691 > !!$!                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
692 > !!$!                  + alphaPi*constEXP*rcuti)) )
693 > !!$             dudr  = -sw*pre11*q_i*q_j * (varERFC*riji*riji &
694 > !!$                  + alphaPi*varEXP*riji)
695 > !!$            
696 > !!$             dudx = dudx + dudr * xhat
697 > !!$             dudy = dudy + dudr * yhat
698 > !!$             dudz = dudz + dudr * zhat
699  
700 <          if (j_is_SplitDipole) then
701 <             BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
702 <             ri = 1.0_dp / BigR
703 <             scale = rij * ri
700 >          elseif (summationMethod .eq. REACTION_FIELD) then
701 >             preVal = pre11 * q_i * q_j
702 >             rfVal = preRF*rij*rij
703 >             vterm = preVal * ( riji + rfVal )
704 >            
705 >             vpair = vpair + vterm
706 >             epot = epot + sw*vterm
707 >            
708 >             dudr  = sw * preVal * ( 2.0d0*rfVal - riji )*riji
709 >            
710 >             dudx = dudx + dudr * xhat
711 >             dudy = dudy + dudr * yhat
712 >             dudz = dudz + dudr * zhat
713 >
714            else
715 <             ri = riji
716 <             scale = 1.0_dp
715 >             vterm = pre11 * q_i * q_j * riji
716 >             vpair = vpair + vterm
717 >             epot = epot + sw*vterm
718 >            
719 >             dudr  = - sw * vterm * riji
720 >            
721 >             dudx = dudx + dudr * xhat
722 >             dudy = dudy + dudr * yhat
723 >             dudz = dudz + dudr * zhat
724 >
725            endif
726  
727 <          ri2 = ri * ri
525 <          ri3 = ri2 * ri
526 <          sc2 = scale * scale
727 >       endif
728  
729 +       if (j_is_Dipole) then
730 +
731            pref = pre12 * q_i * mu_j
529 –          vterm = - pref * ct_j * ri2 * scale
530 –          vpair = vpair + vterm
531 –          epot = epot + sw * vterm
732  
733 <          !! this has a + sign in the () because the rij vector is
734 <          !! r_j - r_i and the charge-dipole potential takes the origin
735 <          !! as the point dipole, which is atom j in this case.
736 <
737 <          dudx = dudx - pref * sw * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
738 <          dudy = dudy - pref * sw * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
739 <          dudz = dudz - pref * sw * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
740 <
741 <          duduz_j(1) = duduz_j(1) - pref * sw * ri2 * xhat * scale
742 <          duduz_j(2) = duduz_j(2) - pref * sw * ri2 * yhat * scale
743 <          duduz_j(3) = duduz_j(3) - pref * sw * ri2 * zhat * scale
733 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
734 > !!$             ri2 = riji * riji
735 > !!$             ri3 = ri2 * riji
736 > !!$
737 > !!$             pref = pre12 * q_i * mu_j
738 > !!$             vterm = - pref * ct_j * (ri2 - rcuti2)
739 > !!$             vpair = vpair + vterm
740 > !!$             epot = epot + sw*vterm
741 > !!$            
742 > !!$             !! this has a + sign in the () because the rij vector is
743 > !!$             !! r_j - r_i and the charge-dipole potential takes the origin
744 > !!$             !! as the point dipole, which is atom j in this case.
745 > !!$            
746 > !!$             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
747 > !!$                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
748 > !!$             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
749 > !!$                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
750 > !!$             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
751 > !!$                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
752 > !!$            
753 > !!$             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
754 > !!$             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
755 > !!$             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
756 > !!$
757 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
758  
759 +          if (summationMethod .eq. REACTION_FIELD) then
760 +             ri2 = riji * riji
761 +             ri3 = ri2 * riji
762 +    
763 +             pref = pre12 * q_i * mu_j
764 +             vterm = - pref * ct_j * ( ri2 - preRF2*rij )
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) - &
777 +                                     preRF2*uz_j(3) )        
778 +             duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
779 +             duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
780 +             duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
781 +
782 +          else
783 +             if (j_is_SplitDipole) then
784 +                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
785 +                ri = 1.0_dp / BigR
786 +                scale = rij * ri
787 +             else
788 +                ri = riji
789 +                scale = 1.0_dp
790 +             endif
791 +            
792 +             ri2 = ri * ri
793 +             ri3 = ri2 * ri
794 +             sc2 = scale * scale
795 +
796 +             pref = pre12 * q_i * mu_j
797 +             vterm = - pref * ct_j * ri2 * scale
798 +             vpair = vpair + vterm
799 +             epot = epot + sw*vterm
800 +            
801 +             !! this has a + sign in the () because the rij vector is
802 +             !! r_j - r_i and the charge-dipole potential takes the origin
803 +             !! as the point dipole, which is atom j in this case.
804 +            
805 +             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
806 +             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
807 +             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
808 +            
809 +             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
810 +             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
811 +             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
812 +
813 +          endif
814         endif
815  
816         if (j_is_Quadrupole) then
# Line 552 | Line 821 | contains
821            cy2 = cy_j * cy_j
822            cz2 = cz_j * cz_j
823  
824 <
825 <          pref =  pre14 * q_i / 3.0_dp
826 <          vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
827 <               qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
828 <               qzz_j * (3.0_dp*cz2 - 1.0_dp))
829 <          vpair = vpair + vterm
830 <          epot = epot + sw * vterm
831 <
832 <          dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( &
833 <               qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
834 <               qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
835 <               qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
836 <          dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + pref * sw * ri4 * ( &
837 <               qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
838 <               qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
839 <               qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
840 <          dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + pref * sw * ri4 * ( &
841 <               qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
842 <               qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
843 <               qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
844 <
845 <          dudux_j(1) = dudux_j(1) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*xhat)
846 <          dudux_j(2) = dudux_j(2) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*yhat)
847 <          dudux_j(3) = dudux_j(3) + pref * sw * ri3 * (qxx_j*6.0_dp*cx_j*zhat)
848 <
849 <          duduy_j(1) = duduy_j(1) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*xhat)
850 <          duduy_j(2) = duduy_j(2) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*yhat)
851 <          duduy_j(3) = duduy_j(3) + pref * sw * ri3 * (qyy_j*6.0_dp*cy_j*zhat)
852 <
853 <          duduz_j(1) = duduz_j(1) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*xhat)
854 <          duduz_j(2) = duduz_j(2) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*yhat)
855 <          duduz_j(3) = duduz_j(3) + pref * sw * ri3 * (qzz_j*6.0_dp*cz_j*zhat)
824 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
825 > !!$             pref =  pre14 * q_i / 3.0_dp
826 > !!$             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
827 > !!$                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
828 > !!$                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
829 > !!$             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
830 > !!$                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
831 > !!$                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
832 > !!$             vpair = vpair + ( vterm1 - vterm2 )
833 > !!$             epot = epot + sw*( vterm1 - vterm2 )
834 > !!$            
835 > !!$             dudx = dudx - (5.0_dp * &
836 > !!$                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
837 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
838 > !!$                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
839 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
840 > !!$                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
841 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
842 > !!$                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
843 > !!$             dudy = dudy - (5.0_dp * &
844 > !!$                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
845 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
846 > !!$                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
847 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
848 > !!$                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
849 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
850 > !!$                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
851 > !!$             dudz = dudz - (5.0_dp * &
852 > !!$                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
853 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
854 > !!$                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
855 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
856 > !!$                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
857 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
858 > !!$                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
859 > !!$            
860 > !!$             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
861 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
862 > !!$             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
863 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
864 > !!$             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
865 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
866 > !!$            
867 > !!$             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
868 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
869 > !!$             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
870 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
871 > !!$             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
872 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
873 > !!$            
874 > !!$             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
875 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
876 > !!$             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
877 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
878 > !!$             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
879 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
880 > !!$        
881 > !!$          else
882 >             pref =  pre14 * q_i / 3.0_dp
883 >             vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
884 >                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
885 >                  qzz_j * (3.0_dp*cz2 - 1.0_dp))
886 >             vpair = vpair + vterm
887 >             epot = epot + sw*vterm
888 >            
889 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
890 >                  qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
891 >                  qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
892 >                  qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
893 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
894 >                  qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
895 >                  qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
896 >                  qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
897 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
898 >                  qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
899 >                  qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
900 >                  qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
901 >            
902 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
903 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
904 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
905 >            
906 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
907 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
908 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
909 >            
910 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
911 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
912 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
913 >          
914 > !!$          endif
915         endif
588 –
916      endif
917  
918      if (i_is_Dipole) then
919  
920         if (j_is_Charge) then
921 <
595 <          if (i_is_SplitDipole) then
596 <             BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
597 <             ri = 1.0_dp / BigR
598 <             scale = rij * ri
599 <          else
600 <             ri = riji
601 <             scale = 1.0_dp
602 <          endif
603 <
604 <          ri2 = ri * ri
605 <          ri3 = ri2 * ri
606 <          sc2 = scale * scale
607 <
921 >          
922            pref = pre12 * q_j * mu_i
923 <          vterm = pref * ct_i * ri2 * scale
924 <          vpair = vpair + vterm
925 <          epot = epot + sw * vterm
923 >          
924 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
925 > !!$             ri2 = riji * riji
926 > !!$             ri3 = ri2 * riji
927 > !!$
928 > !!$             pref = pre12 * q_j * mu_i
929 > !!$             vterm = pref * ct_i * (ri2 - rcuti2)
930 > !!$             vpair = vpair + vterm
931 > !!$             epot = epot + sw*vterm
932 > !!$            
933 > !!$             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
934 > !!$                  - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
935 > !!$             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
936 > !!$                  - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
937 > !!$             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
938 > !!$                  - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
939 > !!$            
940 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*( ri2*xhat - d(1)*rcuti3 )
941 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*( ri2*yhat - d(2)*rcuti3 )
942 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*( ri2*zhat - d(3)*rcuti3 )
943 > !!$
944 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
945 >          if (summationMethod .eq. REACTION_FIELD) then
946 >             ri2 = riji * riji
947 >             ri3 = ri2 * riji
948  
949 <          dudx = dudx + pref * sw * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
950 <          dudy = dudy + pref * sw * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
951 <          dudz = dudz + pref * sw * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
949 >             pref = pre12 * q_j * mu_i
950 >             vterm = pref * ct_i * ( ri2 - preRF2*rij )
951 >             vpair = vpair + vterm
952 >             epot = epot + sw*vterm
953 >            
954 >             dudx = dudx + sw*pref * ( ri3*(uz_i(1) - 3.0d0*ct_i*xhat) - &
955 >                  preRF2*uz_i(1) )
956 >             dudy = dudy + sw*pref * ( ri3*(uz_i(2) - 3.0d0*ct_i*yhat) - &
957 >                  preRF2*uz_i(2) )
958 >             dudz = dudz + sw*pref * ( ri3*(uz_i(3) - 3.0d0*ct_i*zhat) - &
959 >                  preRF2*uz_i(3) )
960 >            
961 >             duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF2*rij )
962 >             duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF2*rij )
963 >             duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF2*rij )
964  
617 –          duduz_i(1) = duduz_i(1) + pref * sw * ri2 * xhat * scale
618 –          duduz_i(2) = duduz_i(2) + pref * sw * ri2 * yhat * scale
619 –          duduz_i(3) = duduz_i(3) + pref * sw * ri2 * zhat * scale
620 –       endif
621 –
622 –       if (j_is_Dipole) then
623 –
624 –          if (i_is_SplitDipole) then
625 –             if (j_is_SplitDipole) then
626 –                BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
627 –             else
628 –                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
629 –             endif
630 –             ri = 1.0_dp / BigR
631 –             scale = rij * ri                
965            else
966 <             if (j_is_SplitDipole) then
967 <                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
966 >             if (i_is_SplitDipole) then
967 >                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
968                  ri = 1.0_dp / BigR
969 <                scale = rij * ri                            
970 <             else                
969 >                scale = rij * ri
970 >             else
971                  ri = riji
972                  scale = 1.0_dp
973               endif
974 +            
975 +             ri2 = ri * ri
976 +             ri3 = ri2 * ri
977 +             sc2 = scale * scale
978 +
979 +             pref = pre12 * q_j * mu_i
980 +             vterm = pref * ct_i * ri2 * scale
981 +             vpair = vpair + vterm
982 +             epot = epot + sw*vterm
983 +            
984 +             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
985 +             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
986 +             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
987 +            
988 +             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
989 +             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
990 +             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
991            endif
992 +       endif
993 +      
994 +       if (j_is_Dipole) then
995  
996 <          ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
996 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
997 > !!$             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
998 > !!$
999 > !!$             ri2 = riji * riji
1000 > !!$             ri3 = ri2 * riji
1001 > !!$             ri4 = ri2 * ri2
1002 > !!$
1003 > !!$             pref = pre22 * mu_i * mu_j
1004 > !!$             vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
1005 > !!$             vpair = vpair + vterm
1006 > !!$             epot = epot + sw*vterm
1007 > !!$            
1008 > !!$             a1 = 5.0d0 * ct_i * ct_j - ct_ij
1009 > !!$            
1010 > !!$             dudx = dudx + sw*pref*3.0d0*( &
1011 > !!$                  ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
1012 > !!$                  - rcuti4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) )
1013 > !!$             dudy = dudy + sw*pref*3.0d0*( &
1014 > !!$                  ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
1015 > !!$                  - rcuti4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) )
1016 > !!$             dudz = dudz + sw*pref*3.0d0*( &
1017 > !!$                  ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
1018 > !!$                  - rcuti4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) )
1019 > !!$            
1020 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1021 > !!$                  - rcuti3*(uz_j(1) - 3.0d0*ct_j*xhat))
1022 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1023 > !!$                  - rcuti3*(uz_j(2) - 3.0d0*ct_j*yhat))
1024 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1025 > !!$                  - rcuti3*(uz_j(3) - 3.0d0*ct_j*zhat))
1026 > !!$             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1027 > !!$                  - rcuti3*(uz_i(1) - 3.0d0*ct_i*xhat))
1028 > !!$             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1029 > !!$                  - rcuti3*(uz_i(2) - 3.0d0*ct_i*yhat))
1030 > !!$             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1031 > !!$                  - rcuti3*(uz_i(3) - 3.0d0*ct_i*zhat))
1032 > !!$          
1033 > !!$          elseif (summationMethod .eq. DAMPED_WOLF) then
1034 > !!$             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1035 > !!$            
1036 > !!$             ri2 = riji * riji
1037 > !!$             ri3 = ri2 * riji
1038 > !!$             ri4 = ri2 * ri2
1039 > !!$             sc2 = scale * scale
1040 > !!$            
1041 > !!$             pref = pre22 * mu_i * mu_j
1042 > !!$             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j)
1043 > !!$             vpair = vpair + vterm
1044 > !!$             epot = epot + sw*vterm
1045 > !!$            
1046 > !!$             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1047 > !!$            
1048 > !!$             dudx = dudx + sw*pref*3.0d0*ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1049 > !!$             dudy = dudy + sw*pref*3.0d0*ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1050 > !!$             dudz = dudz + sw*pref*3.0d0*ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1051 > !!$            
1052 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*ri3 *(uz_j(1) - 3.0d0*ct_j*xhat)
1053 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*ri3 *(uz_j(2) - 3.0d0*ct_j*yhat)
1054 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*ri3 *(uz_j(3) - 3.0d0*ct_j*zhat)
1055 > !!$            
1056 > !!$             duduz_j(1) = duduz_j(1) + sw*pref*ri3 *(uz_i(1) - 3.0d0*ct_i*xhat)
1057 > !!$             duduz_j(2) = duduz_j(2) + sw*pref*ri3 *(uz_i(2) - 3.0d0*ct_i*yhat)
1058 > !!$             duduz_j(3) = duduz_j(3) + sw*pref*ri3 *(uz_i(3) - 3.0d0*ct_i*zhat)
1059 > !!$            
1060 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
1061 >          if (summationMethod .eq. REACTION_FIELD) then
1062 >             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1063  
1064 <          ri2 = ri * ri
1065 <          ri3 = ri2 * ri
1066 <          ri4 = ri2 * ri2
648 <          sc2 = scale * scale
1064 >             ri2 = riji * riji
1065 >             ri3 = ri2 * riji
1066 >             ri4 = ri2 * ri2
1067  
1068 <          pref = pre22 * mu_i * mu_j
1069 <          vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1070 <          vpair = vpair + vterm
1071 <          epot = epot + sw * vterm
1068 >             pref = pre22 * mu_i * mu_j
1069 >              
1070 >             vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
1071 >                  preRF2*ct_ij )
1072 >             vpair = vpair + vterm
1073 >             epot = epot + sw*vterm
1074 >            
1075 >             a1 = 5.0d0 * ct_i * ct_j - ct_ij
1076 >            
1077 >             dudx = dudx + sw*pref*3.0d0*ri4 &
1078 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1079 >             dudy = dudy + sw*pref*3.0d0*ri4 &
1080 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1081 >             dudz = dudz + sw*pref*3.0d0*ri4 &
1082 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1083 >            
1084 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1085 >                  - preRF2*uz_j(1))
1086 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1087 >                  - preRF2*uz_j(2))
1088 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1089 >                  - preRF2*uz_j(3))
1090 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1091 >                  - preRF2*uz_i(1))
1092 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1093 >                  - preRF2*uz_i(2))
1094 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1095 >                  - preRF2*uz_i(3))
1096  
1097 <          a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1098 <
1099 <          dudx=dudx+pref*sw*3.0d0*ri4*scale*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1100 <          dudy=dudy+pref*sw*3.0d0*ri4*scale*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1101 <          dudz=dudz+pref*sw*3.0d0*ri4*scale*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1102 <
1103 <          duduz_i(1) = duduz_i(1) + pref*sw*ri3*(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1104 <          duduz_i(2) = duduz_i(2) + pref*sw*ri3*(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1105 <          duduz_i(3) = duduz_i(3) + pref*sw*ri3*(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1106 <
1107 <          duduz_j(1) = duduz_j(1) + pref*sw*ri3*(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1108 <          duduz_j(2) = duduz_j(2) + pref*sw*ri3*(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1109 <          duduz_j(3) = duduz_j(3) + pref*sw*ri3*(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1097 >          else
1098 >             if (i_is_SplitDipole) then
1099 >                if (j_is_SplitDipole) then
1100 >                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
1101 >                else
1102 >                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
1103 >                endif
1104 >                ri = 1.0_dp / BigR
1105 >                scale = rij * ri                
1106 >             else
1107 >                if (j_is_SplitDipole) then
1108 >                   BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
1109 >                   ri = 1.0_dp / BigR
1110 >                   scale = rij * ri                            
1111 >                else                
1112 >                   ri = riji
1113 >                   scale = 1.0_dp
1114 >                endif
1115 >             endif
1116 >            
1117 >             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1118 >            
1119 >             ri2 = ri * ri
1120 >             ri3 = ri2 * ri
1121 >             ri4 = ri2 * ri2
1122 >             sc2 = scale * scale
1123 >            
1124 >             pref = pre22 * mu_i * mu_j
1125 >             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
1126 >             vpair = vpair + vterm
1127 >             epot = epot + sw*vterm
1128 >            
1129 >             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1130 >            
1131 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
1132 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1133 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
1134 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1135 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
1136 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1137 >            
1138 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
1139 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
1140 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
1141 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
1142 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
1143 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
1144 >            
1145 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
1146 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
1147 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
1148 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
1149 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
1150 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
1151 >          endif
1152         endif
669 –
1153      endif
1154  
1155      if (i_is_Quadrupole) then
# Line 679 | Line 1162 | contains
1162            cy2 = cy_i * cy_i
1163            cz2 = cz_i * cz_i
1164  
1165 <          pref = pre14 * q_j / 3.0_dp
1166 <          vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1167 <               qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1168 <               qzz_i * (3.0_dp*cz2 - 1.0_dp))
1169 <          vpair = vpair + vterm
1170 <          epot = epot + sw * vterm
1171 <
1172 <          dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + pref * sw * ri4 * ( &
1173 <               qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1174 <               qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1175 <               qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1176 <          dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + pref * sw * ri4 * ( &
1177 <               qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1178 <               qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1179 <               qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1180 <          dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + pref * sw * ri4 * ( &
1181 <               qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1182 <               qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1183 <               qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1184 <
1185 <          dudux_i(1) = dudux_i(1) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*xhat)
1186 <          dudux_i(2) = dudux_i(2) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*yhat)
1187 <          dudux_i(3) = dudux_i(3) + pref * sw * ri3 * (qxx_i*6.0_dp*cx_i*zhat)
1188 <
1189 <          duduy_i(1) = duduy_i(1) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*xhat)
1190 <          duduy_i(2) = duduy_i(2) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*yhat)
1191 <          duduy_i(3) = duduy_i(3) + pref * sw * ri3 * (qyy_i*6.0_dp*cy_i*zhat)
1192 <
1193 <          duduz_i(1) = duduz_i(1) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*xhat)
1194 <          duduz_i(2) = duduz_i(2) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*yhat)
1195 <          duduz_i(3) = duduz_i(3) + pref * sw * ri3 * (qzz_i*6.0_dp*cz_i*zhat)
1165 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
1166 > !!$             pref = pre14 * q_j / 3.0_dp
1167 > !!$             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1168 > !!$                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1169 > !!$                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1170 > !!$             vterm2 = pref * rcuti3*( 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 + ( vterm1 - vterm2 )
1174 > !!$             epot = epot + sw*( vterm1 - vterm2 )
1175 > !!$            
1176 > !!$             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1177 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1178 > !!$                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1179 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1180 > !!$                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1181 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1182 > !!$                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1183 > !!$             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1184 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1185 > !!$                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1186 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1187 > !!$                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1188 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1189 > !!$                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1190 > !!$             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1191 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1192 > !!$                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1193 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1194 > !!$                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1195 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1196 > !!$                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1197 > !!$            
1198 > !!$             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1199 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1200 > !!$             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1201 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1202 > !!$             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1203 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1204 > !!$            
1205 > !!$             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1206 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1207 > !!$             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1208 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1209 > !!$             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1210 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1211 > !!$            
1212 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1213 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1214 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1215 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1216 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1217 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1218 > !!$
1219 > !!$          else
1220 >             pref = pre14 * q_j / 3.0_dp
1221 >             vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1222 >                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1223 >                  qzz_i * (3.0_dp*cz2 - 1.0_dp))
1224 >             vpair = vpair + vterm
1225 >             epot = epot + sw*vterm
1226 >            
1227 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1228 >                  qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1229 >                  qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1230 >                  qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1231 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1232 >                  qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1233 >                  qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1234 >                  qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1235 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1236 >                  qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1237 >                  qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1238 >                  qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1239 >            
1240 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1241 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1242 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1243 >            
1244 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1245 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1246 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1247 >            
1248 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1249 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1250 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1251 > !!$          endif
1252         endif
1253      endif
1254  
1255  
1256      if (do_pot) then
1257   #ifdef IS_MPI
1258 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1259 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1258 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1259 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1260   #else
1261         pot = pot + epot
1262   #endif
# Line 822 | Line 1361 | contains
1361      return
1362    end subroutine doElectrostaticPair
1363  
825 –
1364    subroutine destroyElectrostaticTypes()
1365  
1366      if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1367  
1368    end subroutine destroyElectrostaticTypes
1369  
1370 +  subroutine self_self(atom1, eFrame, mypot, t, do_pot)
1371 +    logical, intent(in) :: do_pot
1372 +    integer, intent(in) :: atom1
1373 +    integer :: atid1
1374 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1375 +    real(kind=dp), dimension(3,nLocal) :: t
1376 +    real(kind=dp) :: mu1, c1
1377 +    real(kind=dp) :: preVal, epot, mypot
1378 +    real(kind=dp) :: eix, eiy, eiz
1379 +
1380 +    ! this is a local only array, so we use the local atom type id's:
1381 +    atid1 = atid(atom1)
1382 +
1383 +    if (.not.summationMethodChecked) then
1384 +       call checkSummationMethod()
1385 +    endif
1386 +    
1387 +    if (summationMethod .eq. REACTION_FIELD) then
1388 +       if (ElectrostaticMap(atid1)%is_Dipole) then
1389 +          mu1 = getDipoleMoment(atid1)
1390 +          
1391 +          preVal = pre22 * preRF2 * mu1*mu1
1392 +          mypot = mypot - 0.5d0*preVal
1393 +          
1394 +          ! The self-correction term adds into the reaction field vector
1395 +          
1396 +          eix = preVal * eFrame(3,atom1)
1397 +          eiy = preVal * eFrame(6,atom1)
1398 +          eiz = preVal * eFrame(9,atom1)
1399 +          
1400 +          ! once again, this is self-self, so only the local arrays are needed
1401 +          ! even for MPI jobs:
1402 +          
1403 +          t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1404 +               eFrame(9,atom1)*eiy
1405 +          t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1406 +               eFrame(3,atom1)*eiz
1407 +          t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1408 +               eFrame(6,atom1)*eix
1409 +          
1410 +       endif
1411 +
1412 + !!$    elseif (summationMethod .eq. UNDAMPED_WOLF) then
1413 + !!$       if (ElectrostaticMap(atid1)%is_Charge) then
1414 + !!$          c1 = getCharge(atid1)
1415 + !!$          
1416 + !!$          mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1417 + !!$       endif
1418 + !!$      
1419 + !!$    elseif (summationMethod .eq. DAMPED_WOLF) then
1420 + !!$       if (ElectrostaticMap(atid1)%is_Charge) then
1421 + !!$          c1 = getCharge(atid1)
1422 + !!$          
1423 + !!$          mypot = mypot - (constERFC * rcuti * 0.5_dp + &
1424 + !!$               dampingAlpha*invRootPi) * c1 * c1      
1425 + !!$       endif
1426 +    endif
1427 +    
1428 +    return
1429 +  end subroutine self_self
1430 +
1431 +  subroutine rf_self_excludes(atom1, atom2, sw, eFrame, d, rij, vpair, myPot, &
1432 +       f, t, do_pot)
1433 +    logical, intent(in) :: do_pot
1434 +    integer, intent(in) :: atom1
1435 +    integer, intent(in) :: atom2
1436 +    logical :: i_is_Charge, j_is_Charge
1437 +    logical :: i_is_Dipole, j_is_Dipole
1438 +    integer :: atid1
1439 +    integer :: atid2
1440 +    real(kind=dp), intent(in) :: rij
1441 +    real(kind=dp), intent(in) :: sw
1442 +    real(kind=dp), intent(in), dimension(3) :: d
1443 +    real(kind=dp), intent(inout) :: vpair
1444 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1445 +    real(kind=dp), dimension(3,nLocal) :: f
1446 +    real(kind=dp), dimension(3,nLocal) :: t
1447 +    real (kind = dp), dimension(3) :: duduz_i
1448 +    real (kind = dp), dimension(3) :: duduz_j
1449 +    real (kind = dp), dimension(3) :: uz_i
1450 +    real (kind = dp), dimension(3) :: uz_j
1451 +    real(kind=dp) :: q_i, q_j, mu_i, mu_j
1452 +    real(kind=dp) :: xhat, yhat, zhat
1453 +    real(kind=dp) :: ct_i, ct_j
1454 +    real(kind=dp) :: ri2, ri3, riji, vterm
1455 +    real(kind=dp) :: pref, preVal, rfVal, myPot
1456 +    real(kind=dp) :: dudx, dudy, dudz, dudr
1457 +
1458 +    if (.not.summationMethodChecked) then
1459 +       call checkSummationMethod()
1460 +    endif
1461 +
1462 +    dudx = 0.0d0
1463 +    dudy = 0.0d0
1464 +    dudz = 0.0d0
1465 +
1466 +    riji = 1.0d0/rij
1467 +
1468 +    xhat = d(1) * riji
1469 +    yhat = d(2) * riji
1470 +    zhat = d(3) * riji
1471 +
1472 +    ! this is a local only array, so we use the local atom type id's:
1473 +    atid1 = atid(atom1)
1474 +    atid2 = atid(atom2)
1475 +    i_is_Charge = ElectrostaticMap(atid1)%is_Charge
1476 +    j_is_Charge = ElectrostaticMap(atid2)%is_Charge
1477 +    i_is_Dipole = ElectrostaticMap(atid1)%is_Dipole
1478 +    j_is_Dipole = ElectrostaticMap(atid2)%is_Dipole
1479 +
1480 +    if (i_is_Charge.and.j_is_Charge) then
1481 +       q_i = ElectrostaticMap(atid1)%charge
1482 +       q_j = ElectrostaticMap(atid2)%charge
1483 +      
1484 +       preVal = pre11 * q_i * q_j
1485 +       rfVal = preRF*rij*rij
1486 +       vterm = preVal * rfVal
1487 +      
1488 +       myPot = myPot + sw*vterm
1489 +      
1490 +       dudr  = sw*preVal * 2.0d0*rfVal*riji
1491 +      
1492 +       dudx = dudx + dudr * xhat
1493 +       dudy = dudy + dudr * yhat
1494 +       dudz = dudz + dudr * zhat
1495 +      
1496 +    elseif (i_is_Charge.and.j_is_Dipole) then
1497 +       q_i = ElectrostaticMap(atid1)%charge
1498 +       mu_j = ElectrostaticMap(atid2)%dipole_moment
1499 +       uz_j(1) = eFrame(3,atom2)
1500 +       uz_j(2) = eFrame(6,atom2)
1501 +       uz_j(3) = eFrame(9,atom2)
1502 +       ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
1503 +      
1504 +       ri2 = riji * riji
1505 +       ri3 = ri2 * riji
1506 +      
1507 +       pref = pre12 * q_i * mu_j
1508 +       vterm = - pref * ct_j * ( ri2 - preRF2*rij )
1509 +       myPot = myPot + sw*vterm
1510 +      
1511 +       dudx = dudx - sw*pref*( ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1512 +            - preRF2*uz_j(1) )
1513 +       dudy = dudy - sw*pref*( ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1514 +            - preRF2*uz_j(2) )
1515 +       dudz = dudz - sw*pref*( ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1516 +            - preRF2*uz_j(3) )
1517 +      
1518 +       duduz_j(1) = duduz_j(1) - sw * pref * xhat * ( ri2 - preRF2*rij )
1519 +       duduz_j(2) = duduz_j(2) - sw * pref * yhat * ( ri2 - preRF2*rij )
1520 +       duduz_j(3) = duduz_j(3) - sw * pref * zhat * ( ri2 - preRF2*rij )
1521 +      
1522 +    elseif (i_is_Dipole.and.j_is_Charge) then
1523 +       mu_i = ElectrostaticMap(atid1)%dipole_moment
1524 +       q_j = ElectrostaticMap(atid2)%charge
1525 +       uz_i(1) = eFrame(3,atom1)
1526 +       uz_i(2) = eFrame(6,atom1)
1527 +       uz_i(3) = eFrame(9,atom1)
1528 +       ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
1529 +      
1530 +       ri2 = riji * riji
1531 +       ri3 = ri2 * riji
1532 +      
1533 +       pref = pre12 * q_j * mu_i
1534 +       vterm = pref * ct_i * ( ri2 - preRF2*rij )
1535 +       myPot = myPot + sw*vterm
1536 +      
1537 +       dudx = dudx + sw*pref*( ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1538 +            - preRF2*uz_i(1) )
1539 +       dudy = dudy + sw*pref*( ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1540 +            - preRF2*uz_i(2) )
1541 +       dudz = dudz + sw*pref*( ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1542 +            - preRF2*uz_i(3) )
1543 +      
1544 +       duduz_i(1) = duduz_i(1) + sw * pref * xhat * ( ri2 - preRF2*rij )
1545 +       duduz_i(2) = duduz_i(2) + sw * pref * yhat * ( ri2 - preRF2*rij )
1546 +       duduz_i(3) = duduz_i(3) + sw * pref * zhat * ( ri2 - preRF2*rij )
1547 +      
1548 +    endif
1549 +      
1550 +
1551 +    ! accumulate the forces and torques resulting from the self term
1552 +    f(1,atom1) = f(1,atom1) + dudx
1553 +    f(2,atom1) = f(2,atom1) + dudy
1554 +    f(3,atom1) = f(3,atom1) + dudz
1555 +    
1556 +    f(1,atom2) = f(1,atom2) - dudx
1557 +    f(2,atom2) = f(2,atom2) - dudy
1558 +    f(3,atom2) = f(3,atom2) - dudz
1559 +    
1560 +    if (i_is_Dipole) then
1561 +       t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1562 +       t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1563 +       t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1564 +    elseif (j_is_Dipole) then
1565 +       t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1566 +       t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1567 +       t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1568 +    endif
1569 +
1570 +    return
1571 +  end subroutine rf_self_excludes
1572 +
1573   end module electrostatic_module

Diff Legend

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