ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90 (file contents):
Revision 2301 by gezelter, Thu Sep 15 22:05:21 2005 UTC vs.
Revision 2339 by chrisfen, Wed Sep 28 18:47:17 2005 UTC

# Line 73 | Line 73 | module electrostatic_module
73  
74    !! variables to handle different summation methods for long-range electrostatics:
75    integer, save :: summationMethod = NONE
76 +  logical, save :: summationMethodChecked = .false.
77    real(kind=DP), save :: defaultCutoff = 0.0_DP
78    logical, save :: haveDefaultCutoff = .false.
79    real(kind=DP), save :: dampingAlpha = 0.0_DP
# Line 82 | Line 83 | module electrostatic_module
83    real(kind=DP), save :: constERFC = 0.0_DP
84    real(kind=DP), save :: constEXP = 0.0_DP
85    logical, save :: haveDWAconstants = .false.
86 <
87 <
86 >  real(kind=dp), save :: rcuti = 0.0_dp
87 >  real(kind=dp), save :: rcuti2 = 0.0_dp
88 >  real(kind=dp), save :: rcuti3 = 0.0_dp
89 >  real(kind=dp), save :: rcuti4 = 0.0_dp
90 >  real(kind=dp), save :: alphaPi = 0.0_dp
91 >  real(kind=dp), save :: invRootPi = 0.0_dp
92 >  
93 > #ifdef __IFC
94 > ! error function for ifc version > 7.
95 >  double precision, external :: derfc
96 > #endif
97 >  
98    public :: setElectrostaticSummationMethod
99    public :: setElectrostaticCutoffRadius
100    public :: setDampedWolfAlpha
# Line 117 | Line 128 | contains
128   contains
129  
130    subroutine setElectrostaticSummationMethod(the_ESM)
120 –
131      integer, intent(in) :: the_ESM    
132  
133      if ((the_ESM .le. 0) .or. (the_ESM .gt. REACTION_FIELD)) then
134         call handleError("setElectrostaticSummationMethod", "Unsupported Summation Method")
135      endif
136 +
137 +    summationMethod = the_ESM
138  
139    end subroutine setElectrostaticSummationMethod
140  
# Line 356 | Line 368 | contains
368  
369    subroutine checkSummationMethod()
370  
371 +    if (.not.haveDefaultCutoff) then
372 +       call handleError("checkSummationMethod", "no Default Cutoff set!")
373 +    endif
374 +
375 +    rcuti = 1.0d0 / defaultCutoff
376 +    rcuti2 = rcuti*rcuti
377 +    rcuti3 = rcuti2*rcuti
378 +    rcuti4 = rcuti2*rcuti2
379 +
380      if (summationMethod .eq. DAMPED_WOLF) then
381         if (.not.haveDWAconstants) then
382            
# Line 363 | Line 384 | contains
384               call handleError("checkSummationMethod", "no Damping Alpha set!")
385            endif
386            
387 <          if (.not.have....)
388 <          constEXP =
389 <          constERFC =
387 >          if (.not.haveDefaultCutoff) then
388 >             call handleError("checkSummationMethod", "no Default Cutoff set!")
389 >          endif
390 >
391 >          constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff)
392 >          constERFC = derfc(dampingAlpha*defaultCutoff)
393 >          invRootPi = 0.56418958354775628695d0
394 >          alphaPi = 2*dampingAlpha*invRootPi
395            
396            haveDWAconstants = .true.
397         endif
398      endif
399  
400 +    if (summationMethod .eq. REACTION_FIELD) then
401 +       if (.not.haveDielectric) then
402 +          call handleError("checkSummationMethod", "no reaction field Dielectric set!")
403 +       endif
404 +    endif
405 +
406 +    summationMethodChecked = .true.
407    end subroutine checkSummationMethod
408  
409  
410  
411    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
412 <       vpair, fpair, pot, eFrame, f, t, do_pot, corrMethod, rcuti)
412 >       vpair, fpair, pot, eFrame, f, t, do_pot)
413  
414      logical, intent(in) :: do_pot
415  
416      integer, intent(in) :: atom1, atom2
417      integer :: localError
385 –    integer, intent(in) :: corrMethod
418  
419 <    real(kind=dp), intent(in) :: rij, r2, sw, rcuti
419 >    real(kind=dp), intent(in) :: rij, r2, sw
420      real(kind=dp), intent(in), dimension(3) :: d
421      real(kind=dp), intent(inout) :: vpair
422      real(kind=dp), intent(inout), dimension(3) :: fpair
423  
424 <    real( kind = dp ) :: pot, swi
424 >    real( kind = dp ) :: pot
425      real( kind = dp ), dimension(9,nLocal) :: eFrame
426      real( kind = dp ), dimension(3,nLocal) :: f
427      real( kind = dp ), dimension(3,nLocal) :: t
# Line 414 | Line 446 | contains
446      real (kind=dp) :: pref, vterm, epot, dudr, vterm1, vterm2
447      real (kind=dp) :: xhat, yhat, zhat
448      real (kind=dp) :: dudx, dudy, dudz
449 <    real (kind=dp) :: scale, sc2, bigR, switcher, dswitcher
450 <    real (kind=dp) :: rcuti2, rcuti3, rcuti4
449 >    real (kind=dp) :: scale, sc2, bigR
450 >    real (kind=dp) :: varERFC, varEXP
451  
452      if (.not.allocated(ElectrostaticMap)) then
453         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 424 | Line 456 | contains
456  
457      if (.not.summationMethodChecked) then
458         call checkSummationMethod()
459 +      
460      endif
461  
462  
# Line 443 | Line 476 | contains
476      yhat = d(2) * riji
477      zhat = d(3) * riji
478  
446 –    rcuti2 = rcuti*rcuti
447 –    rcuti3 = rcuti2*rcuti
448 –    rcuti4 = rcuti2*rcuti2
449 –
450 –    swi = 1.0d0 / sw
451 –
479      !! logicals
480      i_is_Charge = ElectrostaticMap(me1)%is_Charge
481      i_is_Dipole = ElectrostaticMap(me1)%is_Dipole
# Line 567 | Line 594 | contains
594         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
595      endif
596    
570 – !!$    switcher = 1.0d0
571 – !!$    dswitcher = 0.0d0
572 – !!$    ebalance = 0.0d0
573 – !!$    ! weaken the dipole interaction at close range for TAP water
574 – !!$    if (j_is_Tap .and. i_is_Tap) then
575 – !!$      call calc_switch(rij, mu_i, switcher, dswitcher)
576 – !!$    endif
577 –
597      epot = 0.0_dp
598      dudx = 0.0_dp
599      dudy = 0.0_dp
# Line 592 | Line 611 | contains
611  
612         if (j_is_Charge) then
613  
614 <          if (corrMethod .eq. 1) then
614 >          if (summationMethod .eq. UNDAMPED_WOLF) then
615 >
616               vterm = pre11 * q_i * q_j * (riji - rcuti)
617 +             vpair = vpair + vterm
618 +             epot = epot + sw*vterm
619 +            
620 +             dudr  = -sw*pre11*q_i*q_j * (riji*riji*riji - rcuti2*rcuti)
621 +            
622 +             dudx = dudx + dudr * d(1)
623 +             dudy = dudy + dudr * d(2)
624 +             dudz = dudz + dudr * d(3)
625  
626 +          elseif (summationMethod .eq. DAMPED_WOLF) then
627 +
628 +             varERFC = derfc(dampingAlpha*rij)
629 +             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
630 +             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
631               vpair = vpair + vterm
632 <             epot = epot + sw * vterm
632 >             epot = epot + sw*vterm
633              
634 <             dudr  = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti)
634 >             dudr  = -sw*pre11*q_i*q_j * ( riji*(varERFC*riji*riji &
635 >                                                 + alphaPi*varEXP) &
636 >                                         - rcuti*(constERFC*rcuti2 &
637 >                                                 + alphaPi*constEXP) )
638              
639               dudx = dudx + dudr * d(1)
640               dudy = dudy + dudr * d(2)
641               dudz = dudz + dudr * d(3)
642  
643            else
608 –             vterm = pre11 * q_i * q_j * riji
644  
645 +             vterm = pre11 * q_i * q_j * riji
646               vpair = vpair + vterm
647 <             epot = epot + sw * vterm
647 >             epot = epot + sw*vterm
648              
649               dudr  = - sw * vterm * riji
650              
# Line 622 | Line 658 | contains
658  
659         if (j_is_Dipole) then
660  
661 <          pref = sw * pre12 * q_i * mu_j
661 >          pref = pre12 * q_i * mu_j
662  
663 <          if (corrMethod .eq. 1) then
663 >          if (summationMethod .eq. UNDAMPED_WOLF) then
664               ri2 = riji * riji
665               ri3 = ri2 * riji
666  
667 +             pref = pre12 * q_i * mu_j
668               vterm = - pref * ct_j * (ri2 - rcuti2)
669 <             vpair = vpair + swi*vterm
670 <             epot = epot + vterm
669 >             vpair = vpair + vterm
670 >             epot = epot + sw*vterm
671              
672               !! this has a + sign in the () because the rij vector is
673               !! r_j - r_i and the charge-dipole potential takes the origin
674               !! as the point dipole, which is atom j in this case.
675              
676 <             dudx = dudx - pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
676 >             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
677                    - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
678 <             dudy = dudy - pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
678 >             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
679                    - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
680 <             dudz = dudz - pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
680 >             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
681                    - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
682              
683 <             duduz_j(1) = duduz_j(1) - pref*( ri2*xhat - d(1)*rcuti3 )
684 <             duduz_j(2) = duduz_j(2) - pref*( ri2*yhat - d(2)*rcuti3 )
685 <             duduz_j(3) = duduz_j(3) - pref*( ri2*zhat - d(3)*rcuti3 )
683 >             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
684 >             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
685 >             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
686  
687            else
688               if (j_is_SplitDipole) then
# Line 660 | Line 697 | contains
697               ri2 = ri * ri
698               ri3 = ri2 * ri
699               sc2 = scale * scale
700 <            
700 >
701 >             pref = pre12 * q_i * mu_j
702               vterm = - pref * ct_j * ri2 * scale
703 <             vpair = vpair + swi * vterm
704 <             epot = epot + vterm
703 >             vpair = vpair + vterm
704 >             epot = epot + sw*vterm
705              
706               !! this has a + sign in the () because the rij vector is
707               !! r_j - r_i and the charge-dipole potential takes the origin
708               !! as the point dipole, which is atom j in this case.
709              
710 <             dudx = dudx - pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
711 <             dudy = dudy - pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
712 <             dudz = dudz - pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
710 >             dudx = dudx - sw*pref * ri3 * ( uz_j(1) - 3.0d0*ct_j*xhat*sc2)
711 >             dudy = dudy - sw*pref * ri3 * ( uz_j(2) - 3.0d0*ct_j*yhat*sc2)
712 >             dudz = dudz - sw*pref * ri3 * ( uz_j(3) - 3.0d0*ct_j*zhat*sc2)
713              
714 <             duduz_j(1) = duduz_j(1) - pref * ri2 * xhat * scale
715 <             duduz_j(2) = duduz_j(2) - pref * ri2 * yhat * scale
716 <             duduz_j(3) = duduz_j(3) - pref * ri2 * zhat * scale
714 >             duduz_j(1) = duduz_j(1) - sw*pref * ri2 * xhat * scale
715 >             duduz_j(2) = duduz_j(2) - sw*pref * ri2 * yhat * scale
716 >             duduz_j(3) = duduz_j(3) - sw*pref * ri2 * zhat * scale
717  
718            endif
719         endif
# Line 688 | Line 726 | contains
726            cy2 = cy_j * cy_j
727            cz2 = cz_j * cz_j
728  
729 <
730 <          pref =  sw * pre14 * q_i / 3.0_dp
693 <
694 <          if (corrMethod .eq. 1) then
729 >          if (summationMethod .eq. UNDAMPED_WOLF) then
730 >             pref =  pre14 * q_i / 3.0_dp
731               vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
732                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
733                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
734               vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
735                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
736                    qzz_j * (3.0_dp*cz2 - 1.0_dp) )
737 <             vpair = vpair + swi*( vterm1 - vterm2 )
738 <             epot = epot + ( vterm1 - vterm2 )
737 >             vpair = vpair + ( vterm1 - vterm2 )
738 >             epot = epot + sw*( vterm1 - vterm2 )
739              
740               dudx = dudx - (5.0_dp * &
741 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + pref * ( &
741 >                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
742                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
743                    qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
744                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
# Line 710 | Line 746 | contains
746                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
747                    qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
748               dudy = dudy - (5.0_dp * &
749 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + pref * ( &
749 >                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
750                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
751                    qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
752                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
# Line 718 | Line 754 | contains
754                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
755                    qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
756               dudz = dudz - (5.0_dp * &
757 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + pref * ( &
757 >                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
758                    (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
759                    qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
760                    (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
# Line 726 | Line 762 | contains
762                    (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
763                    qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
764              
765 <             dudux_j(1) = dudux_j(1) + pref * (ri3*(qxx_j*6.0_dp*cx_j*xhat) - &
765 >             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
766                    rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
767 <             dudux_j(2) = dudux_j(2) + pref * (ri3*(qxx_j*6.0_dp*cx_j*yhat) - &
767 >             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
768                    rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
769 <             dudux_j(3) = dudux_j(3) + pref * (ri3*(qxx_j*6.0_dp*cx_j*zhat) - &
769 >             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
770                    rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
771              
772 <             duduy_j(1) = duduy_j(1) + pref * (ri3*(qyy_j*6.0_dp*cy_j*xhat) - &
772 >             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
773                    rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
774 <             duduy_j(2) = duduy_j(2) + pref * (ri3*(qyy_j*6.0_dp*cy_j*yhat) - &
774 >             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
775                    rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
776 <             duduy_j(3) = duduy_j(3) + pref * (ri3*(qyy_j*6.0_dp*cy_j*zhat) - &
776 >             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
777                    rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
778              
779 <             duduz_j(1) = duduz_j(1) + pref * (ri3*(qzz_j*6.0_dp*cz_j*xhat) - &
779 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
780                    rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
781 <             duduz_j(2) = duduz_j(2) + pref * (ri3*(qzz_j*6.0_dp*cz_j*yhat) - &
781 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
782                    rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
783 <             duduz_j(3) = duduz_j(3) + pref * (ri3*(qzz_j*6.0_dp*cz_j*zhat) - &
783 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
784                    rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
785          
786            else
787 +             pref =  pre14 * q_i / 3.0_dp
788               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
789                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
790                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
791 <             vpair = vpair + swi * vterm
792 <             epot = epot + vterm
791 >             vpair = vpair + vterm
792 >             epot = epot + sw*vterm
793              
794 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
794 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref * ri4 * ( &
795                    qxx_j*(6.0_dp*cx_j*ux_j(1) - 2.0_dp*xhat) + &
796                    qyy_j*(6.0_dp*cy_j*uy_j(1) - 2.0_dp*xhat) + &
797                    qzz_j*(6.0_dp*cz_j*uz_j(1) - 2.0_dp*xhat) )
798 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
798 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref * ri4 * ( &
799                    qxx_j*(6.0_dp*cx_j*ux_j(2) - 2.0_dp*yhat) + &
800                    qyy_j*(6.0_dp*cy_j*uy_j(2) - 2.0_dp*yhat) + &
801                    qzz_j*(6.0_dp*cz_j*uz_j(2) - 2.0_dp*yhat) )
802 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
802 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref * ri4 * ( &
803                    qxx_j*(6.0_dp*cx_j*ux_j(3) - 2.0_dp*zhat) + &
804                    qyy_j*(6.0_dp*cy_j*uy_j(3) - 2.0_dp*zhat) + &
805                    qzz_j*(6.0_dp*cz_j*uz_j(3) - 2.0_dp*zhat) )
806              
807 <             dudux_j(1) = dudux_j(1) + pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
808 <             dudux_j(2) = dudux_j(2) + pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
809 <             dudux_j(3) = dudux_j(3) + pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
807 >             dudux_j(1) = dudux_j(1) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*xhat)
808 >             dudux_j(2) = dudux_j(2) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*yhat)
809 >             dudux_j(3) = dudux_j(3) + sw*pref * ri3*(qxx_j*6.0_dp*cx_j*zhat)
810              
811 <             duduy_j(1) = duduy_j(1) + pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
812 <             duduy_j(2) = duduy_j(2) + pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
813 <             duduy_j(3) = duduy_j(3) + pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
811 >             duduy_j(1) = duduy_j(1) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*xhat)
812 >             duduy_j(2) = duduy_j(2) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*yhat)
813 >             duduy_j(3) = duduy_j(3) + sw*pref * ri3*(qyy_j*6.0_dp*cy_j*zhat)
814              
815 <             duduz_j(1) = duduz_j(1) + pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
816 <             duduz_j(2) = duduz_j(2) + pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
817 <             duduz_j(3) = duduz_j(3) + pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
815 >             duduz_j(1) = duduz_j(1) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*xhat)
816 >             duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
817 >             duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
818            
819            endif
820         endif
# Line 786 | Line 823 | contains
823      if (i_is_Dipole) then
824  
825         if (j_is_Charge) then
826 <
827 <          pref = sw * pre12 * q_j * mu_i
828 <
829 <          if (corrMethod .eq. 1) then
826 >          
827 >          pref = pre12 * q_j * mu_i
828 >          
829 >          if (summationMethod .eq. UNDAMPED_WOLF) then
830               ri2 = riji * riji
831               ri3 = ri2 * riji
832  
833 +             pref = pre12 * q_j * mu_i
834               vterm = pref * ct_i * (ri2 - rcuti2)
835 <             vpair = vpair + swi * vterm
836 <             epot = epot + vterm
835 >             vpair = vpair + vterm
836 >             epot = epot + sw*vterm
837              
838               !! this has a + sign in the () because the rij vector is
839               !! r_j - r_i and the charge-dipole potential takes the origin
840               !! as the point dipole, which is atom j in this case.
841              
842 <             dudx = dudx + pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
842 >             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
843                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
844 <             dudy = dudy + pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
844 >             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
845                    - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
846 <             dudz = dudz + pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
846 >             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
847                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
848              
849 <             duduz_i(1) = duduz_i(1) - pref*( ri2*xhat - d(1)*rcuti3 )
850 <             duduz_i(2) = duduz_i(2) - pref*( ri2*yhat - d(2)*rcuti3 )
851 <             duduz_i(3) = duduz_i(3) - pref*( ri2*zhat - d(3)*rcuti3 )
849 >             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
850 >             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
851 >             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
852  
853            else
854               if (i_is_SplitDipole) then
# Line 825 | Line 863 | contains
863               ri2 = ri * ri
864               ri3 = ri2 * ri
865               sc2 = scale * scale
866 <            
866 >
867 >             pref = pre12 * q_j * mu_i
868               vterm = pref * ct_i * ri2 * scale
869 <             vpair = vpair + swi * vterm
870 <             epot = epot + vterm
869 >             vpair = vpair + vterm
870 >             epot = epot + sw*vterm
871              
872 <             dudx = dudx + pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
873 <             dudy = dudy + pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
874 <             dudz = dudz + pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
872 >             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0 * ct_i * xhat*sc2)
873 >             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0 * ct_i * yhat*sc2)
874 >             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0 * ct_i * zhat*sc2)
875              
876 <             duduz_i(1) = duduz_i(1) + pref * ri2 * xhat * scale
877 <             duduz_i(2) = duduz_i(2) + pref * ri2 * yhat * scale
878 <             duduz_i(3) = duduz_i(3) + pref * ri2 * zhat * scale
876 >             duduz_i(1) = duduz_i(1) + sw*pref * ri2 * xhat * scale
877 >             duduz_i(2) = duduz_i(2) + sw*pref * ri2 * yhat * scale
878 >             duduz_i(3) = duduz_i(3) + sw*pref * ri2 * zhat * scale
879            endif
880         endif
881 <
881 >      
882         if (j_is_Dipole) then
883  
884 <          pref = sw * pre22 * mu_i * mu_j
846 <
847 <          if (corrMethod .eq. 1) then
884 >          if (summationMethod .eq. UNDAMPED_WOLF) then
885               ri2 = riji * riji
886               ri3 = ri2 * riji
887               ri4 = ri2 * ri2
888  
889 +             pref = pre22 * mu_i * mu_j
890               vterm = pref * (ri3 - rcuti3) * (ct_ij - 3.0d0 * ct_i * ct_j)
891 <             vpair = vpair + swi * vterm
892 <             epot = epot + vterm
891 >             vpair = vpair + vterm
892 >             epot = epot + sw*vterm
893              
894               a1 = 5.0d0 * ct_i * ct_j - ct_ij
895              
896 <             dudx = dudx + pref*3.0d0*ri4 &
897 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) - &
898 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
899 <             dudy = dudy + pref*3.0d0*ri4 &
900 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) - &
901 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
902 <             dudz = dudz + pref*3.0d0*ri4 &
903 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) - &
904 <                  pref*3.0d0*rcuti4*(a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
896 >             dudx = dudx + sw*pref*3.0d0*ri4 &
897 >                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1)) &
898 >                         - sw*pref*3.0d0*rcuti4 &
899 >                             * (a1*rcuti*d(1)-ct_i*uz_j(1)-ct_j*uz_i(1))
900 >             dudy = dudy + sw*pref*3.0d0*ri4 &
901 >                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2)) &
902 >                         - sw*pref*3.0d0*rcuti4 &
903 >                             * (a1*rcuti*d(2)-ct_i*uz_j(2)-ct_j*uz_i(2))
904 >             dudz = dudz + sw*pref*3.0d0*ri4 &
905 >                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3)) &
906 >                         - sw*pref*3.0d0*rcuti4 &
907 >                             * (a1*rcuti*d(3)-ct_i*uz_j(3)-ct_j*uz_i(3))
908              
909 <             duduz_i(1) = duduz_i(1) + pref*(ri3*(uz_j(1) - 3.0d0*ct_j*xhat) &
909 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
910                    - rcuti3*(uz_j(1) - 3.0d0*ct_j*d(1)*rcuti))
911 <             duduz_i(2) = duduz_i(2) + pref*(ri3*(uz_j(2) - 3.0d0*ct_j*yhat) &
911 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
912                    - rcuti3*(uz_j(2) - 3.0d0*ct_j*d(2)*rcuti))
913 <             duduz_i(3) = duduz_i(3) + pref*(ri3*(uz_j(3) - 3.0d0*ct_j*zhat) &
913 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
914                    - rcuti3*(uz_j(3) - 3.0d0*ct_j*d(3)*rcuti))
915 <             duduz_j(1) = duduz_j(1) + pref*(ri3*(uz_i(1) - 3.0d0*ct_i*xhat) &
915 >             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
916                    - rcuti3*(uz_i(1) - 3.0d0*ct_i*d(1)*rcuti))
917 <             duduz_j(2) = duduz_j(2) + pref*(ri3*(uz_i(2) - 3.0d0*ct_i*yhat) &
917 >             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
918                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
919 <             duduz_j(3) = duduz_j(3) + pref*(ri3*(uz_i(3) - 3.0d0*ct_i*zhat) &
919 >             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
920                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
921 +
922            else
881 –            
923               if (i_is_SplitDipole) then
924                  if (j_is_SplitDipole) then
925                     BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
# Line 905 | Line 946 | contains
946               ri4 = ri2 * ri2
947               sc2 = scale * scale
948              
949 +             pref = pre22 * mu_i * mu_j
950               vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
951 <             vpair = vpair + swi * vterm
952 <             epot = epot + vterm
951 >             vpair = vpair + vterm
952 >             epot = epot + sw*vterm
953              
954               a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
955              
956 <             dudx = dudx + pref*3.0d0*ri4*scale &
957 <                  *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
958 <             dudy = dudy + pref*3.0d0*ri4*scale &
959 <                  *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
960 <             dudz = dudz + pref*3.0d0*ri4*scale &
961 <                  *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
956 >             dudx = dudx + sw*pref*3.0d0*ri4*scale &
957 >                             *(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
958 >             dudy = dudy + sw*pref*3.0d0*ri4*scale &
959 >                             *(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
960 >             dudz = dudz + sw*pref*3.0d0*ri4*scale &
961 >                             *(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
962              
963 <             duduz_i(1) = duduz_i(1) + pref*ri3 &
964 <                  *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
965 <             duduz_i(2) = duduz_i(2) + pref*ri3 &
966 <                  *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
967 <             duduz_i(3) = duduz_i(3) + pref*ri3 &
968 <                  *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
963 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3 &
964 >                                         *(uz_j(1) - 3.0d0*ct_j*xhat*sc2)
965 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3 &
966 >                                         *(uz_j(2) - 3.0d0*ct_j*yhat*sc2)
967 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3 &
968 >                                         *(uz_j(3) - 3.0d0*ct_j*zhat*sc2)
969              
970 <             duduz_j(1) = duduz_j(1) + pref*ri3 &
971 <                  *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
972 <             duduz_j(2) = duduz_j(2) + pref*ri3 &
973 <                  *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
974 <             duduz_j(3) = duduz_j(3) + pref*ri3 &
975 <                  *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
970 >             duduz_j(1) = duduz_j(1) + sw*pref*ri3 &
971 >                                         *(uz_i(1) - 3.0d0*ct_i*xhat*sc2)
972 >             duduz_j(2) = duduz_j(2) + sw*pref*ri3 &
973 >                                         *(uz_i(2) - 3.0d0*ct_i*yhat*sc2)
974 >             duduz_j(3) = duduz_j(3) + sw*pref*ri3 &
975 >                                         *(uz_i(3) - 3.0d0*ct_i*zhat*sc2)
976            endif
977         endif
978      endif
# Line 945 | Line 987 | contains
987            cy2 = cy_i * cy_i
988            cz2 = cz_i * cz_i
989  
990 <          pref = sw * pre14 * q_j / 3.0_dp
991 <
950 <          if (corrMethod .eq. 1) then
990 >          if (summationMethod .eq. UNDAMPED_WOLF) then
991 >             pref = pre14 * q_j / 3.0_dp
992               vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
993                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
994                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
995               vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
996                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
997                    qzz_i * (3.0_dp*cz2 - 1.0_dp) )
998 <             vpair = vpair + swi * ( vterm1 - vterm2 )
999 <             epot = epot + ( vterm1 - vterm2 )
998 >             vpair = vpair + ( vterm1 - vterm2 )
999 >             epot = epot + sw*( vterm1 - vterm2 )
1000              
1001 <             dudx = dudx - (5.0_dp*(vterm1*riji*xhat - vterm2*rcuti2*d(1))) + &
1002 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1001 >             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1002 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1003                    qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1004                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1005                    qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1006                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1007                    qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1008 <             dudy = dudy - (5.0_dp*(vterm1*riji*yhat - vterm2*rcuti2*d(2))) + &
1009 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1008 >             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1009 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1010                    qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1011                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1012                    qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1013                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1014                    qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1015 <             dudz = dudz - (5.0_dp*(vterm1*riji*zhat - vterm2*rcuti2*d(3))) + &
1016 <                  pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1015 >             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1016 >                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1017                    qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1018                    (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1019                    qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1020                    (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1021                    qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1022              
1023 <             dudux_i(1) = dudux_i(1) + pref * (ri3*(qxx_i*6.0_dp*cx_i*xhat) - &
1023 >             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1024                    rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1025 <             dudux_i(2) = dudux_i(2) + pref * (ri3*(qxx_i*6.0_dp*cx_i*yhat) - &
1025 >             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1026                    rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1027 <             dudux_i(3) = dudux_i(3) + pref * (ri3*(qxx_i*6.0_dp*cx_i*zhat) - &
1027 >             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1028                    rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1029              
1030 <             duduy_i(1) = duduy_i(1) + pref * (ri3*(qyy_i*6.0_dp*cy_i*xhat) - &
1030 >             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1031                    rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1032 <             duduy_i(2) = duduy_i(2) + pref * (ri3*(qyy_i*6.0_dp*cy_i*yhat) - &
1032 >             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1033                    rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1034 <             duduy_i(3) = duduy_i(3) + pref * (ri3*(qyy_i*6.0_dp*cy_i*zhat) - &
1034 >             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1035                    rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1036              
1037 <             duduz_i(1) = duduz_i(1) + pref * (ri3*(qzz_i*6.0_dp*cz_i*xhat) - &
1037 >             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1038                    rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1039 <             duduz_i(2) = duduz_i(2) + pref * (ri3*(qzz_i*6.0_dp*cz_i*yhat) - &
1039 >             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1040                    rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1041 <             duduz_i(3) = duduz_i(3) + pref * (ri3*(qzz_i*6.0_dp*cz_i*zhat) - &
1041 >             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1042                    rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1043  
1044            else
1045 +             pref = pre14 * q_j / 3.0_dp
1046               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1047                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1048                    qzz_i * (3.0_dp*cz2 - 1.0_dp))
1049 <             vpair = vpair + swi * vterm
1050 <             epot = epot + vterm
1049 >             vpair = vpair + vterm
1050 >             epot = epot + sw*vterm
1051              
1052 <             dudx = dudx - 5.0_dp*vterm*riji*xhat + pref * ri4 * ( &
1052 >             dudx = dudx - 5.0_dp*sw*vterm*riji*xhat + sw*pref*ri4 * ( &
1053                    qxx_i*(6.0_dp*cx_i*ux_i(1) - 2.0_dp*xhat) + &
1054                    qyy_i*(6.0_dp*cy_i*uy_i(1) - 2.0_dp*xhat) + &
1055                    qzz_i*(6.0_dp*cz_i*uz_i(1) - 2.0_dp*xhat) )
1056 <             dudy = dudy - 5.0_dp*vterm*riji*yhat + pref * ri4 * ( &
1056 >             dudy = dudy - 5.0_dp*sw*vterm*riji*yhat + sw*pref*ri4 * ( &
1057                    qxx_i*(6.0_dp*cx_i*ux_i(2) - 2.0_dp*yhat) + &
1058                    qyy_i*(6.0_dp*cy_i*uy_i(2) - 2.0_dp*yhat) + &
1059                    qzz_i*(6.0_dp*cz_i*uz_i(2) - 2.0_dp*yhat) )
1060 <             dudz = dudz - 5.0_dp*vterm*riji*zhat + pref * ri4 * ( &
1060 >             dudz = dudz - 5.0_dp*sw*vterm*riji*zhat + sw*pref*ri4 * ( &
1061                    qxx_i*(6.0_dp*cx_i*ux_i(3) - 2.0_dp*zhat) + &
1062                    qyy_i*(6.0_dp*cy_i*uy_i(3) - 2.0_dp*zhat) + &
1063                    qzz_i*(6.0_dp*cz_i*uz_i(3) - 2.0_dp*zhat) )
1064              
1065 <             dudux_i(1) = dudux_i(1) + pref * ri3*(qxx_i*6.0_dp*cx_i*xhat)
1066 <             dudux_i(2) = dudux_i(2) + pref * ri3*(qxx_i*6.0_dp*cx_i*yhat)
1067 <             dudux_i(3) = dudux_i(3) + pref * ri3*(qxx_i*6.0_dp*cx_i*zhat)
1065 >             dudux_i(1) = dudux_i(1) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*xhat)
1066 >             dudux_i(2) = dudux_i(2) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*yhat)
1067 >             dudux_i(3) = dudux_i(3) + sw*pref*ri3*(qxx_i*6.0_dp*cx_i*zhat)
1068              
1069 <             duduy_i(1) = duduy_i(1) + pref * ri3*(qyy_i*6.0_dp*cy_i*xhat)
1070 <             duduy_i(2) = duduy_i(2) + pref * ri3*(qyy_i*6.0_dp*cy_i*yhat)
1071 <             duduy_i(3) = duduy_i(3) + pref * ri3*(qyy_i*6.0_dp*cy_i*zhat)
1069 >             duduy_i(1) = duduy_i(1) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*xhat)
1070 >             duduy_i(2) = duduy_i(2) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*yhat)
1071 >             duduy_i(3) = duduy_i(3) + sw*pref*ri3*(qyy_i*6.0_dp*cy_i*zhat)
1072              
1073 <             duduz_i(1) = duduz_i(1) + pref * ri3*(qzz_i*6.0_dp*cz_i*xhat)
1074 <             duduz_i(2) = duduz_i(2) + pref * ri3*(qzz_i*6.0_dp*cz_i*yhat)
1075 <             duduz_i(3) = duduz_i(3) + pref * ri3*(qzz_i*6.0_dp*cz_i*zhat)
1073 >             duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1074 >             duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1075 >             duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1076            endif
1077         endif
1078      endif

Diff Legend

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