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 2302 by chrisfen, Fri Sep 16 16:07:39 2005 UTC vs.
Revision 2355 by chuckv, Wed Oct 12 18:59:16 2005 UTC

# Line 54 | Line 54 | module electrostatic_module
54  
55    PRIVATE
56  
57 +
58   #define __FORTRAN90
59 + #include "UseTheForce/DarkSide/fInteractionMap.h"
60   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61  
62 +
63    !! these prefactors convert the multipole interactions into kcal / mol
64    !! all were computed assuming distances are measured in angstroms
65    !! Charge-Charge, assuming charges are measured in electrons
# Line 83 | Line 86 | module electrostatic_module
86    real(kind=DP), save :: constERFC = 0.0_DP
87    real(kind=DP), save :: constEXP = 0.0_DP
88    logical, save :: haveDWAconstants = .false.
89 <
90 <
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 >  
96 > #ifdef __IFC
97 > ! error function for ifc version > 7.
98 >  double precision, external :: derfc
99 > #endif
100 >  
101    public :: setElectrostaticSummationMethod
102    public :: setElectrostaticCutoffRadius
103    public :: setDampedWolfAlpha
# Line 118 | Line 131 | contains
131   contains
132  
133    subroutine setElectrostaticSummationMethod(the_ESM)
121 –
134      integer, intent(in) :: the_ESM    
135  
136      if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
137         call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
138      endif
139  
140 +    summationMethod = the_ESM
141 +
142    end subroutine setElectrostaticSummationMethod
143  
144    subroutine setElectrostaticCutoffRadius(thisRcut)
# Line 356 | Line 370 | contains
370    end function getDipoleMoment
371  
372    subroutine checkSummationMethod()
373 +
374 +    if (.not.haveDefaultCutoff) then
375 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
376 +    endif
377 +
378 +    rcuti = 1.0d0 / defaultCutoff
379 +    rcuti2 = rcuti*rcuti
380 +    rcuti3 = rcuti2*rcuti
381 +    rcuti4 = rcuti2*rcuti2
382  
383      if (summationMethod .eq. DAMPED_WOLF) then
384         if (.not.haveDWAconstants) then
# Line 369 | Line 392 | contains
392            endif
393  
394            constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
395 <          constERFC = erfc(dampingAlpha*defaultCutoff)
396 <          
395 >          constERFC = derfc(dampingAlpha*defaultCutoff)
396 >          invRootPi = 0.56418958354775628695d0
397 >          alphaPi = 2*dampingAlpha*invRootPi
398 >  
399            haveDWAconstants = .true.
400         endif
401      endif
# Line 387 | Line 412 | contains
412  
413  
414    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
415 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
415 >       vpair, fpair, pot, eFrame, f, t, do_pot)
416  
417      logical, intent(in) :: do_pot
418  
419      integer, intent(in) :: atom1, atom2
420      integer :: localError
396 –    integer, intent(in) :: corrMethod
421  
422 <    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
422 >    real(kind=dp), intent(in) :: rij, r2, sw
423      real(kind=dp), intent(in), dimension(3) :: d
424      real(kind=dp), intent(inout) :: vpair
425      real(kind=dp), intent(inout), dimension(3) :: fpair
426  
427 <    real( kind = dp ) :: pot, swi
427 >    real( kind = dp ) :: pot
428      real( kind = dp ), dimension(9,nLocal) :: eFrame
429      real( kind = dp ), dimension(3,nLocal) :: f
430      real( kind = dp ), dimension(3,nLocal) :: t
# Line 425 | Line 449 | contains
449      real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
450      real (kind=dp) :: xhat, yhat, zhat
451      real (kind=dp) :: dudx, dudy, dudz
452 <    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
453 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
452 >    real (kind=dp) :: scale, sc2, bigR
453 >    real (kind=dp) :: varERFC, varEXP
454 >    real (kind=dp) :: limScale
455  
456      if (.not.allocated(ElectrostaticMap)) then
457         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 435 | Line 460 | contains
460  
461      if (.not.summationMethodChecked) then
462         call checkSummationMethod()
463 +      
464      endif
465  
466  
# Line 449 | Line 475 | contains
475      !! some variables we'll need independent of electrostatic type:
476  
477      riji = 1.0d0 / rij
478 <
478 >  
479      xhat = d(1) * riji
480      yhat = d(2) * riji
481      zhat = d(3) * riji
482  
457 –    rcuti2 = rcuti*rcuti
458 –    rcuti3 = rcuti2*rcuti
459 –    rcuti4 = rcuti2*rcuti2
460 –
461 –    swi = 1.0d0 / sw
462 –
483      !! logicals
484      i_is_Charge = ElectrostaticMap(me1)%is_Charge
485      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 578 | Line 598 | contains
598         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
599      endif
600    
581 – !!$    switcher = 1.0d0
582 – !!$    dswitcher = 0.0d0
583 – !!$    ebalance = 0.0d0
584 – !!$    ! weaken the dipole interaction at close range for TAP water
585 – !!$    if (j_is_Tap .and. i_is_Tap) then
586 – !!$      call calc_switch(rij, mu_i, switcher, dswitcher)
587 – !!$    endif
588 –
601      epot = 0.0_dp
602      dudx = 0.0_dp
603      dudy = 0.0_dp
# Line 603 | Line 615 | contains
615  
616         if (j_is_Charge) then
617  
618 <          if (corrMethod .eq. 1) then
618 >          if (summationMethod .eq. UNDAMPED_WOLF) then
619 >
620               vterm = pre11 * q_i * q_j * (riji - rcuti)
621 +             vpair = vpair + vterm
622 +             epot = epot + sw*vterm
623 +            
624 +             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
625 +            
626 +             dudx = dudx + dudr * d(1)
627 +             dudy = dudy + dudr * d(2)
628 +             dudz = dudz + dudr * d(3)
629  
630 +          elseif (summationMethod .eq. DAMPED_WOLF) then
631 +
632 +             varERFC = derfc(dampingAlpha*rij)
633 +             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
634 +             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
635               vpair = vpair + vterm
636 <             epot = epot + sw * vterm
636 >             epot = epot + sw*vterm
637              
638 <             dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
638 >             dudr  = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
639 >                                                  + alphaPi*varEXP) &
640 >                                                 - (constERFC*rcuti2 &
641 >                                                    + alphaPi*constEXP)) )
642              
643               dudx = dudx + dudr * d(1)
644               dudy = dudy + dudr * d(2)
645               dudz = dudz + dudr * d(3)
646  
647            else
619 –             vterm = pre11 * q_i * q_j * riji
648  
649 +             vterm = pre11 * q_i * q_j * riji
650               vpair = vpair + vterm
651 <             epot = epot + sw * vterm
651 >             epot = epot + sw*vterm
652              
653               dudr  = - sw * vterm * riji
654              
# Line 633 | Line 662 | contains
662  
663         if (j_is_Dipole) then
664  
665 <          pref = sw * pre12 * q_i * mu_j
665 >          pref = pre12 * q_i * mu_j
666  
667 <          if (corrMethod .eq. 1) then
667 >          if (summationMethod .eq. UNDAMPED_WOLF) then
668               ri2 = riji * riji
669               ri3 = ri2 * riji
670  
671 +             pref = pre12 * q_i * mu_j
672               vterm = - pref * ct_j * (ri2 - rcuti2)
673 <             vpair = vpair + swi*vterm
674 <             epot = epot + vterm
673 >             vpair = vpair + vterm
674 >             epot = epot + sw*vterm
675              
676               !! this has a + sign in the () because the rij vector is
677               !! r_j - r_i and the charge-dipole potential takes the origin
678               !! as the point dipole, which is atom j in this case.
679              
680 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
680 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
681                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
682 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
682 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
683                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
684 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
684 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
685                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
686              
687 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
688 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
689 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
687 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
688 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
689 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
690  
691            else
692               if (j_is_SplitDipole) then
# Line 671 | Line 701 | contains
701               ri2 = ri * ri
702               ri3 = ri2 * ri
703               sc2 = scale * scale
704 <            
704 >
705 >             pref = pre12 * q_i * mu_j
706               vterm = - pref * ct_j * ri2 * scale
707 <             vpair = vpair + swi * vterm
708 <             epot = epot + vterm
707 >             vpair = vpair + vterm
708 >             epot = epot + sw*vterm
709              
710               !! this has a + sign in the () because the rij vector is
711               !! r_j - r_i and the charge-dipole potential takes the origin
712               !! as the point dipole, which is atom j in this case.
713              
714 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
715 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
716 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
714 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
715 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
716 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
717              
718 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
719 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
720 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
718 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
719 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
720 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
721  
722            endif
723         endif
# Line 699 | Line 730 | contains
730            cy2 = cy_j * cy_j
731            cz2 = cz_j * cz_j
732  
733 <
734 <          pref =  sw * pre14 * q_i / 3.0_dp
704 <
705 <          if (corrMethod .eq. 1) then
733 >          if (summationMethod .eq. UNDAMPED_WOLF) then
734 >             pref =  pre14 * q_i / 3.0_dp
735               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
736                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
737                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
738               vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
739                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
740                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
741 <             vpair = vpair + swi*( vterm1 - vterm2 )
742 <             epot = epot + ( vterm1 - vterm2 )
741 >             vpair = vpair + ( vterm1 - vterm2 )
742 >             epot = epot + sw*( vterm1 - vterm2 )
743              
744               dudx = dudx - (5.0_dp * &
745 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
745 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
746                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
747                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
748                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 721 | Line 750 | contains
750                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
751                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
752               dudy = dudy - (5.0_dp * &
753 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
753 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
754                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
755                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
756                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 729 | Line 758 | contains
758                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
759                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
760               dudz = dudz - (5.0_dp * &
761 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
761 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
762                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
763                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
764                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 737 | Line 766 | contains
766                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
767                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
768              
769 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
769 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
770                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
771 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
771 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
772                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
773 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
773 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
774                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
775              
776 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
776 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
777                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
778 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
778 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
779                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
780 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
780 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
781                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
782              
783 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
783 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
784                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
785 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
785 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
786                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
787 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
787 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
788                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
789          
790            else
791 +             pref =  pre14 * q_i / 3.0_dp
792               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
793                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
794                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
795 <             vpair = vpair + swi * vterm
796 <             epot = epot + vterm
795 >             vpair = vpair + vterm
796 >             epot = epot + sw*vterm
797              
798 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
798 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
799                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
800                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
801                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
802 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
802 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
803                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
804                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
805                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
806 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
806 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
807                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
808                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
809                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
810              
811 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
812 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
813 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
811 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
812 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
813 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
814              
815 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
816 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
817 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
815 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
816 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
817 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
818              
819 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
820 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
821 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
819 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
820 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
821 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
822            
823            endif
824         endif
# Line 797 | Line 827 | contains
827      if (i_is_Dipole) then
828  
829         if (j_is_Charge) then
830 <
831 <          pref = sw * pre12 * q_j * mu_i
832 <
833 <          if (corrMethod .eq. 1) then
830 >          
831 >          pref = pre12 * q_j * mu_i
832 >          
833 >          if (summationMethod .eq. UNDAMPED_WOLF) then
834               ri2 = riji * riji
835               ri3 = ri2 * riji
836  
837 +             pref = pre12 * q_j * mu_i
838               vterm = pref * ct_i * (ri2 - rcuti2)
839 <             vpair = vpair + swi * vterm
840 <             epot = epot + vterm
839 >             vpair = vpair + vterm
840 >             epot = epot + sw*vterm
841              
842               !! this has a + sign in the () because the rij vector is
843               !! r_j - r_i and the charge-dipole potential takes the origin
844               !! as the point dipole, which is atom j in this case.
845              
846 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
846 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
847                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
848 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
848 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
849                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
850 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
850 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
851                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
852              
853 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
854 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
855 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
853 >             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
854 >             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
855 >             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
856  
857            else
858               if (i_is_SplitDipole) then
# Line 836 | Line 867 | contains
867               ri2 = ri * ri
868               ri3 = ri2 * ri
869               sc2 = scale * scale
870 <            
870 >
871 >             pref = pre12 * q_j * mu_i
872               vterm = pref * ct_i * ri2 * scale
873 <             vpair = vpair + swi * vterm
874 <             epot = epot + vterm
873 >             vpair = vpair + vterm
874 >             epot = epot + sw*vterm
875              
876 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
877 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
878 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
876 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
877 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
878 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
879              
880 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
881 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
882 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
880 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
881 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
882 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
883            endif
884         endif
885 <
885 >      
886         if (j_is_Dipole) then
887  
888 <          pref = sw * pre22 * mu_i * mu_j
857 <
858 <          if (corrMethod .eq. 1) then
888 >          if (summationMethod .eq. UNDAMPED_WOLF) then
889               ri2 = riji * riji
890               ri3 = ri2 * riji
891               ri4 = ri2 * ri2
892  
893 +             pref = pre22 * mu_i * mu_j
894               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
895 <             vpair = vpair + swi * vterm
896 <             epot = epot + vterm
895 >             vpair = vpair + vterm
896 >             epot = epot + sw*vterm
897              
898               a1 = 5.0d0 * ct_i * ct_j - ct_ij
899              
900 <             dudx = dudx + pref*3.0d0*ri4 &
901 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
902 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
903 <             dudy = dudy + pref*3.0d0*ri4 &
904 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
905 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
906 <             dudz = dudz + pref*3.0d0*ri4 &
907 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
908 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
900 >             dudx = dudx + sw*pref*3.0d0*ri4 &
901 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
902 >                         - sw*pref*3.0d0*rcuti4 &
903 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
904 >             dudy = dudy + sw*pref*3.0d0*ri4 &
905 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
906 >                         - sw*pref*3.0d0*rcuti4 &
907 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
908 >             dudz = dudz + sw*pref*3.0d0*ri4 &
909 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
910 >                         - sw*pref*3.0d0*rcuti4 &
911 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
912              
913 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
913 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
914                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
915 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
915 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
916                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
917 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
917 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
918                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
919 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
919 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
920                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
921 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
921 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
922                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
923 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
923 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
924                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
925 +
926            else
892 –            
927               if (i_is_SplitDipole) then
928                  if (j_is_SplitDipole) then
929                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 916 | Line 950 | contains
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 + swi * vterm
956 <             epot = epot + vterm
955 >             vpair = vpair + vterm
956 >             epot = epot + sw*vterm
957              
958               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
959              
960 <             dudx = dudx + pref*3.0d0*ri4*scale &
961 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
962 <             dudy = dudy + pref*3.0d0*ri4*scale &
963 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
964 <             dudz = dudz + pref*3.0d0*ri4*scale &
965 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
960 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
961 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
962 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
963 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
964 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
965 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
966              
967 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
968 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
969 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
970 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
971 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
972 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
967 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
968 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
969 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
970 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
971 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
972 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
973              
974 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
975 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
976 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
977 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
978 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
979 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
974 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
975 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
976 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
977 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
978 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
979 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
980            endif
981         endif
982      endif
# Line 956 | Line 991 | contains
991            cy2 = cy_i * cy_i
992            cz2 = cz_i * cz_i
993  
994 <          pref = sw * pre14 * q_j / 3.0_dp
995 <
961 <          if (corrMethod .eq. 1) then
994 >          if (summationMethod .eq. UNDAMPED_WOLF) then
995 >             pref = pre14 * q_j / 3.0_dp
996               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
997                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
998                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
999               vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1000                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1001                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1002 <             vpair = vpair + swi * ( vterm1 - vterm2 )
1003 <             epot = epot + ( vterm1 - vterm2 )
1002 >             vpair = vpair + ( vterm1 - vterm2 )
1003 >             epot = epot + sw*( vterm1 - vterm2 )
1004              
1005 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
1006 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1005 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1006 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1007                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1008                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1009                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1010                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1011                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1012 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
1013 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1012 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1013 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1014                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1015                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1016                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1017                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1018                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1019 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
1020 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1019 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1020 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1021                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1022                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1023                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1024                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1025                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1026              
1027 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1027 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1028                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1029 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1029 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1030                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1031 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1031 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1032                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1033              
1034 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1034 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1035                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1036 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1036 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1037                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1038 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1038 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1039                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1040              
1041 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1041 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1042                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1043 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1043 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1044                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1045 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1045 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1046                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1047  
1048            else
1049 +             pref = pre14 * q_j / 3.0_dp
1050               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1051                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1052                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1053 <             vpair = vpair + swi * vterm
1054 <             epot = epot + vterm
1053 >             vpair = vpair + vterm
1054 >             epot = epot + sw*vterm
1055              
1056 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1056 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1057                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1058                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1059                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1060 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1060 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1061                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1062                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1063                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1064 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1064 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1065                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1066                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1067                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1068              
1069 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1070 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1071 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1069 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1070 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1071 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1072              
1073 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1074 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1075 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1073 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1074 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1075 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1076              
1077 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1078 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1079 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1077 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1078 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1079 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1080            endif
1081         endif
1082      endif
# Line 1049 | Line 1084 | contains
1084  
1085      if (do_pot) then
1086   #ifdef IS_MPI
1087 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1088 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1087 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1088 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1089   #else
1090         pot = pot + epot
1091   #endif

Diff Legend

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